Local primordial non-Gaussianity from the large-scale clustering of photometric DESI luminous red galaxies

We use angular clustering of luminous red galaxies from the Dark Energy Spectroscopic Instrument (DESI) imaging surveys to constrain the local primordial non-Gaussianity parameter $\fnl$. Our sample comprises over 12 million targets, covering 14,000 square degrees of the sky, with redshifts in the range $0.2<z<1.35$. We identify Galactic extinction, survey depth, and astronomical seeing as the primary sources of systematic error, and employ linear regression and artificial neural networks to alleviate non-cosmological excess clustering on large scales. Our methods are tested against simulations with and without $\fnl$ and systematics, showing superior performance of the neural network treatment. The neural network with a set of nine imaging property maps passes our systematic null test criteria, and is chosen as the fiducial treatment. Assuming the universality relation, we find $\fnl = 34^{+24(+50)}_{-44(-73)}$ at 68\%(95\%) confidence. We apply a series of robustness tests (e.g., cuts on imaging, declination, or scales used) that show consistency in the obtained constraints. We study how the regression method biases the measured angular power-spectrum and degrades the $\fnl$ constraining power. The use of the nine maps more than doubles the uncertainty compared to using only the three primary maps in the regression. Our results thus motivate the development of more efficient methods that avoid over-correction, protect large-scale clustering information, and preserve constraining power. Additionally, our results encourage further studies of $\fnl$ with DESI spectroscopic samples, where the inclusion of 3D clustering modes should help separate imaging systematics and lessen the degradation in the $\fnl$ uncertainty.


INTRODUCTION
Inflation is a widely accepted paradigm in modern cosmology that explains many important characteristics of our Universe.It predicts that the early Universe underwent a period of accelerated expansion, resulting in the observed homogeneity and isotropy of the Universe on large scales (Guth 1981;Linde 1982;Albrecht & Steinhardt 1982).After the period of inflation, the Universe entered a phase of reheating in which primordial perturbations were generated, setting the initial seeds for structure formation (Kofman et al. 1994;Bassett et al. 2006;Lyth & Liddle 2009).Although inflation is widely accepted as a compelling explanation, the characteristics of the field or fields that drove the inflationary expansion remain largely unknown in cosmology.While early studies of the cosmic microwave background (CMB) and large-scale structure (LSS) suggested that primordial fluctuations are both Gaussian and scaleinvariant (Komatsu et al. 2003;Tegmark et al. 2004;Guth & Kaiser 2005), some alternative classes of inflationary models predict different levels of non-Gaussianities in the primordial gravitational field.Non-Gaussianities are a measure of the degree to which the distribution of matter in the Universe deviates from a Gaussian distribution, which would have important implications for the growth of structure and galaxies in the Universe (see, e.g., Verde 2010;Desjacques & Seljak 2010;Biagetti 2019).
In its simplest form, local primordial non-Gaussianity (PNG) is parameterized by the non-linear coupling constant  NL (Komatsu & Spergel 2001): where Φ is the primordial curvature perturbation and  is assumed to be a Gaussian random field.Local-type PNG generates a primordial bispectrum, which peaks in the squeezed triangle configuration where one of the three wave vectors is much smaller than the other two.This means that one of the modes is on a much larger scale than the other two, and this mode couples with the other two modes to generate a non-Gaussian signal, which then affects the local number density of galaxies.The coupling between the short and long wave-lengths induces a distinct bias in the galaxy distribution, which leads to a  −2 -dependent feature in the two-point clustering of galaxies and quasars (Dalal et al. 2008).Obtaining reliable, accurate, and robust constraints on  NL is crucial in advancing our understanding of the dynamics of the early Universe.For instance, the standard single-field slow-roll inflationary model predicts a small value of  NL ∼ 0.01 (see, e.g., Maldacena 2003).On the other hand, some alternative inflationary scenarios involve multiple scalar fields that can interact with each other during inflation, leading to the generation of larger levels of non-Gaussianities.These models predict considerably larger values of  NL that can reach up to 100 or higher (see, e.g., Chen 2010, for a review).With (  NL ) ∼ 1, we can rule out or confirm specific models of inflation and gain insight into the physics that drove the inflationary expansion (see, e.g., Alvarez et al. 2014;de Putter et al. 2017).
The current tightest bound on  NL comes from Planck's bispectrum measurement of CMB anisotropies,  NL = 0.9 ± 5.1 (Planck Collaboration et al. 2019).Limited by cosmic variance, CMB data cannot enhance the statistical precision of  NL measurements enough to break the degeneracy amongst various inflationary paradigms (see, e.g., Abazajian et al. 2016;Simons Observatory et al. 2019).On the other hand, LSS surveys probe a 3D map of the Universe, and thus provide more modes to limit  NL .However, nonlinearities raised from structure formation pose a serious challenge for measuring  NL with the three-point clustering of galaxies, and these nonlinear effects are non-trivial to model and disentangle from the primordial signal (Baldauf et al. 2011b,a).Currently, the most precise constraints on  NL from LSS reach a level of (  NL ) ∼ 20 − 30, with the majority of the constraining power coming from the two-point clustering statistics that utilize the scaledependent bias effect (Slosar et al. 2008;Ross et al. 2013;Castorina et al. 2019;Mueller et al. 2022;Cabass et al. 2022;D'Amico et al. 2022).Surveying large areas of the sky can unlock more modes and help improve these constraints.
The Dark Energy Spectroscopic Instrument (DESI) is ideally suited to enable excellent constraints on primordial non-Gaussianity from the galaxy distribution.DESI uses 5000 robotically-driven fibers to simultaneously collect spectra of extra-galactic objects (Levi et al. 2013;DESI Collaboration et al. 2016b;Silber et al. 2023).DESI is designed to deliver an unparalleled volume of spectroscopic data covering ∼ 14, 000 square degrees that promises to deepen our understanding of the energy contents of the Universe, neutrino masses, and the nature of gravity (DESI Collaboration et al. 2022).Moreover, DESI alone is expected to improve our constraints on local PNG down to (  NL ) = 5, assuming systematic uncertainties are under control (DESI Collaboration et al. 2016a).With multi-tracer techniques (Seljak 2009), cosmic variance can be further reduced to allow surpassing CMB-like constraints (Alonso et al. 2015).For instance, the distortion of CMB photons around foreground masses, which is referred to as CMB lensing, provides an additional probe of LSS, but from a different vantage point.We can significantly reduce statistical uncertainties below (  NL ) ∼ 1 by cross-correlating LSS data with CMB-lensing, or other tracers of matter, such as 21 cm intensity mapping (see, e.g., Schmittfull & Seljak 2018;Heinrich & Doré 2022;Jolicoeur et al. 2023;Sullivan et al. 2023).
However, further work is needed to fully harness the potential of the scale-dependent bias effect in constraining  NL with LSS.The amplitude of the  NL signal in the galaxy distribution is proportional to the bias parameter   , such that Δ ∝    NL  −2 .Assuming the universality relation,   ∼ ( − ), where  is the linear halo bias and  = 1 is a parameter that describes the response of galaxy formation to primordial potential perturbations in the presence of local PNG (see, e.g., Slosar et al. 2008).The value of  is not very well constrained for other tracers of matter (Barreira et al. 2020;Barreira 2020), and Barreira (2022) showed that marginalizing over  even with wide priors leads to biased  NL constraints because of parameter space projection effects.More simulation-based studies are necessary to investigate the halo-assembly bias and the relationship between   and  for various galaxy samples.For instance, Lazeyras et al. (2023) used N-body simulations to investigate secondary halo properties, such as concentration, spin and sphericity of haloes, and found that halo spin and sphericity preserve the universality of the halo occupation function while halo concentration significantly alters the halo function.Without better-informed priors on , it is argued that the scale-dependent bias effect can only be used to constrain the    NL term (see, e.g., Barreira 2020).However, regardless of the specific value of , a nonzero detection of    NL implies the presence of local PNG, given that   is greater than zero.In this work, we assume the universality relation that links   to  −  and, further, fix the value of .
In addition to the theoretical uncertainties, measuring  NL through the scale-dependent bias effect is a difficult task due to various imaging systematic effects that can modulate the galaxy power spectrum on large scales.The imaging systematic effects often induce wide-angle variations in the density field, and in general, any large-scale variations can translate into an excess signal in the power spectrum (see, e.g., Huterer et al. 2013), that can be misinterpreted as the signature of non-zero local PNG (see, e.g., Thomas et al. 2011).Such spurious variations can be caused by Galactic foregrounds, such as dust extinction and stellar density, or varying imaging conditions, such as astrophysical seeing and survey depth (see, e.g., Ross et al. 2011).The imaging systematic issues have made it challenging to accurately measure  NL , as demonstrated in previous efforts to constrain it using the large-scale clustering of galaxies and quasars (see, e.g., Ross et al. 2013;Pullen & Hirata 2013;Ho et al. 2015), and it is anticipated that they will be particularly problematic for wide-area galaxy surveys that observe regions of the night sky closer to the Galactic plane and that seek to incorporate more lenient selection criteria to accommodate fainter galaxies (see, e.g, Kitanidis et al. 2020).
The primary objective of this paper is to utilize the scaledependent bias signature in the angular power spectrum of galaxies selected from DESI imaging data to constrain the value of  NL .With an emphasis on a careful treatment of imaging systematic effects, we aim to lay the groundwork for subsequent studies of local PNG with DESI spectroscopy.To prepare our sample for measuring such a subtle signal, we employ linear multivariate regression and artificial neural networks to mitigate spurious density fluctuations and ameliorate the excess clustering power caused by imaging systematics.We thoroughly investigate potential sources of systematic error, including survey depth, astronomical seeing, photometric calibration, Galactic extinction, and local stellar density.Our methods and results are validated against simulations, with and without imaging systematics.
This paper is structured as follows.Section 2 describes the galaxy sample from DESI imaging and lognormal simulations with, or without, PNG and synthetic systematic effects.Section 3 outlines the theoretical framework for modelling the angular power spectrum, strategies for handling various observational and theoretical systematic effects, and statistical techniques for measuring the significance of remaining systematics in our sample after mitigation.Our results are presented in Section 4, and Section 5 summarizes our conclusions and directions for future work.

DATA
Luminous red galaxies (LRGs) are massive galaxies that populate massive haloes, lack active star formation, and are highly biased tracers of the dark matter gravitational field (Postman & Geller 1984;Kauffmann et al. 2004).A distinct break around 4000 Å in the LRG spectrum is often utilized to determine their redshifts accurately.LRGs are widely targeted in previous galaxy redshift surveys (see, e.g., Eisenstein et al. 2001;Prakash et al. 2016), and their clustering and redshift properties are well studied (see, e.g., Ross et al. 2020;Gil-Marín et al. 2020;Bautista et al. 2021;Chapman et al. 2022).
DESI is designed to collect spectra of millions of LRGs covering the redshift range 0.2 <  < 1.35.DESI selects its targets for spectroscopy from the DESI Legacy Imaging Surveys, which consist of three ground-based surveys that provide photometry of the sky in the optical , , and  bands.These surveys include the Mayall -band Legacy Survey using the Mayall telescope at Kitt Peak (MzLS; Dey et al. 2018), the Beĳing-Arizona Sky Survey using the Bok telescope at Kitt Peak (BASS; Zou et al. 2017), and the Dark Energy Camera Legacy Survey on the Blanco 4m telescope (DECaLS; Flaugher et al. 2015).As shown in Figure 2, the BASS and MzLS programmes observed the same footprint in the North Galactic Cap (NGC) while the DECaLS programme observed both caps around the galactic plane; the BASS+MzLS footprint is separated from the DECaLS NGC at DEC > 32.375 degrees, although there is an overlap between the two regions for calibration purposes (Dey et al. 2018).Additionally, the DECaLS programme integrates observations executed from the Blanco instrument under the Dark Energy Survey (DES Collaboration et al. 2016), which cover about 1130 deg 2 of the South Galactic Cap (SGC) footprint.The DESI imaging catalogues also integrate the 3.4 (W1) and 4.6  (W2) infrared photometry from the Wide-Field Infrared Explorer (WISE; Wright et al. 2010;Meisner et al. 2018).Collaboration et al. 2023).The redshift evolution of the linear bias is supported by HOD fits to the angular clustering of the DESI LRG targets (Zhou et al. 2021), where  () represents the growth factor.

DESI imaging LRGs
Our sample of LRGs is drawn from the DESI Legacy Imaging Surveys Data Release 9 (DR9; Dey et al. 2018) using the colormagnitude selection criteria designed for the DESI 1% survey (DESI Collaboration et al. 2023), described as the Survey Validation 3 (SV3) selection in more detail in Zhou et al. (2022).The colormagnitude selection cuts are defined in the , ,  bands in the optical and 1 band in the infrared, as summarized in Table 1.The selection cuts vary for each imaging survey, but they are designed to achieve a nearly consistent density of approximately 800 galaxies per square degree across a total area of roughly 14, 000 square degrees.Table 2 summarizes the mean galaxy density and area for each region.This is accomplished despite variations in survey efficiency and photometric calibration between DECaLS and BASS+MzLS.The implementation of these selection cuts in the DESI data processing pipeline is explained in Myers et al. (2022).The redshift distribution of our galaxy sample are inferred respectively from DESI spectroscopy during the Survey Validation phase (DESI Collaboration et al. 2023), and is shown via the solid curve in Figure 1.Zhou et al. (2021) analyzed the DESI LRG targets and found that the redshift evolution of the linear bias for these targets is consistent with a constant clustering amplitude and varies via 1/ (), where  () is the growth factor (as illustrated by the dashed red line in Figure 1).
The LRG sample is masked rigorously for foreground bright stars, bright galaxies, and clusters of galaxies1 to further reduce stellar contamination (Zhou et al. 2022).Then, the sample is binned into HEALPix (Gorski et al. 2005) pixels at nside = 256, corresponding to pixels of about 0.25 degrees on a side, to construct the 2D density map (as shown in the top panel of Figure 2).The LRG density is corrected for the pixel incompleteness and lost areas using a catalogue of random points, hereafter referred to as randoms, uniformly scattered over the footprint with the same cuts and masks applied.Moreover, the density of galaxies is matched to the randoms separately for each of the three data sections (BASS+MzLS, DE-CaLS North / South) so the mean density differences are mitigated (see Table 2).The DESI LRG targets are selected brighter than the imaging survey depth limits, e.g.,  = 24.4, = 23.8, and  = 22.9 for the median 5 detection in AB mag in the DECaLS North region (Table 2); and thus the LRG density map does not exhibit severe spurious fluctuations.

Imaging systematic maps
The effects of observational systematics in the DESI targets have been studied in great detail (see, e.g., Kitanidis et al. 2020;Zhou et al. 2021;Chaussidon et al. 2022).Zhou et al. (2022) has previously identified nine astrophysical properties as potential sources of imaging systematic errors in the DESI LRG targets.These imaging properties are mapped into HEALPix of nside= 256.As illustrated by the 3×3 grid in the bottom panel of Figure 2, the maps include local stellar density constructed from point-like sources with a G-band magnitude in the range 12 ≤  < 17 from the Gaia DR2 (see, Gaia Collaboration et al. 2018;Myers et al. 2022) Schlegel et al. (1998); survey depth (galaxy depth in , , and  and PSF depth in W1) and astronomical seeing (i.e., point spread function, or psfsize) in , , and .The depth maps have been corrected for extinction using the coefficients adapted from Schlafly & Finkbeiner (2011).Table 2 summarizes the median values for the imaging properties in each region.In addition to these nine maps, we consider two external maps for the neutral hydrogen column density (HI) from HI4PI Collaboration et al. (2016) and photometric calibration in the z-band (CALIBZ) from DESI Collaboration et al. (2023) to further test the robustness of our analysis against unknown systematics.
The fluctuations in each imaging map are unique and tend to be correlated with the LRG density map.For instance, large-scale LRG density fluctuations could be caused by stellar density, extinction, or survey depth; while small scale-fluctuations could be caused by psfsize variations.Some regions of the DR9 footprint are removed from our analysis to avoid potential photometric calibration issues.These regions are either disconnected from the main footprint (e.g., the islands in the NGC with DEC < −10) or calibrated using different catalogues of standard stars (e.g., DEC < −30 in the SGC).The potential impact of not imposing these declination cuts on the LRG sample and our  NL constraints is explored in Section 4.
We employ the Pearson correlation coefficient to characterize the correlation between the galaxy density and imaging properties, which for two random variables  and  is given by, where x and ȳ represent the mean estimates of the random variables.Figure 3 shows the Pearson correlation coefficient between the DESI LRG target density map and the imaging systematics maps for the three imaging regions (DECaLS North, DECaLS South, and BASS+MzLS) in the top panel.The horizontal curves represent the 95% confidence regions for no correlation and are constructed by cross-correlating 100 synthetic lognormal density fields, generated with  NL = 0, and the imaging systematic maps.Consistent among the different regions, there are statistically significant correlations between the LRG density and depth, extinction, and stellar density.
There are less significant correlations between the LRG density and the 1-band depth and psfsize.The signs of the correlations imply that there are more targets where extinction is high, and less targets where depth is high.Another interpretation might be that more contaminants are targeted where depth is shallow.Figure 3 (bottom  (Zhou et al. 2022).Magnitudes are corrected for Galactic extinction.The z-band fiber magnitude,  fiber , corresponds to the expected flux within a DESI fiber.

Footprint Criterion Description
fiber < 21.  panel) shows the correlation matrix among the imaging systematic maps for the entire DESI footprint.Significant inner correlations exist among the imaging systematic maps themselves, especially between local stellar density and Galactic extinction; also, the band and -band survey properties are more correlated with each other than with the -band counterpart.Additionally, we compute the Spearman correlation coefficients between the LRG density and imaging systematic maps to assess whether or not the correlations are impacted by outliers in the imaging data, but find no substantial differences from Pearson.

Treatment of imaging systematics
There are several approaches for handling imaging systematic errors, broadly classified into data-driven and simulation-based modeling approaches (see e.g.Ross et al. 2011;Ross et al. 2012Ross et al. , 2017;;Ho et where  0 is a global offset, and a • x  represents the inner product between the parameters, a, and the values for imaging systematics in pixel , x  .The Softplus functional form for   is adapted to force the predicted galaxy counts to be positive (Dugas et al. 2001).Then, Markov Chain Monte Carlo (MCMC) search is performed using the emcee package (Foreman-Mackey et al. 2013) to explore the parameter space by minimizing the negative Poisson log-likelihood between the actual and predicted number counts of galaxies.Spatial coordinates are not included in x  to help avoid overcorrection.As a result, the predicted number counts solely reflect the spurious density fluctuations that arise from varying imaging conditions.The number of pixels is substantially larger than the number of parameters for the linear model, and thus no training-validationtesting split is applied to the data for training the linear model.This aligns with the methodology used for training linear models in previous analyses (see, e.g., Zhou et al. 2022).The predicted galaxy counts are evaluated for each region using the marginalized mean estimates of the parameters, combined with those from other regions to cover the DESI footprint.The linear-based imaging weights are then defined as the inverse of the predicted target density, normalized to a median of unity.
Neural network model: Our neural network-based mitigation approach uses the implementation of fully connected feedforward neural networks from Rezaie et al. (2021).With the neural network approach, a • x  in Equation 3 is replaced with   (x  |a), where   represents the fully connected neural network and a denotes its parameters.The implementation, training, validation, and application of neural networks on galaxy survey data are presented in Rezaie et al. (2021).We briefly summarize the methodology here.
A fully connected feedforward neural network (also called a multi-layer perceptron) is a type of artificial neural network where the neurons are arranged in layers, and each neuron in one layer is connected to every neuron in the next layer.The imaging systematic information flows only in one direction, from input to output.Each neuron applies a non-linear activation function (i.e., transformation) to the weighted sum of its inputs, which are the outputs of the neurons in the previous layer.The output of the last layer is the model prediction for the number counts of galaxies.Our architecture consists of three hidden layers with 20 rectifier activation functions on each layer, and a single neuron in the output layer.The rectifier is defined as max(0, ) to introduce nonlinearities in the neural network (Nair & Hinton 2010).This simple form of nonlinearity is very effective in enabling deep neural networks to learn more complex, non-linear relationships between the input imaging maps and output galaxy counts.
Compared with linear regression, neural networks potentially are more prone to over-fitting, i.e., excellent performance on training data and poor performance on validation (or test) data.Therefore, our analysis uses a training-validation-testing split to avoid over-fitting and ensure that the neural network is well-optimized.Specifically, 60% of the LRG data is used for training, 20% is used for validation, and 20% is used for testing.The split is performed randomly aside from the locations of the pixels.We also test a geometrical split in which neighboring pixels belong to the same set of training, testing, or validation, but no significant performance difference is observed.
The neural networks are trained for up to 70 training epochs with the gradient descent Adam optimizer (Loshchilov & Hutter 2017), which iteratively updates the neural network parameters following the gradient of the negative Poisson log-likelihood.The step size of the parameter updates is controlled via the learning rate hyper-parameter, which is initialized with a grid search and is designed to dynamically vary between two boundary values of 0.001 and 0.1 to avoid local minima (see again, Loshchilov & Hutter 2016).At each training epoch, the neural network model is applied to the validation set, and ultimately the model with the best performance on validation is identified and applied to the test set.The neural network models are tested on the entirety of the LRG sample with the technique of permuting the choice of the training, validation, or testing sets (Arlot & Celisse 2010).With the crossvalidation technique, the model predictions from the different test sets are aggregated together to form the predicted target density map into the DESI footprint.To reduce the error in the predicted number counts, we train an ensemble of 20 neural network models and average over the predictions.The imaging weights are then defined as the inverse of the predicted target density, normalized to a median of unity.

Synthetic lognormal density fields
Density fluctuations of galaxies on large scales can be approximated with lognormal distributions (Coles & Jones 1991;Clerkin et al. 2017).Unlike N-body simulations, simulating lognormal density fields is not computationally intensive, and allows quick and robust validation of data analysis pipelines.Lognormal simulations are therefore considered efficient for our study since the signature of local PNG appears on large-scales and small-scale clustering is not used in our analysis.The package FLASK (Full-sky Lognormal Astro-fields Simulation Kit; Xavier et al. 2016) is employed to generate ensembles of synthetic lognormal density maps that mimic the bias, redshift, and angular distributions of the DESI LRG targets, as illustrated in Figure 1 and 2. Two universes with  NL = 0 and 76.9 are considered.A set of 1000 realizations is produced for every  NL .The mocks are designed to match the clustering signal of the DESI LRG targets on scales insensitive to  NL .The analysis adapts the fiducial BOSS cosmology (BOSS Collaboration et al. 2017) which assumes a flat ΛCDM universe, including one massive neutrino with   = 0.06 eV, Hubble constant ℎ = 0.68, matter density Ω  = 0.31, baryon density Ω  = 0.05, and spectral index   = 0.967.The amplitude of the matter density fluctuations on a scale of 8ℎ −1 Mpc is set as  8 = 0.8225.The same fiducial cosmology is used throughout this paper unless specified otherwise.Our robustness tests show that the none of the cosmological parameters can produce a  NL -like signatures, and therefore, our analysis is not sensitive to the choice of fiducial cosmology.

Contaminated mocks
We employ the linear multivariate model (Equation 3) to introduce synthetic spurious fluctuations in the lognormal density fields, and validate our imaging systematic mitigation methods.The motivation for choosing a linear contamination model is to assess how much of the clustering signal can be removed by applying more flexible models, based on neural networks, for correcting less severe imaging systematic effects.The imaging systematic maps considered for the contamination model are extinction, depth in z, and psfsize in r.As shown in the Pearson correlation (Figure 3) and will be discussed later in Section 3.4, the DESI LRG targets correlate strongly with these three maps.We fit for the parameters of the linear models with the MCMC process, executed separately on each imaging survey (BASS+MzLS, DECaLS North, and DECaLS South).Then, the imaging selection function for contaminating each simulation is uniquely determined by randomly drawing from the parameter space probed by MCMC, and then the results from each imaging survey are combined to form the DESI footprint.The clean density is then multiplied by the contamination model to induce systematics.The same contamination model is used for both the  NL = 0 and 76.9 simulations.
Similar to the imaging systematic treatment analysis for the DESI LRG targets, the neural network methods with various combinations of the imaging systematic maps are applied to each simulation, with and without PNG, and with and without systematics, to derive the imaging weights.Section 3 presents how the simulation results are incorporated to calibrate  NL biases due to overcorrection.We briefly summarize two statistical tests based on the mean galaxy density contrast and the cross power spectrum between the galaxy density and the imaging systematic maps to assess the quality of the data and the significance of the remaining systematic effects (see, also, Rezaie et al. 2021).We calculate these statistics and compare the values to those measured from the clean mocks before looking at the auto power spectrum of the DESI LRG targets.

ANALYSIS TECHNIQUES
We address imaging systematics in DESI data by performing a separate treatment for each imaging region (e.g., DECaLS North) within the DESI footprint to reduce the impact of systematic effects specific to that region.Once the imaging systematic weights are obtained for each imaging region separately, we combine the data from all regions to compute the power spectrum for the entire DESI footprint to increase the overall statistical power and enable more robust measurements of  NL .We then conduct robustness tests on the combined data to assess the significance of any remaining systematic effects.

Power spectrum estimator
We first construct the density contrast field from the LRG density, , where the mean galaxy density  is estimated from the entire LRG sample.As a robustness test, we also analyze the power spectrum from each imaging region individually, in which  is calculated separately for each region.Then, we use the pseudo angular power spectrum estimator (Hivon et al. 2002), where the coefficients  ℓ are obtained by decomposing   into spherical harmonics,  ℓ , where  represents the survey window that is described by the number of randoms normalized to the expected value.
We use the implementation of anafast from the HEALPix package (Gorski et al. 2005) to do fast harmonic transforms (Equation 6) and estimate the pseudo angular power spectrum of the LRG targets and the cross power spectrum between the LRG targets and the imaging systematic maps.

Modelling
The estimator in Equation 5 yields a biased power spectrum when the survey sky coverage is incomplete.Specifically, the survey mask causes correlations between different harmonic modes (Beutler et al. 2014;Wilson et al. 2017), and the measured clustering power is smoothed on scales near the survey size.An additional potential cause of systematic error arises from the fact that the mean galaxy density used to construct the density contrast field (Equation 4) is estimated from the available data, rather than being known a priori.This introduces what is known as an integral constraint effect, which can cause the power spectrum on modes near the size of the survey to be artificially suppressed, effectively pushing it towards zero (Peacock & Nicholson 1991;De Mattia & Ruhlmann-Kleider 2019).Since  NL is highly sensitive to the clustering power on these scales, it is crucial to account for these systematic effects in the model galaxy power spectrum to obtain unbiased  NL constraints (see, also, Riquelme et al. 2022), which we describe below.
The other theoretical systematic issues are however subdominant in the angular power spectrum.For instance, relativistic effects generate PNG-like scale-dependent signatures on large scales, which interfere with measuring  NL with the scale-dependent bias effect using higher order multipoles of the 3D power spectrum (Wang et al. 2020).Similarly, matter density fluctuations with wavelengths larger than survey size, known as super-sample modes, modulate the galaxy 3D power spectrum (Castorina & Moradinezhad Dizgah 2020).In a similar way, the peculiar motion of the observer can mimic a PNG-like scale-dependent signature through aberration, magnification and the Kaiser-Rocket effect, i.e., a systematic dipolar apparent blue-shifting in the direction of the observer's peculiar motion (Bahr-Kalus et al. 2021).

Angular power spectrum
The relationship between the linear matter power spectrum () and the projected angular power spectrum of galaxies is expressed by the following equation: where  shot is a scale-independent shot noise term.The projection kernel ) includes redshift space distortions and magnification bias, and determines the contribution of each wavenumber  to the galaxy power spectrum on mode ℓ.For more details on this estimator, refer to Padmanabhan et al. (2007).The non-linearities in the matter power spectrum are negligible for the scales of interest (see, e.g., Ho et al. 2015).For ℓ = 40, Δ ℓ () peaks at  ∼ 0.02 ℎMpc −1 , which is above the non-linear regime.The FFTLog algorithm and its extension2 as implemented in Fang et al. (2020) are employed to calculate the integrals for the projection kernel Δ ℓ (), which includes the  th order spherical Bessel functions,  ℓ (), and its second derivatives, where  is the linear bias (dashed curve in Figure 1),  represents the linear growth factor normalized as  ( = 0) = 1,  () is the growth rate, and / is the redshift distribution of galaxies normalized to unity and described in terms of comoving distance3 (solid curve in Figure 1).The magnification bias window function where Ω  is the matter density,  0 is the Hubble constant4 ,  is the speed of light, and  represents the slope of the number count function, a metric quantifying the response of the number density of galaxies to achromatic changes in brightness (Loverde et al. 2008).
The estimation of  involves shifting all magnitudes by an infinitesimal amount and re-running the color-magnitude selection.Zhou et al. (2023) developed a strategy to estimate  for a fiber fluxselected sample like the DESI LRG targets, for which the impact of magnification on fiber flux depends on the shape parameters for each morphology type.Following the same strategy, the parameter  is re-calculated for our selection of the DESI LRGs (DESI SV3)5 :  = 0.951 ± 0.011 for BASS+MzLS,  = 0.943 ± 0.007 for DE-CaLS North+DECaLS South, and  = 0.945 ± 0.006 for DESI.For consistency, we fix  to the above central values in our analysis.The PNG-induced scale-dependent shift is given by (Slosar et al. 2008) where  () is the transfer function, and (∞)/(0) ∼ 1.3 with () ≡ (1 + ) () is the growth suppression due to non-zero Λ because of our normalization of  (see, e.g., Reid et al. 2010;Mueller et al. 2019).We assume the universality relation which directly relates   to  via   = 2  ( − ) with   = 1.686 representing the critical density for spherical collapse (Fillmore & Goldreich 1984).We fix  = 1 in our analysis and marginalize over b (see, also, Slosar et al. 2008;Reid et al. 2010;Ross et al. 2013).

Survey geometry and integral constraint
We employ a technique similar to the one proposed by Chon et al. (2004) to account for the impact of the survey geometry on the theoretical power spectrum.The ensemble average for the partial sky power spectrum is related to that of the full sky power spectrum via a mode-mode coupling matrix, M ℓℓ ′ , We convert this convolution in the spherical harmonic space into a multiplication in the correlation function space.Specifically, we first transform the theory power spectrum (Equation 7) to the correlation function, ωmodel .Then, we estimate the survey mask correlation function, ωwindow , and obtain the pseudo-power spectrum, Figure 4 illustrates the survey mask correlation function ωwindow for various masks representing the DESI footprint and its imaging sub-regions.Appendix A2 shows the validation of our method by comparing it to an alternative approach that computes the modemode coupling matrix M ℓℓ ′ and performs the convolution (Equation 13) directly in ℓ-space.We notice as the  NL deviates from zero, our approach introduces a noisy feature in the model, qualitatively in an unbiased manner (Δ  NL < 1.1).Figure B2 indeed demonstrates that our approach can recover the truth  NL in spite of the noisy feature.The integral constraint is another systematic effect which is induced since the mean galaxy density is estimated from the observed galaxy density, and therefore is biased by the limited sky coverage (Peacock & Nicholson 1991).To account for the integral constraint, the survey mask power spectrum is used to introduce a scale-dependent correction factor that needs to be subtracted from the power spectrum as, where Cwindow is the survey mask power spectrum, i.e., the spherical harmonic transform of ωwindow .
The lognormal simulations are used to validate the survey window and integral constraint correction.Figure 5 shows the mean power spectrum of the  NL = 0 simulations (dashed) and the bestfitting theory prediction before and after accounting for the survey mask and integral constraint.The simulations are neither contaminated nor mitigated.The light and dark shades represent the 68% estimated error on the mean and one single realization, respectively.The DESI mask, which covers around 40% of the sky, is applied to the simulations.We find that the survey window effect modulates the clustering power on ℓ < 200 and the integral constraint alters the clustering power on ℓ < 6.

Parameter estimation
Our parameter inference uses standard MCMC sampling.A constant clustering amplitude is assumed to determine the redshift evolution of the linear bias of the DESI LRG targets, () = / (), which is supported by the HOD fits to the angular power spectrum (Zhou et al. 2021).In MCMC, we allow  NL ,  shot , and  to vary, while all other cosmological parameters are fixed at the fiducial values (see §2.2).The galaxy power spectrum is divided into a discrete set of bandpower bins with Δℓ = 2 between ℓ = 2 and 20 and Δℓ = 10 from ℓ = 20 to 300.Each clustering mode is weighted by 2ℓ + 1 when averaging over the modes in each bin.
The expected large-scale power is highly sensitive to the value of  NL such that the the amplitude of the covariance for  ℓ is influenced by the true value of  NL , see also Ross et al. (2013) for a discussion.As illustrated in the top row of Figure 6, we find that the distribution of the power spectrum at the lowest bin, 2 ≤ ℓ < 4, is highly asymmetric and its standard deviation varies significantly from the simulations with  NL = 0 to 76.9.We can make the covariance matrix less sensitive to  NL by taking the log transformation of the power spectrum, log  ℓ .As shown in the bottom panels in Figure 6, the log transformation reduces the asymmetry and the difference in the standard deviations between the  NL = 0 and 76.9 simulations.Therefore, we minimize the negative log likelihood defined as, where Θ represents a container for the parameters  NL , , and  shot ; C (Θ) is the (binned) expected pseudo-power spectrum; C is the (binned) measured pseudo-power spectrum; and C is the covariance on log C constructed from the  NL = 0 log-normal simulations.Log-normal simulations have been commonly used and validated to estimate the covariance matrices for galaxy density fields, and non-linear effects are subdominant on the scales of interest to  NL (see, e.g., Clerkin et al. 2017;Friedrich et al. 2021).We also test for the robustness of our results against an alternative covariance constructed from the  NL = 76.9 mocks.Flat priors are implemented for all parameters:  NL ∈ [−1000, 1000],  shot ∈ [−0.001, 0.001], and  ∈ [0, 5].

Characterization of remaining systematics
One potential problem that can arise in the data-driven mitigation approach is over-correction, which occurs when the corrections applied to the data remove the clustering signal and induce additional biases in the inferred parameter of interest.The neural network approach is more prone to this issue compared to the linear approach due to its increased degrees of freedom.As illustrated in the bottom panel of Figure 3, the significant correlations among the imaging systematic maps may pose additional challenges for modeling the spurious fluctuations in the galaxy density field.Specifically, using highly correlated imaging systematic maps increases the statistical noise in the imaging weights, which elevates the potential for over subtracting the clustering power.These over-correction effects are estimated to have a negligible impact on baryon acoustic oscillations (Merz et al. 2021); however, they can significantly modulate the galaxy power spectrum on large scales, and thus lead to biased  NL constraints (Rezaie et al. 2021;Mueller et al. 2022).Although not explored thoroughly, the over-correction issues could limit the detectability of primordial features in the galaxy power spectrum and that of parity violations in higher order clustering statistics (Beutler et al. 2019;Cahn et al. 2021;Philcox 2022).Therefore, it is crucial to develop, implement, and apply techniques to minimize and control over-correction in the hope of ensuring that the constraints are as accurate and reliable as possible; one such approach is to reduce the dimensionality of the problem.Our goal is to reduce the correlations between the DESI LRG target density and the imaging systematic maps while controlling the over-correction effect.In the following, we describe how we approach this objective, by employing a series of simulations along with the residual systematics that we construct based on the cross power spectrum between the LRG density and imaging maps, and the mean LRG density as a function of imaging.We test different sets of the imaging systematic maps to identify the optimal set of the feature maps: (i) Two maps: Extinction, depth in z.
(ii) Three maps: Extinction, depth in z, psfsize in r.
(iii) Four maps: Extinction, depth in z, psfsize in r, stellar density.
(vi) Eleven maps: same as Nine maps but with two additional maps; Extinction, depth in 1, psfsize in , stellar density, neutral hydrogen density, and photometric calibration in z.
It is imperative to note that these sets are selected prior to examining the auto power spectrum of the LRG sample and unblinding the  NL constraints, and that the auto power spectrum and  NL measurements are unblinded only after our mitigation methods passed our rigorous tests for residual systematics.As detailed in the following, we discover that these tests tend to depend on the  NL value which is used in the mocks for the covariance matrix estimation.

Normalized cross power spectrum
We characterize the cross correlations between the galaxy density and imaging systematic maps by where C  ,ℓ represents the the square of the cross power spectrum between the galaxy density and  th imaging map,   , divided by the auto power spectrum of   : With this normalization, C  ,ℓ estimates the contribution of systematics at every multipole up to the linear order to the galaxy power spectrum.Then, the  2 value for the cross power spectra is calculated via, where the covariance matrix C  =< C,ℓ C,ℓ ′ > is constructed from the lognormal mocks.We consider both sets with  NL = 0 and 76.9, and for each mitigated case, the covariance is from the mocks that have received the same treatment.These  2 values are measured for every clean mock realization with the leave-one-out technique and compared to the values observed in the DESI LRG targets with various imaging systematic corrections.Specifically, we use 999 realizations to estimate a covariance matrix and then apply the covariance matrix from the 999 realizations to measure the  2 for the one remaining realization.This process is repeated for all 1000 realizations to construct a histogram for  2 .We only include the bandpower bins from ℓ = 2 to 20 with Δℓ = 2, which results in a total of 81 bins, and test for the robustness with higher ℓ modes in A1.
Figure 7 shows C,ℓ from the DESI LRG targets before and after applying various corrections for imaging systematics.The dark and light shades show the 97.5 th percentile from the  NL = 0 and 76.9 mocks, respectively, that have had no mitigation applied to them.Without imaging weights (No Weight), the DESI LRG targets have the highest cross-correlations against extinction, stellar density, and depth in z.There are less significant correlations against depth in the g and r bands, and psfsize in the z band, which could be driven because of the inner correlations between the imaging systematic maps.First, we consider cleaning the DESI LRG targets with the linear model using two maps (extinction and depth in z) as identified from the Pearson correlation.Linear two maps is the least aggressive treatment method in terms of both the model flexibility and the number of input maps.With linear two maps, most of the cross correlation signals are reduced below statistical uncertainties, especially against extinction, stellar density, and depth.However, the cross correlations against psfsize in the r and z bands increase slightly on 6 < ℓ < 20 and 6 < ℓ < 14, respectively.This might be indicative of residual trends against psfsize.
The linear three maps approach alleviates the cross correlation against psfsize in r, and it yields similar results to those obtained from linear nine maps, which indicates most of the contaminations can be attributed to these three maps.Therefore, we identify extinction, depth in z, and psfsize in r (three maps) as the primary sources of systematic effects in the DESI LRG targets.Then, we adapt neural network three maps to model non-linear systematic effects.Compared with the linear three maps method, we find that non-linear three maps can reduce the cross correlations against both the r and z-band psfsize maps, which shows the benefit of using a non-linear approach.To further examine the robustness of our cleaning methods, we also show the cross correlations after cleaning the DESI LRG targets using nine imaging property maps (non-linear nine maps).We do not find any significant residuals against the two extra maps for the neutral hydrogen density and photometric calibration in the z band.

Mean galaxy density contrast
We calculate the histogram of the mean galaxy density contrast relative to the  th imaging property,   : where  is the global mean galaxy density,   is the survey window in pixel , and the summations over  are evaluated from the pixels in every bin of   .We compute the histograms against all nine imaging properties (see Figure 2).We use a set of eight equal-width bins for every imaging map, which results in a total of 72 bins.Then, we construct the total mean density contract as, and the total residual error as, where the covariance matrix C  =<     > is constructed from the lognormal mocks, in a consistent manner similar to the normalized cross power spectrum.Figure 8 shows the mean density contrast against the imaging properties for the DESI LRG targets.The dark and light shades represent the 1 level fluctuations observed in 1000 lognormal density fields respectively with  NL = 0 and 76.9 before mitigation.The DESI LRG targets before treatment (No Weight) exhibits a strong trend around 10% against the z-band depth which is consistent with the cross power spectrum.Additionally, there are significant spurious trends against extinction and stellar density at about 5 − 6%.The linear approach is able to mitigate most of the systematic fluctuations with only extinction and depth in the z-band as input; however, a new trend appears against the r-band psfsize map with the linear two maps approach, which is indicative of the psfsize-related systematics in the DESI LRG targets.This finding is in agreement with that from the cross power spectrum test.With The square of the cross power spectra between the DESI LRG targets and imaging systematic maps normalized by the auto power spectrum of the imaging systematic maps; see equation 18.The systematic maps considered are Galactic extinction (EBV), stellar density (nStar), depth in grzw1 (depth  1 ), and seeing in grz (psfsize   ).The dark green curves display the cross spectra before imaging systematic correction (No Weight).The yellow, brown, and light green curves represent the results after applying the imaging weights from the linear models trained with two maps, three maps, and nine maps.The orange and purple curves display the results after applying the imaging weights from the non-linear models trained with three maps and nine maps.The dark and light shades represent the 97.5 percentile from cross correlating the imaging systematic maps and the  NL = 0 and 76.9 lognormal density fields, respectively, without mitigation.
linear three maps, we still observe around 2% residual spurious fluctuations in the low end of depth in z and around 1% in the high end of psfsize in z, which implies the presence of non-linear systematic effects.We find that the imaging weights from the nonlinear model trained with the three identified maps (or four maps including the stellar density) are capable of reducing the fluctuations below 2%.Even with the non-linear three maps, we have about 1% remaining systematic fluctuations against the z-band psfsize.The spurious trends are diminished especially when we adapt non-linear nine maps, especially against the low end of depth in g and r and against the high end of psfsize in z.

Residual error 𝜒 2
We use the  2 statistics to quantitatively assess how significant these mean density and cross power spectrum fluctuations are in comparison to the clean mocks.Figure 9 presents  2 histograms for the mean density contrast (left) and the normalized cross spectrum (right) statistics obtained from the lognormal mocks with different  NL values before and after applying mitigation methods.The mocks with  NL = 0 are shown with the solid curves while the other set with  NL = 76.9 are represented with the dashed curves.The use of the self-consistent covariance matrix (with respect to  NL or mitigation method) results in similar distributions, and therefore the mock histograms are employed as reference to evaluate the significance of residual systematics in the DESI LRG targets.We continue to use the self-consistent covariance, but consider both the  NL = 0 and 76.9 covariance.The DESI LRG target  2 values are compared via the vertical lines and summarized in Table 3.The solid and dashed vertical lines represent the values computed using the covariances based on the  NL = 0 and 76.9 mocks, respectively.Regardless of the covariance used in the  2 calculations, we find that the case without treatment (No Weight) exhibits serious contamination.For instance with the  NL = 0 covariance, the mean density and cross power errors are respectively  2 = 679.8and 20014.8(both with -value < 0.001).These  2 values are significantly high given that the degree  ), and seeing in grz (psfsize   ).The black curves display the results before imaging systematic correction.The red, blue and orange curves represent the relationships after applying the imaging weights from the linear models trained with two maps, three maps, and eight maps, respectively.The green and pink curves display the results after applying the imaging weights from the non-linear models trained with three maps and four maps.The dark and light shades represent the 68% dispersion of 1000 lognormal mocks with  NL = 0 and 76.9, respectively. of freedom is 72 for the mean density and 81 for the cross power spectrum.After cleaning, the  2 values are decreased dramatically for both the mean density and normalized cross spectrum tests.The small impact on  2 from including stellar density suggests that the stellar density trend can be explained by extinction due to the correlation between these properties, such that in regions with high stellar density, there is likely to be a higher concentration of dust, which can cause greater extinction of light.However, neither nonlinear three maps nor non-linear four maps can reduce the mean density  2 enough to be consistent with the mocks, indicating some significant residual error with -values less than 0.005.
The tests conducted here demonstrate the effectiveness of various cleaning approaches for the DESI LRG targets without revealing the measured power spectrum or  NL constraints.Overall, we observe that the non-linear method with the set of nine imaging property maps, successfully passes the mean density test irrespective of the covariance, as indicated by  2 = 71.9(-value = > 0.48) for the  NL = 0 covariance.On the other hand, the non-linear three maps and non-linear four maps methods both fail to sufficiently mitigate systematics in the mean density test, as evidenced by low -values.
Our work shows that it is essential to maintain a consistent covari-ance matrix, involving the same mitigation and ensuring consistency in  NL within the covariance.The sensitivity of the mean density  2 to the  NL assumption in the covariance is notably lower, with greater reliance on the consistent mitigation method.Conversely, the normalized cross spectrum  2 exhibits a higher dependency on the  NL assumption in the covariance.The mean density diagnostic appears to be a less cosmology-sensitive probe of residual systematics.As a robustness test, we also increase the largest ℓ used in the  2 calculation to ℓ = 100, which corresponds to density fluctuations on angles smaller than 2 degrees.But we find no remaining systematic error from higher harmonic modes (see Appendix A1).The conclusion of these tests is that the non-linear method with the set of nine maps passes our null tests for the remaining systematics, and thus is chosen as the default approach for the treatment of imaging systematic effects.In the following subsection, we show how imaging systematic regressions remove clustering modes, with increasing dependence on the number of maps used, and thus bias the best fitting estimates of  NL .Then, we present how we calibrate for the over-correction for our default mitigation method.The values observed in the DESI LRG targets after the non-linear treatments are represented via vertical lines using the  NL = 0 (solid) or 76.9 (dashed) covariance, and the histograms are constructed from 1000 realizations of clean lognormal mocks with  NL = 0 (solid) and 76.9 (dashed).
Table 3. Mean galaxy density contrast  2 and normalized cross power spectrum  2 from the DESI LRG targets and -values that are inferred from the comparison to the  NL = 0 and 76.9 clean mocks that have received the same mitigation.For the case of No Weight, we use the clean mocks without mitigation.

Calibration of over-correction
The template-based mitigation of imaging systematics removes some of the true clustering signal, and mitigating with more maps should remove more modes and thus bias both the  NL estimation and its associated uncertainty.We calibrate the over-correction effect using the mocks presented in §2.Having two sets of mocks with low and high power at large scales (low ℓ) offers a key advantage: it provides a model for mapping the entire posterior distribution, which enables sus to understand how the constraints on  NL degrade as the magnitude of the imaging systematic correction increases.We apply the neural network model to both the  NL = 0 and 76.9 simulations, with and without imaging systematics, using various sets of imaging systematic maps.Specifically, we consider non-linear three maps, non-linear four maps, and non-linear nine maps.Then, we measure the power spectra from the mocks.We fit both the mean power spectrum and each individual power spectrum from the mocks.Appendix B2 outlines the impact of the non-linear methods on the mock power spectra, and here we summarize relevant details for the calibration of over-correction.
Fihgure 10 displays a comparison between the best-fitting estimates of  NL before and after mitigation for the clean mocks.The best-fitting estimates from the mean of the mocks are represented by the solid curves, and the individual spectra results are displayed as the scatter points.The results from fitting the mean power spectrum of the contaminated mocks are also shown via the dashed curves.We find nearly identical results for the biases caused by mitigation, whether or not the mocks have any contamination, which can be seen by observing the solid and dashed curves displayed on Figure 10 (see, also, Figure B4, for a comparison of the mean power spectrum).For clarity, the best-fitting estimates for the individual contaminated data are not shown in the figure .As summarized in Table B2, we observe notable shifts in the best-fitting estimates of  NL obtained from the mean power spectrum of the mocks.Specifically, for the  NL = 0 mocks, we obtain Δ  NL = −12 for non-linear three maps, −20 for non-linear four maps, and −27 for non-linear nine maps.Larger shifts are evident for  NL = 76.9:Δ  NL = −23 for non-linear three maps, −39 for non-linear four maps, and −72 for non-linear nine maps.These factor imply that the effect of systematic mitigation on the inferred  NL depends on the true value of  NL .
To calibrate our methods, we fit a linear curve to the  NL estimates from the mean power spectrum of the mocks,  NL,no mitigation,clean =  1  NL,mitigated +  2 .The  1 and  2 co-  efficients for non-linear three, four, and nine maps are summarized in Table 4.These coefficients represent the impact of the cleaning methods on the likelihood.The uncertainty in  NL after mitigation increases by  1 − 1. Figure 10 also shows that the choice of our cleaning method can have significant implications for the accuracy of the measured  NL , and careful consideration should be given to the selection of the primary imaging systematic maps and the calibration of the cleaning algorithms in order to minimize systematic uncertainties.

RESULTS
We now present our  NL constraints obtained from the power spectrum of the DESI LRG targets.The treatment of the imaging systematic effects is performed on each imaging region (BASS+MzLS, DECaLS North/South) separately.After cleaning, the regions are combined for the measurement of the power spectrum.We unblind the galaxy power spectrum and the  NL values after our cleaning methods are validated and vetted by the cross power spectrum and mean galaxy density diagnostics.As presented in Section 3.4, these tests show that none of the linear methods yields reasonable statistics, and only the nonlinear approach with the nine maps can pass the criteria, which is why we select the nonlinear nine maps as our fiducial method for cleaning systematics.We also conduct additional tests to check the robustness of our constraints against various assumptions, such as analyzing each region separately, applying cuts on imaging conditions, and changing the smallest mode used in fitting for  NL .

DESI imaging LRG sample
We find that the excess clustering signal in the power spectrum of the DESI LRG targets is mitigated after correcting for the imaging systematic effects.Figure 11 shows the measured power spectrum of the DESI LRG targets before and after applying imaging weights and the best-fitting theory curves.The solid grey line and the grey shade represent respectively the mean power spectrum and 1 error, estimated from the  NL = 0 lognormal simulations.The differences between various cleaning methods are significant on large scales (ℓ < 20), but the small scale clustering measurements are consistent.We associate the differences to the over-correction caused by including more maps for the treatment of systematics, which we base upon the validation of the methods on the mocks, or the suppression of excess power from systematics.Comparing non-linear three maps to non-linear four maps, we find that adding stellar density in the non-linear approach (non-linear four maps) further reduces the excess power relative to the mock power spectrum, in particular on modes between 2 ≤ ℓ < 4.However, when calibrated on the lognormal simulations, we find that the over-subtraction due to stellar density is reversed after accounting for over-correction.
Our fiducial approach, non-linear nine maps, yields the lowest (and almost constant) power on large scales among all methods.

Calibrated constraints
All  NL constraints presented here are calibrated for the effect of over-correction using the lognormal simulations.Table 5 describes the best-fitting and marginalized mean estimates of  NL from fitting the power spectrum of the DESI LRG targets before and after cleaning with the non-linear approach given various combinations for the imaging systematic maps.Figure 12 shows the marginalized probability distribution for  NL in the top panel, and the 68% and 95% probability contours for the linear bias parameter and  NL in the bottom panel, from our sample before and after applying various corrections for imaging systematics.Overall, we find the maximum likelihood estimates to be consistent among the various cleaning methods.We obtain 33(21) <  NL < 61(76) at 68%(95%) confidence with  2 = 33.9 for non-linear three maps with 34 degrees of freedom.Accounted for over-correction, we obtain 33(19) <  NL < 62(78) with  2 = 34.4 using nonlinear four maps which includes the additional stellar density map.
With or without stellar density, the confidence intervals are consistent with each other and significantly off from zero PNG; specifically, the probability that  NL is erroneously greater than zero, (  NL > 0) = 99.9 per cent, which we attribute to systematics (see Section 3.4).We also apply a more aggressive systematics treatment that includes regression using the non-linear approach against the set of nine imaging maps we identified, non-linear nine maps, and find that zero  NL is recovered.Specifically, our maximum likelihood value changes to  NL ∼ 34 with  2 = 39.1, and the uncertainty on  NL increases by more than a factor of two, resulting in −10(−39) <  NL < 58(84) at 68%(95%) confidence.This increase is attributed to the aggressive treatment, which removes large-scale clustering information and diminishes the constraining power of the dataset.
Additionally, we explore the sensitivity of the  NL posterior using the non-linear nine maps method while varying the values of  in the range of 0.5 to 1.6 and  in the range of 0.75 to 1.25.Figure 13 illustrates our findings, and Table 6 provides a summary.Regardless of the specific values chosen for  and , we reliably recover  NL = 0 within the 95% confidence interval.The top panel also implies that marginalizing over  can induce projection effects and lead to biased  NL For comparison, we obtain 102(86) <  NL < 140(161) at 68%(95%) confidence with  2 = 45.1 for the no weight case.

Uncalibrated constraints: robustness tests
Figure 14 shows the probability distributions of  NL for various treatments before accounting for the over-correction effect.The method with the largest flexibility and more number of imaging systematic maps is more likely to regress out the clustering signal aggressively and return biased  NL constraints.The non-linear three maps approach returns a best-fitting estimate of  NL = 27 with the 68%(95%) confidence of 17(6) <  NL < 40(53) and  2 = 33.9.With the stellar density map included, non-linear four maps yields a smaller best-fitting estimates of  NL = 14 with the error of 5(−6) <  NL < 26(38).The non-linear nine maps gives an asymmetric posterior with the marginalized mean  NL = −17, and the smallest best-fitting estimate of  NL = −13 with the error of −31(−44) <  NL < −3(9).The disparities in the best-fitting estimates can be linked to over-correction, mirroring the effects observed in the mocks (refer to Figure B5).Consequently, caution is advised when considering the uncalibrated values.Without adjusting for over-correction, non-linear four maps and non-linear nine maps recover zero  NL within 95% and 68% confidence, respectively.However, the non-linear method with three maps tension with  NL = 0 at a confidence level of 99.5 percent.Now we proceed to perform some robustness tests and assess how sensitive the  NL constraints are to the assumptions made in the analysis or the quality cuts applied to the data.For each case, we retrain the cleaning methods and derive new sets of imaging weights.Accordingly, for the cases where a new survey mask is applied to the data, we re-calculate the covariance matrices using the new survey mask to account for the changes in the survey window and integral constraint effects.Calibrating the mitigation biases for all of these experiments is beyond the scope of this work and redundant, as we are only interested in the relative shift in the  NL constraints after changing the assumptions.Therefore, the absolute scaling of the  NL constraints presented here are biased because of the over-correction effect.Table 7 summarizes the uncalibrated  NL constraints from the DESI LRG targets.Our tests are as follows: • Linear methods: Even though the linear methods show remaining systematics (e.g., against depth in z as shown in Figure 8), we obtain identical constraints from linear four maps and linear nine maps, respectively, 20(9) <  NL < 45(58) and 19(9) <  NL < 43(57) at 68%(95%) confidence.For the linear treatment methods, the probability of  NL being greater than zero is erroneously 99.9 per cent.Any attempt to account for the overcorrection would elevate this probability even further.The overestimation of  NL can be attributed to an increase in systematic contamination.
• Imaging regions: We compare how our constraints from fitting the power spectrum of the whole DESI footprint compares to that from the power spectrum of each imaging region individually, namely BASS+MzLS, DECaLS North, and DECaLS South.
Figure 15 shows the 68% and 95% probability contours on  NL and  from each individual region, compared with that from DESI.The cleaning method here is non-linear nine maps, and the covariance matrices are estimated from the  NL = 0 mocks.The bias in DECaLS North is lower than the ones from DECaLS South Table 7.The uncalibrated best-fitting and marginalized mean estimates for  NL from fitting the power spectrum of the DESI LRG targets before and after correcting for systematics.The estimates are not calibrated for over-correction, and thus are subject to mitigation systematics.The number of degrees of freedom is 34 (37 data points -3 parameters) for all cases except the case that combines the data at the likelihood level, 'BASS+MzLS+DECaLS', in which the dof is 104 (3 × 37 − 7).The lowest mode is ℓ = 2 and the covariance matrix is from the  NL = 0 clean mocks (no mitigation) except for the case with '+ Cov' in which the covariance matrix is from the  NL = 76.9 clean mocks (no mitigation).We fix  = 1 for all cases and  = 0.945 for DESI, 0.943 for DECaLS North (and South), and 0.951 for BASS+MzLS.and BASS+MzLS, which might indicate some remaining systematic effects that could not be mitigated with the available imaging systematic maps.This is because given the negative correlation be-tween  and  NL , a larger value of  NL due to excess clustering power needs to be compensated by a smaller value of .Overall, we find that the constraints from analyzing each imaging survey separately are consistent with each other and DESI within 68% confidence.We also consider combining the data at the likelihood level ('BASS+MzLS+DECaLS').In this case the total number of data points is 111 (3 × 37).We allow the bias and shotnoise paramters to vary independently for each sub-region but use a single and common  NL value, which brings the total number of free parameters to 7 and the number of degrees of freedom to 104.We obtain a best-fitting estimate of  NL = −31 with  2 = 114.2and 68% (95%) confidence interval of −41(−53) <  NL < 9(5).Compared with our fiducial analysis which combines the data at the map level, we observe around 13% loss in constraining power.
• Stellar density template (nStar): When not accounting for overcorrection, adding the stellar density map appears to result in significant changes in the  NL constraints, e.g., compare non-linear three maps with non-linear four maps in Table 7.But these changes disappear when we account for the mitigation bias and we find both methods recover the same maximum likelihood estimate for  NL ∼ 46 within 69% confidence, see Table 5, which implies that these changes can be associated with the over-correction issue from the chance correlations between the stellar density map and largescale structure.• Pixel completeness (comp.cut): We discard pixels with fractional completeness less than half to assess the effect of partially complete pixels on  NL .This pixel completeness cut removes 0.6% of the survey area, and no significant changes in the  NL constraints are observed.
• Imaging quality (imag.cut): Pixels with poor photometry are removed from our sample by applying the following cuts on imaging;  [ −] < 0.1,  < 3000, depth  > 23.2, depth  > 22.6, depth  > 22.5, psfsize  < 2.5, psfsize  < 2.5, and psfsize  < 2. Although these cuts remove 8% of the survey mask, there is a neg-ligible impact on the best-fitting estimates of  NL from fitting the DESI power spectrum.However, when each region is fit individually, the BASS+MzLS constraint is more stable than those from DECaLS North and DECaLS South.
• Covariance matrix (cov): We fit the power spectrum of our sample cleaned with non-linear nine maps, but use the covariance matrix constructed from the  NL = 76.9 mocks.With the alternative covariance, a 7% increase in the 68% error on  NL , (  NL ), is observed.We also find that the best-fitting and marginalized mean estimates of  NL increase slightly by Δ  NL = 2. Overall, we find that the differences are not significant in comparison to the statistical precision.
• External maps (CALIBZ+HI): The neural network eleven maps correction includes the additional maps for the neutral column density (HI) and the z-band calibration error (CALIBZ).With this correction, the best-fitting  NL increases from −4 to −2 for DECaLS North and from −43 to −38 for DECaLS South, which might suggest that adding HI and CALIBZ increases the input noise, and thus negatively impacts the performance of the neural network model.This test is not performed on BASS+MzLS due to a lack of coverage from the CALIBZ map.
• Declination mask (no DEC cut): The fiducial mask removes the disconnected islands in DECaLS North and regions with DEC < −30 in DECaLS South, where there is a high likelihood of calibration issues as different standard stars are used for photometric calibrations.We analyze our sample without these cuts, and find that the best-fitting and marginalized  NL mean estimates from DECaLS South shift significantly to higher values of  NL by Δ  NL ∼ 41, which supports the case that there are remaining photometric systematics in the DECaLS South region below DEC = −30.On the other hand, the constraints from DECaLS North do not change significantly, indicating the islands do not induce significant contaminations.
• Scale dependence (varying ℓ min ): We raise the value of the lowest harmonic mode ℓ min used for the likelihood evaluation during MCMC.This is equivalent to utilizing smaller spatial scales in the measurements of the power spectrum.By doing so, we anticipate a reduction in the impact of imaging systematics on  NL inference as lower ℓ modes are more likely to be contaminated.Figure 16 illustrates the power spectra before and after the correction with non-linear nine maps in the top panel.The bottom panel shows the best fitting estimate and 95% error on  NL with non-linear nine maps for the DESI, BASS+MzLS, DECaLS North, and DECaLS South regions.We discover that a slight upward shift in the best fitting estimates of  NL on scales ranging from 10 to 20 for DECaLS North and BASS+MzLS when we utilized a higher ℓ min .This outcome might imply that the imaging systematic maps do not contain enough information to help the cleaning method null out the contaminating signal in the NGC.We also find that the bump is resilient against an alternative correction, in which we apply the neural networks trained on the DECaLS South to the DECaLS North region (see A4). Overall, this result is contrary to what one might predict if a significant systematic-induced spike existed at the very low ℓ, or if we had an extremely large-scale systematic leakage from the ℓ = 1 mode.As a result, it suggests that the underlying issue is more subtle than originally anticipated.

DISCUSSION AND CONCLUSION
We have measured the local PNG parameter  NL using the scaledependent bias in the angular clustering of LRGs selected from the DESI Legacy Imaging Survey DR9.Our sample includes more than 12 million LRG targets covering around 14, 000 square degrees in the redshift range of 0.2 <  < 1.35.We leverage early spectroscopy during DESI Survey Validation (DESI Collaboration et al. 2023) to infer the redshift distribution of our sample (Figure 1).Our power spectrum model accounts for various theoretical and observational effects such as RSD, magnification bias, survey geometry, and integral constraint.Most importantly, we utilize a novel machine learning-method to mitigate the effect of imaging systematics and reduce excess clustering power on large scales (or low ℓ).(Slosar et al. 2008;Ross et al. 2013;Mueller et al. 2022;Cabass et al. 2022), including our analysis with −39 <  NL < 84 (DESI photo LRG) and CMB surveys (Komatsu et al. 2003;Komatsu 2010;Planck Collaboration et al. 2014, 2019).The median  NL value is used in case the maximum likelihood estimate was not reported in the reference.
We use lognormal simulations to estimate the covariance matrices.
As a caveat, this omits the contributions from higher order statistics in the covariance matrix, but we leave that for future as we do not anticipate any major impact on the best fitting estimates of  NL .
In our fiducial analysis, which includes a non-linear treatment of systematics using nine imaging property maps (Galactic extinction, stellar density, depth in 1, and psfsize in ), we obtain  NL = 34 +24(+50) −44(−73) with  = 1 and  = 0.945.This measurement is consistent with recent CMB and LSS measurements, as visualized in Figure 17.The sensitivity of our constraints is explored against  and .The best fitting estimates of  NL decrease as we increase either  or .Specifically, we find that the error on  NL is more sensitive to  than .Compared with the fiducial result, the error increases by more than a factor of two for  = 1.6, and only by 7% for  = 1.25 (Figure 13).The minimum  2 however does not change much, indicating that the impact on the power spectrum fit is negligible.
The signature of local PNG is very sensitive to excess clustering power caused by imaging systematic effects.We have applied a series of robustness tests to investigate the impact of how the galaxy selection function is determined.Specifically, both linear and nonlinear methods are applied using various combinations of imaging systematic maps (including two external maps for the neutral hydrogen column density and photometric calibration error in the z band).We also examine the effect of additional masks based on imaging conditions and survey completeness.Overall, we find that no change in the analysis shifts the maximum likelihood value of  NL to a significantly different value (Figure 15, Figure 16, and Table 7).
Although being essential for the mitigation of imaging systematics, the template-based approach inevitably removes some of the large-scale clustering information.One of the primary highlights of this work is that we present a strategy to calibrate the systematic mitigation's impact on the inferred  NL .As we increase the number of maps for mitigation, more of the power spectrum is removed, introducing a larger bias to the  NL posterior distribution.Our mock tests suggest that this bias is  NL -dependent, such that the mocks with larger  NL experience a more substantial reduction in the low−ℓ power due to systematic mitigation.Therefore, it is crucial to calibrate for this effect using simulations that have gone through the same treatment methods and are subject to the same over-correction effect.
As a greater flexibility in the mitigation increases the overcorrection and decreases the statistical power, we tested if we can reduce the flexibility in our method by using a smaller set of maps, including Galactic extinction, depth in the z band, and astronomical seeing in the r band (nonlinear three maps) to retain some constraining power.Additionally, we consider an additional map for local stellar density (nonlinear four maps).Using three or four maps, we can qualitatively mitigate systematic trends in the mean galaxy density and cross correlations of the galaxy density field and imaging property maps (see Section 3.4).These methods do not degrade the error on  NL as much as the fiducial method which used nine maps.However, when applying our null-tests that are applied in order to detect residual systematic variance (see 3.4), we obtain passing results only for the nonlinear nine map case.In this work, we found updating the covariance matrix for each particular variation (e.g., the mitigation method applied) was important in order to obtain a similar  2 of the null test when applied to the mocks and hence self-consistently obtain a -value for the null test.Another important conclusion from applying our null tests to mocks is that the mean density contrast test is less sensitive to the  NL value for the mocks than the angular cross-power.Given the amount of  NL constraining power that we lose when applying the nine map regression (the uncertainties approximately double), our findings highlight the importance of exploring, developing, and validating alternative mitigation approaches to avoid over-correction for a robust analysis of local PNG.
Our analysis can be considered as the first attempt to identify major systematics in DESI, so we can be ready for constraining  NL with DESI spectroscopy.Internal DESI tests of the photometric calibration were unable to uncover DESI-specific issues, e.g., when comparing to Gaia data.The most significant trends that we find are with the E(B-V) map.The source of such a trend would be a mis-calibration of the E(B-V) map itself or the coefficients applied to obtain Galactic extinction corrected photometry.Such a mis-calibration would plausibly be proportional in amplitude to the estimated E(B-V) map, though it may not have E(B-V)'s spatial distribution.There are ongoing efforts within DESI to obtain improved Galactic extinction information, which will help us address systematics.Additionally, cross-correlations of the DESI LRG density with the CMB lensing map is more stable in terms of systematics and can complement the results presented in this work.We can further avoid the over-fitting issue by combining our neural network-based treatment method with forward-modeling techniques, such as Obiwon (Kong et al. 2020), but we will leave that for future work.

APPENDIX A: EXTRA ROBUSTNESS TESTS A1 Scale dependent systematics
To investigate the statistical significance of the cross power spectrum's  , we examine its dependence on the largest harmonic mode ℓ max .Our fiducial cross power spectrum diagnostic (equation 19) uses harmonic modes up to ℓ = 20, which determines the smallest scale used for characterizing residual systematic errors.We extend ℓ max from 20 to 100, where the latter scale corresponds to density fluctuations on scales smaller than 2 degrees.Figure A1 shows the median of the normalized cross power spectrum's  2 from the clean  NL = 0 mocks after non-linear nine maps as the highest mode ℓ max increases from 20 to 100 (represented by the solid line).The pink circles represent the  2 values for the DESI LRG targets cleaned with the non-linear nine maps method.Overall, we find that for all scales up to ℓ = 100, the nonlinear nine maps approach yields consistent values with the clean mocks.

A2 Survey window convolution
Here we calculate the mode-mode coupling matrix from the DESI mask.This matrix depends only on the survey geometry and can be described in terms of the window power spectrum (Hivon et al. 2002), where the last term in the right hand side represents the Wigner 3-j symbol (or Clebsch-Gordan coefficient), and is calculated using SymPy (Meurer et al. 2017).We benchmark our code against the publicly available software, NaMaster6 (Alonso et al. 2019). NL = 76.9 mocks.When  NL = 0, our config-space convolution of the window aligns with the ℓ-space convolution approach.However, when  NL = 76.9, both the config-space and NaMaster ℓ-space methods yield a convoluted power spectrum with noticeable noise-like numerical artifacts on large scales, therefore possibly in a  NL -dependent manner.To assess the impact of these discrepancies on our  NL constraints, we fit the clustering of the DESI LRG targets, disregarding the integral constraint effect.The best-fitting estimates of  NL will be biased, but our focus is on understanding the relative impact on  NL between the two approaches.For both config-space and ℓ-space methods, we obtain a similar minimum  2 value of 39.6 with 34 degrees of freedom.Notably, the posterior width for the config-space approach is slightly larger than that of the ℓ-space by 10%.The absolute difference in the best-fitting estimates of  NL between the two cases is less than 1.1, considered negligible relative to the statistical precision of our measurements.

A3 Redshift uncertainties
We use the Early Data Assembly Version 1 (EDA V1) to construct the redshift distribution for the DESI LRG targets.We find that the change in the maximum likelihood estimate of  NL is negligible, |Δ  NL | < 1, compared to the statistical precision of our measurements.Figure A3 shows the measured power spectrum of the DESI targets and the corresponding best fit theory curves.The variations in / do not significicantly alter the conclusion of our paper.

A4 Spurious bump in NGC
As shown in Figure 16 (top panel), we realize that the spurious feature at ℓ ∼ 10 − 20 is removed in the DECaLS South region after mitigation, but it remains in the BASS+MzLS and DECaLS North.We use the neural networks trained on the DECaLS South with three and nine maps to mitigate the galaxy density in the DECaLS North region, and then measure the power spectrum.Figure A4 shows the power spectrum before treatment (No Weight) and after the nonlinear three maps and nine maps methods for comparison.We find that whatever causing the bump is different between the DECaLS North and South.The best-fit estimates for  NL from the

APPENDIX B: LOGNORMAL MOCKS
We fit the mean power spectrum of the lognormal mocks to validate the modeling pipeline, and in particular the survey geometry and integral constraint treatments.We investigate the impact of covariance matrix on the inference of  NL .Finally, we show the impact of imaging systematic mitigation and the over-subtraction effect when the cleaning methods are applied to the mocks.

B1 Clean mocks
The 68% and 95% probability contours on the PNG parameter  NL and bias coefficient  are shown in Figure B1 and B2, respectively, for the  NL = 0 and 76.9 mocks.The best-fitting, marginalized mean estimates, as well as the 1 and 2 confidence intervals of  NL are summarized in Table B1.
Measuring the power spectrum from the entire DESI footprint reduces the cosmic variance and thus improves the constraining power.Figure B1 compares the constraints from fitting the log of the mean power spectrum of the mocks when it is measured from the DESI footprint to those obtained from the sub imaging surveys.We find that the underlying true  NL value is recovered within 95% confidence, and that the contours for the DESI region are smaller by a factor of two.
The power spectrum of the mocks at low ℓ is very sensitive to the cosmic variance and the true value of  NL .Consequently, a large value of  NL can induce very large power on low ℓ, and thus significantly change the covariance matrix.We find that applying the log transformation on the power spectrum makes the result more robust against the choice of the covariance matrix.Figure B2 shows the confidence contours when we fit either the power spectrum or its log transform of the  NL = 76.9 mocks, and use different covariance matrices.We consider the  NL = 0 and 76.9 mocks to construct the covariance from one set and use it to fit the mean power spectrum of the other set.When the covariance matrix is constructed from the same set of mocks used for the mean power spectrum, we find that the difference in  NL constraints between fitting the power spectrum and its log transformation is negligible at only 2%.If we use the  NL = 0 mocks to estimate the covariance, and fit the log power spectrum of the  NL = 76.9 mocks, we find that the error on  NL increases only by 7%.However, when the mean power spectrum of the  NL = 76.9 mocks is fit using the covariance matrix estimated from the  NL = 0 mocks, the constraints tighten by a factor of 5 due to a higher signal to noise ratio.Therefore, we argue that fitting the log power spectrum can help mitigate the need for having  NL -dependent covariance matrices and make the constraints less sensitive to covariance construction.
Figure B3 shows the best-fitting estimates for  vs  NL for  NL = 0 and = 76.9 mocks in the top and bottom panels, respectively.Truth values are represented via the dotted lines.The points are color-coded with the minimum  2 from fit for each realization.The histograms of the best-fitting  NL estimates are plotted in the background.For the  NL = 0 mocks, the best-fitting estimates are more symmetric.To understand this behaviour, we consider the first derivative of the likelihood (Equation 16), which is proportional to the first derivative of the log power spectrum.By simplifying the integrals involved in  ℓ , we have  ℓ =  0,ℓ +  1,ℓ  NL +  2,ℓ  2 NL where  123,ℓ are ℓ-dependent terms.Then, the derivative of the likelihood is proportional to For infinitesimal values of  NL , the derivative becomes asymptotically independent from  NL while for large values of  NL it decreases as 2/  NL .This implies that for the  NL = 0 mocks, the likelihood is more likely to be skewed toward negative values.

B2 Contaminated mocks
Our nonlinear neural network-based approach is applied to the  NL = 0 and 76.9 mocks.We only consider the methods that include running the neural network with three, four, and nine imaging systematic maps.The measured mean power spectrum of the mocks are shown in Figure B4 for  NL = 0 (left) and 76.9 (right).The solid and dashed curves show the measurements respectively from the clean and contaminated mocks.We find that the imaging treatment has removed some of the true clustering signal, and the amount of the over-subtraction is almost the same regardless of whether the mocks have systematics.The over-subtraction induces biases in the  NL constraints, as summarized in Table B2.The over-subtraction at low ℓ is so high that we get a poor fit after applying the mitigation with the nonlinear three maps approach, e.g.,  2 = 86.8 for the clean  NL = 0 mocks.
Using the calibration parameters presented in §3.5, we account for the shift in the  NL constraints caused by the imaging systematic mitigation.We show the marginalized probability distributions on  NL before and after accounting for the over-correction in the right and left panels of Figure B5.

Figure 1 .
Figure1.The redshift distribution (solid line and vertical scale on the left) and bias evolution (dashed line and vertical scale on the right) of the DESI LRG targets.The redshift distribution is determined from DESI spectroscopy (DESICollaboration et al. 2023).The redshift evolution of the linear bias is supported by HOD fits to the angular clustering of the DESI LRG targets(Zhou et al. 2021), where  () represents the growth factor.

Figure 2 .
Figure 2. Top: The DESI LRG target density map before correcting for imaging systematic effects in Mollweide projection.The disconnected islands from the North footprint and parts of the South footprint with declination below −30 are removed from the sample for the analysis due to potential calibration issues (see text).Bottom: Mollweide projections of the imaging systematic maps (survey depth, astronomical seeing/psfsize, Galactic extinction, and local stellar density) in celestial coordinates.Not shown here are two external maps for the neutral hydrogen column density and photometric calibration, which are only employed for the robustness tests.The imaging systematic maps are colour-coded to show increasing values from blue to red.

Figure 4 .
Figure 4.The survey mask correlation functions across the imaging regions forming the DESI footprint, plotted against angular separation.The inset focuses on correlations within the angular range of 100 to 180 degrees.

Figure 5 .
Figure5.The mean power spectrum from the  NL = 0 mocks (no contamination) and best-fitting theoretical prediction after accounting for the survey geometry and integral constraint effects.Bottom panel shows the residual power spectrum relative to the mean power spectrum.The dark and light shades represent the 68% error on the mean and one realization, respectively.No imaging systematic cleaning is applied to these mocks.

−Figure 6 .
Figure6.The distribution of the first bin power spectra and its log transformation from the simulations with  NL = 0 (left) and 76.9 (right).The log transformation largely removes the asymmetry in the distributions.
Figure 7.The square of the cross power spectra between the DESI LRG targets and imaging systematic maps normalized by the auto power spectrum of the imaging systematic maps; see equation 18.The systematic maps considered are Galactic extinction (EBV), stellar density (nStar), depth in grzw1 (depth  1 ), and seeing in grz (psfsize   ).The dark green curves display the cross spectra before imaging systematic correction (No Weight).The yellow, brown, and light green curves represent the results after applying the imaging weights from the linear models trained with two maps, three maps, and nine maps.The orange and purple curves display the results after applying the imaging weights from the non-linear models trained with three maps and nine maps.The dark and light shades represent the 97.5 percentile from cross correlating the imaging systematic maps and the  NL = 0 and 76.9 lognormal density fields, respectively, without mitigation.

Figure 8 .
Figure8.The mean density contrast of the DESI LRG targets as a function of the imaging systematic maps: Galactic extinction (EBV), stellar density (nStar), depth in grzw1 (depth  1 ), and seeing in grz (psfsize   ).The black curves display the results before imaging systematic correction.The red, blue and orange curves represent the relationships after applying the imaging weights from the linear models trained with two maps, three maps, and eight maps, respectively.The green and pink curves display the results after applying the imaging weights from the non-linear models trained with three maps and four maps.The dark and light shades represent the 68% dispersion of 1000 lognormal mocks with  NL = 0 and 76.9, respectively.

Figure 9 .
Figure9.The remaining systematic error  2 from the mean galaxy density contrast (left) and the galaxy-imaging normalized cross power spectrum (right).The values observed in the DESI LRG targets after the non-linear treatments are represented via vertical lines using the  NL = 0 (solid) or 76.9 (dashed) covariance, and the histograms are constructed from 1000 realizations of clean lognormal mocks with  NL = 0 (solid) and 76.9 (dashed).

Figure 10 .
Figure10.The No mitigated, clean vs mitigated  NL values from the  NL = 0 and 76.9 mocks.The solid (dashed) lines represent the best-fitting estimates from fitting the mean power spectrum of the clean (contaminated) mocks.The scatter points show the best-fitting estimates from fitting the individual spectra for the clean mocks.

Figure 11 .
Figure11.The power spectrum of the DESI LRG targets before (No weight) and after correcting for imaging systematics using the linear and non-linear methods.The curves represent the corresponding best-fitting theory predictions.The solid curve and grey shade respectively represent the mean power spectrum and 68% error from the  NL = 0 mocks.

Figure 13 .
Figure13.The best-fitting estimates of  NL and their corresponding 68% (95%) errors from the DESI LRG targets using the non-linear nine maps approach given various values of  or .The star symbol represents the fiducial analysis with  = 1 and  = 0.945.

Figure 14 .
Figure 14.Same as Figure 12 but without accouting for over-correction.

Figure 15 .
Figure15.The uncalibrated 2D constraints from the DESI LRG targets using the nonlinear nine maps treatment for each imaging survey compared with that for the whole DESI footprint.The dark and light shades represent the 68% and 95% confidence intervals, respectively.

Figure 16 .
Figure16.Top: The measured power spectrum of the DESI LRG targets before (solid curves) and after non-linear nine maps (scatter points) for the DESI, BASS+MzLS, DECaLS North, and DECaLS South regions.The solid curve and grey shade respectively represent the mean power spectrum and 68% error from the  NL = 0 mocks with the same angular mask for each region.Bottom: The uncalibrated  NL constraints vs the lowest ℓ mode used for fitting  NL .The points represent the best fitting estimates of  NL and error bars represent 95% confidence.The scaling of  NL is not calibrated to account for over-correction caused by mitigation.

Figure 17 .
Figure17.History of constraints on local  NL at 95% confidence from single-tracer LSS(Slosar et al. 2008;Ross et al. 2013;Mueller et al. 2022;Cabass et al. 2022), including our analysis with −39 <  NL < 84 (DESI photo LRG) and CMB surveys(Komatsu et al. 2003;Komatsu 2010; Planck  Collaboration et al. 2014, 2019).The median  NL value is used in case the maximum likelihood estimate was not reported in the reference.

Figure A1 .
Figure A1.The cross power spectrum's  2 between the DESI LRG target density and imaging systematic maps as a function of the highest mode ℓ max when the sample is cleaned with the linear (triangles) and non-linear (squares) three maps.The lowest mode is fixed at ℓ min = 2.The solid curve and dark (light) shade represent the median value and 68% (95%) confidence regions, estimated from the  NL = 0 mocks.

Figure A2 .
Figure A2.The model power spectrum before and after the survey geometry convolution for  NL = 0 and 76.9 using the DESI survey mask.The bottom panel shows the residual error with respect to the NaMaster code.The shade represents the dispersion of the  NL = 76.9 mocks.

Figure A3 .Figure A4 .
Figure A3.Top: The redshift distribution of the DESI LRG targets from the EDA V1 and Denali.Bottom: The power spectrum of the DESI LRG targets and the best fit theory models using different redshift distributions.

Figure B1 .
Figure B1.68% and 95% confidence contours from the mean power spectrum of the  NL = 0 mocks for the DESI footprint and sub-imaging surveys.The truth values are represented by vertical and horizontal lines.

bFigure B3 .
Figure B2.68% and 95% confidence contours of fitting the mean power spectrum or its log transformation from the  NL = 76.9 mocks for the DESI footprint.Using the log  ℓ fitting yield constraints that are insensitive to the covariance used.The truth values are represented by vertical and horizontal lines.

Figure B5 .
Figure B5.Probability distributions of  NL from the mean power spectrum of the  NL = 0 (top) and  NL = 76.9 (bottom) mocks before and after mitigation with the non-linear methods using three, four, and nine maps.The dashed (solid) curves show the distributions for the contaminated (clean) mocks.Left: The posteriors are adjusted to account for the over-correction effect.Right: The posteriors are subject to the over-correction effect, and thus the scaling of  NL values is biased due to mitigation.

Table 1 .
Color-magnitude selection criteria for the DESI LRG targets

Table 2 .
Statistics for DESI imaging data.Median depths are for galaxy/point sources detected at 5.Median psfsize values are computed with a depth-weighted average at each location on the sky.

Table 4 .
Linear parameters employed to de-bias the  NL constraints to account for the over-correction issue.

Table 5 .
The calibrated best-fitting, marginalized mean, and marginalized 68% (95%) confidence estimates for  NL from fitting the power spectrum of the DESI LRG targets before and after correcting for imaging systematic effects.The lowest mode is ℓ min = 2,  = 1, and  = 0.945.68% and 95% probability distribution contours for the bias and  NL from the DESI LRG targets before and after applying the non-linear cleaning methods.The lognormal mocks are used to calibrate these distributions for over-correction.
Figure 12.The calibrated constrains from the DESI LRG targets.Top: probability distribution for  NL marginalized over the shotnoise and bias.Bottom:

Table 6 .
The calibrated best-fitting, marginalized mean, and marginalized 68% (95%) confidence estimates for  NL from the DESI LRG targets cleaned with the non-linear nine maps approach, given various values of  and .The fiducial analysis uses  = 1 and  = 0.945 (DESI footprint).