Constraining cosmology with the Gaia-unWISE Quasar Catalog and CMB lensing: structure growth

We study the angular clustering of Quaia, a Gaia- and unWISE-based catalog of over a million quasars with an exceptionally well-defined selection function. With it, we derive cosmology constraints from the amplitude and growth of structure across cosmic time. We divide the sample into two redshift bins, centered at z = 1.0 and z = 2.1, and measure both overdensity auto-correlations and cross-correlations with maps of the Cosmic Microwave Background convergence measured by Planck. From these data, and including a prior from measurements of the baryon acoustic oscillations scale, we place constraints on the amplitude of the matter power spectrum σ 8 = 0.766 ± 0.034, and on the matter density parameter Ω m = 0.343+0.017 -0.019. These measurements are in reasonable agreement with Planck at the ∼ 1.4σ level, and are found to be robust with respect to observational and theoretical uncertainties. We find that our slightly lower value of σ 8 is driven by the higher-redshift sample, which favours a low amplitude of matter fluctuations. We present plausible arguments showing that this could be driven by contamination of the CMB lensing map by high-redshift extragalactic foregrounds, which should also affect other cross-correlations with tracers of large-scale structure beyond z ∼ 1.5. Our constraints are competitive with those from state-of-the-art 3×2-point analyses, but arise from a range of scales and redshifts that is highly complementary to those covered by cosmic shear data and most galaxy clustering samples. This, coupled with the unprecedented combination of volume and redshift precision achieved by Quaia, allows us to break the usual degeneracy between Ω m and σ 8.


Introduction
Much of the progress in constraining the physical parameters governing the initial conditions and evolution of our Universe is currently driven by the analysis of tracers of the largescale structure.In these analyses, we study the spatial distribution and time evolution of various tracers of the matter density fluctuations, which allows us to constrain the Universe's geometry, as well as the growth of structure within it.Two of the most powerful largescale structure tracers are weak gravitational lensing and galaxy clustering [1][2][3].The former provides largely unbiased maps of the matter fluctuations integrated along the line of sight from the source redshift, while the latter is a biased tracer of these fluctuations at the redshifts of the galaxies being observed.
This complementarity (unbiased vs. biased, cumulative in redshift vs. local) motivates the combined analysis of weak lensing and galaxy clustering data in a technique commonly known as lensing tomography (often also labelled "2×2-point" analysis) [4][5][6][7].In this methodology one measures the two-point auto-correlations of a set of samples of galaxies that are reasonably localised in redshift, together with their cross-correlations with lensing convergence or shear maps from a background source.The cross-correlations effectively "slice" the lensing map into its contributions at different redshifts and, when combined with the auto-correlations, enables us to break the degeneracy between galaxy bias and the amplitude and growth of structure.The same technique can be applied to any projected tracer of the large-scale structure other than weak lensing, and can be used to e.g.measure the evolution in the mean gas pressure [8][9][10][11][12] and star-formation density [13,14], or to reconstruct the redshift distribution of different probes [15,16].
The lensing of the Cosmic Microwave Background (CMB) has become a treasure trove for lensing tomography.At z ∼ 1100, the CMB is a lensing source that sits behind all tracers of the large-scale structure, and thus can be used to reconstruct the growth history at all redshifts.Various analyses have made use of lensing tomography to constrain cosmology, using both photometric [17][18][19][20][21][22] and spectroscopic samples [7,23,24].Similar studies can be carried out in cross-correlation with galaxy weak lensing (cosmic shear), as well as combining all three tracers of structure [25,26].Constraints can be strengthened significantly by adding the auto-correlation of the weak lensing data, in what's commonly called a "3 × 2-point" analysis (sometimes called 5 × 2-point when combining galaxy clustering, cosmic shear, and CMB lensing) [2,3,22,[25][26][27][28][29][30][31].It is often not possible to fully break the degeneracy between growth and geometry with these datasets, and thus results are often reported in terms of the combination S 8 ≡ σ 8 (Ω m /0.3) 0.5 , where σ 8 is defined as the standard deviation of the linear matter overdensity on spheres of radius 8 Mpc/h, and Ω m is the fractional energy density of non-relativistic matter.A growing number of these analyses have consistently recovered measurements of the amplitude of matter fluctuations that, qualitatively, seem to be in tension with the value inferred by CMB experiments such as Planck [32], or the Atacama Cosmology Telescope (ACT) [33,34] (the so-called "S 8 tension").There are hints that the source of this tension may lie in the impact of non-linear baryonic effects, particularly in the case of cosmic shear data, which is more sensitive to small-scale structures (k ∼ 0.1−1 Mpc −1 ) [35][36][37].Nevertheless, similar levels of tension have been found from pure CMB lensing tomography analyses [19,21], and thus the origin and actual significance of this potential tension remains unclear.
In this context, quasars constitute an interesting probe of the large-scale structure, allowing us to trace the matter fluctuations at significantly larger redshifts than standard optical galaxy surveys, covering commensurately large volumes.Their study has been carried out through three-dimensional analysis of their auto-correlation [38][39][40], as well as their crosscorrelation with other tracers, such as the Lyman-α forest [41], with the aim of studying both cosmology and the astrophysics of quasars.Since its first detection [42], the cross-correlation between quasars and CMB lensing maps has increasingly gained interest, particularly using quasar samples selected from the Sloan Digital Sky Survey [43][44][45], and as a way to improve our understanding of quasar properties.The use of this cross-correlation for cosmology, however, has so far been limited, although some works have exploited it to constrain the growth history at high redshifts [25], as well the level of primordial non-Gaussianity [46].The availability of a full-sky quasar sample with a well-understood selection function would facilitate this kind of studies, as well as increase the sensitivity of these measurements by increasing the area overlap with existing CMB lensing maps from experiments such as Planck [47,48].
In this paper we will study CMB lensing tomography making use of Quaia [49], a catalog of high-redshift quasars constructed based on the third data release of the Gaia satellite [50] as well as unWISE infrared observations [51][52][53].Covering the redshift range z 4 with good redshift accuracy (|∆z/(1 + z)| < 0.01 (0.1) for 62% (83%) of the sources), and spanning the full sky, Quaia allows us to constrain the growth of structure over a range of cosmic times that is highly complementary to those accessed by previous analyses, as well as giving us access to significantly larger scales, where the impact of non-linearities in the matter power spectrum is less important.This complementarity is highly relevant to ascertain whether the source of the S 8 tension lies in small-scale physics at low redshifts, an overall lack of power on all scales throughout the Universe's evolution, or a combination of statistical and systematic uncertainties.Furthermore, as we will show, the large-scale reach of the sample will allow us to measure σ 8 and Ω m independently (instead of through the S 8 combination).
The paper is structured as follows.The Quaia sample and the Planck CMB lensing map used in the analysis are described in Section 2. Section 3 presents the theoretical framework underpinning CMB lensing tomography.The data analysis methods used are outlined in Section 4. Section 5 presents the main results of this analysis, as well as a large battery of tests aimed at confirming their robustness.Our conclusions are then summarised in Section 6.

The Gaia-unWISE Quasar Catalog
We use the Gaia-unWISE Quasar Catalog, Quaia [49], as our quasar sample.The catalog is all-sky and highly homogeneous, with 1,295,502 sources brighter than its magnitude limit of G < 20.5 and a median redshift of z = 1.47.It is derived from the quasar candidates identified in the third data release of the space-based Gaia mission, which classified ∼6.6 million sources as potential quasars and estimated their redshifts Gaia's low-resolution BP/RP spectra.Quaia combines Gaia data with mid-infrared information from the unWISE reprocessing of the Wide-field Infrared Survey Explorer (WISE) observations [51][52][53] to improve the sample selection and redshift estimation.The redshift estimation uses a k-Nearest Neighbors model with Gaia and unWISE photometry and the Gaia-estimated spectral redshift as features, and cross-matched SDSS DR16Q quasars with high-precision spectroscopic redshifts [54] as labels.The Quaia redshift estimate, z Quaia , is taken to be the median of the spectroscopic redshifts of the K = 27 nearest neighbors, and the uncertainty σ z is their standard deviation.Of the Quaia redshifts, an estimated 62% (83%) have |∆z/(1 + z)| < 0.01 (0.1) with respect to the SDSS redshifts, and 91% having |∆z/(1 + z)| < 0.2.
In this analysis we use a modified version of the Quaia-modelled selection function, re-fit for each of the two redshift bins used for our fiducial results (see Section 4.1).While Quaia already has relatively minimal selection effects, partly thanks to its space-based observations, the selection function model includes the effects of dust extinction, stellar crowding, and the Gaia scanning law (as detailed in [49]).In Section 5.2.1 we will also explore a brighter sample with G < 20.0, for which the selection function was also refit.

Planck CMB lensing
Our baseline analysis uses the CMB lensing map reconstructed with the Planck PR4 data release based on the NPIPE processing pipeline [55].The release improved over various assumptions of the PR3 2018 analysis [47] that made the analysis sub-optimal in terms of noise treatment and exploited the improved low-level data processing of NPIPE (encompassing additional data and more accurate simulations of the end-to-end data processing) to reach a ∼ 20% improvement in signal-to-noise at all scales.The details of the analysis are summarized in [48,56].The major pipeline improvements included the adoption of a more optimal generalized minimum-variance estimator (GMV) that performs a joint inversevariance Wiener filtering of the temperature and polarization CMB maps at the same time instead of treating both observables separately [57], as well as an additional post-processing Wiener filtering of the lensing maps according to the local reconstruction noise level on the sky that improves the sensitivity of signal-dominated scales.The PR4 lensing release also includes an estimate of the lensing potential based on the 2018 pipeline and using CMB temperature-only and polarization-only data [58].We use these CMB lensing maps to assess the impact of foregrounds in our results, given that the extralagactic foregrounds such as thermal Sunyaev-Zeldovich (tSZ) effect and the Cosmic Infrared Background (CIB) are largely unpolarized at this level of sensitivity.We used the PR3 and PR4 associated lensing simulations to evaluate the consistency between the measurements carried out with the different lensing maps, and to evaluate mask-dependent correction to the lensing map normalization.This is discussed in Section 4.2, and in the PR4 lensing documentation1 .All the data products are publicly available through the Planck legacy archive2 .

Projected galaxy clustering
The first probe we will use in our analysis is the projected distribution of quasars.The core observable in this case is the angular quasar overdensity3 δ g (n) ≡ n g (n)/n g − 1, where n g (n) is the surface density of quasars in the sky direction defined by unit vector n, and ng is its ensemble average.The angular overdensity is related to the three-dimensional overdensity ∆ g (x, t) (defined analogously to the angular quasar overdensity), where x is the tracer position, through the line-of-sight (LOS) projection where p(z) is the redshift distribution of quasars in the sample, z is the redshift to radial comoving distance χ, and H(z) is the expansion rate at redshift z4 .
The relation between the quasar overdensity and the underlying matter overdensity ∆ m is, in general, non-linear, non-local, and stochastic.On sufficiently large-scales (r 10 Mpc), such as those used in this work, a linear bias relation can be used, so that where b g is the linear quasar bias and is an uncorrelated stochastic field representing the non-deterministic component of the ∆ g -∆ m relation.In our fiducial analysis, we will assume that this component is dominated by Poisson shot noise, in which case the power spectrum of is simply P (k) = 1/ Ng , where Ng is the mean 3-dimensional quasar density.In general, may receive additional contributions due to physics on scales below the Lagrangian size of the protohalo hosting the galaxies composing our samples (e.g.halo exclusion effects).The low density of sources in the quasar sample used here makes it unlikely for these contributions to dominate over the pure Poisson term, but we will study the potential impact of stochastic bias in Section 5.2.2.
The observed projected overdensity of sources is also affected by modifications to the properties of the light they emit caused by the large-scale structure.These contributions may be local (e.g.additional red-or blue-shifting caused by peculiar velocities or gravitational potentials), or LOS-integrated (e.g.lensing, integrated Sachs-Wolfe effect, Shapiro time delay) [59,60].Of these, the most relevant effect for distant sources with broad radial projection kernels, such as the samples used here, is the contribution from lensing magnification [61,62].Gravitational lensing causes a coherent distortion in the observed positions of background sources away from foreground overdensities, as well as modifies their observed flux, which can add sources to or remove sources from the sample (depending on the slope of the flux distribution, and on the sign of the lensing convergence).The combination of both effects leads to an extra additive contribution to δ g in Eq. 3.1 of the form where the lensing magnification kernel is and s(z) is the logarithmic slope of the cumulative number counts of sources at the magnitude limit of the sample: where m is the source apparent magnitude, and m lim is the sample's limiting magnitude.Section 4.1.2describes the methods used to calculate s in practice.Finally, ∇ 2 ⊥ and ∇ −2 are the transverse Laplacian and the inverse Laplacian respectively.In the Limber approximation (see Eq. 3.11), they partially cancel each other up to a factor which can be safely ignored on scales > 10 (although our theoretical calculations will account for this factor).We will therefore ignore the combination ∇ 2 ⊥ ∇ −2 in what follows for simplicity.
Finally, it is worth noting that Eq. 3.3 assumes that the matter overdensity is related to the Newtonian potential Φ via Poisson's equation, and that the two scalar metric potentials in the Newtonian gauge, Φ and Ψ are equal (which is valid for general relativity in the absence of anisotropic stress).

CMB weak lensing
Gravitational lensing perturbs the trajectories of CMB photons, thus distorting the observed pattern of CMB anisotropies, and inducing a statistical coupling between different harmonic scales.This distortion can be used to reconstruct the lensing displacement that causes it, and thus produce maps of the lensing convergence κ(n), defined as Here, χ LSS is the comoving distance to the last-scattering surface (LSS).As in the case of lensing magnification, CMB lensing is subject to the same K correction factor, which can be neglected on the scales used here.

Statistics of projected tracers
Consider a projected tracer u(n) of a three-dimensional quantity U (x, z): where q u (χ) is the associated projection kernel.Assuming statistical homogeneity and isotropy, let us define the angular power spectrum C uv , and the three-dimensional power spectrum P U V (k, z) between two pairs of such quantities ((u, U ) and (v, V )) as: Here • • • denotes ensemble averaging, δ K is the Kronecker delta, δ D is the Dirac delta, and u m and U (k, z) are the spherical harmonic transform and the Fourier transform, respectively, of the corresponding configuration-space fields: The angular power spectrum can be connected with the three-dimensional P U V (k) via This equation is valid in the Limber approximation [63,64], which is valid when the kernels under study overlap over a range of χ significantly larger than the correlation length of U and V .This is an excellent approximation for the broad CMB lensing kernel, as well as the redshift distribution of the quasar samples we will study.
Within the linear bias model we use here, both the projected quasar overdensity δ g and the lensing convergence κ can be treated as tracers of the 3D matter overdensities ∆ m , and therefore P U V (k, z) above is simply the matter power spectrum for all correlations studied.The radial kernels for both tracers are: with the magnification kernel q mag given in Eq. 3.4.(1 + z) The matter power spectrum as observed by a hypothetical probe sensitive to different redshifts, as a function of the corresponding angular scale ≡ kχ − 1/2, where χ is the comoving distance to that redshift.The central white band roughly corresponds to the range of scales used in cosmic shear analyses.Within these scales, at low redshifts we are mostly sensitive to the monotonic power-law-like part of the power spectrum, which is also affected by non-linear effects (the scale of non-linearities k NL , defined in Eq. 3.13, is shown by the black solid line).At higher redshifts we are able to better reconstruct the overall shape of the power spectrum at the accessible scales, and to use scales that are less affected by non-linearities.For reference, the plot also shows the position of k = 0.1 Mpc −1 and k = 0.01 Mpc −1 , roughly corresponding to the location of the BAO wiggles and the horizon scale at matter-radiation equality, respectively.

The importance of high-redshift data
With the theoretical model described above, we can now explore the benefits of studying the growth of structure at high redshifts, as enabled by the Quaia quasar sample we use here.
One of the most compelling reasons to study the history of structure growth, is the ongoing debate surrounding the so-called "S 8 tension": the fact that a number of low-redshift probes of structure recover an amplitude of the matter fluctuations that is in slight tension with the value inferred from CMB anisotropies.The evidence of this tension from individual experiments, and from some combinations of them, is mild, although not negligible (between 2σ and 3σ), and the fact that several different probes have consistently found low values of the S 8 parameter has motivated a widespread search for both further evidence of this tension from other datasets, as well as potential explanations for it in terms of both astrophysical effects and potential new physics.
The S 8 tension is largely dominated by cosmic shear measurements [25,[65][66][67][68][69].The weak lensing kernel for sources at z 1 peaks at z ∼ 0.2 − 0.5, and extends down to z ∼ 0 [70].Thus, cosmic shear is mostly sensitive to the matter inhomogeneities at low redshifts and, due to projection effects, on relatively small scales k ∼ 0.1 − 1 Mpc.On the other hand, probes such as the CMB lensing convergence auto-correlation, which access significantly higher redshifts (the CMB lensing kernel peaks at z ∼ 2), and are sensitive to larger, more linear scales, find no evidence of a low S 8 [34].This has motivated the exploration of the mis-modelling of the small-scale matter power spectrum due to non-linearities and baryonic effects as a plausible solution to the S 8 tension [35,37].Additional, complementary analyses of structure formation at high redshifts and on linear scales are important in order to investigate this hypothesis; the clustering of high-redshift quasars offers such an opportunity.
Regardless of its potential ability to shed light on the ongoing S 8 tension, the study of large-scale structure at redshifts z ∼ 2 − 4 is interesting for cosmology in its own right.Although we expect the Universe's expansion to be dominated by non-relativistic matter beyond z ∼ 1, direct observations of this epoch are few and far between.The matter fluctuations at z > 2 could therefore bear the imprint of new early-Universe physics that we are not sensitive to at lower redshift.At the same time, the comoving volume available at high-z gives us access to significantly larger scales than those probed by low-redshift tracers.Fig. 1 illustrates this point.The figure shows the matter power spectrum as it would be observed by an imaginary probe of structure at different redshifts 5 , as a function of the angular multipole, related to the comoving wavenumber via kχ ≡ + 1/2, where χ is the distance to the redshift being probed.The coloured lines show the power spectrum at different redshifts in the range 0.3 < z < 4. The solid black line shows the angular multipole corresponding to the non-linear scale k NL , defined as the comoving scale at which the overdensity variance reaches unity: (3.13) The power spectra are shown in transparent colours above this scale, which marks the regime in which non-linear and baryonic effects become relevant.The dashed black line shows the position of the wavenumber k = 0.1 Mpc −1 , roughly corresponding to the range of scales where baryon acoustic oscillations are most relevant.The dotted black line shows the wavenumber k = 0.01 Mpc −1 , which corresponds approximately to the horizon scale at the matterradiation equality, where the power spectrum reaches its maximum.Finally, the shaded gray bands mask out the scales 30, and 1000 (corresponding to 6 • and 10 respectively), usually discarded in cosmic shear analyses.As we can see, at low redshifts z 0.5, which dominate current cosmic shear data, we mostly have access to the range of scales where the power spectrum is approximately a monotonic power law (projection effects in weak lensing wash out the BAO wiggles), which also receives significant contribution from non-linear scales.In turn, at higher redshifts (z 1) we gain access to larger scales, and are able to better resolve the shape of the power spectrum, reducing also the contribution from non-linear scales.Measurements at these high redshifts thus can allow us to break degeneracies between different cosmological parameters, which modify the power spectrum in complementary ways.
To illustrate this further, let us briefly jump ahead and examine some of the data we will analyse later on.The left panel of Fig. 2 shows the angular auto-correlation of the high-redshift bin we use in our analysis, centered at z 2 (see Fig. 4 and Section 4.1), and the right panel shows the cross-correlation of the same bin with the Planck CMB lensing map.The measured power spectra are shown as black points with error bars, and the red line shows the best-fit model describing these data.The blue solid and dotted lines show the result of changing the cosmic matter fraction Ω m by δΩ m = ±0.1 with respect to its best fit value, while the analogous yellow lines show the result of perturbing σ 8 by the same amount.We can see that, although on small scales ( 400) both parameters modify the measured Quasar auto-spectrum (left panel) and cross-spectrum with the CMB lensing convergence (right panel) for the high-redshift bin used in our analysis, centered at redshift z = 2.1 (see Fig. 4 and Section 4.1).The measurements are shown as black points with error bars, with the best-fit model shown in red.The blue and yellow lines show the theoretical predictions after changing Ω m and σ 8 by 0.1 in both directions, respectively.On the redshifts and scales our data are sensitive to, the effect of both parameters is markedly different, allowing us to break the degeneracy between them.
power spectra in a similar manner (scaling it up and down), when including the full range of scales used in our analysis (marked by the gray dashed bands), their impact is markedly different.While σ 8 modifies the overall normalisation of the power spectrum, Ω m changes its spectral tilt 6 .As we will see, this will allow us to break the degeneracy between these parameters, which affects cosmic shear and lower-redshift data sets in general, and measure them independently of one another.

Quasar overdensity maps
Our fiducial quasar sample will be the full G < 20.5 Quaia catalog, which we divide into two redshift bins (z Quaia < 1.47 and z Quaia > 1.47) each containing approximately the same number of sources (647,749 and 647,753 respectively).Splitting the sample into two bins allows us to explore any trends with redshift (both astrophysical and due to systematics), but we determined that further dividing the sample did not lead to significantly tighter cosmological constraints, In what follows, we will refer to these two bins as "Low-z" and "High-z" respectively.The following subsections describe the methodology used to produce quasar overdensity maps, calibrate their contamination from observational systematics, estimate the slope of the cumulative flux distribution to quantify the impact of magnification bias, and calibrate the redshift distribution of each redshift bin.

From catalog to maps
We construct a map of the projected quasar overdensity in each redshift bin using the HEALPix pixelization scheme [71], with resolution parameter N side = 512 (corresponding to an angular resolution δθ ∼ 0.11 • ).In each pixel p we calculate the quasar overdensity as Here N p is the number of quasars in p, N is the mean number density of quasars in the sample, and w p is the selection function value in p.The selection function is described below in more detail, and can be interpreted as the mean fraction of observed objects (out of those that should have been observed).The mean number density is calculated correcting for the selection function as: 2) The method used to construct the selection function is described in detail in [49].In short, a non-linear relation, modelled as a Gaussian process, is found between the number of sources observed in each pixel and the value of a number of systematics templates in those pixels.This relation was found using pixels of resolution N side = 64 (δθ 1 • ), and considering 4 systematic templates: Milky Way extinction as measured by [72], a star map constructed from Gaia, the specific contribution to the star map from the LMC and SMC, and a measure of depth in Gaia (the so-called M 10 quantity [73]) that maps the Gaia scanning law and source completeness.It is worth noting that the method intrinsically assumes that the underlying galaxy distribution is statistically isotropic.
To avoid instabilities when the selection function is close to zero, we mask all pixels where more than half of the sources are missed on average (w p < 0.5), and only the resulting unmasked pixels are used in Eq. 4.2.Elsewhere, the quasar overdensity is set to zero, and the impact of this masking is accounted for via the pseudo-C estimator (see Section 4.3).The resulting footprint covers 74% and 67% of the sky in the Low-z and High-z bins, respectively.The results presented here were found to be stable against changes in the selection function threshold.In order to downweight regions of the map where a smaller fraction of objects were observed, we again use the selection (truncated at w p = 0.5 as above) as a weighting map/mask when estimating power spectra involving δ g .
To further remove any residual correlation between the quasar overdensity map and the sky systematics listed above, we make use of linear deprojection, as described in [74,75].Effectively, we fit for, and project out, any linear relation between δ g,p and the fluctuations of the different systematic templates about their mean 7 .Since deprojecting only 4 templates leads to the loss of only a small number of modes, we do not correct for the associated deprojection bias when computing angular power spectra.
The use of linear deprojection after having calibrated the selection function using a Gaussian process (which to some extent should also reconstruct linear trends in the data), may seem redundant.However, it is worth emphasizing that, in the presence of small levels of contamination, where a linear expansion is sufficiently accurate, linear deprojection is exact, non-degenerate with other potential trends in the data, and further guarantees no correlation between the cleaned maps and the different systematic templates.We thus employ it to further reduce the impact of any residual contamination.As we will see in Section 5.2.4,although we find that deprojection is indeed able to remove the presence of residual systematics in the data, within the scale cuts used in our analysis, their impact is negligible.

Magnification bias
In order to quantify the impact of magnification on the measured power spectra, we need to estimate the slope of the cumulative source number counts with limiting magnitude s (see Eq. 3.5).The dotted lines show the value of s found for the full sample (i.e.not separated in redshift).The value of s in our fiducial sample, s ∼ 0.4, implies that the impact of lensing magnification on our analysis is negligible.
The most relevant cuts in the Quaia sample are in G-band magnitude, and in colour space.Magnification leads to a coherent shift in magnitude in all bands, and thus the G-band cut is sensitive to it, but not the colour-space cuts.Therefore, as a first method to estimate s, we treat Quaia as a purely magnitude-limited sample, and simply calculate the slope of the number counts at the limiting magnitude 8 .The resulting value of s for the two redshift bins, and for a magnitude limit G < 20.5 is s(z < 1.47) = 0.388 ± 0.004, s(z ≥ 1.47) = 0.420 ± 0.005. (4. 3) The statistical uncertainties were estimated via bootstrap resampling.The result for a brighter sample with G < 20 (which we will explore in Section 5.2.1) is More in detail, the Quaia sample is not exactly magnitude-limited, as its definition also involves a cut in the joint space of magnitude and proper motion.Thus, as an alternative estimate of s, we generate two additional catalogs by first perturbing all magnitudes in the raw catalog by δG = ±0.05,and then imposing the cuts defining our sample.s is then estimated from the difference in the number of objects between these two catalogues.The result is: So far we have treated s as a constant in each redshift bin, whereas in reality it likely evolves with redshift.To explore this, we have divided the sample into 6 redshift bins covering the range z Quaia ∈ [0, 4], and calculated s from the number counts slope in each of them.The resulting s(z) for the G < 20.5 sample is shown in Fig. 3.We see that s only departs from its global value (s = 0.404 ± 0.004 and s = 0.515 ± 0.006 for the G < 20.5 and G < 20 samples respectively) significantly at low redshifts, where the contribution from magnification is negligible.We thus assume s to be constant in our fiducial analysis.
Since, for our sample, the counts slope is close to s ∼ 0.4, the impact of lensing magnification is small, and likely negligible (although we account for it in our theory prediction).This is in contrast with studies carried out with previous optical quasar samples, such as [76,77], where lower values were found (e.g.s = 0.2764 for the DESI-selected quasar sample [46]).This may be due to differences in the selection of the various samples, or to the fact that Quaia targets brighter quasars on average, and hence is sensitive to a different part of the quasar luminosity function.A value of s < 0.4 corresponds to a negative contribution from magnification, which mostly affects the CMB lensing cross-correlation.To make up for this deficit, we would then infer a larger value of σ 8 , hence the importance of accurately accounting for lensing magnification.In our case, however, we find that this effect can be largely neglected.
In Section 5.2.5, we will present an alternative test, using cross-correlations with a lowredshift sample, to validate our inferred value of s in comparison with the value found for other quasar samples.

The redshift distribution
The redshift distribution p(z) of the samples under study is a crucial element of the theory prediction, and we must both estimate it and propagate the uncertainties of that estimate to the final inferred cosmological parameters.In this work we make use of two different estimates of the redshift distribution.
Our fiducial estimate of p(z) is based on the direct calibration (DIR) method of [78].In this approach, a cross-matched spectroscopic sample is re-weighted to reproduce the magnitude/colour distribution of the target sample.The p(z) is then simply estimated as a weighted histogram of the spectroscopic redshifts.In our case, we make use of the crossmatched SDSS spectroscopic sample used to determine the Quaia redshift estimate z Quaia , and re-weight it in the 5-dimensional space of magnitudes in the (G, BP, RP, W 1, W 2) bands using a k-nearest-neighbours algorithm.This spectroscopic sample contains a total of 243,206 sources.The weight of each spectroscopic quasar is proportional to the number of objects in Quaia within a sphere in magnitude space containing the 10 nearest spectroscopic neighbors (for a more detailed description of the method see e.g.Section 2.2 of [79]).We verified that the re-weighted spectroscopic catalog follows the same magnitude distribution as Quaia in all four bands.We estimate the uncertainties in the resulting redshift distribution using jackknife resampling with 100 samples.The estimated p(z)s are shown in Fig. 4 as points with error bars for the two redshift bins used in our analysis.We propagate the uncertainties in the DIR redshift distribution making use of the approximate analytical marginalisation method proposed in [80,81] (see also [82]).To account for any potential mis-calibration of the p(z) uncertainties, we verified that our final results did not change significantly when multiplying the jackknife errors by up to a factor of 4.
One potential drawback of of the DIR-calibrated redshift distributions is the presence of stochastic noise associated with the specifics of the spectroscopic sample used (sample variance, shot noise, selection effects).If this is not correctly captured by the calibrated uncertainties (estimated via jackknife resampling above), this noise may lead to spurious shifts in the final inferred parameters.As an alternative to the DIR p(z), we estimate the redshift distributions of both bins via PDF stacking.The Quaia photometric redshifts z Quaia carry an associated uncertainty σ z .Assuming a Gaussian error model, we estimate the sample redshift distribution by adding one Gaussian probability distribution function (PDF) N (z Quaia , σ 2 z ) for each quasar in the sample.The resulting redshift distribution is shown as solid lines in Fig. 4. We will only use this alternative estimate of the p(z) to quantify the sensitivity of our final results to mis-modelling of the true underlying redshift distribution (which, we will show, is small), and therefore, unlike in our fiducial case, we will not associate any uncertainties to it or marginalise over them.
As Fig. 4 shows, the redshift distributions exhibit non-negligible tails beyond the threshold z Quaia , although both estimates of the p(z) are in qualitative agreement as to their extent.The final cosmological constraints are sensitive to the width of these tails, as they also control the amplitude of the p(z), and hence it is important to ensure that they are correctly modelled.In Section 5.2.3 we will make use of the cross-correlation between quasars in both bins as a sensibility check for the presence and extent of these tails.

Bias evolution
Although calibrating the true redshift distribution is important in order to understand the samples under study, the power spectra (C gg , C gκ ) are in fact only sensitive to the combi-nation b(z)p(z).In other words, the final cosmological constraints depend not only on the redshift distribution, but also on the evolution of the sample's bias with redshift.
The redshift evolution of the quasar clustering amplitude has been the subject of several works in the literature, including spectroscopic quasars in 2dF [83][84][85], SDSS [38,[43][44][45][86][87][88][89][90] and DESI [46,91], and photometrically-selected quasars [92][93][94][95].Although there is significant scatter, especially at high redshifts, the general consensus is that quasar bias evolves strongly with redshift.Thus, although assuming a constant bias within each of the redshift bins studied here would potentially be a more agnostic approach, ignoring the potential evolution within each bin could lead to significant shifts in the final parameter constraints (although in Section 5.2.2we will show that our final constraints are rather resilient to this effect).Thus, in our fiducial analysis we will assume that the quasar bias within each redshift bin i takes the form with b i g a free parameter (one in each bin).b i g = 1 would recover the best-fit bias model found by [38] from the eBOSS quasar sample.To quantify the impact of the assumed level of bias evolution encoded in this model, in Section 5.2.2 we will explore the constant-bias case, and a steeper evolution model, used in [93].
Another method to account for the redshift evolution of the sample bias is the so-called clustering redshifts approach.This technique [15,[96][97][98] is able to precisely constrain the combination b(z)p(z) from the cross-correlations of the target galaxy sample with a series of spectroscopic (or spectrophotometric) samples with well-determined redshifts.In Section 5.2.3 we will make use of clustering redshifts to validate the robustness of our constraints to uncertainties in the combination b(z)p(z) in a model-independent way.

The CMB κ map
We rotate the harmonic coefficients of the GMV CMB lensing convergence of the PR4 release into Equatorial coordinates, and truncate them to a maximum multipole max = m max = 3N side − 1, before transforming it into a HEALPix map of N side = 512 through an inverse spherical harmonic transform.We use the sky mask provided with the same release, rotate it in real space to Equatorial coordinates, apodize it with a C1 kernel (see Section IIIB of [99]), and downgrade it to the same angular resolution.
Compared to previous releases, the PR4 lensing map has a non-trivial normalization that accounts, in a more optimal patch-dependent way, for the effect of anisotropic noise in the CMB maps, as well as its propagation in lensing reconstruction.However, simplifications in the response to the true sky signal and noise can introduce small biases to the normalization, in particular close to the mask boundaries.We compute the multiplicative correction to the estimator normalization (the so-called Monte Carlo (MC) correction) using the simulations associated with the PR4 release (which include realistic sky and noise modelling) measuring the ratio of the power spectrum of the input lensing map spectra to the cross-correlation between the input and the reconstructed lensing maps (see also discussion in [46,100]).We compute the power spectra involved in this ratio on the Planck DR4 lensing mask with a binning of δ = 5 and interpolate to the relevant effective multipoles of our chosen binning scheme.Such MC correction is of the order of ∼5% on small scales and only becomes significant on the largest scales 10, which have a minor statistical weight in our analysis.A mismatch in the normalization might also affect the accuracy of the subtraction of the mean-field anisotropy induced by the mask.This mean field becomes important in the CMB lensing autospectrum analysis, and dominates on the largest scales.However, it only contributes to the covariance of the spectrum in the cross-correlations with external tracers and does not introduce any bias.The use of 600 end-to-end Monte-Carlo simulations released with the NPIPE maps to estimate the mean field of the PR4 lensing estimator allowed us to reduce the uncertainty connected to this effect compared to the PR3 release.As such we do not propagate further uncertainties connected to the mean field.

Angular power spectra and covariances
To calculate all auto-and cross-power spectra we make use of the pseudo-C approach as implemented in NaMaster9 [75].The estimator is described in detail in [75,101], and we only provide a brief description here.The pseudo-C estimator starts by assuming that all observed maps ã are simply given by a locally-weighted version of the true underlying signal maps a where the weights map w a is commonly called the "mask".In the simplest case, we may consider a binary mask, where w a = 1 for observed pixels, and w a = 0 for unobserved ones.In practice, the mask can be used as an inverse-variance weight, which down-weights regions of the sky where the uncertainty on the true signal is high.This motivates us using the selection function itself as the mask for the quasar overdensity map, since it tracks the observed number density of the sample which, in the Poisson limit, corresponds to the inverse local variance of δ g .A product in real space (Eq.4.7) corresponds to a convolution in harmonic space, and the harmonic coefficients of ã and a are related to one another via The pseudo-C of two observed maps (a, b) is defined as It is then straightforward to show that the ensemble average of Cab is related to the true undelying power spectrum via where Ñ ab is the pseudo-C of the noise components of both fields.M ab is the so-called mode-coupling matrix, and encodes the statistical coupling between different harmonic modes caused by the presence of a mask.M ab depends solely on the pseudo-C of the masks of both fields, and can be calculated analytically through an O( 3 max ) operation.An unbiased estimator for the true power spectrum can thus be constructed by estimating and subtracting the noise bias from Cab , and multiplying the result by the inverse of the mode-coupling matrix.To improve the condition of M ab before inversion, and to reduce the size of the resulting data vector, it is often common to bin the pseudo-C into finite bands of , commonly called bandpowers.The estimated bandpowers Ĉab q (where q labels each bandpower), can then be related to the true underlying power spectrum via where the bandpower window functions F ab q encode the effects of mode-coupling, binning, and inversion of the binned mode-coupling matrix.We bin all power spectra used here into bandpowers of width ∆ =30 between = 2 and = 3N side − 1.
Out of all power spectra used in our analysis, only the quasar auto-correlations have a non-zero noise bias due to shot noise.Assuming Poisson statistics, this bias can be estimated analytically as where w p is the sky average of the selection function, and NΩ ≡ N /Ω pix is the angular number density (in units of srad −1 ), with the mean number of objects per pixel given by Eq. 4.2, and Ω pix the solid angle of each HEALPix pixel.In Section 5.2.2 we will explore the impact of residual non-Poissonian shot noise in the data.All power spectra are corrected for the effects of the pixel window function as described in [75], although we find the impact of this correction to be negligible on the scales used here.Besides this, as described in the previous section, the quasar-κ cross-correlation is corrected for the CMB lensing transfer function.Finally, as discussed in Section 4.1.1,we make use of linear deprojection as implemented in NaMaster, ignoring the negligible deprojection bias it incurs in our case.
We calculate the covariance matrix of the estimated power spectra analytically, following the methodology described in [102].The method assumes that all fields involved are Gaussian, which is sufficiently accurate in our case, since all fields under study are largely noise-dominated.The approach described in [102] has been thoroughly validated, and accurately accounts for the effects of mode-coupling caused by the mask.To achieve this, the method only assumes that the true input power spectra are smooth functions of that do not vary significantly compared to the width of the mode-coupling kernel caused by the masks involved (the so-called Narrow-Kernel Approximation, NKA).Given the broad sky coverage of the probes we study, this approximation is only mildly inaccurate on the largest scales, which we discard in our analyses in any case.

Likelihood and priors
To infer free parameters of our model θ, given our data d, we make use of a Gaussian (multivariate normal) likelihood: where t( θ) is the theoretical prediction for d, C is the covariance matrix of d, and K is an arbitrary normalisation constant that does not depend on θ.The posterior distribution p( θ|d) is calculated from the likelihood above using Bayes' theorem log p( θ|d) = log p(d| θ) + log p( θ) + K , where p( θ) is the prior.

Parameter Prior
Parameter Prior Table 1.Free parameters and priors used in our model.We use BAO data from BOSS and eBOSS [103,104] to put a joint prior on Ω m and h.( * )The shot noise amplitudes A i SN are fixed to their Poisson level in our fiducial analysis, and we vary them with a 10% prior to test the dependence of our results on this assumption.( † )In Section 5.2.7 we explore the impact of imposing a prior on the horizon scale at recombination from the positions of the acoustic peaks in the CMB power spectrum.Effectively, this puts a tight constraint on the combination Ω m h 3 as shown on the table.
In our case, d is a collection of power spectra, comprising the quasar auto-correlations of the redshift bins under study, and their cross-correlations with the CMB lensing convergence.C is the covariance matrix of these measurement, calculated as described in the previous section, and the model used to calculate the theoretical prediction t was presented in Section 3. We only include bandpowers in d satisfying the following scale cuts: • We discard the first bandpower ( < 32) in all power spectra, where the impact of systematic contamination was found to be most relevant (see Section 5.2.4).
• We impose a high-cut at max = k max χ(z), where z is the mean redshift of the bin under study, and k max is the 3D wavenumber parametrizing our small-scale cut.In our fiducial results we choose k max = 0.15 Mpc −1 , corresponding to max = (497, 827) in the two redshift bins under study.The main motivation for this cut is ensuring the validity of the linear bias model used in our theoretical predictions.We will study the impact of this choice in Section 5.2.1).
With these scale cuts, we are left with 16 and 27 bandpowers for power spectra involving the Low-z and High-z quasar redshift bins, respectively, and thus the fiducial size of the data vector combining all auto-and cross-correlations is N d = 86.It is worth noting that the Gaussian likelihood of Eq. 4.13, with a parameter-independent covariance, has been found to be an excellent approximation to the true likelihood of power spectra on scales 10 [105,106].
Within the ΛCDM model, the main cosmological parameters that affect the power spectra considered here are the total non-relativistic matter density parameter Ω m , the amplitude of linear matter fluctuations σ 8 , and the expansion rate h, and we will vary all three of them in our likelihood.We will fix the baryon density parameter, and the scalar spectral index to the best-fit values found by Planck: Ω b h 2 = 0.0224, n s = 0.9665 [32].Our data is not able to fully break the degeneracy between Ω m and h and hence we will impose a prior on both background parameters using low-redshift baryon acoustic oscillations (BAO) data.In particular, we include three measurements of the angular and radial BAO scales from Luminous Red Galaxies (LRGs) at redshifts z = 0.38, 0.51, 0.698 obtained by the BOSS and eBOSS surveys in the DR12 and DR16 of the Sloan Digital Sky Survey [103,104].In Section 5.2.7 we will explore the sensitivity of our final constraints to this choice of prior.
Besides the cosmological parameters, we consider the two linear bias parameters {b 1 g , b 2 g } with flat priors in the range b i g ∈ [0.1, 3.0].We assume that the quasar bias evolves as in Eq. 4.6 in our fiducial analysis, and explore the impact of this assumption on our final results in Section 5.2.2.We will also quantify the impact of non-Poissonian shot noise on  the final parameter constraints by freeing-up the amplitude of the shot-noise contribution to the quasar auto-correlations.This will add two new free parameters (A 1 SN , A 2 SN ), one for each redshift bin, with a 10% Gaussian prior centered around 1. Since both parameters enter linearly in the theoretical prediction, we will marginalise over them analytically.Table 1 lists all the parameters of the model, and their priors.
As described in Section 4.1.3we marginalise over the uncertainties in the calibrated redshift distribution.We do so following the analytical approach of [80], which results in a modification of the power spectrum covariance matrix, without adding any new free parameters to the model.
All theoretical calculations were carried out using the Core Cosmology Library10 (CCL [107]).The linear matter power spectrum was computed with CLASS [108], and we used the HALOFit prescription of [109] to connect it with the non-linear power spectrum.We sample the posterior distribution using a the Metropolis-Hastings Monte-Carlo Markov Chain (MCMC) algorithm as implemented in Cobaya [110].We determine chains to be converged when the Gelman-Rubin statistic reaches R < 0.01.We will present parameter constraints by quoting the peak of the marginalised 1D distribution, together with its 68% asymmetric confidence-level interval.

Fiducial cosmological constraints
In our fiducial setup, we consider the G < 20.5 Quaia sample split into two redshift bins (see Fig. 4), with DIR-calibrated redshift distributions, fixed magnification bias slope parameters as in Eq. 4.3, assuming a linear bias evolving as in Eq. 4.6, and purely Poisson shot noise.We use only scales with k < k max = 0.15 Mpc −1 , and impose a BAO prior on (Ω m , h).
Fig. 5 shows the angular power spectra measured in both bins (C gg in black, C gκ in red), together with the best-fit theoretical prediction (blue lines).The shaded bands show the angular scales discarded in our analysis.The bottom panels of both figures show the residuals with respect to the best-fit prediction normalised by the 1σ uncertainties.We find that the best-fit model provides a reasonable fit to the data, with a χ 2 = 96.9.With N d = 86 data points and N p = 5 free parameters, the associated probability to exceed (PTE) is p = 0.11, which corresponds to an adequate fit.
The cross-correlation between Quaia and the Planck CMB lensing map is detected, within the scales used in this analysis, with a signal-to-noise ratio of 16.9 and 13.6 in the Low-z and High-z bins respectively, for a total detection significance of 21.7σ.This is actually stronger than the significance of the clustering auto-correlation signal itself (9σ and 8.4σ respectively), due to the impact of shot noise.Hence the cross-correlation carries a larger statistical weight in our final cosmological constraints.
The constraints on the cosmological parameters (σ 8 , Ω m , h, S 8 ) found with this fiducial setup are shown in blue in Fig. 6.As explained in Section 3.4, we are able to break the degeneracy between σ 8 and Ω m , measuring them independently: σ 8 = 0.766 ± 0.034, Ω m = 0.343 +0.017 −0.019 . (5.1) The figure also shows, in red and black, the constraints on the same parameters found by the DES Y3 3×2-point analysis [3] MagLim sample, and by Planck [32], respectively.Our results are in reasonable agreement with the Planck measurements, with only a mild tension at the 1.3σ and 1.4σ-level for σ 8 and Ω m respectively.The constraints are also fully compatible with those found by DES Y3.It is worth placing our constraints in the context of measurements obtained by other experiments.Fig. 7 shows our constraints on σ 8 , Ω m , and S 8 , together with those found in some of the latest analyses of galaxy clustering, cosmic shear, and CMB data.Our constraints are in broad agreement with those found by most of these experiments.We find that we are able to place a competitive constraint on σ 8 , unmatched by cosmic-shear experiments due to its degeneracy with Ω m , and comparable to that found by 3 × 2pt, 5 × 2pt, and lensing tomography analyses, although with a significantly smaller sample size in the case of Quaia.Our constraint on S 8 , in turn, is less stringent than weak lensing and 3 × 2pt analyses, but comparable with those found by full-shape power spectrum analyses in spectroscopic galaxy surveys [111,112]., Ω m , S 8 found in this analysis (salmon band), and found in other datasets.The points shown correspond to cosmic-shear analyses (marked γγ [66][67][68]), 3 × 2pt and 5 × 2pt data [2,3,25,26], lensing tomography [19,21], full-shape power spectrum analyses [111,112], the CMB lensing power spectrum + BAO data [34,48], and CMB data [32,33].The Planck constraints are shown as a blue vertical band.
As highlighted in Section 3.4, this is one of the highest-redshift constraints on the amplitude of matter fluctuations.That it agrees reasonably with probes at both higher and lower redshifts is therefore a non-trivial result that confirms the ability of the simple ΛCDM model to describe the large-scale properties of the Universe across most of its history.

Robustness tests
In this section we study the robustness of these results against the various analysis choices made, and potential sources of systematic uncertainty.Table 2 lists the final parameter constraints found for the different alternatives explored here, together with our fiducial measurements.The constraints on (σ 8 , Ω m , S 8 ) from these cases are also shown in Fig. 8.In all these cases we find a reasonably good fit to the data with χ 2 PTE values ranging between 0.04 and 0.5.

Internal consistency
We begin by exploring whether we obtain consistent results from different sub-sections of our data: • We explore the impact of our small-scale cut (k max = 0.15 Mpc −1 ) on the final constraints by repeating our analysis imposing both more and less conservative cuts.Rows 2 and 3 in Table 2 and Fig. 8 show the constraints found for k max = 0.1 Mpc −1 and k max = 0.2 Mpc −1 .In both cases we recover constraints that are compatible with our fiducial measurements within ∼ 0.5σ, with a ∼ 6% increase and reduction in the final uncertainties for the more and less conservative cuts respectively.Our results are thus robust to this choice.This also reinforces the validity of the simple linear bias model used in our analysis, since the presence of large non-linear corrections would lead to significant changes in the final constraints when adding or removing smaller scales.
• We repeat the analysis for a brighter quasar sample, cutting at G < 20.This brighter sample is potentially less sensitive to sky contamination and photo-z misestimation than the fiducial catalog used in our analysis.For this sample, the magnification bias slope is given by Eq. 4.4.As shown in row 4 of Table 2 and Fig. 8, the resulting parameter constraints are fully compatible with the fiducial case, with a ∼ 30% increase in the final uncertainties, caused by the higher shot noise of this sample.
• We obtain cosmological constraints from all Quaia sources combined into a single redshift bin (with the redshift distribution shown in black in Fig. 4).This test allows us to check for the impact of any potential mis-estimation of the intermediateredshift tails of the redshift distributions in the 2-bin case, as well as increased sky contamination caused by the additional cut in redshift.The results, shown in row 5 of Table 2 and Fig. 8, are largely compatible with our fiducial constraints, with a shift of less than 1σ in both σ 8 and Ω m , and barely any change in the final uncertainties.
• Finally, we repeat our analysis for each of the two redshift bins separately.The   2. Together with the results found in the different robustness described in Section 5.2, we show the constraints found in our fiducial analysis (first row, and salmon band), and those found by Planck (last row and blue band).
constraints are shown in rows 6 and 7 of Table 2 and Fig. 8, and the corresponding 2D contours are also shown separately in Fig. 9.The cosmological parameters inferred from each redshift bin are compatible with one another, although the High-z bin favours a value of σ 8 that is 2.2σ lower than the Low-z bin, 1.5σ lower than our fiducial constraints, and 2.5σ lower than Planck.We thus find that the slightly lower value of σ 8 found in our analysis in comparison with Planck is mostly driven by the higher redshift bin.We will therefore scrutinise the results obtained from this bin further in the following sections.

Bias model
Our fiducial constraints assume the specific redshift dependence of the quasar bias in Eq. 4.6.
For concreteness, we will label this Model B1 in this section.To quantify the dependence of our results on the assumed evolution of the quasar bias, we repeat our analysis for two other models: .Constraints on Ω m , σ 8 , and S 8 found in our fiducial analysis (blue contours), in comparison with those found using only each of the two Quaia redshift bins independently (solid orange and red contours for the Low-z and High-z bins respectively).The Planck constraints are shown in black.The higher redshift bin recovers a value of σ 8 that is smaller than that found by Planck at the 2.5σ level, and is hence the main driver of the slightly lower value of σ 8 found in our fiducial analysis.The dashed red contours show the constraints found with the High-z when cross-correlating with the polarizationonly CMB lensing map (which is more reliable against extragalactic foreground contamination in the CMB map), as presented in Section 5.2.6.The ∼1.5σ upwards shift in σ 8 observed is a hint that the low σ 8 found with the fiducial CMB κ map may be caused by extragalactic foreground contamination in that map.
• Model B2: the model of [93], which is significantly steeper, and likely ruled out by current data: • Model B3: a constant bias model (b(z) = b i g ), also ruled out by current data if applied at all redshifts.
The constraints found with these two alternative models are shown in rows 8 and 9 of Table 2 and Fig. 8.We observe that models with a steeper bias evolution lead to lower values of both σ 8 , and Ω m .Nevertheless, for both models, B2 and B3, which encompass the degree of  [83] (blue), BOSS [89] (green), and eBOSS [38] (red).The grey line shows the fiducial bias evolution model used in our analysis, with the amplitude parameter b g fixed to 1, to reproduce the best-fit of [38] for eBOSS.The bias values measured in this analysis for the Quaia sample are shown in black for the three different bias evolution models considered (models B1, B2, and B3, described in Section 5.2.2).The value found when using Tomographer to constrain the redshift dependence of the product b(z)p(z) is shown in orange.Our constraints are in rough agreement with previous measurements.The systematically higher value of the quasar bias found for Quaia is consistent with it being slightly brighter than the spectroscopic samples shown in the figure .evolution exhibited by existing measurements, the parameter shifts are mild (less than 1σ).Hence, our results do not depend strongly on this choice.
It is useful to place our constraints on quasar bias in the context of previous measurements.For the three bias evolution models studied, we calculate the effective bias in each of the two redshift bins as and assign to it an effective redshift given by the mean z in each bin weighted by the product b(z)p(z).The results are shown in Fig. 10 for all three models (black symbols), together with the measurements of [38,83,89] (blue, green, and red points), and our fiducial model (Eq.4.6, with b i g = 1).Reassuringly, our measurements are in rough agreement with other quasar samples.In particular, we find no strong evidence of the lower quasar bias at high redshifts (z 2) reported by previous analyses using CMB lensing cross-correlations [43,46,90,113].As discussed in Section 4.1.2,the impact of magnification bias, an effect that the CMB lensing cross-correlation is particularly sensitive to, is likely negligible in our sample.Finally, it is worth noting that, although our bias constraints lie somewhat above most previous measurements in both redshift bins, this is not entirely surprising, as our sample is different (brighter, in fact) from the spectroscopic samples targeted by previous measurements.
Our fiducial bias model is purely linear, and ignores contributions from higher-order terms.This assumption is likely valid given the relatively conservative scale cuts used [114][115][116], and we are further reassured by the robustness of our constraints to changes in these -25 - and prediction from our best-fit model.This cross-correlation is sensitive to the overlap between the redshift distributions of both bins, shown in Fig. 4. Our best-fit model provides a good fit to this cross-correlation, which was not used in that analysis.
scale cuts (as shown in Section 5.2.1).However, we have assumed that the stochastic term in Eq. 3.2 is completely described by Poisson statistics, and can be subtracted analytically.Since a stochastic term would affect the quasar auto-correlation on all harmonic scales, we test this assumption by repeating our analysis varying the amplitude of the shot-noise contribution to C gg with a 10% prior.This also allows us to test for any misestimation of the purely-Poisson noise level.The result is shown in row 10 of Table 2 and Fig. 8.The mean value of all parameters remains virtually unchanged, the only effect being a ∼ 10 − 20% increase in the final uncertainties.The data thus shows no preference for additional flat contribution to C gg , and its inclusion causes no shift on the final parameter constraints.

Redshift uncertainties
To test the robustness of our constraints to a mis-calibration of the quasar redshift distributions, we repeat our analysis assuming the p(z) calculated via PDF stacking, as described in Section 4.1.3.The results are shown in row 11 of Table 2 and Fig. 8.The constraints are virtually unchanged, which shows that our results are insensitive to mild changes in the redshift distributions (such as the height of the bump11 at z ∼ 2.25 shown in Fig. 4).Another important systematic related to the p(z) is a potential mis-estimate of the inter-bin tails of both samples.Driven by photometric redshift outliers, they control the amplitude of the redshift distribution, and hence the amplitude of the quasar auto-correlation.Fortunately, the cross-correlation between both bins is particularly sensitive to the details of the redshift overlap between them [20,117].This cross-correlation is shown in Fig. 11, together with the theoretical prediction for our fiducial best-fit model.On scales 30 ≤ ≤ 800, the cross-correlation is detected at the 3.2σ level, and our best-fit model provides a good fit to it (PTE = 0.51), without any additional parameters.This shows that the combination of best-fit bias parameters found in both bins, and the intermediate-redshift tails of their redshift distributions, are compatible with a complementary piece of the data vector that was not used in the likelihood.As shown in Section 5.2.1, the lower value of σ 8 we recover compared to Planck is largely driven by the auto-and cross-correlation of the High-z quasar redshift bin.Since both the shape of the redshift distribution, and the evolution of the linear bias may depart from the fiducial models used here, particularly at high redshifts.(e.g.due to photo-z outliers or a steeper-than-expected b(z) relation), we perform one further check to test the robustness of the low σ 8 value measured for that sample.We made use of the clustering redshifts technique, as implemented in Tomographer12 , to measure the product b(z)p(z) up to an unknown normalization factor.Tomographer uses a spectroscopic reference sample from SDSS comprising 2 million sources at redshifts z 3. Through its online interface, the code transforms an input catalog into a HEALPix map at resolution N side = 2048, imposes a high-pass filter removing scales larger than ∼ 2 • to mitigate the impact of sky systematics, and calculates the correlations between the map and the reference sample in thin redshift bins on scales r ∈ [2, 8] Mpc.In combination with the reference sample auto-correlations, this is then used to provide constraints (mean and standard deviation) of b(z)p(z) in each spectroscopic redshift bin.
We use Tomographer to estimate b(z)p(z) for the High-z bin in our sample.The redshift coverage of the reference spectroscopic sample in Tomographer prevents a reliable determination of b(z)p(z) beyond z ∼ 3, and thus, for this test, we impose an additional cut in the catalog at z Quaia < 3.As with the DIR redshift distribution, we marginalise over the uncertainties in b(z)p(z) obtained by Tomographer using the analytical marginalisation method of [80].The cosmological constraints obtained from this test are shown in line 12 of Table 2 and Fig. 8. Compared with the constraints found with the High-z bin using the DIR redshift distribution and our fiducial bias evolution template, we observe only a small upwards shift in σ 8 (δσ 8 0.02, less than 0.5σ).Hence, the data still favours a lower value of σ 8 compared to Planck, although at a lower significance of 2.1σ.This slightly better agreement with Planck could be simply a statistical fluctuation, it could be caused by the slightly lower effective redshift of the sample, or could be evidence of a departure from the quasar bias evolution model assumed here.The effective bias obtained from our MCMC chains for this bin using clustering redshifts is shown in orange in Fig. 10, and agrees well with the values found for the three other bias evolution models explored.Fig. 12 shows the values of b(z)p(z) recovered from Tomographer in orange, together with the DIR-calibrated p(z) multiplied by the B1 bias evolution model (Eq.4.6).Both predictions for the product b(z)p(z) are qualitatively similar in terms of the redshift support and high-redshift tails.The small difference observed in the final constraints might be caused by the enhancement in the bump at z ∼ 2.25 in the DIR distribution caused by the strong bias evolution at that redshift.Overall, however, we find that both predictions yield highly compatible constraints for the High-z bin, and hence, when combined with the results presented in this section and Section 5.2.2, it is unlikely for the lower value of σ 8 found in this bin to be caused by mis-modelling of the sample's redshift distribution or its bias evolution.

Sky contamination
The presence of residual sky contaminants in the quasar overdensity maps could cause a bias in the final cosmological constraints.For instance, if purely additive, contamination would raise the amplitude of the quasar auto-correlation, leaving the cross-correlation intact.When fitting for bias and σ 8 , this would then lead to an upwards shift in the quasar bias (to accommodate the higher auto-correlation) and, correspondingly, a downwards shift in σ 8 (to compensate for the larger bias in the cross-correlation).To quantify the robustness of our results to the presence of residual contamination, we repeat the analysis skipping the final linear deprojection step when computing all power spectra.As a reminder, this step was meant to ensure that any small residual contamination from the known sky systematics was removed from the maps.The top panel of Fig. 13 shows the four angular power spectra used in our analysis with (points) and without (crosses) performing the final linear deprojection step.We see that linear deprojection has a significant effect on the first bandpower, which is discarded in our analysis, and less of an effect on smaller scales.The constraints found without linear deprojection are shown in row 13 of Table 2 and Fig. 8.We observe only negligible shifts in the best-fit parameters ( 0.1σ), and thus conclude that our results are likely robust against residual contamination from the known sky systematics.
To further understand this result, we compute the cross-correlation between the quasar overdensity maps and the 4 contaminant maps described in Section 4.1.1.The detection of a non-zero correlation would constitute evidence of residual contamination in the maps.We find that, without linear deprojection, the quasar overdensity maps exhibit large-scale residual contamination, particularly from dust extinction and stars.This contamination disappears after deprojection (unsurprisingly, since linear deprojection is precisely intended to cancel these correlations at the pixel level).An example of this is shown in the lower panel of Fig. 13 for the two Quaia bin-systematics combinations exhibiting the most contamination.Although linear deprojection modifies some of the cross-correlations on all scales, Figure 13.Top: quasar-quasar (black) and quasar-κ (red) angular power spectra for the two redshift bins used in this analysis (top and bottom subpanels).The measurements made using linear deprojection to remove residual sky contamination are shown as points with error bars, while the crosses show the results found without deprojection.The impact of residual contamination is mostly relevant for the first bandpower of the auto-correlation, which we discard in our analysis.Bottom: example cross-correlations between the quasar overdensity maps and the most relevant sky contaminants for each of the Quaia redshift bins.Results are shown with and without linear deprojection (points and crosses respectively), divided by their statistical uncertainties to better highlight significant departures from zero correlation.Evidence of contamination is mostly present in the first bandpower, in agreement with the results of the top panel.Left: redshift distribution of the High-z Quaia bin (points with error bars), and of the WISE-SuperCOSMOS sample we cross-correlate it with to validate our estimate of the magnification bias slope.Right: cross-correlation between the High-z Quaia bin and the low-redshift WISE-SuperCOSMOS sample.The solid lines show the predictions for the best-fit value of the magnification bias slope s (red), and for the value found by [46] in the DESI quasar sample (blue).The cross-correlation is able to rule out a value of s as low as that found in the DESI sample at ∼ 3.5σ.they are only significantly different from zero in the first bandpower without deprojection, in agreement with the results found above.

Magnification bias
As discussed in Section 4.1.2,the magnification bias slope of our fiducial quasar sample, measured with two different methods, is close to s 0.4, and hence our measurements are largely unaffected by magnification bias.To quantify the sensitivity of our final constraints to a potential mis-estimation of this quantity, we repeat our analysis fixing the number counts slope to the value found for the DESI quasar sample, s = 0.2764 [46].This value is in rough agreement with previous estimates for similar spectroscopic quasar samples (e.g. in BOSS and eBOSS).The result from this exercise, shown in row 14 of Table 2 and Fig. 8, leads to a < 1σ upwards shift in σ 8 , and almost no change in the Ω m constraint.It is worth noting that this value of s is 3σ away from our estimate for Quaia (see Section 4.1.2),and is therefore not supported by the data.Nevertheless, the impact on our constraints of such a significant misestimate of s is small.
As an additional test to validate our estimate of s for the High-z bin (which shows the largest tension with Planck), we measure the cross-correlation between this bin and a lowredshift sample of galaxies selected from the WISE-SuperCOSMOS catalog (WI-SC, [118]).In particular, we select 14.7 million galaxies with photometric redshifts z p ∈ [0.1, 0.3].The methods used to analyse this sample are described in detail in [9], and we will not repeat them here for brevity.As shown in the left panel of Fig. 14, both samples have negligible overlap in redshift, and thus, any detection of the cross-correlation between them would be largely driven by magnfication bias with s = 0.4.
Fixing cosmological parameters to the best-fit values found by Planck [32], we determine the linear bias of the WI-SC sample from its auto-correlation on large scales to be b g = 1.1, in good agreement with [5].Fixing the galaxy bias to this value, and the quasar bias to the value found for the High-z bin in our fiducial analysis, we use the cross-correlation between both for the polarization-only map (orange triangles).When the power spectrum is negative, we show its absolute value as empty circles with dashed error bars.The polarization-only data is systematically above the GMV estimator on scales 300, although with larger uncertainties, with all deviations at the level of ∼1σ.The blue solid line shows the best-fit theoretical prediction found in our fiducial analysis.In turn, the dashed blue line shows the theoretical prediction resulting from multiplying the fiducial best-fit curve by a free amplitude and fitting it to the difference between both power spectra.The best-fit curve found this way is ∼30% higher than our fiducial best fit, which is significant at the ∼2σ level.This translates into a ∼1.5σ upwards shift in σ 8 for this bin when using the polarization-only map, making it fully compatible with the Planck value.
samples to constrain the value of s.The right panel of Fig. 14 shows this cross-correlation, together with the prediction for the best-fit value of the number counts slope (see Eq. 5.4) in red, and that for the value of s found by [46] for DESI in blue.The value of s obtained from this cross-correlation is s(z > 1.47) = 0.440 ± 0.047, (5.4) and provides a good fit to the data (χ 2 = 26.7 for 26 degrees of freedom).This value is in good agreement with that found with two different methods in Section 4.1.2,whereas the DESI value lies ∼ 3.5σ away from it.It is thus unlikely for magnification bias to be the cause of the low σ 8 value found from the higher-redshift quasar sample.

Foreground contamination in CMB lensing
So far we have mostly explored the impact of various systematics affecting the Quaia sample.One of the most important sources of systematic uncertainty for cross-correlation studies, however, is contamination in the CMB lensing map, particularly from extragalactic foregrounds, tracing the same large-scale structure as the probe being cross-correlated.
To quantify the robustness of our results against this potential systematic, we computed the cross-correlation of the Quaia sample with an alternative map of the CMB convergence constructed using only polarization information (and therefore largely devoid of any extragalactic foreground contamination that is mainly unpolarized) and a PR3-like pipeline run on NPIPE maps.We note that lensing reconstruction based on polarization can be more prone to biases generated by Galactic foregrounds, such as Galactic dust emission, which might in turn correlate with uncorrected dust extinction in the Quaia sample.However, [119] showed that dust biases are mitigated to a negligible level if lensing reconstruction is performed on a CMB map cleaned of foregrounds through component separation (such as those used for PR3 and PR4 lensing analyses), even for data much deeper than Planck.The systematic deprojection performed on the Quaia sample would reduce this effect further.Using the methods described below we find no indication of contamination in the cross-correlation with the Low-z bin, and therefore we will limit our discussion to the High-z bin, which displays the largest deviation in σ 8 with respect to Planck (see Section 5.2.1).
Figure 15 shows the cross-correlation of the High-z bin with our fiducial generalised minimum variance (GMV) κ map in red, and with the polarization-only map κ Pol in orange.Although, in each bandpower, both power spectra are compatible within 1σ, we find that the polarization-only coss-correlation is systematically ∼ 30% higher on large scales ( 300).At the Planck sensitivity, the polarization-only reconstructed κ map is significantly noisier than the GMV estimator (the errors in the cross-correlation are between ∼2 and ∼4 times larger on the scales explored here), and hence we find only mild evidence that this shift is indeed caused by contamination in the κ map.To quantify this in a model-independent way, we first calculate the χ 2 of the difference between both power spectra with respect to a null model 13 .The PTE of the resulting χ 2 value is 0.37, and hence, by this metric, both power spectra are compatible with each other within 1σ.To quantify the evidence for a coherent deviation between both spectra, we instead take the theoretical prediction for the best-fit parameters found in our fiducial analysis for this cross-correlations, and multiply it by a free amplitude A Pol .We fit this 1-parameter model to the difference between both power spectra on scales 30 ≤ ≤ 600, finding: 1 + A pol = 1.296 ± 0.125. (5.5) We thus find weak evidence, at the 2.3σ level, that the power spectrum computed with the polarization-only κ map is ∼ 30% higher than that for the minimum-variance map.A more in-depth investigation of this issue, making use of the Planck CMB lensing simulations in the PR3 release, is described in Appendix A. Finally, we propagate the impact of this potential systematics to our final constraints.For the High-z Quaia bin, using the polarization-only κ map, we find the cosmological constraints shown in row 15 of Table 2 and Fig. 8.These constraints are also shown as dashed blue lines in Fig. 9.We observe a ∼1.4σ upwards shift in σ 8 , which is now fully compatible with Planck within 1σ.Although the evidence is relatively weak (at the ∼2σ level), this provides an attractive explanation for the low amplitude of our fiducial cross-correlation between CMB lensing and the High-z bin with respect to the Planck prediction.The most likely source of contamination behind this effect would be the Cosmic Infrared Background (CIB), given its relevance at redshifts 2, and the fact that the resulting bias to the CMB lensing map and power spectrum is expected to be negative (see Fig. 23 of [47]).The cosmological constraints found using using the GMV κ map with the Low-z bin, and the polarization-only map with the High-z bin, are shown in row 16 of Table 2 and Fig. 8.The main effect is an upwards shift of σ 8 by ∼1σ, and a ∼30% increase in the final uncertainties caused by the significantly higher noise of the polarization-only map. ) from our analysis using a BAO prior (our fiducial case, in blue), and imposing a prior on the expansion rate today, from Planck (red) and from the SH0ES collaboration (orange).Although the constraints on σ 8 are only mildly affected by the choice of prior, the effect on Ω m is larger, shifting it by up to 1.3σ.

Impact of priors
As explained in Section 4.4, our data is not sufficient to break the degeneracy between Ω m and h, and we have relied on low-redshift BAO data to pin down Ω m .To explore the dependence of our results on this choice, we have repeated our analysis with three alternative prior choices.First, we consider the option of adding, to our fiducial BAO prior, a constraint from the angular sound horizon scale at recombination, θ LSS , from the positions of the peaks in the CMB power spectrum.This is one of the cleanest CMB observables, purely geometrical by nature, and is remarkably well measured (θ LSS = 1.04109 ± 0.00030 [32]).Following [19], we include this constraint by imposing a Gaussian prior on the combination Ω m h 3 = 0.09633 ± 0.00029 (and varying this quantity instead of h in our MCMC chains).The result is shown in row 17 of Table 2 and Fig. 8.The addition of this prior leads to a significant tightening of the constraint on Ω m , which now agrees fully with Planck with almost matching uncertainties.The constraint on σ 8 , on the other hand, remains almost unchanged, with only a 0.3σ downwards shift (away from Planck), and a 5% tightening of its uncertainty.It is worth noting that, in spite of the relatively large change in the posterior distribution of Ω m , the best-fit model in this case is still a good fit to the data (the χ 2 only increases by ∼2 with respect to our fiducial best fit).
We then study the impact of replacing the BAO prior by a Gaussian prior directly on h.We consider two possibilities: a Planck prior of the form h = 0.674 ± 0.005, and a supernova prior from the SH0ES team [120]: h = 0.7304 ± 0.0104.The results are shown in rows 18 and 19 of Table 2 and Fig. 8.The 2D constraints in the (σ 8 , Ω m , S 8 ) plane are additionally shown in Fig. 16.Using the Planck prior does not lead to significant changes in the final constraints beyond a broadening of the Ω m uncertainty.Interestingly, using the SH0ES prior instead, shifts Ω m and σ 8 closer to the values preferred by Planck (by ∼ 1.3σ in the case of Ω m ).The origin of the H 0 tension thus has an impact on the interpretation of our constraints (as would be the case for any data set displaying a degeneracy with h).

Conclusions
In this work, we have presented a study of the tomographic clustering of quasars in the Quaia sample, a quasar catalog selected from the third Gaia data release in combination with unWISE [49].The analysis is based on the angular auto-power spectrum of Quaia in two redshift bins, centered at z = 1.0 and z = 2.1, as well as the cross-correlation with CMB lensing data from Planck.
We find that, thanks to the large range of scales made available by Quaia's large cosmological volume, we are able to break the degeneracy between Ω m and σ 8 that is usually present in tomographic large-scale structure analyses with lower-redshift samples.We are able to constrain σ 8 at a level that is competitive with (and in some cases better than) the precision obtained by state-of-the-art weak lensing and galaxy clustering experiments (see Fig. 7).Furthermore, our measurement is highly complementary to those datasets, since it makes use of significantly higher-redshift data, and mostly linear scales, where the modelling of astrophysical systematics is, arguably, cleaner.Our constraints on the amplitude of matter fluctuations, σ 8 = 0.766 ± 0.034, and on the fractional abundance of non-relativistic matter, Ω m = 0.343 +0.017 −0.019 , are compatible with Planck at the 1.4σ level.We also find no evidence of the so-called "S 8 tension", although our data is less sensitive to that parameter combination than weak lensing experiments.We have shown that our measurement is robust against a large suite of potential systematic effects and analysis choices, including the assumed scale cuts, the assumed redshift evolution of the quasar bias, the impact of non-Poissonian shot noise, uncertainties in the redshift distribution of our samples, contamination from residual sky systematics in the galaxy overdensity map, a mis-estimation of the magnification bias slope, and the choice of priors on other cosmological parameters.Furthermore, we find that our measurements of the quasar bias is in reasonable agreement with previous studies.
When extracting constraints from each of the two redshift bins of the catalog independently, we find that the higher redshift sample favours a value of σ 8 that is lower than the Planck measurement at the ∼2.5σ level.This is caused by the relatively low amplitude of the quasar-lensing cross-correlation in that bin, and drives the lower value of σ 8 found in our fiducial analysis.This result is robust against all sources of uncertainty in the modelling of the Quaia sample we have explored (quasar bias evolution, redshift distribution calibration, magnification bias).The explanation for this low amplitude thus most likely lies in the CMB lensing map itself.In particular, by reanalysing the data using an alternative CMB κ map, constructed using only polarization information, we have found compelling (although not definitive) evidence of potential contamination from extragalactic foregrounds, which would lower the amplitude of this particular cross-correlation.The most likely foreground source responsible for this effect would be the Cosmic Infrared Background, given its relevance at z 2. We find that, using the polarization-only κ map, the amplitude of the quasar-κ crosscorrelation can grow by up to ∼20-30%, although the significantly higher noise level of this map prevents us from detecting this shift at more than ∼2σ.This correction is substantially larger than the nominal level of contamination due to CIB leakage in the CMB lensing power spectrum determined by Planck (at the level of ∼1%), but it is able to bring our results into complete agreement with the Planck constraints.
Fully understanding the presence and origin of this potential contamination would be highly relevant, since it could have a substantial effect on other cross-correlation studies involving high-redshift datasets (such as [43,46,90,121,122]).A better characterisation of the effect could be achieved using more sensitive CMB lensing estimators that are still robust against CIB contamination, such as including T E information, or using hardening methods [123].The use of high-resolution ground-based CMB lensing data [7,34] will also be able to shed light on this effect in the near future.
Another caveat of our analysis is the degeneracy with the value of the Hubble constant h, which we mitigated by making use of low-redshift BAO data.In the future, other spectroscopic quasar catalogs, such as the sample targeted by DESI [124], may be able to break this degeneracy independently, making use of their own measurements of the redshift-distance relation.We also showed that our final constraints, especially in the case of Ω m , are sensitive to the choice of prior on h (if one is used), and thus may be affected by the ongoing tension between high redshift CMB and large-scale structure data, and direct measurements from supernovae.
Our analysis has demonstrated the significant and complementary constraining power of high-redshift large-scale structure data sets covering large volumes.The characteristics of the Quaia catalog (comoving volume, redshift precision, homogeneity, reduced observational systematic effects enabled by space-quality data), make it ideal to carry out a number of other cosmological studies.These include constraining the level of large-scale homogeneity [125] and isotropy [126], and detecting ultra-large scale effects, such as primordial non-Gaussianity [127].The main limiting factor of this sample is its relatively low density, which significantly reduces the amount of information that can be extracted from small (but still only mildly non-linear) scales.This will likely be improved upon by future quasar catalogs from nextgeneration experiments such as DESI [124] or the Vera C. Rubin Observatory [128], which will be able to improve on the results obtained here.Future Gaia data releases may also enable improvements to the catalog in terms of redshift precision, quasar identification, and decontamination, which could allow for improved measurements.
We then assessed the significance of the observed shifts in the data bandpowers C gκ MV −C gκ Pol comparing the value obtained for the χ 2 performed on data with the χ 2 distribution obtained from simulations using only the multipoles in the 30 ≤ ≤ 600 range.We adopted the same approach for the χ 2 test for each bandpower observed for the data.We found a PTE of 6% for the full spectrum, supporting the conclusions derived with our data-based model for PR4, although with worse PTE approaching a 2σ-level tension as we found fitting the A pol parameter.When considering the χ 2 test per bandpower, we found that almost all the observed shifts in the data are consistent with simulations (PTEs above 5%).However, the bin at ≈ 100 and the one at ≈ 370 have a PTE of about 2%, hinting to some mild evidence of inconsistency with the hypothesis of bandpower shifts only due to statistical fluctuations.The feature at ≈ 100 in particular is the one driving the low (yet acceptable) PTE for the full χ 2 test.We also note that lowering the covariance estimated from our Monte Carlo approach by a factor of 20% to mimic the gain in noise obtained in PR4 GMV analysis compared to PR3 did not change these conclusions.We also repeated a similar test comparing the MV with TT-only reconstruction, the MV with the TT lensing reconstruction obtained explicitly nulling the tSZ contribution when performing component separation, and the TT with the TT-only tSZ-deprojected reconstruction.In all these cases we found that the observed shifts in the gκ bandpowers are consistent with the differences induced by statistical fluctuations, reinforcing the conclusion that, among the extragalactic foregrounds, the CIB emission rather than tSZ is more likely to drive the low amplitude of the gκ signal at high redshift.

Figure 2 .
Figure2.Quasar auto-spectrum (left panel) and cross-spectrum with the CMB lensing convergence (right panel) for the high-redshift bin used in our analysis, centered at redshift z = 2.1 (see Fig.4and Section 4.1).The measurements are shown as black points with error bars, with the best-fit model shown in red.The blue and yellow lines show the theoretical predictions after changing Ω m and σ 8 by 0.1 in both directions, respectively.On the redshifts and scales our data are sensitive to, the effect of both parameters is markedly different, allowing us to break the degeneracy between them.

Figure 3 .
Figure 3. Number counts slope (see Eq. 3.5) governing the impact of lensing magnification, for our fiducial G < 20.5 sample (black) and for the brighter G < 20 sample (red), as a function of redshift.The dotted lines show the value of s found for the full sample (i.e.not separated in redshift).The value of s in our fiducial sample, s ∼ 0.4, implies that the impact of lensing magnification on our analysis is negligible.

. 5 )Figure 4 .
Figure 4. Redshift distribution of the quasar samples studied in this work.The points with error bars show the DIR estimate of the redshift distributions for the two redshift bins used in our fiducial analysis (red and blue respectively).The solid lines show the same redshift distributions estimated via PDF stacking, as described in the main text.The PDF-stacked redshift distribution of the full sample is shown in dashed black (as used in Section 5.2.3).

Figure 5 .
Figure 5. Angular auto-power spectra of the quasar overdensity field (black) and cross-correlation with the Planck CMB lensing convergence map (red) in the two Quaia redshift bins (top and bottom panels for the Low-z and High-z bins respectively).The blue lines show the best-fit theoretical prediction found in our fiducial configuration.The bottom sub-plots in each panel show the residuals with respect to the best-fit prediction normalised by the 1σ uncertainties.The central white band shows the angular scales used in our analysis for each redshift bin.

Figure 9
Figure 9. Constraints on Ω m , σ 8 , and S 8 found in our fiducial analysis (blue contours), in comparison with those found using only each of the two Quaia redshift bins independently (solid orange and red contours for the Low-z and High-z bins respectively).The Planck constraints are shown in black.The higher redshift bin recovers a value of σ 8 that is smaller than that found by Planck at the 2.5σ level, and is hence the main driver of the slightly lower value of σ 8 found in our fiducial analysis.The dashed red contours show the constraints found with the High-z when cross-correlating with the polarizationonly CMB lensing map (which is more reliable against extragalactic foreground contamination in the CMB map), as presented in Section 5.2.6.The ∼1.5σ upwards shift in σ 8 observed is a hint that the low σ 8 found with the fiducial CMB κ map may be caused by extragalactic foreground contamination in that map.

Figure 10 .
Figure 10.Measurements of the bias b(z) for various quasar samples: the 2dF QSO Redshift Survey[83] (blue), BOSS[89] (green), and eBOSS[38] (red).The grey line shows the fiducial bias evolution model used in our analysis, with the amplitude parameter b g fixed to 1, to reproduce the best-fit of[38] for eBOSS.The bias values measured in this analysis for the Quaia sample are shown in black for the three different bias evolution models considered (models B1, B2, and B3, described in Section 5.2.2).The value found when using Tomographer to constrain the redshift dependence of the product b(z)p(z) is shown in orange.Our constraints are in rough agreement with previous measurements.The systematically higher value of the quasar bias found for Quaia is consistent with it being slightly brighter than the spectroscopic samples shown in the figure.

2 ×10 − 7 Figure 11 .
Figure 11.Cross-correlation between the two Quaia redshift bins (black points with error bars), and prediction from our best-fit model.This cross-correlation is sensitive to the overlap between the redshift distributions of both bins, shown in Fig.4.Our best-fit model provides a good fit to this cross-correlation, which was not used in that analysis.

Figure 12 .
Figure 12.Product of the redshift distribution and linear bias of the High-z Quaia redshift bin.The blue points show the DIR estimate of the redshift distribution multiplied by our fiducial bias evolution model (Eq.4.6).The orange points show the measurement of b(z)p(z) made using the clustering redshifts approach of Tomographer.
Figure 14.Left: redshift distribution of the High-z Quaia bin (points with error bars), and of the WISE-SuperCOSMOS sample we cross-correlate it with to validate our estimate of the magnification bias slope.Right: cross-correlation between the High-z Quaia bin and the low-redshift WISE-SuperCOSMOS sample.The solid lines show the predictions for the best-fit value of the magnification bias slope s (red), and for the value found by[46] in the DESI quasar sample (blue).The cross-correlation is able to rule out a value of s as low as that found in the DESI sample at ∼ 3.5σ.

Figure 15 .
Figure 15.Cross-correlation between the High-z Quaia redshift bin, and maps of the CMB convergence.Results are shown for the fiducial generalised minimum variance lensing map (red points), and for the polarization-only map (orange triangles).When the power spectrum is negative, we show its absolute value as empty circles with dashed error bars.The polarization-only data is systematically above the GMV estimator on scales 300, although with larger uncertainties, with all deviations at the level of ∼1σ.The blue solid line shows the best-fit theoretical prediction found in our fiducial analysis.In turn, the dashed blue line shows the theoretical prediction resulting from multiplying the fiducial best-fit curve by a free amplitude and fitting it to the difference between both power spectra.The best-fit curve found this way is ∼30% higher than our fiducial best fit, which is significant at the ∼2σ level.This translates into a ∼1.5σ upwards shift in σ 8 for this bin when using the polarization-only map, making it fully compatible with the Planck value.

Table 2 .
Constraints (marginalised mean and 68% confidence level interval) on (σ 8 , Ω m , S 8 ) for the different cases studied here.Our fiducial result, described in Section 5.1, is shown in the first row.All other rows show the results found in the various robustness tests described in Section 5.2.