Prepared for submission to JCAP

Constraining primordial
non-Gaussianity from DESI DR1
arXiv:2512.17865v1 [[Link]] 19 Dec 2025

quasars and Planck PR4 CMB

Lensing

S. Chiarenza 1,2,3 A. Krolewski1,2,3 M. Bonici1,2,3

E. Chaussidon 4 R. de Belsunce 4 W. J. Percival 1,2,3
J. Aguilar4 S. Ahlen 6 A. Baleato Lizancos 4,5 D. Bianchi 7,8
D. Brooks9 T. Claybaugh4 A. Cuceu 4 K. S. Dawson 10 A. de
la Macorra 11 P. Doel9 S. Ferraro 4,5 A. Font-Ribera 13
J. E. Forero-Romero 14,15 E. Gaztañaga 16,17,18 S. Gontcho A
Gontcho 4,19 G. Gutierrez20 H. K. Herrera-Alcantar 21,23
K. Honscheid 23,24,25 D. Huterer 26,27 M. Ishak 28 R. Joyce 29
D. Kirkby 30 A. Kremin 4 O. Lahav9 C. Lamman 25
M. Landriau 4 L. Le Guillou 31 M. E. Levi 4 M. Manera 32,13
P. Martini 23,25,33 A. Meisner 29 R. Miquel13,34 S. Nadathur 17
J. A. Newman 35 G. Niz 36,37 N. Palanque-Delabrouille 22,4
C. Poppett4,38,5 F. Prada 39 I. Pérez-Ràfols 40 G. Rossi41
E. Sanchez 42 D. Schlegel4 M. Schubnell26,27 H. Seo 43
J. Silber 4 D. Sprayberry29 G. Tarlé 27 B. A. Weaver29
C. Yèche 22 R. Zhou 4 H. Zou 44
1 Department of Physics and Astronomy, University of Waterloo, 200 University Ave W,
Waterloo, ON N2L 3G1, Canada
2 Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5,

Canada
3 Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo,

ON N2L 3G1, Canada

For other aliations, see Appendix C.
 E-mail: schiaren@[Link]

Abstract. We present the rst measurement of local-type primordial non-Gaussianity from

the cross-correlation between 12 million spectroscopically conrmed quasars from the rst
data release (DR1) of the Dark Energy Spectroscopic Instrument (DESI) and the Planck PR4
CMB lensing reconstructions. The analysis is performed in three tomographic redshift bins
covering 08 < z < 35, covering a sky fraction of ∼ 20%. We adopt a catalog-based pseudo-
Cℓ estimator and apply linear imaging weights validated on noiseless mocks. Compared to
previous analyses using photometric quasar samples, our results benet from the high purity
of the DESI spectroscopic sample, the reduced noise of PR4 lensing, and the absence of excess
large-scale power in the spectroscopic quasar auto-correlation. Fitting simultaneously for the
non-Gaussianity parameter fNL and the linear bias amplitude in each redshift bin, we obtain
fNL = 2+28 +20
−34 for a response parameter p = 16, and fNL = 6−24 for p = 10. These results
improve the constraints on fNL by ∼ 35% compared to the previous analysis based on the
Legacy Imaging Survey DR9. Our results demonstrate the statistical power of DESI quasars
for probing inationary physics, and highlight the promise of future DESI data releases.

Keywords: cosmological parameters from LSS – power spectrum – CMB – galaxy clustering
Contents

1 Introduction 1

2 Theory 3

3 Data 7
3.1 DESI DR1 Spectroscopically conrmed quasars 7
3.2 Planck PR4 CMB lensing map 9

4 Measuring Angular Power Spectra 9

4.1 Monte Carlo normalization correction for the Planck lensing maps 11

5 Imaging systematic weights validation 11

5.1 The DESI systematic weights 11
5.2 Correcting for excess mode removal 12
5.2.1 Creating the mocks 12
5.2.2 Test on systematics-free mocks 15

6 Optimal Redshift Weighting for fNL 16

7 Results and Discussion 19

7.1 The optimal weights results 22
7.2 Fitting a single bias amplitude for all three bins 23
7.3 Impact of the bias relation 24
7.4 The impact of the primordial power spectrum 25

8 Conclusions 27

A Derivation of the optimal fNL weights for Cℓκg 37

B Restricting to quasars with z < 31 38

C Affiliations 39

1 Introduction

Ination is the leading framework for explaining the emergence of cosmic structure. A key
prediction of this theory is that the initial curvature perturbations are almost perfectly Gaus-
sian [1–3]. This prediction is supported by precise measurements of the Cosmic Microwave
Background (CMB) temperature and polarization [4]. Still, many ination models allow
small departures from Gaussianity, constraining these deviations, known as primordial non-
Gaussianity (PNG), is a powerful direct test of early-Universe physics. For example, in
single-eld slow-roll (SFSR) ination the expected level of local PNG is small, fNLloc ∼ O(10−2 )

[5, 6]. Multi-eld scenarios, on the other hand, can naturally reach fNL ∼ O(1) [7]. Achieving
loc

order-unity precision on fNLloc is thus an important goal in ruling out signicant parts of the

parameter space of ination.

–1–
 The best current limit comes from the Planck bispectrum, fNL = −09 ± 51 [8]. This is
already close to the cosmic-variance limit for the CMB, and future CMB experiments will be
able to improve it by only a factor of two [9]. Stronger constraints must therefore come from
large-scale structure (LSS).
Local PNG introduces a scale-dependent bias in the clustering of galaxies [10], which
manifests as an enhancement of the quasars (QSOs) power spectrum, P (k), on the very
large scales. The eect scales as ∝ k −2 and grows with redshift, making quasars particularly
powerful tracers: their high-redshift distribution provides access to the large cosmic volumes
and low-k modes where the signal is most pronounced. In fact, the best constraint to date
from LSS comes from the 3D power spectrum analysis done on the combination of DESI DR1
QSO and LRG samples, yielding fNL = −36+9.0 −9.1 [11]. However, measuring galaxy clustering
on very large scales is inherently challenging. The signal uctuations are small in amplitude,
and systematic eects can introduce signicant spurious power, particularly on the largest
scales [12–19]. As a result, all large-scale structure constraints on fNL to date have been
limited by systematic uncertainties rather than fundamental statistical noise [20].
Cross-correlations between galaxy surveys and CMB lensing provide a robust alternative
to galaxy auto-correlations for probing primordial non-Gaussianity. The main advantage
comes from the fact that the dominant sources of noise and systematics in galaxy surveys are
generally uncorrelated with those in CMB experiments [21–24]. This property makes cross-
correlations especially powerful for controlling large-scale systematics and isolating the true
cosmological signal [24, 25]. In fact, while systematics can bias the galaxy auto-correlation
signal, in the cross-correlation they primarily act to increase the statistical variance without
introducing a systematic shift. We note that, in principle, correlated systematics could exist:
galactic emission can bias lensing reconstruction and is correlated with extinction, but no
signicant eects have been detected so far. Moreover, the CMB lensing reconstruction can
suer from a potential contamination from extragalactic foregrounds, the main one being the
Cosmic Infrared Background (CIB) and thermal Sunyaev-Zel’dovich (tSZ) eect, particularly
when the reconstruction is derived from temperature data [26]. However, such eects are
expected to be small for our measurement, and no signicant bias has been reported in
similar cross-correlation analyses [27]. As a supporting test, we performed the measurement
using the Planck polarization-only lensing reconstruction (available within the PR4 data
release), which is strongly unaected by CIB and tSZ contamination. While the polarization-
only data yield much larger statistical uncertainties, the recovered signal shows no signicant
deviation from our baseline result. Another benet of CMB lensing cross-correlations is their
sensitivity to high redshifts, where the fNL signature is strongest. The lensing kernel for
CMB lensing peaks at z ≈ 2, overlapping well with the redshift distribution of quasars and
other high-redshift tracers. There have been several recent applications of this idea [28–30].
In particular the constraint from the Quaia photometric quasar sample cross-correlated with
Planck PR4 lensing [31] yielded σfNL ≈ 25 [29]. The reported constraint using DESI LRGs
is of σfNL ≈ 40 [28]. Finally, the analysis of the photometric DESI Legacy Survey quasar
sample cross-correlated with Planck PR4 lensing presented in [30], which forms the basis
of the approach used in this work, obtained σfNL ≈ 45. Looking ahead, forecasts suggest
that ongoing spectroscopic surveys like DESI, and future Stage 4 surveys like SPHEREx and
LSST, combined with CMB lensing, are advancing toward the era of σfNL < 1 [32–39].
Driven by this motivation, we present the rst measurement of local-type non-Gaussianity
from the cross-correlation between DESI DR1 quasars and Planck PR4 CMB lensing maps
[4]. Our analysis focuses exclusively on the cross-correlation to minimize contamination from

–2–
auto-correlation systematics, providing a clean and complementary test of the standard LSS
analysis which employs the 3D power spectrum of high redshift tracers. The analysis builds
upon that by [30]: DESI DR1 spectroscopically conrmed quasars are employed instead of
the photometric targets in the DESI Legacy Survey [40], so the sample is aected by fewer
imaging systematics. Planck PR4 lensing maps are employed, which have 10 − 20% less noise
than the Planck 2018 maps [41] due to use of additional data and more optimal ltering and
analysis methods [4, 31]. Our constraining power now reaches σfNL ∼ 20, showing that the
purer quasar sample and the new pixel-free pipeline improve the constraining power by a
factor  2 over [30], who nd σfNL ∼ 45 despite the ∼ 2× larger sky fraction.
In Section 2, we describe the theory necessary to compute angular correlation functions
at low ℓ. In Section 3, we describe the quasar and CMB lensing data, and in Section 4,
we describe our angular power spectrum pipeline. In Section 5, we validate our pipeline on
mocks, to verify that we do not over-correct when mitigating the eect of imaging systematics.
In Sec. 6 the implementation of optimal redshift weights to maximize the fNL signal in the
cross-correlation is discussed. Finally, in Section 7, we present the results, and in Section 8 we
compare them to previous fNL constraints. Throughout this paper, we x the cosmological
parameters to the Planck 2018 at ΛCDM model [42] with h = 06766, As = 2105 × 10−9 ,
ns = 09665, Ωm = 03096, Ωb = 0049, one neutrino with mass 006 eV, and σ8 = 08102,
but we also test the impact of freeing the primordial power spectrum (i.e., As and ns ) on the
constraints we obtain on fNL .

2 Theory
Local PNG is characterized by the dimensionless parameter fNL :
 
Φ(x) = Φ(x) + fNL Φ2 (x) − ⟨Φ2 ⟩ , (2.1)
where Φ is a Gaussian random eld representing the primordial potential [43]. This local
transformation introduces mode-coupling between long and short wavelengths, which man-
ifests as a scale-dependent bias in large-scale structure tracers [10, 12, 44]. The physical
origin of PNG-induced scale-dependent bias lies in the modulation of halo formation by long-
wavelength potential uctuations. Taking the Laplacian of Eq. 2.1, we see that the presence
of PNG increases the density in the peaks of the density eld:
δNG = δG (1 + 2fNL ϕG ) (2.2)
This shifts the eective collapse threshold:
δc → δc (1 − 2fNL ϕG ) (2.3)
And, following [10], one can show that the resulting scale-dependent correction to the linear
bias b1 (z) is:
bΦ (z)
b(k, z) = b1 (z) + fNL , (2.4)
TΦ→δ (k, z)
where b1 (z) is the linear bias of the tracer, bΦ is the PNG bias, giving the response to the
presence of local PNG of the tracer, and the transfer function TΦ→δ (k, z) is the transfer
function between the primordial gravitational eld ΦG and the matter density perturbation,
computed as:
  ns −1
Plin (k, z) 9 2π 2 k
TΦ→δ (k, z) = with PΦ (k) = As , (2.5)
PΦ (k) 25 k 3 kpivot

–3–
where PΦ (k) is the primordial potential power spectrum, ns is the spectral index, and As the
amplitude of the initial power spectrum at kpivot = 005 Mpc−1 . Hence, through the Poisson
equation, TΦ→δ (k, z) has the well-known scale dependence:

TΦ→δ (k, z) ∝ k 2 TΦ→Φ (k, z) (2.6)

where TΦ→Φ (k, z) is the usual total matter transfer function. Since fNL is always paired with
bΦ , a measurement of scale-dependent bias measures the product bΦ fNL . Nonetheless, this
degeneracy does not aect the signicance of a potential non-zero detection, which would
remain a robust indication of the presence of PNG [45, 46]. The PNG bias bΦ quanties the
logarithmic response of the galaxy number density to a change in amplitude of the matter
clustering. It is dened as:
∂ log n̄
bΦ =  (2.7)
∂ log σ8
However, the theoretical modeling of bΦ is a widely discussed topic [45–48] that goes beyond
the scope of this present work. Therefore, we will follow the standard prescription, also known
as universality relation [49]:
bΦ (z) = 2δc (b1 (z) − p) (2.8)
where δc ≈ 1686 is the critical spherical collapse density threshold and b1 (z) is the redshift-
dependent linear bias. The most typical choices for the response parameter p is 1 for a mass-
selected sample, while for a sample dominated by recent mergers, such as quasars, p = 16 is
a more appropriate choice [49, 50]. We will report constraints for both prescriptions.
In this work, we probe ∆b(k) using the matter-galaxy cross-power spectrum, Pgm (k).
Specically, since CMB lensing is a 2-dimensional projected eld, our observable is the angular
cross-power spectrum:
  
κδ 2 Wκ (z2 )
Cℓ = dz1 Wδ (z1 ) dz2 2 dk Pmm (k, z1 , z2 )jℓ (kχ(z1 ))jℓ (kχ(z2 )) (2.9)
π χ (z2 )

where Wδ (z) is the number counts kernel:

H(z)
Wδ (z) = b1 (z)n(z), (2.10)
c
n(z) being the quasar target redshift distribution, and Wκ (z) is the CMB lensing kernel
 
3 H02 χ(z)
Wκ (z) = Ωm χ(z)(1 + z) 1 − (2.11)
2 c2 χ(z⋆ )

with χ(z⋆ ) the comoving distance to the surface of last scattering, at z⋆ = 1090. When
including the scale-dependent bias induced by fNL , we replace the linear galaxy bias b(z)
with a scale- and redshift-dependent term b(k, z), as dened in Eq. 2.4. This assumes a
linear galaxy bias, valid over the scales used here (4 < ℓ < 300, i.e. k  01 h Mpc−1 for
08 < z < 35).
In accordance with other PNG analysis within DESI, we x the redshift evolution of
the bias to follow the functional form of [11], which captures the expected increase of quasar
bias with redshift. In practice, this redshift dependence is integrated over in the model,
and therefore modies the eective redshift at which the measurement is interpreted. We
introduce a free normalization parameter bi0 for each tomographic bin, which preserves the

–4–
redshift evolution but allows the overall amplitude to vary, resulting in discontinuities between
bins. In Sec. 7.3 we demonstrate that fNL depends only on the eective bias amplitude of
each bin, as the constraints remain consistent across all alternative bias evolution models
considered.
Eq. 2.9 is a numerically challenging three-dimensional integral over oscillatory spherical
Bessel functions. The Limber approximation, typically employed to simplify the problem,
approximates the spherical Bessel functions as Dirac delta functions in their rst peak [51]:

π
jℓ (kχ) → δD (ℓ + 12 − kχ) (2.12)
2ℓ + 1
The three-dimensional integral is reduced to a one-dimensional integral in k or χ, but this
approximation is not valid on the very large angular scales, which happen to contain the
most fNL information. Therefore, we evaluate our theory model using Blast1 [52, 53], an
algorithm for calculating angular power spectra without employing the Limber approximation
or assuming a scale independent growth rate, based on the use of Chebyshev polynomials.
The code assumes that the power spectrum can be factorized as

P (k) = Plin (k) + (Pnl − Plin )(k) (2.13)

The non-linear part is negligible until ℓ > 200, where the Limber approximation works very
well. Therefore, for multipoles ℓ < 200, we compute the signal using the full non-Limber ex-
pression evaluated using Plin (k), and then add the non-linear correction (Pnl − Plin ) evaluated
with Limber approximation. For ℓ > 200, Blast uses the Limber approximation with Pnl (k).
The modeling also includes the contributions from lensing magnication and redshift-
space distortions. Lensing magnication accounts for the changes on the background density
caused by foreground structures. This alters the observed number counts and is eectively
correlated with CMB lensing, leading to an extra contribution to the cross-correlation:
  
κµ 2 Wκ (z1 ) Wµ (z2 )
Cℓ = dz1 2 dz2 2 dk Pmm (k, z1 , z2 )jℓ (kχ(z1 ))jℓ (kχ(z2 )) (2.14)
π χ (z1 ) χ (z2 )
where the magnication bias kernel Wµ (z) is
 +∞  
3 H02 ′ ′ ′ χ(z)
Wµ (z) = Ωm χ(z)(1 + z) dz n(z )(5s(z ) − 2) 1 − (2.15)
2 c2 z χ(z ′ )

The inner integral runs over the distribution of galaxies, and s ≡ d log
dm
10 n
[54] is the response
of the number density n to achromatic changes in the brightness dm. In Sec. 3, we describe
how s is estimated for our sample.
Redshift space distortions also contribute to the observed number counts by causing
galaxies to be observed in a dierent redshift shell due to their peculiar velocities:
  
2
CℓκRSD
= dz1 Wκ (z1 ) dz2 WRSD (z2 ) dkPmm (k, z1 , z2 )jℓ (kχ(z1 ))jℓ′′ (kχ(z2 )) (2.16)
π
The RSD kernel can be expressed as:
H(z)
WRSD (z) = f (z)n(z) (2.17)
c
1
[Link]

–5–
 ×10−7

3.5 100

3.0
80
2.5

100 · Cκi/Cκg
2.0 60 Number counts
Magnification bias
Cκg

1.5 RSD
40 PNG
1.0

0.5 20

Total (s = 0.099, b0 = 1, fNL = 0)

0.0
Total, fNL = 50 0
0 50 100 150 200 250 300 101 102
 

Figure 1: Left: Total quasar-CMB lensing cross-correlation in the rst tomographic bin
(08 < z < 21). The ducial model with fNL = 0 is presented in blue, and with fNL = 50 in
magenta. Right: Contributions of each term to Eq. 2.22 as a fraction of the ducial model
with fNL = 50. Negative terms are shown as dashed lines.

where f (z) is the logarithmic derivative of the growth rate with respect to scale factor,
f ≡ ddln D
ln a , and the spherical Bessel function is replaced with its second derivative. The only
sensitivity to fNL is through Cℓκδ , as neither Cℓκµ nor CℓκRSD depend on the galaxy bias. As
discussed, fNL leaves an imprint in the clustering signal through the scale dependent bias,
dened in Eq. 2.4. In the presence of PNG, the number counts kernel, dened in Eq. 2.10,
becomes:
H(z)
Wδ+PNG (z) = (b1 (z) + ∆b(k, z)) n(z) (2.18)
c
This can be factorized:
Wδ+PNG (z) = Wδ (z) + WPNG (z), (2.19)
with
H(z)
WPNG (z) =
fNL bΦ (z)n(z) (2.20)
c
Keeping in mind that Blast performs the inner k-integral rst, and that this step is aected
by the presence of the transfer function in the scale-dependent bias, it can be worked out that
the PNG contribution to the cross-correlation is:
  
κPNG 2 Wκ (z1 )
Cℓ = dz1 2 dz2 WPNG (z2 ) dk PΦ (k) TΦ→δ (k, z1 ) jℓ (kχ1 ) jℓ (kχ2 ), (2.21)
π χ (z1 )

where the matter power spectrum has been written as Pmm (k, z1 , z2 ) = PΦ (k) TΦ→δ (k, z1 ) TΦ→δ (k, z2 ).
The reason why Eq. 2.21 contains only one transfer function is that the scale-dependent bias
scales as ∆b(k, z) ∝ 1TΦ→δ (k, z), which cancels the second transfer-function factor from
Pmm . Eectively, we can write our nal model as:

Cℓκg = Cℓκδ + Cℓκµ + CℓκRSD + CℓκPNG (2.22)

–6–
with CℓκPNG ∝ fNL . For a more technical description of the algorithm, see [53].
Fig. 1 shows the fractional contributions of magnication, RSD and PNG to the total
Cℓκg signal with fNL = 50 in the rst tomographic bin (08 < z < 21). The ducial values of
bi0 = 1 and s = 0099 are assumed. The magnication term is a considerable fraction (∼ 15%)
of the clustering term, and rises in a scale-dependent way that is approximately degenerate
with fNL at ℓ > 30, emphasizing the importance of low multipoles. Despite this degeneracy,
the magnication bias slope is measured suciently accurately that it does not worsen the
fNL constraints. The RSD term is very subdominant, only rising to 1% of the ducial model
at ℓ < 6. Finally, the fNL contribution dominates for ℓ < 10 and constitutes a non negligible
part of the signal up to ℓ ≈ 100. Although fNL = 50 is a very high and unrealistic value for
this parameter, it shows how a non-Limber evaluation of the theory model is crucial to probe
the PNG signal.

-0.60 0.52 0.00 0.26 0.00 3.51

(a) Planck PR4 lensing map. (b) DESI DR1 quasar sample. (c) Completeness mask.

Figure 2: Data used for the cross-correlation (same as in Ref. [55]). Panel (a) shows the
Planck PR4 lensing convergence map, κ, together with its smoothed mask, obtained by
applying a Gaussian lter with 1◦ FWHM. Panel (b) displays the DESI DR1 quasar number
counts for the full sample spanning the redshift interval 08 ≤ z ≤ 35, while the corresponding
completeness mask is shown in panel (c). For visualization, we adopt a HEALPix resolution of
Nside = 128, whereas all computations are performed at Nside = 2048. All maps are displayed
in a Mollweide projection and are presented in the Galactic coordinate system.

3 Data

We cross-correlate quasars from the DESI survey with lensing mass maps obtained by the
Planck satellite, namely the PR4 convergence maps [4]. The data are visualized in Fig. 2.
In Sec. 3.1 we summarize the quasar sample and in Sec. 3.2 we briey present the employed
CMB lensing data.

3.1 DESI DR1 Spectroscopically confirmed quasars

Our analysis employs a spectroscopic catalogue of 1 223 391 quasars in the redshift interval
08 ≤ z ≤ 35 taken from the DESI DR1 quasar sample [56] released with the survey’s rst
public data release [57]. DESI is a highly multiplexed ber-fed spectrograph located on the
4m Mayall telescope at the Kitt Peak National Observatory [58, 59]. It has the capability to
record up to 5 000 spectra in a single exposure [60–64]. Throughout a eight-year campaign,
DESI aims to map up to 17 000 deg2 , which will provide precise measurements of baryon
acoustic oscillations (BAO) and redshift-space distortions (RSD) from both galaxies and
quasars across 0  z  35 [60, 66]. Following successful validation of the survey [67] and an
early data release [68], the data from the rst year has already provided competitive BAO

–7–
 0.6 Bin 1: 0.8 < z < 2.1
Bin 2: 2.1 < z < 2.5
0.5
Bin 3: 2.5 < z < 3.5
0.4
n(z)

0.3

0.2

0.1

0.0
1.0 1.5 2.0 2.5 3.0 3.5
z

Figure 3: Normalized redshift distribution n(z) of the DESI DR1 spectroscopic quasar
sample. The dierent colors identify the three non-overlapping redshift bins.

PR4 × DR1
Label z-range Nqso Shot Noise z ze sµ fsky
g1 08 ≤ z < 21 856 831 26 × 10−6 1.49 1.44 0.099 19.8
g2 21 ≤ z < 25 194 754 112 × 10−6 2.28 2.27 0.185 18.9
g3 25 ≤ z ≤ 35 171 806 127 × 10−6 2.85 2.75 0.244 18.6

Table 1: Overview of DESI DR1 quasar samples used for the cross-correlation with Planck
CMB lensing. We list the three redshift bins {g1 , g2 , g3 }, their redshift ranges, quasar counts,
shot-noise levels, mean redshift z, eective redshift ze [65], magnication bias sµ , and the
PR4×DR1
sky-fraction overlap with the PR4 CMB-lensing convergence maps, fsky (in percent).

constraints [69–71], accurate RSD measurements [72], and the currently strongest large-scale
structure limit on primordial non-Gaussianity, derived from the quasar power spectrum [11].
Table 1 summarizes the DR1 quasar subsamples, while Fig. 3 illustrates their redshift
distributions alongside the tomographic binning scheme used in this analysis. This scheme is
consistent with the accompanying study of DESI DR1 quasars provided in [55]. The sample
is divided into three redshift bins to retain radial information. The rst bin corresponds to
the DESI ducial range, 08 < z < 21, while the remaining quasars, extending up to z = 35,
are split into two approximately equally populated bins, 21 < z < 25 and 25 < z < 35.
We utilize the clustering catalogues described in [73, 74] together with redshifts produced by
the DESI spectroscopic pipeline [75–78].
Table 1 also reports the measured values of magnication bias for the three tomographic
bins. We measure the magnication bias slope, s ≡ d log10 Ndm, by perturbing the Legacy
Imaging Survey DR9 photometry [40] uniformly by ±005 mag, re-running quasar target se-
lection, and measuring the resulting change in number density. The values obtained from
brightening and dimming are consistent, so we adopt their mean as our ducial slope. Since
quasars are dominated by point sources and selected using total magnitudes rather than
ber magnitudes, it is an excellent approximation to treat lensing magnication as a simple
brightening or dimming of their ux. This procedure measures the slope using the observed

–8–
(post-magnication) magnitude distribution, but the induced bias is negligible, as magni-
cation only scatters a small fraction of quasars across the ux limit. We nd that the slope
varies with redshift, and therefore use the value appropriate for each bin: s = 0099 for
08 < z < 21, s = 0185 for 21 < z < 25, and s = 0244 for 25 < z < 35.

3.2 Planck PR4 CMB lensing map

Gravitational lensing by large–scale structure deects CMB photons, imprinting subtle dis-
tortions in the temperature and polarization elds that can be exploited to reconstruct the
lensing potential [see, e.g., 79, for a review]. For this work we employ the Planck PR4 CMB
lensing convergence maps, κ, which are publicly available2 [4]. Starting from the global
minimum–variance (GMV) spherical-harmonic coecients κℓm , we resample the data onto a
HEALPix grid with NSIDE=2048 [80].
The PR4 release benets from the NPIPE processing pipeline [4], which provides im-
proved, uniformly processed CMB maps that serve as inputs to the lensing reconstruction.
The reconstruction itself was upgraded in PR4: the collaboration adopted a GMV estimator
that incorporates temperature–polarization correlations [81] and implemented an anisotropic
ltering and local-noise weighting scheme [31]. Combined with the ∼ 8% additional CMB
observations obtained during satellite repointing, these improvements yield an overall ∼ 20%
increase in signal–to–noise. As a result, the PR4 lensing map remains signal–dominated up to
multipoles ℓ ≃ 70, compared with ℓ ≃ 40 for PR3. The eective reconstruction noise power
spectrum, Nℓκκ , supplied by the Planck collaboration, is incorporated in the analytical covari-
ance matrix, described in Sec. 4. The lensing mask retains a sky fraction of fsky = 671%3 .
Its overlap with the tomographic DESI DR1 quasar catalogue is summarized in Table 1.

4 Measuring Angular Power Spectra

We measure clustering directly from the un-pixelised DESI quasar catalogue using the Na-
Master4 implementation [82] of the catalog-based pseudo–Cℓ estimator [82, 83], which builds
upon the algorithm rst described in [84] and relies on the DUCC library5 , used, i.e., in [85]. In
the standard pixelized approach, galaxies and randoms are rst binned onto a high–resolution
HEALPix grid, and the overdensity eld is formed as 1 + δgal by dividing the two maps.
While conceptually straightforward, this division can lead to numerical instabilities in pixels
where the survey completeness is very low. In early tests with this pixel-based procedure, we
observed a signicant excess of power on large scales in the measured angular power spectra,
which could not be accounted for by known imaging systematic eects. By carefully tracking
the origin of this excess, we found it was caused by instabilities in the overdensity maps
calculation in low-completeness pixels. Switching to the catalog-based pseudo–Cℓ estimator
completely resolved this issue, yielding robust measurements across all scales. We analyze
angular scales in the range 4 ≤ ℓ ≤ 300. The binning scheme uses ∆ℓ = 10 for ℓ ≤ 100, where
most of the signal resides, and ∆ℓ = 50 for ℓ ≥ 100, yielding 14 bins in total. Finally, we
stress that the NaMaster estimator automatically subtracts the Poisson shot noise from the
measured auto–power spectra. We have veried that the shot–noise levels estimated internally
2
[Link] and [Link]
3 Npix
The sky fraction is defined as fsky = i=1 tot
xi /Npix , where Npix
tot 2
= 12 Nside ≈ 5 × 107 and xi is the value
of the mask in pixel i.
4
[Link]
5
[Link]

–9–
by the code agree with the theoretical expectations obtained from the number density in each
redshift bin (Table 1).
Furthermore, the pixelized approach is subject to aliasing and pixel–window function
eects. Mitigating these problems usually requires ultra-ne grids (with an arbitrary choice
of NSIDE) and very large random catalogues, driving up both memory and computational
costs. The catalog-based pipeline instead adopts the “FKP” approach from three-dimensional

MC normalization correction
1.15
Correction
1.12
Mean value = 1.08
1.10
C

1.08

1.05

1.03

1.00
50 100 150 200 250


Figure 4: Monte Carlo normalization correction applied to the cross-spectra Cℓκgi which
amounts to an approximately constant amplitude change of ≈ 8%.

analyses by working with the linear dierence ng − α nr instead of the ratio of data over
randoms catalogs. Because this operation is linear, we can simply subtract the direct spher-
ical–harmonic transforms of data and scaled randoms, avoiding numerical instabilities and
pixel-related problems. In the limit that the objects are point-like, the spherical harmonic
transform of these elds is simply a sum of spherical harmonics evaluated at the positions
of the objects, Yℓm (θi , ϕi ), making the evaluation of the estimator straightforward. For a
detailed description of the algorithms, the reader is referred to [83, 84].
Pixel-free estimators oer a promising approach for future analyses involving cross-
correlations with discrete tracers, however the corresponding analytic covariance calculations
for catalog-based elds are still under development [86]. Therefore, following the approach
used, for example, in [29, 87], we adopt the standard analytic Gaussian, window-convolved
covariance approximation implemented in the gaussian_covariance method of the NaMaster
package to estimate the errors in our analysis. Since this covariance routine requires a pix-
elized map of the data eld, we construct such a map from the data aℓm ’s obtained with
the catalog-based estimator, generating HEALPix maps at NSIDE = 2048. Moreover, this ap-
proach also requires some theory spectra, which, as discussed in Sec. 2, we generate using
g g
Blast [52, 53]. The code gives predictions for Cℓκgi and Cℓ i j (where we add the shot noise
gi gi
contribution as in Table 1 to Cℓ ) in the ducial cosmology. For the lensing power spectrum
Cℓκκ we use the publicly available PR4 lensing spectra with the provided lensing reconstruc-
tion noise Nκκ . As in [55], we use an iterative approach: we t the cosmological and nuisance
parameters with the covariance matrix evaluated with the theory curves just described, and
use the resulting best-t parameters to recalculate the ducial spectra and tune the covariance
matrices.

– 10 –
4.1 Monte Carlo normalization correction for the Planck lensing maps
A well-known eect in CMB lensing reconstruction is a misnormalization of the reconstructed
eld caused by the masks and anisotropic ltering applied in order to perform the reconstruc-
tion. The mode coupling introduced in this procedure is not modeled into the cross-spectra
obtained with the NaMaster algorithm. As a result, the reconstructed convergence κ̂ does not
have the same normalization as the true convergence κ, and a correction is required [88–90].
In practice, this correction is determined through simulations: the basic idea is to com-
pare the cross-correlation between the input convergence (appropriately masked) and the
reconstruction with the known input power spectrum. This yields a Monte Carlo normaliza-
tion factor that is dependent on the mask and the footprint:
κ κin,g−mask
Cℓ in,κ−mask
AMC
ℓ = κ̂ κin,g−mask
, (4.1)
Cℓ

where κ̂ is the masked CMB lensing reconstruction, κin,κ−mask is the input lensing convergence
masked with the lensing mask, and κin,g−mask is the input convergence masked using the
galaxy mask. An unbiased estimate of the CMB lensing cross-spectrum is then obtained as

Cℓκ̂MC g = Cℓκ̂g AMC

ℓ  (4.2)

Following [55], we implement this correction in a simulation-based, mode-by-mode fashion.

Specically, we generate a suite of CMB lensing reconstructions using the appropriate mask:
 Nsim ℓ  κ i   q i ∗
ℓ WLℓ i m=−ℓ M κ ℓm M κ ℓm
TL ≡  Nsim ℓ , (4.3)
κ i q i ∗
ℓ WLℓ i m=−ℓ {M κ̂ }ℓm {M κ }ℓm

where M denotes the masks, κ  the input (and κ̂ the reconstructed) convergence, i indexes
the simulations, and XY ℓm ≡ d n̂Yℓm2 ∗ (n̂)X(n̂)Y (n̂) with Y
ℓm the spherical harmonics. For
the Planck PR4 maps over the DESI DR1 footprint, this correction is at the level of ≈ 8%
on the scales of interest (see Fig. 4).

5 Imaging systematic weights validation

Spatial variations in imaging quality and foregrounds imprint spurious modulations in the
angular number density of DESI targets. If left uncorrected, these modulations contaminate
the largest-scale Fourier modes that are most sensitive to primordial non-Gaussianity.

5.1 The DESI systematic weights

DESI adopts the template-tting approach introduced in recent large-scale structure surveys.
The technique derives a per-object weight, wsys , by quantifying how the target surface den-
sity varies with each imaging attribute and then correcting for those dependencies. More
specically, DESI attaches three multiplicative weights to every object:

wtot = wcomp wsys wzfail , (5.1)

ensuring that the ratio of weighted data to weighted randoms remains at across the footprint
once all known observational selection eects are removed. Completeness weights wcomp
correct for the probability that a target is allocated a ber, and redshift failure weights,

– 11 –
wzfail , mitigate spatial uctuations in redshift success produced by variations in exposure
time, focal-plane position, and hardware status. The imaging systematics weights wsys , the
focus of this section, remove spurious target-density uctuations that track imaging conditions
such as depth, seeing, dust, and stellar density. Full implementation details for the rst two
weights are given in [74].
Below, we discuss only the imaging systematic weights, since this is the element we test
extensively for its impact on the quasar angular power spectrum. The other two weights,
completeness and redshift-failure, are not a concern for our analysis. Systematic weights have
many possible choices, including dierent templates and implementations, and can have a
large impact on the results. By contrast, it is well understood how to correctly account for
completeness. Redshift-failure weights have a much smaller impact than the systematic errors
and therefore do not signicantly aect the analysis [74, 91]. The DESI catalog in the latest
released version (v1.5), comes with three possible sets of systematic weights: the neural net
method sysnet, denoted WEIGHT_SN in the catalogs; the random forest method regressis,
denoted WEIGHT_RF; the linear method applied to eBOSS [91] LSS catalogs, identied as
WEIGHT_IMLIN. In the present work, we only consider two dierent sets of linear weights as
it is known that both machine learning based weights, namely the sysnet neural network
[92] and the regressis6 random forest weights [93], over-correct for excess large scale power
introduced by systematics, especially in the case of the quasar sample [11, 74, 93]. The second
set of linear weights is also computed with regressis at HEALPix map level using a linear
regression algorithm. Note that, compared to the catalog linear weights, the linear regression
here does not t the data to binned statistics but rather the uctuation at HEALPix map level.
The rst are identied as WEIGHT_IMLIN, the latter as WEIGHT_Linear.

5.2 Correcting for excess mode removal

As described above, we test two sets of linear systematic weights that are derived from the
same imaging templates but dier in their implementation. Accurately assessing the perfor-
mance of such weights is crucial, as they can over-correct and remove large-scale clustering
signal, thereby weakening the fNL constraints. The correction for this mode removal by
reweighting modes is related to the angular integral constraint (AIC) correction [73]. This
eect has been shown to be particularly severe for neural network–based weighting schemes
[11], which we therefore do not consider in this analysis. Nevertheless, it remains important
to quantify any residual bias for the two linear weighting schemes tested here. To quantify the
impact of our systematic–mitigation weighting scheme, we use custom-generated mocks with
no imaging systematic contamination, so that tracer densities are, on average, uncorrelated
with the imaging features. Any deviation in the angular power spectra after applying the
systematic weights, will be due to the mitigation method itself and must be corrected on the
measurement on the data.

5.2.1 Creating the mocks

We have created 500 uncontaminated, catalog-based lognormal mocks that accurately repro-
duce the DESI DR1 footprint and clustering properties. This number of realizations is su-
cient to ensure that the uncertainty associated with the mock-based correction remains well
below the statistical errors across the relevant scales. The rst step is to generate correlated
Gaussian maps for the three redshift bins used in the analysis (08 < z < 21, 21 < z < 25,
6
[Link]

– 12 –
25 < z < 35) and for the CMB lensing map. This is done using the healpy package [94],
specically the synfast function, which synthesizes maps at a chosen resolution from input
angular power spectra Cℓ . The theoretical Cℓ are computed with the GLASS7 package [95],
which allows the generation of lognormal elds with prescribed two-point statistics. Because
the lognormal transformation, dened as:
 
f (X; λ) = λ eX − 1 , (5.2)

is nonlinear and modies the underlying statistics of the eld, GLASS determines, through an
iterative algorithm, the Gaussian input power spectrum Gℓ that will yield the target Cℓ after
applying the transformation. In our case, we set the shift-parameter λ to 1 and transform
the three Gaussian elds representing the tomographic bins. The target angular spectra Cℓ
are computed with Blast [52, 53] and the bias is adjusted to match the data, ensuring that
the nal mocks reproduce the same large-scale clustering statistics as the data. We convert

×10−7 Mocks Cκg bin 1

Theory
4
Mocks

2
C

−2

−4
0 50 100 150 200 250


Figure 5: Left: A zoom into the survey footprint. In darker blue, the DESI DR1 quasars,
while in light blue are the galaxies making up one of the 500 mocks. The agreement between
the two footprint is very good, conrming the success of our procedure to replicate the DESI
survey geometry. Right: Average angular power spectrum Cℓκg over 500 mocks with the
corresponding 1σ error bars. As desired, the mocks have a power spectrum that matches the
input theory curve.

the pixelized lognormal maps into mock catalogs by Poisson-sampling the number of galaxies
in each pixel according to its overdensity. Each mock galaxy is then assigned sky coordinates
(RA, DEC) by randomly perturbing the position around the pixel center using an eective
pixel radius θpix = Apix 2. During this process, the maps are normalized so that the
resulting mock catalogs contain approximately the same number of objects per redshift bin
as the data. As this process is very approximate on small scales, we generate the maps at
a very high resolution (NSIDE = 8192) to ensure that, on the relevant scales (ℓ < 300), the
resulting angular clustering would not be aected by this sampling procedure.
7
[Link]

– 13 –
 The nal step to obtain realistic mocks is to apply angular cuts and downsampling to
reproduce the survey footprint and completeness. The DESI DR1 footprint exhibits substan-
tial small-scale structure, as the survey completeness remains low at this stage. To reproduce
these variations as accurately as possible, we perform four successive downsampling steps, de-
signed to account for the main observational eects impacting the real survey, while avoiding
the computational cost of running each mock through the fiberassign software [96] and de-
tailed veto masks. In particular, we select areas of the sky covered during DR1, remove areas
around bright stars and regions aected by bad imaging systematic conditions. Finally, we
account for the survey completeness by constructing completeness maps directly from existing
survey random catalogs. Specically, we use the ocial DESI DR1 quasar random les, which
encode the spatial completeness pattern of the survey, together with the Legacy Survey DR9
randoms [40], which are uniform over the sky. Both sets of randoms are pixelized at NSIDE
= 2048, and the ratio of the two maps provides an estimate of the angular completeness.
Each mock catalog is then downsampled according to this completeness map, ensuring that
the mocks reproduce the same spatial variations in target completeness as the data. This
procedure is correct up to the resolution of the maps, which is enough for our analysis, re-
stricted to scales ℓ < 300. The results of these footprint and clustering validations are shown
in Fig. 5. The left panel demonstrates that the mocks accurately reproduce the DESI DR1
survey geometry, as evidenced by their excellent agreement with the ocial quasar randoms.
This conrms that the completeness correction procedure described above eectively transfers
the angular selection function to the mock catalogs. The right panel shows that the aver-
age cross-correlation Cℓκg measured from 500 mocks is consistent with the input theoretical
prediction, within the expected statistical uncertainties. Together, these tests verify that our
catalog-based lognormal mocks correctly reproduce both the survey footprint and the target
clustering statistics.

×10−6 Cgg bin 1: AIC test

1.2
4

1.0
2
C/Cunw
C

0.8
unweighted
0
regressis 0.6 regressis
catalog catalog
0.4
0 50 100 150 200 250 300 0 50 100 150 200 250 300
 

Figure 6: Left: Average angular auto power spectrum Cℓgg of 500 mocks with dierent
linear weighting schemes, in the rst redshift bin considered in the analysis (08 < z < 21).
Right: Ratio of the weighted to unweighted power spectra. The ratio deviates from 1 on large
scales, showing that the weights are removing signal from the analysis. The shaded grey area
highlights the 1% region around unity.

– 14 –
 ×10−7 Cκg bin 1: AIC test
1.2
4
1.0
2

C/Cunw
C

0.8
0
unweighted
−2 regressis 0.6 regressis
catalog catalog
−4 0.4
0 50 100 150 200 250 300 0 50 100 150 200 250 300
 

Figure 7: Left: Average angular cross-power spectrum Cℓκg of 500 mocks with dierent linear
weighting schemes (pixel-based, labeled as “regressis,” or object-based, labeled as “catalog”), in
the rst redshift bin considered in the analysis (08 < z < 21). Right: Ratio of the weighted
to unweighted power spectra. In this case, no signicant deviation from 1 is observed. The
shaded grey area highlights the 1% region around unity.

5.2.2 Test on systematics-free mocks

To evaluate the signicance of mode removal using either the linear systematic correction
weight or those creating using the linear regressor within the regressis package, we measure
the angular power spectra of the mocks using each set of weights and compare them to the
unweighted case (where all weights are set to unity). As the mocks do not have systematic
contamination, a suppression of large-scale power in the weighted cases would indicate that
the weights are spuriously removing cosmological signal instead of systematics, requiring a
corresponding correction. We examine how the systematic weights impact both the galaxy
auto-spectrum, Cℓgg , which is most sensitive to observational systematics and contributes to
the covariance, and the cross-spectrum with CMB lensing, Cℓκg , which is used in our main
analysis. Fig. 6 shows the results for the galaxy autocorrelation Cℓgg , indicating that both
sets of weights are removing power on very large scales (ℓ < 20). We use the ratio, displayed
on the right side of the plot, as a multiplicative correction to the power spectra measured on
the data. The impact is more pronounced for the “catalog” linear weights, and the underlying
cause is currently under investigation, as it must reside in the implementation of the weights.
Fig. 8 shows the corrected autocorrelation measurement: we note that, even after accounting
for angular mode removal, the auto power spectrum does not show signs of excess large scale
power, a huge improvement with respect to the rst analysis using the DESI DR9 Legacy
Survey [30]. This is can be attributed to the fact that the DESI DR1 spectroscopic catalog
has a much higher purity compared to the photometric targets in the Legacy Survey and
to the employment of the new catalog based estimator, which is more stable and robust to
numerical instabilities and pixel related issues, as discussed in Sec. 4. In this analysis we
will adopt the more conservative and robust approach of using only the cross-correlation to
constrain fNL . Nonetheless, this result is encouraging as there is no need to model excess
large scale power in the covariance, as in the previous analysis. Moreover, it indicates that
we are going in the direction of being able to safely use the auto-correlation information for
the fNL analysis in future DESI data releases.

– 15 –
 ×10−6 Cgg bin 1 with AIC correction
5 Theory
catalog
4 regressis

3
C

0 50 100 150 200 250 300



Figure 8: Auto-correlation angular power spectra measured from the catalog with the two
possible sets of weights and corrected for mode removal. The error bars are computed from an
analytical gaussian covariance computed using the NaMaster package, as discussed in Sec 4.

Fig. 7 shows no sign of over-subtraction of angular modes in the CMB lensing cross-
correlation, which we employ in this work. Hence, we decided to not correct the measurement,
an example of which is shown in Fig. 9.

6 Optimal Redshift Weighting for fNL

In the present work, we tested the use of optimal redshift weights for the angular two-point
statistics, with the aim of improving constraining power of the analysis. The idea of opti-
mal redshift weighting was rst developed in [97, 98] and is now routinely applied in three-
dimensional power spectrum analyses. For reference, the optimal fNL weight for the galaxy
monopole is given by [97, 99, 100]:
 
2 f (z)
wP̂0 (z) = wFKP bΦ (z) D(z) b(z) + , (6.1)
3

where wFKP is the standard FKP weight [101]. Since this is a weight for the power spectrum,
the corresponding per-galaxy weight is its square root.
In practice, this corresponds to cross-correlating
 two weighted elds, respectively propor-
f (z)
tional to wFKP (b(z) − p) and wFKP b(z) + 3 , as in [100]. This weighting scheme enhances
the constraining power on fNL by explicitly accounting for the redshift dependence of the
PNG signal. We extend this framework to the angular cross-correlation between CMB lens-
ing and galaxy overdensity, Cℓκg , and derive an optimal estimator for this observable. The
details of the derivation are given in App. A. We nd that the optimal redshift weight that

– 16 –
 ×10−7 Cκg bin 1
Theory
5
catalog
regressis
4

3
C

−1
0 50 100 150 200 250 300


Figure 9: Cross-correlation angular power spectra measured from the catalog with the two
possible sets of weights. In this case, no correction is applied as there is no sign of angular
mode removal introduced in the cross correlation. The error bars come from the analytic
estimation of the covariance matrix performed using the NaMaster package, described in
Sec. 4.

Bin 1 Bin 2 Bin 3

3.0 Fiducial Fiducial Fiducial
Optimal weights Optimal weights Optimal weights

2.0
n(z)

1.0

0.0
1.0 1.5 2.0 2.0 2.2 2.4 2.6 2.5 2.8 3.0 3.2 3.5
z z z

Figure 10: Redshift distribution of the DESI DR1 quasars. The plot shows the n(z) of the
three tomographic bins in which the data is divided for this analysis, in order of increasing
redshift. The solid lines represent the ducial n(z), and the dashed lines show the redshift
distribution after weighting the galaxies with the weights in Eq. 6.2. The displayed redshift
distributions are normalized in each bin.

maximizes the fNL sensitivity of Cℓκg is:

Wκ (χ(z))
wopt (z) = bΦ (z) D(z) wFKP , (6.2)
χ2 (z)

normalized such that wopt (z) dz = 1. As in the power spectrum case, the resulting weights

– 17 –
eectively up-weight higher-redshift galaxies, where the fNL signature is strongest. As one
might expect, we only get one power of the FKP weight here, as opposed to the two powers
in the power spectrum case.

Measured Cκg bin 1

×10−7
7
Theory
Theory with optimal weights
6
Default
Default with optimal weights
5

3
C

−1

0 50 100 150 200 250 300



Figure 11: Cℓκg in the rst tomographic bin. In navy, the measurement and theory curve
in the default conguration (i.e., using the regressis linear weights, see Sec. 5). In light
green, the case where the optimal weights are applied. This plot showcases the impact of the
optimal weights on the pipeline.

The implementation of these weights is straightforward. We evaluate the optimal weights

for each galaxy in the catalog using the expressions derived above. The CMB lensing kernel
Wκ is dened in Eq. 2.11, while bΦ is given in Eq. 2.8, with the linear bias b(z) following the
model of [11]. The growth factor D(z) is computed using the Boltzmann solver CAMB [102],
and the FKP weights [101] are already provided as an attribute of the galaxy catalog. In
the default computation of the redshift distribution n(z), each galaxy is weighted according
to Eq. 5.1 before binning. When including the optimal weighting scheme, the total weight
assigned to each galaxy becomes w = wtot wopt = wcomp wsys wzfail wopt , eectively modifying
the inferred n(z). As shown in Fig. 10, the optimal weights up-weight the high-redshift
part of the sample, where the sensitivity to fNL is greatest. This occurs because galaxies at
higher redshift are more aected by a potential PNG signal, which enhances their large scale
clustering. On the other hand, at the lowest redshifts the weights become negative. In the
rst redshift bin (08 < z < 09), some galaxies receive negative weights. This is expected and
does not pose any issue for the analysis: the negative values simply indicate that galaxies in
this range would respond to a stronger PNG signal with a reduction in clustering amplitude.
Measuring fNL from such an absence of clustering is intrinsically more dicult, and the
optimal weighting appropriately accounts for this eect.
To apply the optimal weights consistently throughout the analysis, we must consider all
components of the pipeline that are aected. A change in the galaxy redshift distribution
modies the theoretical predictions, since n(z) enters the number counts kernel Wδ (χ) dened
in Eq. 2.10. Consequently, the theory curves used in the analytic covariance computation (see

– 18 –
Sec. 4) must also be updated. The analysis performed with optimal weights therefore employs
a covariance matrix recomputed following the same procedure described in Sec. 4, but using
the weighted n(z). Finally, the weighting scheme must also be applied to the data when
measuring the angular power spectra. Fig. 11 illustrates the eect of the optimal weights on
the measured Cℓκg and on the corresponding theoretical predictions. The comparison between
the default conguration (dened in Sec. 5) and the weighted case highlights how both the
measurement and the theoretical model respond to the modied redshift distribution.

7 Results and Discussion

In this section we present tomographic constraints on primordial non-Gaussianity obtained

from the cross-correlation between the Planck PR4 CMB lensing convergence maps and spec-
troscopically conrmed quasars from DESI DR1. Our analysis covers the redshift range
08 ≤ z ≤ 35, divided into three tomographic bins, and includes a total of 1,223,391 quasars
across 7,200 deg2 of sky. We use the angular power spectrum estimator described in Sec. 4
to measure the cross-spectra Cℓκg . Our parameter inference is performed with Turing.jl8
[103], which provides an interface to dene probabilistic models in terms of explicit priors
and our dierentiable [Link] likelihood. Posterior exploration is carried out using the
No-U-Turn Sampler (NUTS) algorithm [104], a self-tuning Hamiltonian Monte Carlo sampler
that eciently leverages gradient information and has been widely adopted for cosmological
inference. The main results is presented in Table 2, visualized in Fig. 12 and discussed in the
following text. Throughout this section we also present results for various analysis variations
and tests that we performed, like the inclusion of optimal fNL weights and dierent redshift
cuts and bins. Everything is summarized in Fig. 17, showcasing the robustness of the anal-
ysis. As discussed in Sec. 1, uncertainties in the theoretical prediction of the non-Gaussian

×10−7 Cκg bin 1 ×10−6 Cκg bin 2 ×10−6 Cκg bin 3

0.5 regressis weights
2
4 catalog weights
Best fit (p = 1.6)
0.0
C

2 1

−0.5
0 0

−1.0
101 102 101 102 101 102
  

Figure 12: As in Fig. 9, the scatter points represent the measurements on the data per-
formed with two possible weighting schemes: the default choice regressis in green, and the
alternative linear weighting scheme implemented in the DESI catalog in purple. The dashed
lines represent the corresponding best t curves when assuming p = 16.

bias coecient bΦ create a degeneracy between fNL and bΦ . In general, one could constrain
the product fNL bΦ directly, without assuming a specic relation between bΦ and b1 [see e.g.
45], since a detection of the product would remain physically meaningful. Our goal here,
however, is not to explore this broader parameterization, but to assess how the inferred fNL
shifts under reasonable choices for the bΦ (b1 ) relation. We therefore consider two benchmark
8
[Link]

– 19 –
cases commonly used in the literature [see e.g. 11, 29, 30]: p = 16, corresponding to a recent-
merger scenario, and p = 10, corresponding to the universality assumption. We impose at
priors on all parameters, with fNL ∼ U (−200, 200) and bi ∼ U (0, 2).
To assess the robustness of our results, we performed a series of analysis variations
and validation tests, described below. Fig. 12 presents the measured angular power spectra
together with the best-t model for the two sets of systematic weights considered. In the
third tomographic bin, the lowest multipole (ℓ = 65) shows a deviation of ≈ 2σ from the
best-t prediction. This feature is consistent with statistical noise and does not dominate
the total χ2 . However, we know that the third bin is the one where the quasar sample has
lower purity, so it is possible that some remaining systematic in the CMB lensing maps could
correlate with uncorrected systematics in the sample resulting in the large scale extra power.
For this reason, we run the analysis in the default conguration removing that point, which
is located at scales very relevant to constrain fNL . The result in this case is fNL = −22+27
−34 as
opposed to fNL = 2+28 −34 when all the points are considered. The error bars stay consistent, but
there is a ≈ 08σ shift in the best t value. We decided to keep the analysis as it is and not
exclude that point from the pipeline as the evidence for doing that is not compelling enough
to avoid falling into conrmation bias.
The marginalized constraints corresponding to each conguration are summarized in
Table 2. The best-t values are consistent within 1σ across all weighting schemes, and the
reduced χ2 values conrm the overall quality of the ts. If we use p = 1 instead of p = 16,
the constraints become tighter: this is to be expected as, in that case, the quasars are more
sensitive to variation in fNL . The same results can be visualized in Fig. 13, which shows the
contours for the 4 free parameters. The constraints on the bias scaling amplitudes bi0 are
not aected by variations in the analysis: Table 3 shows the marginalized posteriors for the
bias amplitudes for the two possible choices of systematic weights. Those results have been
obtained with p = 16. Interestingly, the bias amplitude that we t for the third bin comes
out to be ≈ 20 − 30% lower than expected, with a best t value of b30 = 068, very robust to
variations in the analysis. This discrepancy has emerged in other DESI analyses [55], however
the reason is still unknown and is currently being investigated. Nonetheless, since this fact
can raise doubts on the validity of the assumed bias relation [11], we test the impact of that
choice in Sec. 7.3.

p regressis linear weights catalog linear weights

16 fNL = 2+28 2
−34 (χ = 56) fNL = −14+27 2
−33 (χ = 53)
10 fNL = 6+20 2
−24 (χ = 54) fNL = −4+19 2
−22 (χ = 53)

Table 2: Constraints on fNL for dierent choices of p and weighting schemes. The results
from the ducial scenario are reported in the top left corner. The eective number of degrees
of freedom in the analysis is 38, giving reduced χ2 values of 147 and 139 respectively.

To test the impact of the tomographic redshift bins on the constraints, we also performed
the analysis using a single redshift bin 08 < z < 35. The measurement on the data in
this case is displayed in Fig. 14, while the marginalized constraints for fNL are reported in
Table 4. The constraining power in this case deteriorates, indicating how the redshift binning
or the inclusion of optimal weights, as discussed in Sec. 7.1, is important in extracting fNL
information.

– 20 –
 regressis linear weights (default) catalog linear weights
p = 1.6 p = 1.6
p=1 p=1

1.2 1.2
b10

b10
1.0 1.0
1.5 1.5
b20

b20
1.0 1.0

1.0 1.0
b30

b30
0.5 0.5

0 100 1.0 1.2 1.0 1.4 0.5 1.0 −50 50 1.0 1.2 1.0 1.5 0.5 1.0
fNL b10 b20 b30 fNL b10 b20 b30

Figure 13: Corner plots showing the constraints for the 4 free parameters in the analysis,
namely fNL and a bias amplitude for each of the 3 tomographic bins. Left: constraints in
the ducial conguration, using the regressis linear weights. Right: results of the analysis
performed with a dierent choice of systematic weights, i.e. the ones implemented in the
ocial DESI catalog.

regressis linear weights catalog linear weights

b10 113 ± 005 113 ± 005
b20 114 ± 013 117 ± 012
b30 068 ± 012 068 ± 012

Table 3: Bias amplitude constraints for the two possible choices of linear weights tested
throughout this paper. The default choice are the regressis linear weights, but the weight
choice does not have an impact on the bias amplitudes as shown in this table.

p regressis linear weights catalog linear weights

16 fNL = 17+30 2
−40 (χ = 11) fNL = 4+33 2
−41 (χ = 11)
10 fNL = 12+24 2
−28 (χ = 11) fNL = 3+23 2
−27 (χ = 11)

Table 4: Constraints on fNL for the analysis without tomographic bins. Again, dierent
choices of p and weighting schemes have been tested and the results from the ducial scenario
are reported in the top left corner. The eective number of degrees of freedom in the analysis
is 12, giving reduced χ2 values of 091. This corresponds to a p-value of 05, indicating that
the model provides a good t to the data.

– 21 –
 ×10−7 Cκg single bin
Theory
catalog
6 regressis

4
C

0 50 100 150 200 250 300



Figure 14: Cℓκg when the quasars are not divided into tomographic bins.

7.1 The optimal weights results

In existing literature [11, 100, 105] about optimal weighting schemes to maximize the fNL
signal, it was found that the improvement of constraining power coming from this technique
is around 8 − 10% in the case of the 3D power spectrum. However, as the marginalized
constraints in Table 5 show, we found a limited gain (≈ 3%) when using the regressis
weights, and no gain at all when using the alternative choice of linear weights. We also
tested the impact of optimal weights in the single bin analysis, whose results are displayed in
Table 6: while the error bars remain unchanged, the best t values for fNL suer from quite
signicant shifts (from 0 to 20, ≈ 07σ, and from −16 to 10, ≈ 08σ). Despite the substantial
shifts in the posterior means, the mean obtained with optimal weighting falls within the 1σ
range of our default analysis. Nonetheless, these shifts suggest that the data is quite noisy
(as illustrated in Fig. 12) and characterized by large error bars. This noise could potentially
obscure the advantages of employing optimal weights. Furthermore, given that the fNL signal
is dependent on redshift and that we are dividing our dataset into three tomographic bins, it
is possible we are already extracting the majority of the available information.

Optimal weights regressis linear weights catalog linear weights

No fNL = 2+28 2
−34 (χ = 56) fNL = −14+27 2
−33 (χ = 53)
Yes fNL = 20+30 2
−30 (χ = 62) fNL = 10+27 2
−32 (χ = 59)

Table 5: Comparison of the results with and without applying the optimal weighting scheme.
The ducial value of p = 16 is assumed.

To investigate these hypotheses regarding the under-performance of the optimal weights,

we conducted the analysis using a synthetic noiseless data vector created with fNL = 0 and
bi0 = 1. This was done in two dierent setups: the default case with 3 tomographic redshift

– 22 –
 Optimal weights regressis linear weights catalog linear weights
No fNL = 17+30 2
−40 (χ = 14) fNL = 4+33 2
−41 (χ = 14)
Yes fNL = 9+30 2
−38 (χ = 14) fNL = −3+29 2
−36 (χ = 14)

Table 6: Comparison of the results with and without applying the optimal weighting scheme
for the analysis with a single tomographic bin. The ducial value of p = 16 is assumed.

bins, and a simplied case without tomography, which utilized a single redshift bin ranging
from z = 08 to z = 35. We tested the impact of the optimal weights in these two scenarios,
obtaining the constraints reported in Table 7 and visualized in Fig. 15. The results for the
bias amplitudes are not reported because, as expected, the optimal weights for fNL do not
aect the bias constraints, which are completely consistent to that reported in Table 3 for
the analysis without optimal weights.

Optimal weights 3 redshift bins (default) Single bin (08 < z < 35)
No fNL = 0+28
−35 fNL = 0+40
−54
Yes fNL = 0+27
−33 fNL = 0+35
−44

Table 7: Marginalized constraints for fNL when assessing the impact of the tomographic bins
and noisy data on the optimal weights eectiveness.

The results of this test are very interesting: in fact, in the case of the tomographic
analysis, when using noiseless data the optimal weights are providing a ≈ 2% improvement
in the constraining power, similar to what we see in the analysis with real data. On the
other hand, when we employ only one redshift bin, then using optimal weights tightens the
constraints by ≈ 15% in this noiseless data vector test. This suggests that, by binning
our data, we are already extracting most of the redshift information that helps in better
constraining fNL . Overall, the optimal weighting scheme is straightforward to implement
and computationally inexpensive, making it a promising extension of the analysis. However,
since the tomographic approach already achieves comparable performance and the current
measurements are still limited by systematic uncertainties, the practical gain from applying
optimal weights is modest. For this reason, we adopt the results without optimal weighting as
our ducial choice. Looking ahead, the method is expected to become increasingly valuable
for future DESI data releases and other Stage-IV surveys, where a ∼ 10% improvement in
constraining power could play a signicant role in reaching the target precision of σfNL < 1.

7.2 Fitting a single bias amplitude for all three bins

As an alternative to tting a bias amplitude per tomographic bin, we tried to t the same
amplitude b0 to all three bins. The resulting marginalized fNL constraints are presented in the
second row of Table 8. These constraints are compared with the ducial results. The resulting
error bars on fNL are tighter because the analysis has more degrees of freedom (40 instead of
38); however, as indicated by the increase in the reduced χ2 from 14 to 17, indicating that
this t is not as good as the ducial one. This issue arises because the third tomographic
bin has been observed to prefer a much lower amplitude (b30 = 068 compared to b1,2 0 = 11),
making it dicult to accurately represent this with a single amplitude t. In fact the best

– 23 –
 3 redshift bins (default) Single bin (0.8 < z < 3.5)
Default Default
With optimal weights With optimal weights

1.1
b10

1.0
0.9

1.3

1.1
b20

1.0

0.7
1.0

b0
1.3
0.9
b30

1.0

0.7 0.8
0 100 0.9 1.1 0.7 1.0 1.3 0.7 1.0 1.3 −100 0 100 200 0.9 1.0 1.1
fNL b10 b20 b30 fNL b0

Figure 15: Results of the test using a noiseless data vector to assess the performance of the
optimal weights. Left: Analysis in the default conguration, where the dataset is split into
three redshift bins. In this case, the optimal weights (pink contours) do not improve the fNL
constraint relative to the baseline (blue). Right: Results of the analysis performed without
tomography, grouping the full quasar sample into a single bin. Here, the error bars on fNL
are ≈ 10% tighter when applying the optimal weights (pink) compared to the unweighted
case (blue). In both panels, regressis linear weights with p = 16 were used.

Bias regressis linear weights catalog linear weights

3 biases fNL = 2+28 2
−34 (χ = 56) fNL = −14+27 2
−33 (χ = 53)
1 bias fNL = −3+25 2
−28 (χ = 70) fNL = −20 ± 25 (χ2 = 66)

Table 8: Results for the test of tting a single bias amplitude b0 for all three bins, compared
to the baseline conguration where we allow for each redshift bin to have a free amplitude bi0 .
When tting a single bias amplitude the reduced χ2 rises from 147 to 170, indicative of a
worst t. The results reported here were obtained using p = 16 and without employing the
optimal weights.

t values for b0 in this scenario are b0 = 107 ± 005 and b0 = 108 ± 005 for the regressis
and catalog weights respectively. The decision of having a free amplitude per tomographic
bin as our ducial conguration was driven by this goodness-of-t test, and for consistency
with other similar analyses [29, 55]. Nevertheless, the notably low amplitude in the third bin
demonstrates that our baseline bias model does not fully capture the clustering properties of
the quasars, thereby motivating the additional tests presented in Sec. 7.3.

7.3 Impact of the bias relation

The baseline analysis reports some inconsistencies emerging with the third tomographic bin:
rst of all, the lowest multipole ℓ = 65 shows a > 2σ deviation from the best t model

– 24 –
which might be a sign of some remaining uncorrected systematic. Moreover, the marginalized
posterior for the bias amplitude b30 is not compatible with 1 at ≈ 3σ level, indicating that
the bias we assumed is incorrect for high redshift quasars. Despite this discrepancy emerging
in other DESI analysis, it is important to test its impact on the fNL marginalized posterior,
the focus of this work. As a rst test, we tried to completely remove the third bin from
the analysis, performing it with only the rst two tomographic bins (08 < z < 21, and
21 < z < 25). The resulting marginalized constraint is fNL = −15+29−35 . Noticeably, the best
+28
t value shifts by ≈ 05σ compared to the baseline result (fNL = 2−34 ), but the constraining is
almost unchanged, suggesting that the third bin could be removed from the analysis without
any substantial loss of information. As this was discovered later in the analysis, we do not
change our baseline to this case to avoid conrmation bias and to stay consistent with the
analysis presented in [55].
The marginalized constraints for the bias amplitudes bi0 , reported in Table 3, are in-
dicative of the fact that the bias evolution we choose does not represent well the bias of our
QSO sample. In order to better describe our sample, we tested two alternative bias rela-
tions: we nd the best t eective bias of each tomographic bin by evaluating our bias model
at the eective redshift of each bin (those values are reported in Table 1), with bi0 being
the marginalized posteriors (reported in Table 3). The resulting values are: b(ze 1 ) = 246,
2 ) = 377, b(z 3 ) = 279. Starting from these values, we construct a constant bias model:
b(ze e
in each bin, the bias is constant and set to the values just reported. Alternatively, we dene
a redshift dependent bias model b(z) by linearly interpolating between the three values. In

Parameter Constant b(z) Linear Interpolation b(z)

fNL −3+28
−35 2+28
−33
b10 101 ± 005 102 ± 005
b20 112 ± 012 116 ± 013
b30 107 ± 017 116 ± 018

Table 9: Constraints on the parameters obtained from the analyses performed with dierent
prescriptions for the bias model b(z). In both cases, the reduced χ2 of the t is 15

both cases, we still allow for a free bias amplitude in each tomographic bin. Table 9 reports
the marginalized posteriors for the 4 free parameters of the analysis. The fNL constraints are
completely consistent with the baseline choice for the bias, showing deviations in the best
st values < 01σ and the same constraining power, showing how fNL constraints are mostly
insensitive to the choice of bias evolution. As expected for how the alternative bias models
are set up, the bias amplitudes bi0 in those cases are compatible with 1.

7.4 The impact of the primordial power spectrum

The parameter fNL primarily aects the largest cosmic scales, where the observed power spec-
trum is most sensitive to primordial uctuations. Consequently, it can be partially degenerate
with the fundamental parameters of the primordial power spectrum, namely the amplitude
As and the spectral index ns . Allowing these parameters to vary freely improves the prop-
agation of uncertainties and provides a more robust characterization of the constraints on
fNL . In this extended analysis, we introduced As and ns as additional free parameters within
our dierentiable likelihood framework implemented in Blast. This exibility made the joint

– 25 –
 Results with free As and ns

Free As, ns
log(1010As)
Baseline
3.06 Prior
3.00
ns

0.96

1.2
b10

1.0

1.2
b20

0.8
1.0
b30

0.5
0 3.06 0.96 1.0 1.2 0.8 1.4 0.5 1.0
fNL log(1010As) ns b10 b20 b30

Figure 16: Contours for the default analysis, but performed with two extra free parameters:
As and ns .

sampling of all parameters computationally ecient and straightforward. While at priors
were assumed for the standard cosmological and nuisance parameters, we applied Gaussian
priors on As and ns following the Planck 2018 results [42]:

As = (2101 ± 0034) × 10−9 , ns = 09649 ± 00042

The analysis employed the default conguration: regressis linear weights with p = 16, one
bias amplitude per tomographic bin, and no application of the optimal weighting scheme.
The resulting posterior distributions are shown in Fig. 16, and the corresponding parameter
constraints are reported in Table 10. As expected, enlarging the parameter space to include
As and ns slightly degrades the precision on fNL , increasing its uncertainty by approximately
10%. The best-t value, however, remains consistent with that obtained in the ducial
analysis (Table 2). The same trend is observed for the bias parameters: their uncertainties
broaden, but the central values agree well with those reported in Table 3. This conrms
that the inclusion of primordial parameters does not introduce signicant tension, and that
the treatment of degeneracies between fNL and the primordial power spectrum is under good
control.

– 26 –
 Parameter 68% limits

fNL 1+33
−36
As (2098 ± 0032) · 10−9
ns 09648 ± 00040
b10 1101 ± 0074
b20 112 ± 013
b30 066 ± 014

Table 10: Constraints on the parameters obtained in our ducial analysis, but with two
extra free parameters: As and ns .

8 Conclusions

In this work, we have established constraints on primordial non-Gaussianity by cross-correlating

the Planck PR4 CMB lensing maps with quasars from DESI DR1 that have been spectroscop-
ically conrmed. Using scales 4 < ℓ < 300, we measure the parameter fNL , which encodes
information about local-type primordial non-Gaussianity. We focus on the cross-correlation
signal, as it is more robust to systematic errors, whose presence only increases the noise of
the measurement. While the quasar autocorrelation does not exhibit signicant signs of con-
tamination (Fig. 8), such as the large-scale excess power previously observed in the DESI
Legacy Survey photometric quasar sample [30], for this work we restrict the analysis to Cℓκg ,
with plans to incorporate Cℓgg in future data releases. To ensure that the systematic weights
applied to the quasars do not remove true clustering signal, we tested them extensively on
noiseless mocks (Sec. 5).
In our ducial setup, we use the regressis linear weights (Sec. 5), adopt p = 16, and
t simultaneously for fNL and a free bias amplitude bi0 for each of the three tomographic
redshift bins. This analysis yields fNL = 2+28 +20
−34 , while for p = 10, we nd fNL = 6−24  In
comparison to the analysis of photometric quasars from the DESI Legacy Survey [30, 40],
our constraints are ∼ 35% tighter. The improvement can be attributed to three key factors:
the use of a spectroscopic quasar sample, which is purer and less aected by systematics; the
employment of Planck PR4 lensing maps that exhibit ≈ 20% lower noise compared to the
2018 release; and the application of a new catalog-based estimator (refer to Sec. 4), which
eectively avoids numerical instabilities and pixelization complications [83, 84]. Furthermore,
while the quasar auto-correlation Cℓgg in [30] was signicantly contaminated by systematics,
we detect no excess large-scale power in the DESI DR1 sample. As a result, we do not need
any additional noise modeling in the covariance matrix.
A similar analysis to the one presented here was previously performed using the Quaia
quasar sample [29] and combining Cℓκg and Cℓgg to obtain fNL = −28+26 −24 for p = 16 and
+19 +26.7
fNL = −20−18 for p = 10, and fNL = −138−25 when using the cross-correlation alone (and
p = 1). Quaia is a photometric catalog of 13 million quasars covering nearly the full sky.
Notably, our DESI DR1 spectroscopic sample, which only covers ∼ 17% of the sky and relies
exclusively on Cℓκg , achieves comparable constraining power. This highlights the impressive
constraining power of DESI, which will only grow with upcoming data releases. At the same
time, our constraints are weaker than the tightest bounds from LSS, which come from the 3D

– 27 –
power spectrum of DESI DR1 QSOs and LRGs, which yield fNL = −3 ± 9 [11].

50
fNL

−50

−100
1

3
al

ht s

as
s
n
ht

1 ght

ht

ig ht

3.
bi
ci

bi

on
s,
ig

ig

we eig

Si s

bi

bi

bi ola r
du

=
le

A
p= ei

rp ea
p

nt
we

we

as ti
ng
+ gw

pt w

in

o
Fi

te in
ax

ta
N
zm
g

al

O log

In L
5

ns
lo

lo

6.
im

Co
ta

ta

+ ta

=
pt
Ca

Ca

Ca
O

o
N
Figure 17: Constraints on fNL obtained in our ducial analysis (shaded horizontal band),
and adopting alternative analysis strategies that test the robustness of our results. The
numerical values shown here are listed in Table 2-8 and described in Sec. 7.

In Sec. 6, we derived the optimal weights for Cℓκg to maximize the fNL signal. This
optimal weighting scheme improves the constraints by only 4–5%, smaller than the 10–15%
gain typically achieved in three-dimensional power spectrum analyses. Tests with a synthetic
data vector and alternative tomographic binning (Sec. 7.1) suggest that most of the additional
information from the redshift evolution of the PNG signal is already captured by splitting the
sample into tomographic bins, which is the baseline conguration of this analysis. Finally, we
test the impact of freeing the primordial power spectrum parameters As and ns , nding that
our fNL constraints remain completely consistent.
We presented a cosmological analysis of the large-scale clustering of DESI DR1 spectro-
scopically conrmed quasars, aimed at constraining the fNL parameter, while a full cosmo-
logical interpretation of the sample is provided in the companion paper [55]. Our constraints
on PNG, which proved to be very stable against many analysis variations (see Fig. 17), are
statistically weaker than those from three-dimensional power-spectrum analyses, as the mea-
surement remains noise-dominated on the large scales that carry most of the information.
However, this situation will improve signicantly with future DESI data releases, which will
oer greater statistical power, cover a larger sky area, and access additional scales where
the fNL signal is strongest. The inclusion of the quasar auto-correlation will further enhance
constraining power, provided that imaging systematics are well understood and under control.

Data Availability
The data used in this work are publicly available as part of DESI Data Release 1 (see https:
//[Link]/doc/releases/). The data points corresponding to the gures, as
well as the chains required to reproduce the corner plots, are available on Zenodo at https:
//[Link]/uploads/17965252. All codes and packages employed in this analysis are
publicly accessible and are referenced throughout the paper.

Acknowledgments
WP acknowledges support from the Natural Sciences and Engineering Research Council of
Canada (NSERC), [funding reference number RGPIN-2025-03931] and from the Canadian

– 28 –
Space Agency. Research at Perimeter Institute is supported in part by the Government of
Canada through the Department of Innovation, Science and Economic Development Canada
and by the Province of Ontario through the Ministry of Colleges and Universities. This
research was enabled in part by support provided by Compute Ontario ([Link])
and the Digital Research Alliance of Canada ([Link]).
This material is based upon work supported by the U.S. Department of Energy (DOE),
Oce of Science, Oce of High-Energy Physics, under Contract No. DE–AC02–05CH11231,
and by the National Energy Research Scientic Computing Center, a DOE Oce of Sci-
ence 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 Con-
tract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Labo-
ratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon
and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative En-
ergies and Atomic Energy Commission (CEA); the National Council of Humanities, Science
and Technology of Mexico (CONAHCYT); the Ministry of Science, Innovation and Universi-
ties of Spain (MICIU/AEI/10.13039/501100011033), and by the DESI Member Institutions:
[Link]
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 Observa-
tory, 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 agree-
ment 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-eld Infrared Survey Explorer (NEO-
WISE), 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, Oce of Science, Oce of High Energy Physics of the U.S. Department of
Energy; the National Energy Research Scientic Computing Center, a DOE Oce of Sci-
ence 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 ac-
knowledgments can be found at [Link] Any opinions, ndings,
and conclusions or recommendations expressed in this material are those of the author(s)
and do not necessarily reect 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 scientic research on I’oligam Du’ag (Kitt Peak), a mountain with
particular signicance to the Tohono O’odham Nation.

– 29 –
References
[1] J. M. Bardeen, P. J. Steinhardt and M. S. Turner, Spontaneous creation of almost scale-free
density perturbations in an inationary universe, Physical Review D 28 (1983) 679.
[2] T. Falk, R. Rangarajan and M. Srednicki, Dependence of density perturbations on the coupling
constant in a simple model of ination, Physical Review D 46 (1992) 4232.
[3] V. Acquaviva, N. Bartolo, S. Matarrese and A. Riotto, Second-order cosmological
perturbations from ination, arXiv preprint astro-ph/0209156 (2002) .
[4] Y. Akrami, K. J. Andersen, M. Ashdown, C. Baccigalupi, M. Ballardini, A. J. Banday et al.,
Planck intermediate results-lvii. joint planck l and h data processing, Astronomy &
Astrophysics 643 (2020) A42.
[5] J. Maldacena, Non-gaussian features of primordial uctuations in single eld inationary
models, Journal of High Energy Physics 2003 (2003) 013 [astro-ph/0210603].
[6] V. Acquaviva, N. Bartolo, S. Matarrese and A. Riotto, Gauge-invariant second-order
perturbations and non-Gaussianity from ination, Nuclear Physics B 667 (2003) 119
[astro-ph/0209156].
[7] N. Bartolo, E. Komatsu, S. Matarrese and A. Riotto, Non-Gaussianity from ination: theory
and observations, PhysRep 402 (2004) 103 [astro-ph/0406398].
[8] Planck Collaboration, Y. Akrami, F. Arroja, M. Ashdown, J. Aumont, C. Baccigalupi et al.,
Planck 2018 results. IX. Constraints on primordial non-Gaussianity, A&A 641 (2020) A9
[1905.05697].
[9] K. N. Abazajian, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, K. S. Arnold et al., CMB-S4
Science Book, First Edition, arXiv e-prints (2016) arXiv:1610.02743 [1610.02743].
[10] N. Dalal, O. Doré, D. Huterer and A. Shirokov, Imprints of primordial non-Gaussianities on
large-scale structure: Scale-dependent bias and abundance of virialized objects, PRD 77 (2008)
123514 [0710.4560].
[11] E. Chaussidon, C. Yèche, A. de Mattia, C. Payerne, P. McDonald, A. Ross et al.,
Constraining primordial non-Gaussianity with DESI 2024 LRG and QSO samples, arXiv
preprint arXiv:2411.17623 (2024) .
[12] A. Slosar, C. Hirata, U. Seljak, S. Ho and N. Padmanabhan, Constraints on local primordial
non-Gaussianity from large scale structure, JCAP 2008 (2008) 031 [0805.3580].
[13] J.-Q. Xia, C. Baccigalupi, S. Matarrese, L. Verde and M. Viel, Constraints on primordial
non-Gaussianity from large scale structure probes, JCAP 2011 (2011) 033 [1104.5015].
[14] N. Nikoloudakis, T. Shanks and U. Sawangwit, Clustering analysis of high-redshift luminous
red galaxies in Stripe 82, MNRAS 429 (2013) 2032 [1204.3609].
[15] A. R. Pullen and C. M. Hirata, Systematic Eects in Large-Scale Angular Power Spectra of
Photometric Quasars and Implications for Constraining Primordial Non-Gaussianity, Publ.
Astron. Soc. Pac. 125 (2013) 705 [1212.4500].
[16] B. Leistedt and H. V. Peiris, Exploiting the full potential of photometric quasar surveys:
optimal power spectra through blind mitigation of systematics, MNRAS 444 (2014) 2
[1404.6530].
[17] T. Giannantonio, A. J. Ross, W. J. Percival, R. Crittenden, D. Bacher, M. Kilbinger et al.,
Improved primordial non-Gaussianity constraints from measurements of galaxy clustering and
the integrated Sachs-Wolfe eect, PRD 89 (2014) 023511 [1303.1349].
[18] B. Leistedt, H. V. Peiris and N. Roth, Constraints on Primordial Non-Gaussianity from 800
000 Photometric Quasars, Phys. Rev. Lett. 113 (2014) 221301 [1405.4315].

– 30 –
[19] S. Ho, N. Agarwal, A. D. Myers, R. Lyons, A. Disbrow, H.-J. Seo et al., Sloan Digital Sky
Survey III photometric quasar clustering: probing the initial conditions of the Universe, JCAP
2015 (2015) 040 [1311.2597].
[20] M. Rezaie, A. J. Ross, H.-J. Seo, H. Kong, A. Porredon, L. Samushia et al., Local primordial
non-Gaussianity from the large-scale clustering of photometric DESI luminous red galaxies,
Monthly Notices of the Royal Astronomical Society 532 (2024) 1902.
[21] K. M. Smith, O. Zahn and O. Dore, Detection of Gravitational Lensing in the Cosmic
Microwave Background, Phys. Rev. D76 (2007) 043510 [0705.3980].
[22] C. M. Hirata, S. Ho, N. Padmanabhan, U. Seljak and N. A. Bahcall, Correlation of CMB with
large-scale structure: II. Weak lensing, Phys. Rev. D78 (2008) 043520 [0801.0644].
[23] T.-C. Chang, U.-L. Pen, K. Bandura and J. B. Peterson, Hydrogen 21-cm Intensity Mapping
at redshift 0.8, arXiv e-prints (2010) arXiv:1007.3709 [1007.3709].
[24] J. Rhodes, S. Allen, B. A. Benson, T. Chang, R. de Putter, S. Dodelson et al., Exploiting
Cross Correlations and Joint Analyses, arXiv e-prints (2013) arXiv:1309.5388 [1309.5388].
[25] T. Giannantonio and W. J. Percival, Using correlations between cosmic microwave background
lensing and large-scale structure to measure primordial non-Gaussianity., MNRAS 441 (2014)
L16 [1312.5154].
[26] A. B. Lizancos, W. Coulton, A. Challinor, B. Sherwin and Y. Mehta, A halo model of
extragalactic contamination to CMB lensing, delensing, and cross-correlations, Journal of
Cosmology and Astroparticle Physics 2025 (2025) 031.
[27] G. Piccirilli, G. Fabbian, D. Alonso, K. Storey-Fisher, J. Carron, A. Lewis et al., Growth
history and quasar bias evolution at z< 3 from Quaia, Journal of Cosmology and Astroparticle
Physics 2024 (2024) 012.
[28] J. Bermejo-Climent, R. Demina, A. Krolewski, E. Chaussidon, M. Rezaie, S. Ahlen et al.,
Constraints on primordial non-gaussianity from the cross-correlation of desi luminous red
galaxies and planck cmb lensing, Astronomy & Astrophysics 698 (2025) A177.
[29] G. Fabbian, D. Alonso, K. Storey-Fisher and T. Cornish, Constraints on primordial
non-gaussianity from quaia, 2504.20992.
[30] A. Krolewski, W. J. Percival, S. Ferraro, E. Chaussidon, M. Rezaie, J. N. Aguilar et al.,
Constraining primordial non-Gaussianity from DESI quasar targets and Planck CMB lensing,
Journal of Cosmology and Astroparticle Physics 2024 (2024) 021.
[31] J. Carron, M. Mirmelstein and A. Lewis, Cmb lensing from planck pr4 maps, Journal of
Cosmology and Astroparticle Physics 2022 (2022) 039.
[32] U. Seljak, Extracting Primordial Non-Gaussianity without Cosmic Variance, Phys. Rev. Lett.
102 (2009) 021302 [0807.1770].
[33] O. Doré, J. Bock, M. Ashby, P. Capak, A. Cooray, R. de Putter et al., Cosmology with the
SPHEREX All-Sky Spectral Survey, arXiv e-prints (2014) arXiv:1412.4872 [1412.4872].
[34] D. Yamauchi, K. Takahashi and M. Oguri, Constraining primordial non-Gaussianity via a
multitracer technique with surveys by Euclid and the Square Kilometre Array, PRD 90 (2014)
083520 [1407.5453].
[35] D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo and L. Verde, Constraining
primordial non-Gaussianity with bispectrum and power spectrum from upcoming optical and
radio surveys, MNRAS 478 (2018) 1341 [1801.09280].
[36] M. Schmittfull and U. Seljak, Parameter constraints from cross-correlation of CMB lensing
with galaxy clustering, PRD 97 (2018) 123540 [1710.09465].

– 31 –
[37] S. Ferraro and M. J. Wilson, Ination and Dark Energy from spectroscopy at z > 2, Bull.
AAS 51 (2019) 72 [1903.09208].
[38] D. Gualdi, H. Gil-Marín and L. Verde, Joint analysis of anisotropic power spectrum,
bispectrum and trispectrum: application to N-body simulations, JCAP 2021 (2021) 008
[2104.03976].
[39] D. J. Schlegel, S. Ferraro, G. Aldering, C. Baltay, S. BenZvi, R. Besuner et al., A
Spectroscopic Road Map for Cosmic Frontier: DESI, DESI-II, Stage-5, arXiv e-prints (2022)
arXiv:2209.03585 [2209.03585].
[40] A. Dey, D. J. Schlegel, D. Lang, R. Blum, K. Burleigh, X. Fan et al., Overview of the DESI
Legacy Imaging Surveys, AJ 157 (2019) 168 [1804.08657].
[41] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini et al.,
Planck 2018 results-viii. gravitational lensing, Astronomy & Astrophysics 641 (2020) A8.
[42] Planck Collaboration, Y. Akrami, F. Arroja, M. Ashdown, J. Aumont, C. Baccigalupi et al.,
Planck 2018 results. I. Overview and the cosmological legacy of Planck, ArXiv e-prints (2018)
[1807.06205].
[43] E. Komatsu and D. N. Spergel, Acoustic signatures in the primary microwave background
bispectrum, PRD 63 (2001) 063002 [astro-ph/0005036].
[44] V. Desjacques and U. Seljak, Primordial non-Gaussianity from the large-scale structure,
Classical and Quantum Gravity 27 (2010) 124011 [1003.5020].
[45] A. Barreira, Can we actually constrain fNL using the scale-dependent bias eect? An
illustration of the impact of galaxy bias uncertainties using the BOSS DR12 galaxy power
spectrum, Journal of Cosmology and Astroparticle Physics 2022 (2022) 013.
[46] A. Barreira, On the impact of galaxy bias uncertainties on primordial non-gaussianity
constraints, Journal of Cosmology and Astroparticle Physics 2020 (2020) 031.
[47] M. Biagetti, T. Lazeyras, T. Baldauf, V. Desjacques and F. Schmidt, Verifying the
consistency relation for the scale-dependent bias from local primordial non-gaussianity,
Monthly Notices of the Royal Astronomical Society 468 (2017) 3277.
[48] J. M. Sullivan, T. Prijon and U. Seljak, Learning to concentrate: multi-tracer forecasts on
local primordial non-gaussianity with machine-learned bias, Journal of Cosmology and
Astroparticle Physics 2023 (2023) 004.
[49] A. Slosar, C. Hirata, U. Seljak, S. Ho and N. Padmanabhan, Constraints on local primordial
non-gaussianity from large scale structure, Journal of Cosmology and Astroparticle Physics
2008 (2008) 031.
[50] P. Breiding, M. Chiaberge, E. Lambrides, E. T. Meyer, S. Willner, B. Hilbert et al., Powerful
radio-loud quasars are triggered by galaxy mergers in the cosmic bright ages, The
Astrophysical Journal 963 (2024) 91.
[51] D. N. Limber, The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating
Density Field., ApJ 117 (1953) 134.
[52] S. Chiarenza, M. Bonici, W. Percival and M. White, BLAST: Beyond Limber Angular power
Spectra Toolkit. A fast and ecient algorithm for 3x2pt analysis, arXiv preprint
arXiv:2410.03632 (2024) .
[53] S. Chiarenza, M. Bonici et al., “Evolving [Link] into a fast and dierentiable toolkit for
full-physics 6×2pt cosmology.” in preparation, 2026.
[54] R. Scranton, B. Ménard, G. T. Richards, R. C. Nichol, A. D. Myers, B. Jain et al., Detection
of Cosmic Magnication with the Sloan Digital Sky Survey, ApJ 633 (2005) 589
[astro-ph/0504510].

– 32 –
[55] R. de Belsunce, A. Krolewski, S. Chiarenza, E. Chaussidon, S. Ferraro, B. Hadzhiyska et al.,
Cosmology from Planck CMB Lensing and DESI DR1 Quasar Tomography, arXiv preprint
arXiv:2506.22416 (2025) .
[56] E. Chaussidon, C. Yèche, N. Palanque-Delabrouille, D. M. Alexander, J. Yang, S. Ahlen et al.,
Target selection and validation of DESI quasars, The Astrophysical Journal 944 (2023) 107.
[57] M. Abdul-Karim, A. Adame, D. Aguado, J. Aguilar, S. Ahlen, S. Alam et al., Data Release 1
of the Dark Energy Spectroscopic Instrument, arXiv preprint arXiv:2503.14745 (2025) .
[58] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen et al., The
DESI Experiment Part II: Instrument Design, arXiv e-prints (2016) arXiv:1611.00037
[1611.00037].
[59] B. Abareshi, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander, R. Alfarsy et al., Overview of
the instrumentation for the Dark Energy Spectroscopic Instrument, The Astronomical Journal
164 (2022) 207.
[60] A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen, C. A. Prieto et al., The DESI
experiment part I: science, targeting, and survey design, arXiv preprint arXiv:1611.00036
(2016) .
[61] J. H. Silber, P. Fagrelius and K. Fanning, The Robotic Multi-Object Focal Plane System of the
Dark Energy Spectroscopic Instrument, .
[62] T. N. Miller, P. Doel, G. Gutierrez, R. Besuner, D. Brooks, G. Gallo et al., The Optical
Corrector for the Dark Energy Spectroscopic Instrument, arXiv preprint arXiv:2306.06310
(2023) .
[63] C. Poppett, L. Tyas, J. Aguilar, C. Bebek, D. Bramall, T. Claybaugh et al., Overview of the
Fiber System for the Dark Energy Spectroscopic Instrument, The Astronomical Journal 168
(2024) 245.
[64] E. F. Schlay, D. Kirkby, D. J. Schlegel, A. D. Myers, A. Raichoor, K. Dawson et al., Survey
Operations for the Dark Energy Spectroscopic Instrument, The Astronomical Journal 166
(2023) 259.
[65] N. Sailer, G. Farren, S. Ferraro and M. White, A cookbook for cmb lensing cross-correlations,
2024.
[66] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein et al., The DESI Experiment, a
whitepaper for Snowmass 2013, arXiv preprint arXiv:1308.0847 (2013) .
[67] A. Adame, J. Aguilar, S. Ahlen, S. Alam, G. Aldering, D. Alexander et al., Validation of the
scientic program for the Dark Energy Spectroscopic Instrument, The Astronomical Journal
167 (2024) 62.
[68] A. Adame, J. Aguilar, S. Ahlen, S. Alam, G. Aldering, D. Alexander et al., The early data
release of the Dark Energy Spectroscopic Instrument, The Astronomical Journal 168 (2024)
58.
[69] A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez et al., DESI 2024 VI:
Cosmological constraints from the measurements of baryon acoustic oscillations, Journal of
Cosmology and Astroparticle Physics 2025 (2025) 021.
[70] A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez et al., DESI 2024 IV:
Baryon acoustic oscillations from the Lyman alpha forest, Journal of Cosmology and
Astroparticle Physics 2025 (2025) 124.
[71] A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez et al., DESI 2024 III:
baryon acoustic oscillations from galaxies and quasars, Journal of Cosmology and
Astroparticle Physics 2025 (2025) 012.

– 33 –
[72] A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, C. A. Prieto et al., DESI 2024 VII:
Cosmological Constraints from the Full-Shape Modeling of Clustering Measurements, arXiv
preprint arXiv:2411.12022 (2024) .
[73] A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez et al., DESI 2024 II:
sample denitions, characteristics, and two-point clustering statistics, Journal of Cosmology
and Astroparticle Physics 2025 (2025) 017.
[74] A. Ross, J. Aguilar, S. Ahlen, S. Alam, A. Anand, S. Bailey et al., The construction of
large-scale structure catalogs for the Dark Energy Spectroscopic Instrument, Journal of
Cosmology and Astroparticle Physics 2025 (2025) 125.
[75] J. Guy, S. Bailey, A. Kremin, S. Alam, D. Alexander, C. A. Prieto et al., The spectroscopic
data processing pipeline for the Dark Energy Spectroscopic Instrument, The Astronomical
Journal 165 (2023) 144.
[76] S. Bailey et al., Redrock: Spectroscopic Classication and Redshift Fitting for the Dark Energy
Spectroscopic Instrument, 2024.
[77] A. Brodzeller, K. Dawson, S. Bailey, J. Yu, A. J. Ross, A. Bault et al., Performance of the
quasar spectral templates for the Dark Energy Spectroscopic Instrument, The Astronomical
Journal 166 (2023) 66.
[78] A. Anand, J. Guy, S. Bailey, J. Moustakas, J. Aguilar, S. Ahlen et al., Archetype-based redshift
estimation for the Dark Energy Spectroscopic Instrument survey, The Astronomical Journal
168 (2024) 124.
[79] A. Lewis and A. Challinor, Weak gravitational lensing of the cmb, Physics Reports 429 (2006)
1.
[80] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke et al.,
Healpix: A framework for high-resolution discretization and fast analysis of data distributed
on the sphere, The Astrophysical Journal 622 (2005) 759.
[81] A. S. Maniyar, Y. Ali-Haïmoud, J. Carron, A. Lewis and M. S. Madhavacheril, Quadratic
estimators for cmb weak lensing, Physical Review D 103 (2021) 083524.
[82] D. Alonso, J. Sanchez, A. Slosar and L. D. E. S. Collaboration, A unied pseudo-Cℓ
framework, Monthly Notices of the Royal Astronomical Society 484 (2019) 4127.
[83] K. Wolz, D. Alonso and A. Nicola, Catalog-based pseudo-cℓ ’s, Journal of Cosmology and
Astroparticle Physics 2025 (2025) 028.
[84] A. Baleato Lizancos and M. White, Harmonic analysis of discrete tracers of large-scale
structure, Journal of Cosmology and Astroparticle Physics 2024 (2024) 010.
[85] M. Reinecke, S. Belkner and J. Carron, Improved cosmic microwave background (de-) lensing
using general spherical harmonic transforms, Astronomy & Astrophysics 678 (2023) A165.
[86] N. Tessore and A. Hall, Shot noise in clustering power spectra, 2025.
[87] M. Maus, M. White, N. Sailer, A. B. Lizancos, S. Ferraro, S. Chen et al., A joint analysis of
3D clustering and galaxy\times CMB-lensing cross-correlations with DESI DR1 galaxies,
arXiv preprint arXiv:2505.20656 (2025) .
[88] G. S. Farren, A. Krolewski, N. MacCrann, S. Ferraro, I. Abril-Cabezas, R. An et al., The
Atacama Cosmology Telescope: Cosmology from cross-correlations of unWISE galaxies and
ACT DR6 CMB lensing, The Astrophysical Journal 966 (2024) 157.
[89] A. Benoit-Lévy, T. Dechelette, K. Benabed, J.-F. Cardoso, D. Hanson and S. Prunet, Full-sky
CMB lensing reconstruction in presence of sky-cuts, Astronomy & Astrophysics 555 (2013)
A37.

– 34 –
 [90] J. Carron, Real-world cmb lensing quadratic estimator power spectrum response, Journal of
Cosmology and Astroparticle Physics 2023 (2023) 057.
[91] A. J. Ross, J. Bautista, R. Tojeiro, S. Alam, S. Bailey, E. Burtin et al., The Completed
SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Large-scale structure catalogues
for cosmological analysis, Monthly Notices of the Royal Astronomical Society 498 (2020) 2354.
[92] M. Rezaie, H.-J. Seo, A. J. Ross and R. C. Bunescu, Improving galaxy clustering
measurements with deep learning: analysis of the DECaLS DR7 data, Monthly Notices of the
Royal Astronomical Society 495 (2020) 1613.
[93] E. Chaussidon, C. Yèche, N. Palanque-Delabrouille, A. de Mattia, A. D. Myers, M. Rezaie
et al., Angular clustering properties of the DESI QSO target selection using DR9 Legacy
Imaging Surveys, MNRAS 509 (2022) 3904 [2108.03640].
[94] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon et al., healpy: equal area
pixelization and spherical harmonics transforms for data on the sphere in Python, Journal of
Open Source Software 4 (2019) 1298.
[95] N. Tessore, A. Loureiro, B. Joachimi, M. von Wietersheim-Kramsta and N. Jerey, GLASS:
Generator for Large Scale Structure, arXiv preprint arXiv:2302.01942 (2023) .
[96] D. Bianchi, M. Hanif, A. C. Rosell, J. Lasker, A. Ross, M. Pinon et al., Characterization of
DESI ber assignment incompleteness eect on 2-point clustering and mitigation methods for
DR1 analysis, Journal of Cosmology and Astroparticle Physics 2025 (2025) 074.
[97] E. Castorina, N. Hand, U. Seljak, F. Beutler, C.-H. Chuang, C. Zhao et al., Redshift-weighted
constraints on primordial non-Gaussianity from the clustering of the eBOSS DR14 quasars in
Fourier space, JCAP 2019 (2019) 010 [1904.08859].
[98] E.-M. Mueller, W. J. Percival and R. Ruggeri, Optimizing primordial non-Gaussianity
measurements from galaxy surveys, MNRAS 485 (2019) 4160 [1702.05088].
[99] E.-M. Mueller, W. Percival, E. Linder, S. Alam, G.-B. Zhao, A. G. Sánchez et al., The
clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey:
constraining modied gravity, MNRAS 475 (2018) 2122 [1612.00812].
[100] M. S. Cagliari, E. Castorina, M. Bonici and D. Bianchi, Optimal constraints on Primordial
non-Gaussianity with the eBOSS DR16 quasars in Fourier space, JCAP 2024 (2024) 036
[2309.15814].
[101] H. A. Feldman, N. Kaiser and J. A. Peacock, Power spectrum analysis of three-dimensional
redshift surveys, arXiv preprint astro-ph/9304022 (1993) .
[102] A. Lewis, A. Challinor and A. Lasenby, Ecient Computation of Cosmic Microwave
Background Anisotropies inClosed Friedmann-Robertson-Walker Models, The Astrophysical
Journal 538 (2000) 473.
[103] T. E. Fjelde, K. Xu, D. Widmann, M. Tarek, C. Per, M. Trapp et al., [Link]: a
general-purpose probabilistic programming language, ACM Trans. Probab. Mach. Learn. (2025)
.
[104] M. D. Homan, A. Gelman et al., The no-u-turn sampler: adaptively setting path lengths in
hamiltonian monte carlo., J. Mach. Learn. Res. 15 (2014) 1593.
[105] E. Castorina, N. Hand, U. Seljak, F. Beutler, C.-H. Chuang, C. Zhao et al., Redshift-weighted
constraints on primordial non-Gaussianity from the clustering of the eBOSS DR14 quasars in
Fourier space, JCAP 2019 (2019) 010 [1904.08859].
[106] C. Modi, M. White and Z. Vlah, Modeling CMB lensing cross correlations with CLEFT,
JCAP 2017 (2017) 009 [1706.03173].

– 35 –
[107] S.-F. Chen, M. White, J. DeRose and N. Kokron, Cosmological analysis of three-dimensional
BOSS galaxy clustering and Planck CMB lensing cross correlations via Lagrangian
perturbation theory, JCAP 2022 (2022) 041 [2204.10392].

– 36 –
A Derivation of the optimal fNL weights for Cℓκg

As stated in Sec. 6, we follow the derivation of the optimal fNL weights in [99] as well as the
denitions of the angular power spectra in narrow redshift bins as described in [106, 107].
Following [99], the optimal weight for a generic observable Ô is dened as:

∂O(z)
wÔ (z) = dW, (A.1)
∂fNL

where dW = C −1 accounts for the statistical uncertainty of the observable, while ∂O(z) ∂fNL
captures the redshift evolution of the theoretical response to fNL . Note that we distinguish
between the measured, redshift-averaged, observable Ô and its redshift-dependent theoretical
prediction O(z). In [99], for example, the observable is the monopole P̂0 , while its redshift
evolution is described by P0 (z), which depends on b(z), D(z), and bΦ (z), introduced by
the scale dependent bias (Eq. 2.8). For simplicity of the derivation, we assume the Limber
approximation is valid and therefore use:
  
κ(δ+PNG) Wκ (χ) ℓ + 12
Cℓ = dχ Wδ+PNG (χ) Pmm kℓ = ,z , (A.2)
χ2 χ

where
H(z(χ)) dN
Wδ+PNG (χ) = (b(z(χ)) + ∆b(k, z(χ)))  (A.3)
c dχ

Here, we expressed n(z) = dNdχ to make it evident how we can replace the integral dχ dN
dχ
with a sum over galaxies. Starting from here, we want to apply the denition of optimal
weights in Eq. A.1 to the observable Ĉℓκg . The rst component is ∂O(z)
∂fNL :

κ(δ+PNG)  
∂Cℓ Wκ (χ) bΦ (z) ℓ + 12
= D(z)Pmm k =  (A.4)
∂fNL χ2 TΦ→δ (k) χ

Where we have factored out the redshift dependence of the power spectrum using the standard
scalings: Pmm (k, z) = D 2 (z)Pmm (k) and TΦ→δ (k, z) = D(z) TΦ→δ (k). Note that, when
normalizing the weights, the k-dependent factor Pmm (k)TΦ→δ (k) cancels in each k-mode, as
well as constant factors like 2δc . Thus, the rst piece of the weight is:

Wκ (χ)
D(z) (b(z) − p) (A.5)
χ2

For the second piece of the weight, we consider the covariance of the cross-correlation between
a narrow redshift bin and CMB lensing, given by:

(Cℓgg + Nℓgg )(Cℓκκ + Nℓκκ )

Cov(Cℓκg ) = (A.6)
(2ℓ + 1) ∆ℓ fsky
 2
where we are neglecting the subdominant contribution from the cross-correlation Cℓκg . For
galaxies in a single narrow bin of width ∆χ, we get:
 
4π 1 (Cℓκκ + Nℓκκ )
Cov(Cℓκδ ) = Pmm + (A.7)
dV n̄ (2ℓ + 1) ∆ℓ fsky

– 37 –
where dV = 4πχ2 dχ is the volume of each shell. The term entering the optimal weights is
then the inverse covariance:
 
κg −1 dV 1 (2ℓ + 1)∆ℓfsky
(Cov(Cℓ )) = n̄ (A.8)
4π 1 + n̄P (Cℓκκ + Nℓκκ )

The rst term will give the usual FKP weight:

1
wFKP = , (A.9)
1 + n̄P
while the rest is redshift-independent, and thus cancels when we normalize the weights. The
volume integral in dV = 4πχ2 dχ can be transformed in a sum over galaxies through n̄, and
after performing it to normalize the weights, we nd:
Wκ (χ(z))
wopt (z) = bΦ (z) D(z) wFKP  (A.10)
χ2 (z)

B Restricting to quasars with z < 3.1

The best fNL measurement from LSS available to date was performed in [11]. In that analysis,
only the DESI quasars in the range 08 < z < 31 were considered. This excludes ≈ 32500
quasars from the sample. That choice was motivated by the fact that the covariance matrix
is estimated from simulations, the EZMocks, which have a simulation box that is sucient to
emulate the quasars up to zmax = 31 without repeating the box. We test the impact of this

Measured Cκg bin 3

×10−6 regressis weights Catalog weights
2.5
Theory Theory
zmax = 3.5 zmax = 3.5
2.0
zmax = 3.1 zmax = 3.1

1.5

1.0
C

0.5

0.0

−0.5

0 50 100 150 200 250 300 0 50 100 150 200 250 300
 

Figure 18: Measured angular cross-power spectrum Cℓκg in the third tomographic bin, chang-
ing the maximum redshift. Left: measurement performed with the default weighting scheme,
the regressis linear weights. Right: Same measurement, but the systematic weights are the
linear weights implemented in the catalog. See Sec. 5 for a more detailed explanation of the
weights.

choice to our analysis: we keep the default conguration of 3 tomographic bins, t a bias
amplitude for each bin and we do not apply optimal weights. This analysis change aects
the measurement, theory curve and therefore the covariance matrix, which is re-computed
and tuned following the procedure outlined in Sec. 4. As in the previous cases, results are

– 38 –
 p regressis linear weights catalog linear weights
16 fNL = −10+29 2
−35 (χ = 48) fNL = −23+28 2
−34 (χ = 48)
1 fNL = −1+20 2
−24 (χ = 48) fNL = −9+20 2
−23 (χ = 48)

Table 11: Results for the analysis performed with a value of zmax = 31 instead of zmax = 35.

presented for variation of systematic weights choice and values of the parameter p in Table 11.
Removing the higher redshift part of the sample does not aect the error bars, but has an
eect on the best t values, which move around of ≈ 13 σ in the ducial scenario. It is
interesting to look at the eect on the bias amplitudes: the rst two bins do not show a
relevant shift, with the best t values being completely compatible with the ones showed in
Table 3. However, the best t amplitude in the third bin is b3 = 062 ± 013, a value shift
of ≈ 05 σ. Finally, the χ2 in this case is lower than in the default scenario (48 vs 56). The
dierence is driven by the rst multipole (ℓ = 65) in the third tomographic bin, as shown
in Fig. 18. Excluding the higher redshift part of the quasar sample causes the measurement
to shift by ≈ 1 σ, with a signicant impact of the χ2 value. The reasons for this shift are
unknown and currently under investigation. From the gure, it is also possible to see why a
lower bias is preferred, in fact, the measurement in the last multipole (ℓ = 275) is lower for
both weighting schemes, explaining the preference for a lower bias value.

C Affiliations
4 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA
5 Department of Physics, University of California, Berkeley, 366 LeConte Hall MC 7300, Berke-
ley, CA 94720-7300, USA
6 Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215

USA
7 Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, Via Celoria 16,

I-20133 Milano, Italy

8 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy
9 Department of Physics & Astronomy, University College London, Gower Street, London,

WC1E 6BT, UK
10 Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt

Lake City, UT 84112, USA

11 Instituto de Física, Universidad Nacional Autónoma de México, Circuito de la Investigación

Cientíca, Ciudad Universitaria, Cd. de México C. P. 04510, México

12 University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, USA
13 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Tech-

nology, Edici Cn, Campus UAB, 08193, Bellaterra (Barcelona), Spain

14 Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edicio Ip, CP

111711, Bogotá, Colombia

15 Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edicio H, CP

111711 Bogotá, Colombia

16 Institut d’Estudis Espacials de Catalunya (IEEC), c/ Esteve Terradas 1, Edici RDIT,

Campus PMT-UPC, 08860 Castelldefels, Spain

17 Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Build-

– 39 –
ing, Portsmouth, PO1 3FX, UK
18 Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913

Bellaterra, Barcelona, Spain

19 University of Virginia, Department of Astronomy, Charlottesville, VA 22904, USA
20 Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA
21 Institut d’Astrophysique de Paris. 98 bis boulevard Arago. 75014 Paris, France
22 IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
23 Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West

Woodru Avenue, Columbus, OH 43210, USA

24 Department of Physics, The Ohio State University, 191 West Woodru Avenue, Columbus,

OH 43210, USA
25 The Ohio State University, Columbus, 43210 OH, USA
26 Department of Physics, University of Michigan, 450 Church Street, Ann Arbor, MI 48109,

USA
27 University of Michigan, 500 S. State Street, Ann Arbor, MI 48109, USA
28 Department of Physics, The University of Texas at Dallas, 800 W. Campbell Rd., Richard-

son, TX 75080, USA

29 NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA
30 Department of Physics and Astronomy, University of California, Irvine, 92697, USA
31 Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Ener-

gies (LPNHE), FR-75005 Paris, France

32 Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra

(Barcelona), Spain
33 Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W

18th Avenue, Columbus, OH 43210, USA

34 Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010

Barcelona, Spain
35 Department of Physics & Astronomy and Pittsburgh Particle Physics, Astrophysics, and

Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara Street, Pittsburgh,
PA 15260, USA
36 Departamento de Física, DCI-Campus León, Universidad de Guanajuato, Loma del Bosque

103, León, Guanajuato C. P. 37150, México

37 Instituto Avanzado de Cosmología A. C., San Marcos 11 - Atenas 202. Magdalena Contr-

eras. Ciudad de México C. P. 10720, México

38 Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA

94720, USA
39 Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008

Granada, Spain
40 Departament de Física, EEBE, Universitat Politècnica de Catalunya, c/Eduard Maristany

10, 08930 Barcelona, Spain

41 Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro, Gwangjin-gu,

Seoul 05006, Republic of Korea

42 CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain
43 Department of Physics & Astronomy, Ohio University, 139 University Terrace, Athens, OH

45701, USA
44 National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road,

Chaoyang District, Beijing, 100101, P. R. China

– 40 –