© The Authors 2025

1. Introduction

Cosmic inflation was proposed as a theory in the early 1980s (Guth 1981; Starobinsky 1980). The inflation framework was initially formulated to solve big bang problems such as the horizon, flatness, and magnetic monopole problems; however, inflation is also able to explain the formation of primordial density perturbations (Starobinsky 1982; Guth & Pi 1985; Bardeen et al. 1983). Inflation is defined as a phase in which the Universe is expanding exponentially, driven by a scalar field, ϕ. Several models of inflation have been proposed in the literature (see e.g., Langlois 2010; Vazquez Gonzalez et al. 2020 for a review). The model of inflation and its predictions are defined by choosing the form of the potential, V(ϕ). The simplest inflationary models predict Gaussian initial conditions; however, alternative inflationary models predict different levels of non-Gaussianity in the primordial density perturbations (Chen 2010; Takahashi 2014). The level of non-Gaussianity has usually been characterized in the literature with the f_NL non-Gaussianity parameter, such that detecting f_NL ≠ 0 is a signature of having non-Gaussian initial conditions.

The tightest constraint on f_NL is currently provided by the measurements from the cosmic microwave background (CMB) bispectrum. Using Planck 2018 data, a value f_NL = −0.9 ± 5.1 at the 68% confidence level is found (Planck Collaboration IX 2020). However, Dalal et al. (2008) first noticed that local non-Gaussian initial conditions lead to a characteristic scale-dependent signature in the galaxy bias, following a 1/k² scale-dependence in the ratio between the total matter to observed density of galaxies. During the last decade, many works have performed measurements from the large-scale structure using the scale-dependent bias effect (Ross et al. 2012; Castorina et al. 2019; Mueller et al. 2022; Cabass et al. 2022; D’Amico et al. 2025, among others). In the last years, some measurements from large scale structure (LSS) using the scale-dependent galaxy bias have been achieved using the quasars from the data release DR16 of the Extended Baryon Oscillation Spectroscopic Survey (eBOSS): Mueller et al. (2022) measure f_NL = −12 ± 21 and Cagliari et al. (2024) obtain −4 < f_NL < 27 using different methodologies. More recently, Chaussidon et al. (2024) have improved this constraint to $f_{\mathrm{NL}}=-3.6_{-9.0}^{+9.1}$ using the 3D power spectrum of the Dark Energy Spectroscopic Instrument (DESI) data release DR1 galaxies and quasars. This is a challenging measurement because it requires very accurate control of the largest scales where the scale-dependent bias effect due to f_NL arises. Further, it is also important to mention that, when performing f_NL measurements from LSS, unless using certain assumptions we actually measure the product of f_NL times an unknown bias (see Sect. 2.1 for more details).

DESI (Levi et al. 2013) is a spectroscopic survey that is currently being carried out from the 4 m Mayall telescope at Kitt Peak National Observatory (Arizona, USA). Its unique design with 5000 fibers with robotic positioners allows it to take thousands of spectra in a single exposure (DESI Collaboration 2016a, 2022; Silber et al. 2023; Miller et al. 2024; Guy et al. 2023; Schlafly et al. 2023). Theoretical forecasts (DESI Collaboration 2016b) expect that the full five-year survey will have the ability to achieve a sensitivity σ(f_NL)∼5, similar to the best current CMB bispectrum constraint, if there is a good control of systematic effects. Before the first spectroscopic data releases (DESI Collaboration 2024a, 2024b) and science results (DESI Collaboration 2025a, 2025b, 2025c; DESI Collaboration 2024c, 2024d) came out, a full imaging survey was performed in order to select the spectroscopic targets. This targeting survey is called the DESI Legacy Survey (Dey et al. 2019) and covers a broad area (≳20 000 deg²), making it useful for measuring f_NL using the scale-dependent galaxy bias. Two previous works have already used the DESI Legacy Survey information to put a constraint on f_NL: Rezaie et al. (2024) used the angular power spectra of the louminous red galaxy (LRG) targets, and Krolewski et al. (2024) used the cross-correlation between quasar targets and the Planck CMB lensing.

CMB lensing describes the remapping of the CMB anisotropies due to gravitational lensing by structures along the line of slight. The CMB lensing potential can be easily measured from the observations of the lensed sky (Hu & Okamoto 2002) and was first detected by Smith et al. (2007). Since it contains information about the large-scale structure geometry, its cross-correlation with galaxy tracers can be useful to constrain cosmology. Although CMB lensing and galaxy tracers probe the same structures, they are affected by different systematics, making the cross-correlation between the two a powerful additional tool for measurements into the systematics-dominated regime. Several authors have stressed, using theoretical forecasts, the capabilities of the cross-correlation between CMB lensing and galaxy matter tracers to better constrain f_NL (e.g., Jeong et al. 2009; Schmittfull & Seljak 2018; Giusarma et al. 2018; Ballardini et al. 2019; Bermejo-Climent et al. 2021; Perna et al. 2023). More recently, Krolewski et al. (2024) find $f_{\mathrm{NL}}=-{26}_{-40}^{+45}$ using the cross-correlation between Planck lensing and DESI quasar targets. Additionally, recent data analysis works have been performed to constrain other cosmological parameters such as the amplitude of matter density fluctuations, commonly parametrized in terms of σ₈ (the RMS density contrast smoothed on a scale of 8 h/Mpc), and matter density, Ω_m, using CMB cross-correlations with the DESI Legacy Survey (e.g., White et al. 2022a, Sailer et al. 2024; Kim et al. 2024).

In this paper, we intend to extend the analysis done by Rezaie et al. (2024) with the DESI LRG sample to the inclusion of the CMB lensing cross-correlation as additional information and an additional technique to limit the impact of observational systematics. In Rezaie et al. (2024), an extensive and detailed effort was made to remove observational systematics that could bias the primordial non-Gaussianity (PNG) measurement. Nonetheless, they conclude that their results motivate further studies of PNG with samples less sensitive to systematics like LRG spectroscopic data. Here, our aim is to explore the ability of CMB lensing – LRG cross-correlation to constrain f_NL and its stability in terms of systematics, alone and in combination with the LRG autocorrelation.

In Sect. 2 we review the theoretical framework for the imprints of f_NL in a scale-dependent galaxy bias and for the cosmological observables we study in the angular domain. In Sect. 3 we present the DESI LRG and Planck lensing datasets used for our analysis. In Sect. 4 we describe the methodology followed for treating our datasets, including a systematics mitigation, computation of observables, and parameter inference pipeline. In Sect. 5 we discuss a validation of our analysis pipeline with mock LRG and CMB lensing fields. In Sect. 6 we present our results for the measurement of PNG and some robustness tests, and in Sect. 7 we summarize our conclusions.

2. Theory

In this section we first provide a description of the physical model that originates a scale-dependent galaxy bias due to a local PNG. We then describe the cosmological observables in the 2D harmonic space involved in our analysis.

2.1. Primordial non-Gaussianity and scale-dependent bias

If we assume a type of non-Gaussianity that depends only on the local value of the potential, the parametrization of the primordial potential can be written as (Komatsu & Spergel 2001)

$$\begin{array}{c}\hfill \mathrm{\Phi }=\varphi +f_{\mathrm{NL}}({\varphi }^2-{⟨\varphi ⟩}^2),\end{array}$$(1)

where f_NL is the parameter that describes the amplitude of the quadratic non-Gaussian term and ϕ is a random Gaussian field. We studied f_NL through its impact on the scale-dependent galaxy bias, as introduced in Dalal et al. (2008). If we assume the so-called universality relation (Slosar et al. 2008), the contribution to the galaxy bias is expressed as

$$\begin{array}{c}\hfill \mathrm{\Delta }b(k,z)=2(b_g-p)f_{\mathrm{NL}}\frac{{\delta }_{\mathrm{crit}}}{\alpha (k)},\end{array}$$(2)

where δ_crit = 1.686 is the threshold overdensity for spherical collapse, b_g is the galaxy bias as a function of the redshift z, p is a parameter characterizing the tracers’ response to PNG, assumed as p ≃ 1 for the case of LRG, and α(k) is the relation between potential and density field, such that δ(k) = α(k)Φ(k). The value of α(k) is given by

$$\begin{array}{c}\hfill \alpha (k)=\frac{2k^{2}T(k)D(z)}{3{\mathrm{\Omega }}_{\mathrm{m}}}\frac{c^2}{H_0^2}\frac{g(0)}{g(\infty )},\end{array}$$(3)

where T(k) is the transfer function, D(z) is the growth factor (normalized to be 1 at z = 0), Ω_m the matter density, and the factor g(∞)/g(0)≃1.3 accounts for the different normalizations of D(z) in the CMB and LSS literature. This definition of f_NL is therefore the so-called CMB convention. Other authors (e.g., Carbone et al. 2008; Afshordi & Tolley 2008; Grossi et al. 2009) refer to the use of the ‘LSS convention’, where the g(∞)/g(0) factor is absorbed into the definition of f_NL, such that $f_{\mathrm{NL}}^{\mathrm{LSS}}\simeq 1.3f_{\mathrm{NL}}^{\mathrm{CMB}}$. Throughout this paper we use p = 1 as a baseline; however, many works based on dark-matter-only simulations (e.g., Adame et al. 2024) find significant deviations from p = 1. In this direction, other authors (e.g., Barreira 2020, 2022) have stressed that we can only constrain the product b_ϕf_NL through the scale-dependent bias effect, where b_ϕ is a parameter usually defined as b_ϕ = 2δ_crit(b_g − p), in order to account for the uncertainties on p.

2.2. Cosmological observables

In this work, we focus on the angular power spectrum of the galaxy – CMB lensing cross-correlation, C_ℓ^κG, as well as the galaxy autocorrelation, C_ℓ^GG. The angular power spectrum can be calculated as

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{XY}}=4\pi {\displaystyle \int \frac{\mathrm{d}k}{k}}\mathcal{P}(k)I_{\ell }^X(k)I_{\ell }^Y(k),\end{array}$$(4)

where 𝒫(k) ≡ k³P(k)/(2π²) is the dimensionless primordial power spectrum and I_ℓ^X(k) is the kernel for the X field for unit of primordial power spectrum.

All the weak lensing quantities can be defined from the lensing potential

$$\begin{array}{c}\hfill \varphi (\widehat{\mathbf{n}},\chi )=\frac{2}{c^2}{\displaystyle {\int }_0^{\chi }}\mathrm{d}\chi \prime \frac{\chi -\chi \prime }{\chi \chi \prime }\mathrm{\Phi }(\chi \prime \widehat{\mathbf{n}},\chi \prime ),\end{array}$$(5)

where $\mathrm{\Phi }(\widehat{\mathbf{n}},\chi )$ is the gravitational potential. The comoving distance is

$$\begin{array}{c}\hfill \chi (z)={\displaystyle {\int }_0^z}\frac{c\mathrm{d}z\prime }{H(z\prime )}·\end{array}$$(6)

The observable two-dimensional lensing potential, averaged over background sources with a redshift distribution W_b(χ), is given by

$$\begin{array}{c}\hfill \varphi \left(\widehat{\mathbf{n}}\right)=\frac{2}{c^2}{\displaystyle {\int }_0^{\chi }}\frac{\mathrm{d}\chi \prime }{\chi \prime }{q}_b(\chi \prime )\mathrm{\Phi }(\chi \prime \widehat{\mathbf{n}},\chi \prime ),\end{array}$$(7)

where q_b(χ) is the lensing efficiency for a given background distribution, W_b, defined as

$$\begin{array}{c}\hfill q_b\left(\chi \right)={\displaystyle {\int }_0^{\chi }}\mathrm{d}\chi \prime \frac{\chi \prime -\chi }{\chi \prime }{W}_b(\chi \prime ).\end{array}$$(8)

By expanding the gravitational potential in Fourier space and using the plane-wave expansion, we can define the lensing potential kernel as

$$\begin{array}{c}\hfill I_{\ell }^{\varphi }(k)=2\left(\frac{3{\mathrm{\Omega }}_m{H}_0^2}{2k^2{c}^2}\right){\displaystyle \int \frac{\mathrm{d}\chi }{{(2\pi )}^{3/2}}}\frac{q_b\left(\chi \right)}{\chi a(\chi )}{j}_{\ell }\left(k\chi \right)\delta (k,\chi ),\end{array}$$(9)

where Ω_m is the present-day matter density, a(χ) is the scale factor, H₀ is the Hubble constant, δ(k, χ) is the comoving-gauge linear matter density perturbation, and j_ℓ the spherical Bessel functions. In case of CMB lensing, the source distribution can be approximated by W_CMB(χ) ≃ δ_D(χ−χ_*) and the lensing efficiency by

$$\begin{array}{c}\hfill q_{\mathrm{CMB}}\left(\chi \right)\simeq \frac{{\chi }_{\ast }-\chi }{{\chi }_{\ast }},\end{array}$$(10)

where χ_* is the comoving distance at the surface of last scattering, and Eq. (9) reduces to

$$\begin{array}{c}\hfill I_{\ell }^{{\varphi }_{\mathrm{CMB}}}(k)=2\left(\frac{3{\mathrm{\Omega }}_m{H}_0^2}{2k^2{c}^2}\right){\displaystyle \int \frac{\mathrm{d}\chi }{{(2\pi )}^{3/2}}}\frac{{\chi }_{\ast }-\chi }{{\chi }_{\ast }\chi }\frac{1}{a(\chi )}{j}_{\ell }\left(k\chi \right)\delta (k,\chi ).\end{array}$$(11)

Finally, the convergence κ = ∇²ϕ/2 can be expanded in spherical harmonics as

$$\begin{array}{c}\hfill \kappa \left(\widehat{\mathbf{n}}\right)=-\frac{1}{2}{\displaystyle \sum _{\ell ,m}}\ell (\ell +1){\varphi }_{\ell m}{Y}_{\ell }^m\left(\widehat{\mathbf{n}}\right).\end{array}$$(12)

We then relate the two kernels by

$$\begin{array}{c}\hfill I_{\ell }^{\kappa }(k)=\frac{\ell (\ell +1)}{2}{I}_{\ell }^{\varphi }(k).\end{array}$$(13)

The two-dimensional integrated window function for the galaxy number counts is

$$\begin{array}{c}\hfill I_{\ell }^G(k)={\displaystyle \int \frac{\mathrm{d}\chi }{{(2\pi )}^{3/2}}}W(\chi ){\mathrm{\Delta }}_{\ell }(k,\chi ),\end{array}$$(14)

where Δ_ℓ(k, χ) is the observed number count and W(χ) is a window function given by the redshift distribution of galaxies. At first order, the most important contribution to Δ_ℓ(k, χ) is given by the synchronous gauge source counts Fourier transformed and expanded into multipoles, Δ_ℓ^s(k, χ). We assumed that Δ_ℓ^s(k, χ) is related to the underlying matter density field through a scale- and redshift-dependent galaxy bias, b_g, as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\ell }^s(k,\chi )=b_g(k,\chi )\delta (k,\chi )j_{\ell }\left(k\chi \right),\end{array}$$(15)

where b_g(k, χ) is given by the sum of a linear bias, which is not scale dependent, plus the scale-dependent contribution given by Eq. (2). We also considered nonlinear contributions to the power spectrum using halofit (Takahashi et al. 2012). In this paper, we also consider two important contributions to the total observed number counts: the effects of redshift space dirtortions (RSD) and lensing magnification (see Fig. 4 for more details). The RSD term is given by

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\ell }^{\mathrm{RSD}}(k,\chi )=\frac{k{v}_k}{\mathcal{H}}{j}_{\ell }″(k\chi ),\end{array}$$(16)

where v_k is the velocity of the sources and ℋ is the conformal Hubble parameter. The lensing convergence contribution is given by

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\ell }^{\mathrm{lensing}}(k,\chi )=\ell (\ell +1)(2-5s){\displaystyle {\int }_0^{\chi }}\mathrm{d}\chi \prime \frac{\chi -\chi \prime }{\chi \chi \prime }\mathrm{\Phi }(\chi \prime \widehat{\mathbf{n}},\chi \prime ),\end{array}$$(17)

where s is the magnification bias that accounts for the fact that observed galaxies are magnified by gravitational lensing. We note that there are other contributions from general relativity to the number counts, but we consider them of second order since the most important contribution to C_ℓ^κG is the lensing magnification (see e.g., Appendix A of Bermejo-Climent et al. 2021).

In practice, to numerically evaluate the observables described in this section, we used the Limber approximation (Limber 1953) for ℓ ≥ 100 and computed the full quantities for ℓ < 100, following the implementation in the CAMBSOURCES module (Challinor & Lewis 2011). This was done in order to improve the computational efficiency while having a good estimation of the theory model at the lower multipoles, where the Limber approximation is less accurate. More details on the code used for the computation of theoretical observables can be found in Sect. 4.3.

3. Datasets

We describe in this section the datasets involved in our analysis. The two main ingredients were an LRG photometric catalog from the DR9 Legacy Survey (Zhou et al. 2023a) and the Planck PR4 public CMB lensing maps (Carron et al. 2022). We also used an LRG spectroscopic sample from the DESI Survey Validation data to calibrate the redshift distribution of the photometric DESI LRGs.

3.1. Luminous red galaxies

Our galaxy sample consisted of an LRG catalog obtained from the DESI Imaging Legacy Surveys¹ DR9 (Dey et al. 2019). These surveys were a combination of three projects: the Dark Energy Camera Legacy Survey using the Blanco 4m telescope in Chile (DECaLS, Flaugher et al. 2015), the Beijing-Arizona Sky Survey using the Bok telescope at Kitt Peak (BASS, Zou et al. 2017), and the Mayall z-band Legacy Survey (MzLS, Dey et al. 2019) using the Mayall telescope at Kitt Peak. BASS and MzLS observed the same region in the north galactic cap (NGC), while DECaLS observed in both the NGC and south galactic cap (SGC). The combination of the three projects covered ∼19 000 deg² in the sky to select the spectroscopic targets that are currently being observed with DESI.

We calibrated the redshift distribution, dN/dz, of the sample using the actual LRG spectra measured with the DESI Survey Validation (DESI Collaboration 2024a, 2024b). In Fig. 1 we show the redshift distribution of the sample as obtained from the LRG spectra. These spectra are not available for declinations lower than −30°, hence we removed the Dec <  − 30° region from the photometric LRG footprint. Zhou et al. (2023a) also describe the presence of a photometric zero-point systematic effect at low declinations. The resulting final sample contains around 9 million galaxies covering a ∼16 000 deg² area. Then, we applied to the LRG catalog the mask designed by Zhou et al. (2023b) to reduce contamination from effects such as stars and foregrounds. We pixelized the LRG catalog, converting it into a HEALPix (Gorski et al. 2005) galaxy counts map at N_side = 256. This map is corrected for the pixel incompleteness effect, which accounts for area losses on scales smaller than a N_side = 256 HEALpix pixel, such as cutouts around bright stars. Lastly, the galaxy counts map can be easily converted into an overdensity field by normalizing and subtracting the mean density. We show the LRG overdensity field in Fig. 2.

Fig. 1. Normalized redshift distribution of the LRG sample, directly measured using the spectroscopic redshifts from DESI Y1 data.

Fig. 2. Left panel: LRG overdensity field in galactic coordinates, after applying the δ < −30° cut. Right panel: CMB lensing field in galactic coordinates, obtained from the Planck PR4 maps.

We note that, since this sample contains photometric redshifts, one could design optimal weights in order to enhance the f_NL signal by emphasizing the higher redshift part of the LRG sample. We do not perform this kind of analysis here since we consider it to be beyond the scope of this paper.

3.2. CMB lensing

The other ingredient in our analysis was the Planck CMB lensing potential map. We used the Planck PR4 reconstruction of the CMB lensing potential (Carron et al. 2022), which was obtained from the Planck NPIPE temperature and polarization maps (Akrami et al. 2020). In particular, we used the minimum-variance estimate from temperature and polarization, after mean-field subtraction of the lensing convergence. This latest release of CMB lensing maps improves the noise with respect to the previous Planck PR3 maps (Planck Collaboration VI 2020); in particular, the large-scale noise is lower and the mean-field is better understood thanks to the larger number of simulations. The maps and the mask are publicly available². We show the CMB lensing field in Fig. 2.

We note that the CMB lensing map does not have the Monte Carlo multiplicative correction applied in Carron et al. (2022). We computed this correction based on simulations as in Krolewski et al. (2024), using mode-decoupled pseudo-C_ℓ, and applied the result as a multiplicative factor to the measured cross-correlation angular power spectra, C_ℓ^κG. This correction cannot be applied in a general way, since it depends on the footprint mask for each tracer involved in the analysis, due to local variations of the normalization. The order of this correction is generally ≲5%, but it becomes larger (∼10–12%) for the largest scales, having thus a non-negligible impact on the f_NL measurements.

4. Analysis pipeline

In this section, we describe the pipeline implemented to analyze the LRG and CMB lensing HEALPix maps. The first step was to apply an imaging systematics mitigation code to the LRG maps. This mitigation treatment operates at the map level and returns a series of systematic weights for each pixel that are applied to the raw LRG maps. Then, we computed the angular power spectrum, C_ℓ, of the LRGs, their cross-correlation with the Planck lensing, and the covariance matrices. The final step was to perform a Monte Carlo Markov chain (MCMC) based parameter inference to constrain f_NL.

4.1. Imaging systematics mitigation

Imaging systematics due to effects such as extinction, stellar contamination, or changes in the observational conditions usually generate an excess of power at large scales (low multipoles), where the f_NL signal arises. Thus, an efficient imaging systematics treatment is key for measuring an unbiased f_NL. Rezaie et al. (2020) present a neural network approach for systematics mitigation as a way to model the relation between the observed galaxy density fields and the imaging systematics templates. This pipeline is implemented in the SYSnet code, which is publicly available³. In Rezaie et al. (2024), a detailed study of the performance of SYSnet was done using this LRG sample, with the aim of measuring f_NL from the C_ℓ^GG autospectra. A different number of approaches were explored, given that SYSnet can return different results depending on the selection of features (imaging systematics template maps) used to perform the regression.

In this work, we used a mitigation recipe optimized for measurements from the cross-correlation between the DESI LRG and Planck lensing, C_ℓ^κG. At first order, the cross-correlation itself as a technique would be enough to remove the effects of the systematics if we assumed there were no correlated systematics between the CMB lensing and LRG maps. However, there could still be potential correlated systematics due to galactic foregrounds that affect both probes. Furthermore, the noise in C_ℓ^GG contributes to the C_ℓ^κG covariance matrix: systematics affecting the C_ℓ^GG power spectrum would also lead to an increase in the C_ℓ^κG covariance, and as a result, larger uncertainties on f_NL. At the same time, a very aggressive mitigation recipe could remove real clustering signal and overfit the C_ℓ^κG power spectrum, biasing the f_NL measurements toward lower values. For this reason, in order to find a compromise between removing the excess of power spectrum due to systematics for C_ℓ^GG and avoiding a strong overfit of C_ℓ^κG, we selected the “nonlinear three maps” recipe from Rezaie et al. (2024). This choice is the most conservative one among the recipes applied in Rezaie et al. (2024) with the SYSnet code, and relies on selecting three features or systematics templates for the regression: extinction, galactic depth in the z-band, and point spread function (PSF) size in the r-band. We refer the reader to Rezaie et al. (2024) for more details on the different possibilities of feature selection. In Sect. 5 we also show that this prescription provides unbiased measurements of the angular power spectrum on contaminated mocks.

4.2. Angular power spectra and covariance matrix

In order to estimate the angular power spectra of the LRG autocorrelation and their cross-correlation with the Planck lensing, we used the pseudo-C_ℓ approach implemented in the publicly available NaMaster code by Alonso et al. (2019). The pseudo-C_ℓ of a pair of fields can be defined as

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_{\ell }^{\mathit{XY}}=\frac{1}{2\ell +1}{\displaystyle \sum _m}{X}_{\ell m}{Y}_{\ell m}^{\ast },\end{array}$$(18)

where X and Y are the observed fields. Then, the difference between the true and measured C_ℓ due to the effects of the mask is accounted for through the mode-coupling matrix M_ℓℓ′ as

$$\begin{array}{c}\hfill ⟨{\stackrel{\sim }{C}}_{\ell }⟩={\displaystyle \sum _{\ell \prime }}{M}_{\ell \ell \prime }{C}_{\ell \prime }.\end{array}$$(19)

In our case, we directly applied a completeness mask to the observed LRG density, so we just deconvolved the binary footprint mask. In practice, the inversion of the M_ℓℓ′ matrix is done using the MASTER algorithm (Hivon et al. 2002), which requires a discrete binning of the angular power spectrum. We used the implementation in the compute_full_master function from the NaMaster code to calculate C_ℓ^κG and C_ℓ^GG. We binned the theory curves using the same bandpower window functions.

As scale cuts, for C_ℓ^κGwe adopted ℓ_min = 2 and ℓ_max = 300, and we binned the power spectra using Δℓ = 5 for ℓ < 100 and Δℓ = 50 for ℓ > 100. The reason for using this scheme was to obtain a good sampling of the angular power spectra at large scales where the f_NL signal arises and ℓ_max = 300 is safe enough to minimize the limitations of modeling the nonlinear scales, which could potentially affect the measurement of the linear bias. We checked, with mock fields, that the results are stable with both N_side and ℓ_max (see Sect. 5 for more details on the mocks). For C_ℓ^GG, we implemented the same multipole binning scheme, but we dropped the first multipole bin of the analysis and set ℓ_min = 7. The motivation for this choice was to use a scale cut for C_ℓ^GG that limits the impact of remaining systematics according to the mitigation pipeline tests on mocks (see Sect. 5 for more details).

We show in Fig. 3 the measured power spectra of the LRG autocorrelation, C_ℓ^GG, and its cross-correlation with the Planck lensing, C_ℓ^κG, obtained using uncorrected (raw) LRG maps and LRG maps mitigated with the neural network pipeline. For the covariance matrices, we used the analytic gaussian_covariance function in NaMaster (García-García et al. 2019) to compute the full Gaussian covariance for a masked field. We accounted for the potential extra power in C_ℓ^GG due to systematics by smoothing the measured angular power spectra from the data and using this as input for the covariance matrix computation. As a test, we also computed the covariance using mock fields, obtaining compatible results. More details of the mock fields used can be found in Sect. 5. We show in Fig. 5 the computed joint correlation matrix for C_ℓ^GG and C_ℓ^κG.

Fig. 3. Angular power spectra of the LRG autocorrelation, CℓGG, and CMB lensing – LRG cross-correlation, CℓκG, obtained from the raw (uncorrected) data and after applying the systematics mitigation pipeline to the LRG maps.

4.3. Likelihood and parameter inference

We defined the likelihood, ℒ, as

$$\begin{array}{c}\hfill -2log\mathcal{L}\equiv {\chi }^2={\displaystyle \sum _{\ell ,\ell \prime }}(C_{\ell }^{\mathrm{obs}}-{\stackrel{\sim }{C}}_{\ell }(\theta )){\mathrm{Cov}}_{\ell \ell \prime }^{-1}(C_{\ell \prime }^{\mathrm{obs}}-{\stackrel{\sim }{C}}_{\ell \prime }(\theta )),\end{array}$$(20)

where C_ℓ^obs are the elements of the data vector, ${\stackrel{\sim }{C}}_{\ell }(\theta )$ is the theoretical model of the angular power spectrum for a given set of parameters, θ, and Cov⁻¹ is the inverse of the covariance matrix.

Our theoretical model is based on the Code for Anisotropies in the Microwave Background, CAMB (Lewis & Challinor 2011). This code is publicly available⁴ and allows the computation of the angular power spectra of the CMB fields and also galaxy fields, given a window function determined by the galaxy bias, b_g, and the redshift distribution, dN/dz. We slightly modified the code to include the scale dependence of the galaxy bias induced by f_NL. The dN/dz was directly measured using the DESI LRG spectra from the Survey Validation (DESI Collaboration 2024a, 2024b). The rest of the fiducial cosmological parameters in CAMB are fixed to the Planck 2018 best-fit estimates (Planck Collaboration VI 2020), except for σ₈, which we set to a fiducial value σ₈ = 0.77 to be in agreement with the measurements from the cross-correlation between the Atacama Cosmology Telescope (ACT) lensing and this LRG sample (Sailer et al. 2024; Kim et al. 2024).

In the computation of the theoretical angular power spectra with the CAMB code, we included the lensing magnification contribution to the galaxy number counts. This effect accounts for the fact that the light from distant galaxies is affected by structures along the line of sight, resulting in an increased flux. The lensing magnification contribution to the angular power spectra could be included in CAMB once we specified the magnification bias, s, as an input parameter. This parameter depends on the galaxy tracer, and for this LRG sample we fixed the magnification bias value to s = 0.999, as determined in Kitanidis & White (2021). More recent measurements using this sample split into four redshift bins (White 2022b; Zhou et al. 2023a) find compatible results within the uncertainty for the various redshift bins; hence, we considered as safe enough the assumption of a z-independent value for s. In Fig. 4 we show the relative importance of the lensing magnification contribution to the C_ℓ^GG and C_ℓ^κG theoretical angular power spectra, which can be expressed as $C_{\ell }^{XX,\mathrm{density}+\mathrm{lensing}}/C_{\ell }^{XX,\mathrm{density}}-1$. This was calculated with CAMB by adding or neglecting the contribution to the galaxy number counts described in Eq. (17). The impact at the lowest multipoles can reach up to ∼15% for the LRG autocorrelation and ∼35% for the LRG – CMB lensing cross-correlation; hence, this effect is not negligible in our theoretical model for this tracer.

Fig. 4. Relative impact of the lensing magnification contributions for the CℓGG (blue line) and CℓκG (red line) theoretical angular power spectra obtained with CAMB.

Fig. 5. Correlation matrix as obtained from the joint covariance matrix for CℓGG and CℓκG computed with the NaMaster code.

To compute the constraints on the cosmological parameters, we implemented our likelihood using the MCMC sampler emcee⁵ (Foreman-Mackey et al. 2013). Our analysis includes f_NL and the galaxy bias at z = 0, b₀, as the two main cosmological parameters of interest to constrain. We assumed a fiducial redshift evolution of the galaxy bias following

$$\begin{array}{c}\hfill b_g=b_0\times D{(z)}^{-1},\end{array}$$(21)

where D(z) is the growth factor normalized to be 1 at z = 0. This choice was motivated by the analysis in Zhou et al. (2021), where a bias evolution compatible with b(z)≃1.5/D(z) was found for the DESI LRG targets. For C_ℓ^GG, we also considered the shot noise, N_shot, as a nuisance parameter. We did not impose any prior knowledge on the b₀ and N_shot parameters, but restricted their sampling to positive values to avoid nonphysical results.

5. Validation with mocks

Before applying the analysis pipeline described in the previous section to the real LRG and CMB lensing data, we tested it with mock Gaussian fields that simulated the DESI Legacy Survey LRG sample and the Planck lensing observations. For the LRG, we used a sample of 100 Gaussian fields for f_NL = 0, 50, and −50. For the CMB lensing, we used a set of 100 correlated maps. These correlated fields were generated using the healpy⁶ Python package.

We first used CAMB to compute the theoretical angular power spectra C_ℓ^GG, C_ℓ^κG, and C_ℓ^κκ, given the redshift window function dN/dz obtained from the spectroscopic LRG redshift distribution and the different values of f_NL. The galaxy bias parameter at z = 0 was set to a fiducial value of b₀ = 1.5, based on Zhou et al. (2021). The fiducial cosmological parameters were set according to the Planck 2018 measurements (Planck Collaboration VI 2020).

The theoretical angular power spectra were used as input for the healpy.synfast function. Following the recipe in Giannantonio et al. (2008), Serra et al. (2014) for simulating correlated fields, we first built a CMB lensing map with a given seed and power spectrum, C_ℓ^κκ. Then, we generated another map with the same seed and power spectrum, (C_ℓ^κG)²/C_ℓ^κκ, and we added to this last map another component with a different seed and power spectrum, C_ℓ^GG − (C_ℓ^κG)²/C_ℓ^κκ. These maps will have amplitudes given by

$$\begin{array}{cc}\hfill a_{\ell m}^{\kappa \kappa }& ={\xi }_1{(C_{\ell }^{\kappa \kappa })}^{1/2},\hfill \\ \hfill a_{\ell m}^{\mathit{GG}}& ={\xi }_1{C}_{\ell }^{\kappa G}/{(C_{\ell }^{\kappa \kappa })}^{1/2}+{\xi }_2{(C_{\ell }^{\mathit{GG}}-{(C_{\ell }^{\kappa G})}^2/C_{\ell }^{\kappa \kappa })}^{1/2},\hfill \end{array}$$(22)

where ξ₁ and ξ₂ are random amplitudes, that is, complex numbers with zero mean and unity variance. It is easy to show that

$$\begin{array}{cc}\hfill ⟨a_{\ell m}^{\kappa \kappa }{a}_{\ell m}^{\kappa \kappa \ast }⟩& =C_{\ell }^{\kappa \kappa },\hfill \\ \hfill ⟨a_{\ell m}^{\kappa \kappa }{a}_{\ell m}^{GG\ast }⟩& =C_{\ell }^{\kappa G},\hfill \\ \hfill ⟨a_{\ell m}^{\mathit{GG}}{a}_{\ell m}^{GG\ast }⟩& =C_{\ell }^{\mathit{GG}}.\hfill \end{array}$$(23)

The output products are a subsample of 100 mock CMB lensing fields and 100 correlated LRG overdensity fields, with N_side = 256 for each value of f_NL. Lastly, we applied to the mocks the same masks we used for the LRG and CMB lensing data.

In order to test the performance of the systematics mitigation pipeline, we contaminated the LRG maps. For this purpose, after generating the mock LRG and CMB lensing fields, we used the regressis⁷ code (Chaussidon et al. 2021). From the LRG density fields, we simulated LRG discrete number count maps using Poisson sampling, based on the expected LRG density for each pixel. This procedure was safe enough due to the low map resolution, which made it unlikely we would get pixels with predicted density lower than 0 for the Gaussian fields. We also checked the angular power spectra were in agreement with the input after Poisson sampling. We then generated “high density” simulated LRG maps by duplicating the number of objects per pixel. Then, we created a mock LRG catalog by assigning random coordinates within a given pixel to each galaxy. Finally, the “high density” catalog was used as input for regressis. Based on the systematics weights, which were estimated from running SYSnet on the real data, the code created a contaminated mock catalog by removing galaxies and matching the final catalog to the expected LRG density. The result was converted back into a contaminated map.

We note that we did not add any contamination to the CMB lensing mock fields. This was equivalent to assuming there is no correlation in systematics between the two probes, and the main purpose of this test was to check whether the systematics mitigation removed real cross-correlation signal in the C_ℓ^κG power spectra.

We show in Fig. 6 an example of a clean mock map and the same mock map contaminated after applying regressis. These contaminated maps were used as input for SYSNet, the neural network code for systematics mitigation described in Sect. 4.1. For each realization, we applied SYSNet separately in three different regions of the sky: one corresponding to the BASS and MzLS footprints and two corresponding to DECaLS in the NGC and SGC (see Fig. 2 of Rezaie et al. 2024). After recombining the output in a single map, we computed the angular power spectra of the lensing – LRG cross-correlation, C_ℓ^κG, and the LRG autocorrelation, C_ℓ^GG, for the mocks without the contamination (clean power spectrum), then we added the contamination and finally applied the systematics weights obtained with SYSNet.

Fig. 6. Single Gaussian realization of an LRG map before adding contamination (left panel) and after applying regressis to contaminate the mock (right panel).

The comparison of angular power spectra is shown in Fig. 7. For C_ℓ^κG, we find no significant differences between the contaminated and true power spectra. While the mitigated spectra are systematically lower than the truth, this difference is smaller than the dispersion of the clean mock power spectra (i.e., the $C_{\ell }^{\kappa g}$ uncertainty). In relative terms, at the lowest multipoles where the f_NL signal arises, we find differences between the mitigated and true spectra of up to ∼15%, ∼25%, and ∼30–35% for the f_NL = −50, 0, and 50 mocks, respectively. Nonetheless, constraining f_NL from C_ℓ^κG using uncorrected maps would lead to a much larger $C_{\ell }^{\kappa g}$ covariance from the large amount of extra power in C_ℓ^GG, significantly increasing the f_NL uncertainty. For this reason, we used the mitigated maps as the baseline in our cross-correlation pipeline.

Fig. 7. Mean angular power spectra of the full sample LRG autocorrelation (left panels) and LRG-CMB lensing cross-correlation (right panels), computed from the 100 mock realizations used to test our pipeline for various values of fNL (−50 top panels, 0 middle panels, and 50 bottom panels). The blue lines show the true angular power spectra before adding contamination, the shaded blue area corresponds to the dispersion on the true spectra, the orange lines correspond to the angular power spectra of the contaminated mocks, and the green lines to the contaminated mocks after imaging systematics mitigation.

For C_ℓ^GG, the systematics mitigation presents a more complicated picture: we find that for the f_NL = 50 mocks the mitigated spectra are compatible with the dispersion of the true spectra, but the performance of the mitigation pipeline depends on f_NL: for f_NL = −50 and f_NL = 0, the first multipole bin of the mitigated spectra is clearly still too high, lying outside the 1σ range of the clean mocks. Thus, for C_ℓ^GG we imposed a scale cut of the first five multipoles contained in this bin and adopted ℓ_min = 7 for the parameter inference with the mocks and real data. Having a larger ℓ_min for C_ℓ^GG than for C_ℓ^κG also allowed us to put a joint constraint on f_NL, while having better control of the systematics.

We then applied the MCMC parameter inference pipeline on the mean of the clean and mitigated angular power spectra of the 100 mock realizations for each value of f_NL. We analyzed C_ℓ^κG and C_ℓ^GG separately and jointly in order to understand the impact of each observable on the constraints, as well as the effects of possible remaining systematics. We list the median likelihood values with 68% confidence intervals in Table 1 for the mocks without systematics (clean) and in Table 2 for the mocks with systematics (mitigated). For the clean mocks, we find almost no intrinsic bias on the recovered f_NL values, always having a ≲0.3σ agreement with the input. For the constraints from C_ℓ^κG on mocks after contamination and systematics mitigation, the input f_NL is always within the 1σ error bars of the measurements on mocks, and the bias on the measured f_NL is also lower than 0.5σ for the f_NL = 0 and f_NL = −50 mocks. The larger disagreement between the measured and true values for the f_NL = 50 mocks can be interpreted as a consequence of an overfitted power spectrum after applying the systematics mitigation pipeline (with the bottom right panel of Fig. 7 showing the largest offset between mitigated and clean in the lowest ℓ bin, in units of the standard deviation). For C_ℓ^GG, we also find a level of agreement between the measured and true f_NL values within 1σ; however, in this case, the highest bias on the measurement happens for the f_NL = −50 mocks. For the joint C_ℓ^κG + C_ℓ^GG tests, using ℓ_min = 7 for C_ℓ^GG, we find a good agreement for the zero and positive f_NL cases, and a ≲1σ bias on the f_NL = −50 mocks. Additionally, due to the systematics uncertainties on C_ℓ^GG, we measured the constraints from mocks using a joint approach in which we kept the C_ℓ^GG information for ℓ > 32 only, in such a way that it constrained the galaxy bias but did not affect the scales where the f_NL signal arose.

The results of the tests on mocks with systematics reported in Table 2 are driven by the impact of the systematics mitigation on the angular power spectra shown in Fig. 7: for the negative f_NL mocks, there is a quite limited (∼15%) overfit on C_ℓ^κG and some remaining power for C_ℓ^GG, while for positive f_NL mocks there is a stronger overfit (≳30%) on C_ℓ^κG and the true power spectrum is almost perfectly recovered for C_ℓ^GG. In consequence, for the negative f_NL mocks, the constraints from C_ℓ^κG present almost no bias with respect to the true f_NL, while the constraints from C_ℓ^GG are biased by ≲1σ toward higher f_NL values, likely due to remaining systematics. The opposite behavior is found for the positive f_NL mocks: the constraints from C_ℓ^κG are biased by ≲1σ, while the constraints from C_ℓ^GG present a better accuracy. For the f_NL = 0 mocks, we find almost no bias for C_ℓ^GG and a 0.5σ bias for C_ℓ^κG, also driven by the overfit of systematics mitigation. In the case of the joint constraints, when adopting the baseline ℓ_min = 7 for C_ℓ^GG, the constraints behave in a similar way to those from the autocorrelation, performing accurately for zero and positive f_NL. Our results show that, even if there could be possible biases on f_NL due to remaining or overcorrected systematics, the input f_NL is within the 1σ uncertainties of the measurement for every single case.

6. Results

We list in Table 3 the median likelihood values with 68% confidence intervals for f_NL and b₀, obtained from the MCMC analysis applied to C_ℓ^κG and C_ℓ^GG separately and jointly. We obtain an uncertainty σ(f_NL)≲ 40 when using the C_ℓ^κG cross-correlation only. Using a scale cut at ℓ_min = 7, this uncertainty is reduced to σ(f_NL)∼ 25 from C_ℓ^GG only and to σ(f_NL)∼ 20 from the combination of both observables. The median likelihood values of f_NL suggest a preference of this LRG sample for positive values of f_NL, which could be interpreted as a consequence of remaining systematics. However, taking into account the error bars, our results are consistent with a ΛCDM universe with Gaussian initial conditions at a ∼1σ level for C_ℓ^κG and ≳1σ when combining this information with the LRG autospectra, C_ℓ^GG. For the galaxy bias parameter, b₀, we obtain a compatible result with previous measurements using this LRG sample: for example, Zhou et al. (2021) measured b₀ ∼ 1.5 assuming a redshift evolution proportional to D(z)⁻¹. Fig. 8 shows the 68% and 95% confidence ellipses for the f_NL – b₀ plane obtained from the two observables and their combination.

Fig. 8. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation (grey contours), the CℓGG autocorrelation (red contours), and both observables jointly (blue contours).

We assumed as a baseline fiducial cosmology the σ₈ ∼ 0.77 measurement from the cross-correlation between ACT lensing and this LRG sample (Sailer et al. 2024; Kim et al. 2024). As a robustness test, we recomputed the constraints assuming other σ₈ values in the literature: the Planck 2018 best-fit value for this parameter, σ₈ = 0.811 (Planck Collaboration VI 2020), as well as the DESI full-shape analysis determination, σ₈ = 0.842 (DESI Collaboration 2024d) and the measurement from BOSS × Planck CMB lensing, σ₈ = 0.688 (Chen et al. 2024). We show in Fig. 9 the comparison between the C_ℓ^κG constraints on the f_NL – b₀ plane assuming the baseline and the alternative values from the literature for σ₈. Assuming the Planck 2018, DESI 2024, and BOSS × LRG σ₈ best-fit values, we find $f_{\mathrm{NL}}={42}_{-40}^{+44}$, $f_{\mathrm{NL}}={43}_{-47}^{+49}$, and f_NL = 40 ± 32, respectively. The median likelihood f_NL value is unaffected in every case with respect to the σ₈ = 0.77 case listed in Table 3, while the 68% confidence interval is ∼10% larger for the Planck 2018 case (σ₈ ∼ 0.81), ∼20% larger for the DESI 2024 case (σ₈ ∼ 0.84), and ∼20% smaller for the BOSS × LRG case (σ₈ ∼ 0.69). The trend shows larger uncertainties on f_NL as a consequence of assuming larger σ₈ values. This is essentially due to the degeneracy between the galaxy bias and σ₈: higher σ₈ values result in lower b₀ measurements, and this affects the amplitude of the (b − p) term that enhances the f_NL signal (see Eq. (2)).

Fig. 9. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The green contours assume the baseline σ8 value preferred by the LRG sample used in this paper, while the grey, red, and blue contours assume other σ8 determinations in the literature using Planck 2018, DESI 2024, and BOSS × Planck data, respectively.

We also tested the robustness of our results to the choice of the fiducial galaxy bias redshift evolution assumed in Eq. (21). For this, we assumed a model in which the galaxy bias evolution was given by b_g = b₀ × [D(z)]^a, and we left the a parameter free, together with f_NL, fixing b₀ = 1.51 to the best fit in Table 3. We find $f_{\mathrm{NL}}={40}_{-36}^{+38}$ and a = −1.01 ± 0.09, obtaining full consistency with the baseline results and the fiducial galaxy bias evolution.

Our baseline analysis assumed the universality relation for the scale-dependent galaxy bias induced by f_NL, with a value p = 1 in Eq. (2). However, works on simulations, such as Barreira (2020), suggest a more realistic value for this parameter could be p = 0.55. We recomputed the constraints on f_NL from C_ℓ^κG, assuming p = 0.55, and compared them to the baseline p = 1 case in Fig. 10. When assuming p = 0.55, we find $f_{\mathrm{NL}}={29}_{-27}^{+25}$. The median likelihood value of f_NL is shifted as predicted by Barreira (2020), but the compatibility with ΛCDM remains at the ∼1σ level. Further, we also computed the constraints for a model-independent case in which, instead of assuming the universality relation where b_ϕ = 2δ_c(b_g − p), no specific form for b_ϕ is assumed and we constrained the product b_ϕf_NL instead of f_NL. In this case, we find $b_{\varphi }{f}_{\mathrm{NL}}={146}_{-142}^{+154}$, a measurement that is still compatible with ΛCDM at the ∼1σ level.

Fig. 10. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The red contours correspond to the baseline constraints assuming p = 1, while the blue contours assume p = 0.55, as suggested in Barreira (2020).

In order to understand the stability of our f_NL measurement from C_ℓ^κG in terms of systematics, we compare in Fig. 11 the constraints from the Planck lensing – DESI LRG cross-correlation with the case in which no systematics mitigation was performed. In the latter case, we find $f_{\mathrm{NL}}={64}_{-73}^{+80}$. The result is compatible with the baseline case using systematics weights, but the uncertainties are larger due to the impact of the extra C_ℓ^GG power spectrum in the covariance matrix. This also emphasizes the importance of performing an accurate systematics mitigation, even for cross-correlation analysis.

Fig. 11. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The blue contours correspond to the baseline constraints after applying the SYSnet mitigation weights to the data and the red contours to the case in which no systematics mitigation was performed.

We also tested the robustness of our results to the systematics looking at the stability of the results as a function of the minimum multipole, ℓ_min. In Fig. 12 we represent the constraints on f_NL from C_ℓ^κG and C_ℓ^GG for ℓ_min = 2, 7, and 12. For C_ℓ^κG, the constraints are compatible at a ≲1σ level with f_NL = 0 for the three ℓ_min values. Instead, for C_ℓ^GG, we find a ∼1σ deviation for ℓ_min = 7, while this deviation rises to ∼2σ for ℓ_min = 2 and ℓ_min = 12. This indicates that the lower multipoles on $C_{\ell }^{\mathit{GG}}$ might be dominated by observational systematics, as suggested by the test on mocks, while for C_ℓ^κG it is safer to include all the scales in the analysis.

Fig. 12. Constraints on fNL from CℓκG (orange dots) and CℓGG (blue dots) with 68% confidence interval as a function of the minimum multipole ℓmin included in the analysis.

In addition, we tested the compatibility of our results among different regions of the sky. We computed our constraints on f_NL by applying masks of the NGC and SGC and list the results in Table 4. We find compatible results between the two galactic caps, and even if the best-fit f_NL values we obtain from C_ℓ^κG are consistently higher than those measured from the full footprint, the measurements are compatible with zero PNG, considering the larger error bars.

Our constraint on f_NL from the DESI LRG – Planck lensing cross-correlation is also compatible with the main result from Rezaie et al. (2024) using the LRG autocorrelation only. Rezaie et al. (2024) discuss in detail the impact of the mitigation recipe on the f_NL error bars and best-fit values. They highlight $f_{\mathrm{NL}}={34}_{-44}^{+24}$ when applying a more aggressive nonlinear regression for systematics mitigation using nine features maps, while we selected a more conservative treatment in order to avoid a strong overfit in the cross-correlation, $C_{\ell }^{\kappa G}$. When applying the same nonlinear treatment with three maps that we chose for this study, they find $f_{\mathrm{NL}}={28}_{-11}^{+12}$ (before calibration) and f_NL = 47 ± 14 (after calibration). We note that, in Rezaie et al. (2024), a calibration of the constraints is proposed, based on the results of the tests on mocks, and that all multipoles were included in the baseline analysis, while we did not apply any calibration and we set ℓ_min = 7 for $C_{\ell }^{\mathit{GG}}$. Another important difference in our analysis is the contamination model on the mocks: we used nonlinear weights to implement the contamination on the mocks, while in Rezaie et al. (2024), linear weights were applied to generate this contamination. A calibration of the f_NL constraints based on the mocks would be model dependent in any case, and our aim is instead to show the capabilities of $C_{\ell }^{\kappa G}$ to produce more stable measurements of f_NL, either alone or in combination with $C_{\ell }^{\mathit{GG}}$.

7. Conclusions

In this work, we used an LRG catalog from the DESI Legacy imaging surveys, calibrated with the spectroscopic redshifts that have been observed for the DESI Survey Validation, in combination with the Planck PR4 CMB lensing map, to put a constraint on the primordial local non-Gaussianity parameter, f_NL. We measured f_NL through the scale-dependent bias effect using as observables the cross-correlation between the LRG and CMB lensing maps in the angular domain, C_ℓ^κG, and the autocorrelation of the LRG field, C_ℓ^GG.

In order to limit the impact of imaging systematics on large scales where the f_NL signal is present, we used a neural network code for imaging systematics mitigation. Our measurement was performed without blinding, but the full analysis methodology was tested on mock fields including imaging systematics for different f_NL values. Our end-to-end pipeline works at a reasonable level of agreement with the input f_NL values, being especially robust when combining both C_ℓ^κG and C_ℓ^GG for the positive and zero f_NL mocks. From the systematics tests on mocks, we find residual contamination in the behavior of the five first multipoles (ℓ = 2 − 6) of C_ℓ^GG; hence, we imposed a cut at ℓ_min = 7 for the autocorrelation power spectrum.

From the CMB – LRG cross-correlation, C_ℓ^κG, alone we find $f_{\mathrm{NL}}={39}_{-38}^{+40}$ at the 68% confidence level. We performed some robustness tests on this result, such as changing the fiducial σ₈ to other values in the literature, different from the DESI LRG best fit, changing the value of the p parameter, modifying the fiducial bias evolution, or varying the minimum multipole, ℓ_min, to evaluate the possible impact of systematics, and we find consistent results for every test. If we combine our result with the information from the angular LRG autocorrelation, C_ℓ^GG, adopting the ℓ_min = 7 cut, we find $f_{\mathrm{NL}}={24}_{-21}^{+20}$, although C_ℓ^GG presents a larger statistical fluctuation as a function of ℓ_min, suggesting that C_ℓ^κG is more stable and less sensitive to the effect of imaging systematics. Our results are consistent with the f_NL measurements from the DESI LRG by Rezaie et al. (2024) and motivate the future use of CMB cross-correlations, together with the tracers’ autocorrelations and an appropriate systematics mitigation technique, in order to perform accurate f_NL measurements.

Data availability

All data points, maps, covariances, and MCMC chains shown in the figures of this paper are publicly available in https://doi.org/10.5281/zenodo.14401463

Acknowledgments

J. R. Bermejo-Climent and R. Demina acknowledge support from the U.S. Department of Energy under the grant DE-SC0008475.0. JRBC acknowledges the support of the Spanish Ministry of Science and Innovation under the grants PID2021-126616NB-I00 and “DarkMaps” PID2022-142142NB-I00, and from the European Union through the grant “UNDARK” of the Widening participation and spreading excellence program (project number 101159929). AK was supported as a CITA National Fellow by the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference #DIS-2022-568580. MR is supported by the U.S. Department of Energy grants DE-SC0021165 and DE-SC0011840. This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science, Innovation and Universities of Spain (MICIU/AEI/10.13039/501100011033), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. The DESI Legacy Imaging Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall z-band Legacy Survey (MzLS). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory. Legacy Surveys also uses data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Legacy Surveys was supported by: the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility; the U.S. National Science Foundation, Division of Astronomical Sciences; the National Astronomical Observatories of China, the Chinese Academy of Sciences and the Chinese National Natural Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. The complete acknowledgments can be found at https://www.legacysurvey.org/. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies. The authors are honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

References

All Tables

All Figures

	Fig. 1. Normalized redshift distribution of the LRG sample, directly measured using the spectroscopic redshifts from DESI Y1 data.
In the text

	Fig. 2. Left panel: LRG overdensity field in galactic coordinates, after applying the δ < −30° cut. Right panel: CMB lensing field in galactic coordinates, obtained from the Planck PR4 maps.
In the text

	Fig. 3. Angular power spectra of the LRG autocorrelation, CℓGG, and CMB lensing – LRG cross-correlation, CℓκG, obtained from the raw (uncorrected) data and after applying the systematics mitigation pipeline to the LRG maps.
In the text

	Fig. 4. Relative impact of the lensing magnification contributions for the CℓGG (blue line) and CℓκG (red line) theoretical angular power spectra obtained with CAMB.
In the text

	Fig. 5. Correlation matrix as obtained from the joint covariance matrix for CℓGG and CℓκG computed with the NaMaster code.
In the text

	Fig. 6. Single Gaussian realization of an LRG map before adding contamination (left panel) and after applying regressis to contaminate the mock (right panel).
In the text

Fig. 7. Mean angular power spectra of the full sample LRG autocorrelation (left panels) and LRG-CMB lensing cross-correlation (right panels), computed from the 100 mock realizations used to test our pipeline for various values of fNL (−50 top panels, 0 middle panels, and 50 bottom panels). The blue lines show the true angular power spectra before adding contamination, the shaded blue area corresponds to the dispersion on the true spectra, the orange lines correspond to the angular power spectra of the contaminated mocks, and the green lines to the contaminated mocks after imaging systematics mitigation.

In the text

	Fig. 8. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation (grey contours), the CℓGG autocorrelation (red contours), and both observables jointly (blue contours).
In the text

	Fig. 9. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The green contours assume the baseline σ8 value preferred by the LRG sample used in this paper, while the grey, red, and blue contours assume other σ8 determinations in the literature using Planck 2018, DESI 2024, and BOSS × Planck data, respectively.
In the text

	Fig. 10. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The red contours correspond to the baseline constraints assuming p = 1, while the blue contours assume p = 0.55, as suggested in Barreira (2020).
In the text

	Fig. 11. 1σ and 2σ confidence ellipses for the joint posterior distribution of the fNL and b0 parameters obtained from the CℓκG cross-correlation. The blue contours correspond to the baseline constraints after applying the SYSnet mitigation weights to the data and the red contours to the case in which no systematics mitigation was performed.
In the text

	Fig. 12. Constraints on fNL from CℓκG (orange dots) and CℓGG (blue dots) with 68% confidence interval as a function of the minimum multipole ℓmin included in the analysis.
In the text