Growth history and quasar bias evolution at z < 3 from Quaia

We make use of the Gaia-unWISE quasar catalogue, Quaia, to constrain the growth history out to high redshifts from the clustering of quasars and their cross-correlation with maps of the Cosmic Microwave Background (CMB) lensing convergence. Considering three tomographic bins, centred at redshifts z̅i = [0.69, 1.59, 2.72], we reconstruct the evolution of the amplitude of matter fluctuations σ 8(z) over the last ∼ 12 billion years of cosmic history. In particular, we make one of the highest-redshift measurements of σ 8 (σ 8(z = 2.72) = 0.22 ± 0.06), finding it to be in good agreement (at the ∼ 1σ level) with the value predicted by ΛCDM using CMB data from Planck. We also used the data to study the evolution of the linear quasar bias for this sample, finding values similar to those of other quasar samples, although with a less steep evolution at high redshifts. Finally, we study the potential impact of foreground contamination in the CMB lensing maps and, although we find evidence of contamination in cross-correlations at z ∼ 1.7 we are not able to clearly pinpoint its origin as being Galactic or extragalactic. Nevertheless, we determine that the impact of this contamination on our results is negligible.


Introduction
The evolution of the amplitude of density fluctuations across cosmic time is a key cosmological observable, sensitive to the energy content of the Universe, as well as the laws of gravity that govern the growth of structure.The increase in the quantity and variety of cosmological datasets in the last two decades has provided us with a wide range of probes able to reconstruct this growth [1].On the one hand, probes sensitive to the peculiar velocity field, such as redshift-space distortions [2,3] or velocity surveys [4], are able to measure the growth rate with remarkable precision, constraining the speed at which the structure grows over time.On the other hand, gravitational lensing directly probes the amplitude of matter inhomogeneities.For instance, cosmic shear data and the lensing of the Cosmic Microwave Background (CMB) are sensitive to the cumulative distribution of matter along the line of sight from the source, and thus constrain the amplitude of inhomogeneities at intermediate redshifts, providing also some information about the evolution of these fluctuations, if sources at different redshifts are combined.The combined analysis of galaxy clustering and cosmic shear data (the now commonplace "3×2-point" analysis [5]), provides additional sensitivity to the redshift dependence of growth.This is possible by exploiting the fact that the galaxy overdensity largely probes matter structures locally, and thus cross-correlations of cosmic shear with galaxies at different redshifts can be used to reconstruct the growth history with finer granularity.The technique is particularly powerful when employing CMB lensing (a methodology now commonly labelled "CMB lensing tomography" [6,7]).With a lensing source (i.e. the CMB) located behind any other tracer of the large-scale structure (LSS), CMB lensing tomography can be used to reconstruct the growth history up to high redshifts.
Previous attempts at reconstructing the growth history, using redshift-space distortion, 3 × 2 point data, and/or CMB lensing tomography, have revealed hints of a potentially lower growth at low redshifts with respect to the predictions from ΛCDM favoured by primary CMB data from e.g.Planck [8][9][10][11].These outcomes align with the evidence for the so-called S 8 tension, the mild disagreement in the value of S 8 ≡ σ 8 Ω M /0.3 between low-redshift LSS data and CMB experiments.Since this result is not borne out by constraints based on the CMB lensing power spectrum in combination with Baryon Acoustic Oscillations (BAO) data [12,13], it is therefore unclear whether the source of this tension is a potential astrophysical systematic affecting low-redshift tracers on small scales.Alternative measurements of growth at high redshifts are, therefore, highly valuable to pin down the origin of this tension.However, previous growth reconstruction efforts have been somewhat limited by the lack of local LSS tracers at high redshifts (z ≳ 1, [14][15][16]).Although infrared data from, e.g.unWISE are able to reach redshifts z ∼ 1.5 [17,18], the z ≳ 2 regime is largely unexplored.Radio continuum galaxy surveys are sensitive to this regime, but the complexity of obtaining a reliable redshift distribution for high-density continuum samples makes it difficult to exploit them for precision cosmology [19,20].Other probes of the high-redshift LSS will become available with nextgeneration surveys, such as large Lyman-break galaxy samples (LBGs -see [21] for the first cosmological analysis of an LBG sample), or intensity mapping [22].Until these data are available at scale, perhaps our best avenue to explore the high-redshift Universe are quasar catalogues.
Quasars, together with the Lyman-α forest, have provided our highest-redshift measurements of the redshift-distance relation via BAO [23], and they have been exploited to constrain large-scale cosmological observables, such as primordial non-Gaussianity, due to the enormous volumes they are able to probe.Recently, [24] presented Quaia, a quasar catalogue that covers almost the entire celestial sphere, constructed by combining the Gaia catalogue of quasar candidates with infrared photometry from unWISE.The availability of accurate and precise spectrophotometric redshifts allowed [25] to derive relatively tight constraints on ΛCDM from the auto-correlation of Quaia and its cross-correlation with CMB lensing maps from Planck.Our aim in this paper is to extend the analysis of [25], which focused on ΛCDM constraints, considering only two broad redshift bins.Here, we will take advantage of the wide range of redshifts covered by Quaia to reconstruct the growth history in a model-independent way out to z ∼ 3, beyond the range explored by previous analyses (e.g.[17,18,[26][27][28]).In doing so, we will also study the redshift evolution of the quasar bias in the Quaia sample and attempt to isolate potential sources of systematic contamination in the CMB lensing maps, first identified in [25].
The paper is structured as follows.Section 2 presents the Quaia catalogue and the various CMB lensing maps used in our analysis.The theoretical model used in this study, as well as the various data analysis and parameter inference methods employed here, are described in Section 3. The results of our analysis regarding growth and quasar bias evolution, as well as CMB lensing systematics, are presented in Section 4. We then summarise our findings and conclude in Section 5.

Quaia
The Gaia-unWISE Quasar catalogue (Quaia hereafter) is a quasar sample that covers the entire sky and contains almost 1.3 million sources with magnitude G < 20.5.The catalogue, described in detail in [24], was obtained by combining the sample of quasar candidates in the third Gaia data release [29] together with infrared photometry from the unWISE reprocessing of the Wide-Field Infrared Survey Explorer (WISE) [30] to improve the purity of the sample and the quality of the redshift estimates.Quaia sources are assigned a spectrophotometric redshift estimated using a k-nearest neighbours algorithm that combines the Gaia-estimated spectroscopic redshift and photometry from both Gaia and unWISE.The algorithm is trained on spectroscopic data from the Sloan Digital Sky Survey DR16Q [31].Note that we use the latest version of the Quaia catalogue1 , which is slightly different from that used in [25].The latest version uses an updated dust map, which improves the redshift uncertainties slightly; it also includes an updated selection function taking into account unWISE survey properties, as well as improved modelling of high-completeness regions.In our analysis, we divided the total redshift distribution of the sources in three redshift bins, defined by the bin boundaries z Q < 1.0, 1.0 < z Q < 2.3, and z Q > 2.3, where z Q is the Quaia spectro-photometric redshift estimate.The resulting bins are centred at zi = [0.69,1.59, 2.7].This provides us with a better handle on growth evolution than the two bins used in [25], allowing us to also isolate the signal from the structure at z ∼ 3. The angular selection function for each of these bins was generated using the methodology described in [24].We generate a projected ovedensity 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3. .map HEALPix2 [32] for each redshift bin, defined as where n is the unit vector in the pixel direction, N (n) is the number of sources in the pixel, w(n) is the value of the selection function in the pixel, and N the mean number of quasars per pixel, computed as: We mask all areas of the sky in which the selection function takes values w < 0.5.This cut leaves 57%, 59% and 38% available in each redshift bin.The sky mask for the first redshift bin is shown in Figure 1.
The redshift distribution for each bin was calculated by stacking the individual redshift probability density functions (PDFs) of each source (parameterized as normal distributions determined by the estimated redshift and its uncertainty) in the bin.This was found in [25] to be a reasonably accurate estimate of the redshift distribution when compared with direct calibration estimates.The normalised redshift distributions of the three bins are shown in Figure 2, together with the CMB lensing kernel (see Section 3.1).

CMB lensing maps
We use the CMB lensing data from the Planck satellite that measured the angular power spectrum of the CMB lensing convergence on a range of multipoles 8 ≤ ℓ ≤ 2048 [33,34], detecting the signal at 40σ significance from the minimum variance estimator combining temperature and polarization information.Ref. [25] observed hints of potential residual systematics in the Planck CMB lensing maps when cross-correlated with Quaia sources at z ∼ 2. Contamination from the Cosmic Infrared Background (CIB) is an attractive candidate to explain this effect (since it appears in a cross-correlation and is mostly evident at high redshifts), but a more thorough analysis is needed to fully determine its origin and statistical significance.To explore this in detail, we have made use of eight different convergence maps reconstructed from the Planck data, resulting from different data combinations and lens reconstruction techniques.In particular, the maps were obtained from Planck PR43 data generated through the NPIPE pipeline [12].The maps considered are the following.
i.) Generalized Minimum-Variance ("GMV") map.Based on PR4 data, this lensing reconstruction has an improved signal-to-noise ratio of up to ∼ 20% with respect to the previous Planck data release (PR34 ) results due to a joint inverse-variance Wiener filtering of the CMB temperature and polarisation maps accounting for inhomogeneous noise with an optimal weighting.
ii.) Polarisation-only map ("Pol-only") map, constructed through a quadratic estimator employing only polarisation data, and hence largely immune to contamination from unpolarized Galactic or extragalactic foregrounds, such as the CIB or the Sunyaev-Zeldovich effect.
iii.) No-temperature ("No-TT") map.This corresponds to a minimum-variance quadratic estimator that down-weights T T correlations but preserves T E information.Specifically, this map is built in the same way described in the Planck PR3 paper [34] (with the PR4 CMB maps), but leaving out the T T -block in the calculation of the minimum-variance combination.This should still be robust to unpolarized foregrounds while preserving more sensitivity than Pol-only.
iv.) Temperature alone ("TT-only") map.A quadratic estimator that uses only temperature information.
v.) TE-only map: a CMB reconstruction making use only of the correlation between temperature and E-mode polarization, in order to establish whether any differences observed between the GMV/TT maps, and the Pol-only and No-TT maps are driven only by the polarization data.
vi.) SZ-deprojected map ("TT-noSZ").A temperature-only map based on the SMICA foreground-cleaned CMB maps, where the thermal Sunyaev-Zeldovich effect is deprojected [34].The response of this map to CIB contamination should be different from the GMV or TT-only maps (potentially suffering more from it, since the frequency weighting prioritises removing SZ at the expense of increasing contamination from any other component).
vii.) Minimum-variance, source-hardened map ("MVh"): Bias-hardening [35,36] effectively nulls out the contribution to the reconstructed κ map from any contaminant with a known shape of its power spectrum.Source-hardening corresponds to the application of bias-hardening to the case of Poisson-sampled point sources.The resulting κ map is therefore largely immune to any unclustered point-source-like contaminant.Here, the minimum variance consists of an inverse variance weighted combination of all the different reconstruction channels, where the CMB temperature and polarisation maps are filtered separately prior to lensing reconstruction.The bias-hardening operation was performed only for the temperature estimator.
The processing of all CMB lensing maps involves a standardised procedure.First, the spherical harmonic coefficients are rotated to the equatorial coordinates and then transformed into a HEALPix map with a resolution parameter of N side = 512.To avoid aliasing, all multipoles with ℓ > 3N side are filtered out before creating the map.The same angular mask, publicly provided by Planck and shown in Figure 1, is applied to all CMB lensing maps after an analogous rotation and degradation.This mask excludes areas susceptible to contamination by Galactic foregrounds, and removes SZ-clusters detected at high significance as well as a series of bright extragalactic sources, leaving approximately 67% of the total sky.Additionally, in order to test for potential contamination from Galactic foregrounds (either in the κ maps or in Quaia), we will also make use of the 40% Galactic mask made available by Planck.Results using this mask will be labelled Gal.mask.

Theoretical model
In our analysis, we consider two projected probes: the angular overdensity of Quaia quasars δ g (n), and the CMB lensing convergence κ(n).
On the one hand, the quasars are biased tracers of the underlying matter density field, δ m (x).Assuming a linear bias model (valid on the scales used in this analysis), δ g (n) is related to δ m (x) via where N (z) is the normalized redshift distribution of the quasars, and b(z) is the linear bias.On the other hand, CMB lensing is related to the distortion of the CMB photon trajectories due to the presence of the gravitational potential of the LSS of the Universe, and hence it is an unbiased tracer of the matter fluctuations that source this gravitational potential.The relation between κ(n) and δ m (x), assuming ΛCDM, is: where Ω m,0 is the matter density today, H(z) is the Hubble parameter, and χ * is the comoving distance to the last-scattering surface at z * ≃ 1100.Now, let X(n) be a projected field related to the matter fluctuations through the kernel W x (χ) via The angular power spectrum between two such quantities, X and Y , is related to the matter power spectrum P (k, z) via Note that the above expression holds for the Limber approximation [37], valid on the scales used in this analysis.The power spectrum in Equation 3.4 depends upon the overlap between the radial kernels for the quasars overdensity and the lensing convergence, which, based on the above discussion, are: Our analysis will make use of the galaxy auto-correlation, C gg ℓ , and its cross-correlation with CMB lensing, C κg ℓ .To compute these spectra, we need four main ingredients: a cosmological model (e.g. the ΛCDM model), a non-linear model for the matter power spectrum P (k, z) for which we will use the HaloFit model as described in [38], the redshift distribution N (z), and a model for galaxy bias b(z).As discussed in the previous section, we calculate the redshift distributions by stacking the individual redshift PDFs of all sources in each bin, which was found to be sufficiently accurate in [25].
We assume the redshift evolution of b(z) to be the same in all three redshift bins, although with potentially different amplitudes.Specifically, we assume that the redshift dependence of the bias within each redshift bin is given by the fitting function of [39], and thus the bias in bin i is given by: b where b i g is a free parameter of the model that controls the amplitude of b(z).The specific redshift dependence of Equation 3.6 was obtained by [39] by fitting a quadratic polynomial to the auto-correlation of eBOSS quasars at different redshifts.The best-fit of [39] would be recovered for b i g = 1.Note that the specific form of the redshift dependence is only relevant within each of the redshift bins explored, since we allow for a different overall amplitude in each bin.The detailed form of the bias evolution was found by [25] to have little impact on the final ΛCDM constraints.
As described in [25], magnification bias has a negligible effect for Quaia, since the slope of the cumulative flux distribution is close to s = 0.4.We nevertheless account for the impact of this effect in our analysis.

Growth evolution
The main aim of this paper is to use Quaia to reconstruct the growth history, parametrized in terms of σ 8 (z), in a model-independent way.This is particularly interesting at high redshifts (z ≳ 1.5) where the availability of LSS probes is limited, and which Quaia is ideally suited to explore.
Within ΛCDM, and neglecting the presence of massive neutrinos (i.e.ν m ν = 0), the time dependence of the linear matter power spectrum is factorizable as: where D(z) is the linear growth factor.Here, we will parametrize deviations from the ΛCDM prediction for D(z) defining the growth factor as: where D fid (z) is the linear growth factor for a fixed fiducial cosmology, and ∆(z) is a linear interpolation of the three free parameters ∆(z) i i ∈ {1, 2, 3}, corresponding to the value of ∆(z) at the mean redshift of the three redshift bins.Note that, when reconstructing D(z), we do not enforce the normalization condition D(z = 0) = 1, and thus ∆(z) parametrizes both deviations in the redshift dependence of the growth factor and in the overall amplitude of the linear matter power spectrum at z = 0.For this reason, in this case, we fix the value of σ 8 to σ fid 8 = 0.8102 (the best-fit value found by Planck) when calculating P lin (k, z = 0).We will present our results in terms of the redshift-dependent σ 8 (z), which we calculate as σ 8 (z) = σ fid 8 D(z), with D(z) given by Equation 3.8.

Angular power spectra
The auto and cross correlation spectra used in the analysis, are computed through the pseudo-C ℓ formalism [40] as implemented in NaMaster5 [41], which we describe below.Let us assume a generic field X, defined on a pixelized 2D map such as: , where the weight function w X (n) represents the effect of the mask.The pseudo angular power spectrum of two generic masked fields {X, Y } is defined as where ãX ℓm are the harmonic coefficients of map X.This estimator is biased due to the mode coupling induced by the presence of the mask, but it can be related to the true underlying power spectrum via: The term Ñ XY ℓ is the noise pseudo angular power spectrum and it is different from zero for X = Y .In Equation 3.10, M XY ℓℓ ′ is the mode-coupling matrix, which depends exclusively on the pseudo cross-spectrum of the masks W XY ℓ : The unbiased power spectrum estimator is then: To further account for the loss of information induced by the masked sky, the pseudo-C ℓ is divided into band powers.In our analysis, we binned all the cross-and auto-spectra using the scheme proposed by [26]: equally spaced bins with ∆ℓ = 30 in the range ℓ = [2,240] and logarithmic bins for greater ℓs, with ∆ log 10 ℓ = 0.055.In our analysis, we used the NaMaster code for the computation of power spectra and the analytical derivation of their covariance matrices, following the implementation described in [42][43][44].The covariance matrix is computed by considering each field (i.e.κ and g in our analysis) to be a Gaussian random field and correcting for the mode-coupling induced by the presence of the cut sky.This correction is performed using the narrow-kernel approximation (NKA), where the width of the mode-coupling kernel is much smaller than the range of multipoles over which power spectra vary significantly.The mode-coupled noise bias for the quasar auto-correlation Ñ gg ℓ , caused by Poisson shot-noise, can be computed analytically as: In the previous equation, the mask (i.e. the selection function in this analysis) is averaged over the sky, Ω p is the solid angle of each HEALPix pixel and N is the mean number of quasars per pixel.We applied scale cuts to both cross-and auto-spectra.The minimum value of multipole, ℓ min = 30, is set to exclude the range of multipoles in which the galaxy auto-correlation may be dominated by systematics due to dust, stellar contamination and the Gaia scanning pattern, as described in [25].This cut effectively removes the first band power in the power spectrum, and was selected in [25] by observing that it was the only ℓ-bin that changed significantly after linearly deprojecting the contaminants listed above at the pixel level.While this demonstrates that the selection function probably accounts for systematic contamination on smaller scales, it is worth noting that it may be possible to extract information from angular scales larger than that associated with ℓ = 30 by, e.g., estimating the power spectrum on narrower bandpowers below that ℓ.This was not necessary in this analysis, as most of the information is concentrated on smaller angular scales, but other science cases in which large-scale information is critical (e.g.primordial non-Gaussianity) may benefit from doing so.The maximum multipole used in the analysis ℓ max (i.e. the smallest angular scale included) is given, for each redshift bin, by ℓ max = k max χ(z), where k max = 0.15 Mpc −1 , and the comoving distance χ is computed at the mean redshift of each bin.The resulting values of ℓ max for the three redshift bins used in this work are ℓ max = [388, 709, 950].These scale cuts leave a total of 11, 16, and 18 bandpowers in the power spectra involving the first, second, and third redshift bins, respectively.

Likelihood
To constrain the free parameters of our model from the measured cross-and auto-spectra, we use a Gaussian likelihood: where K is an arbitrary constant, d is the data vector storing all the estimated C ℓ s, t(θ θ θ) is the theoretical prediction for d given a set of parameters θ θ θ and Cov −1 is the inverse of the covariance matrix described in Section 3.3.Concretely, in our fiducial analysis, the data vector contains the galaxy autocorrelation and cross-correlation with κ for the three redshift bins described in Section 3.3 (i.e. a total of six power spectra), with the scale cuts described in Section 3. In all the cases studied here, the free parameters of our model will, at least, contain the bias parameters b i g for all the redshift bins that are being analysed.In addition to these, we will also consider variations in different sets of cosmological parameters, depending on the model being studied.The first two models are labelled respectively CosmoFix and CosmoMarg.The CosmoFix model is defined by fixing the six ΛCDM cosmological parameters to the bestfit values found by Planck [33].The CosmoMarg model has three free ΛCDM parameters, {σ 8 , Ω M , h}, fixing the physical baryon density Ω b h 2 , and the scalar spectral index n s to the Planck best-fit values while, as mentioned in Section 3.2, we set the mass of neutrinos to zero.In the third model explored here, labelled CosmoGrowth, we aim to reconstruct deviations from ΛCDM in the growth history.As described in Section 3.2, we quantify this with three parameters, ∆(z i ), which quantify the relative deviation with respect to the growth factor predicted by the Planck best-fit ΛCDM model at the mean redshifts of the three redshift bins considered here.Since we do not enforce any normalisation on the linear growth factor in this case, the value of σ 8 is absorbed by these parameters (in other words, we fix the value of σ 8 used in the template linear power spectrum of Equation 3.7).In addition to b i g , the free parameters of the CosmoGrowth model are therefore {Ω M , h, ∆(z 1 ), ∆(z 2 ), ∆(z 3 )}.The different priors used in our analyses are listed in Table 1 and, in particular, the priors on ∆(z i ) were designed to ensure that D(z) remains non-negative throughout.As in [25], we include a BAO prior, using data from BOSS and eBOSS [2,45].
To explore the multidimensional likelihood, we use a Monte-Carlo Markov Chain algorithm as implemented in Cobaya6 [46], with a convergence criterion of R − 1 < 0.01, where R is the Gelman-Rubin parameter [47].In what follows, when referring to "best-fit" parameters from our data, these will correspond to the MCMC sample with the highest log-posterior.All theoretical calculations were carried out using the Core Cosmology Library (CCL7 , [48]).

Measured power spectra
We measure the quasar auto-correlations in the three redshift bins described in Section 2.1 and their cross-correlation with the different κ maps described in Section 2.2.The resulting auto-spectra are shown in Figure 3 in blue, orange and magenta for the first, second, and third bins, respectively, while their cross-spectra with five of the lensing reconstruction maps used here are shown in Figure 4.In both power spectra, the statistical uncertainties in these measurements are significantly larger for bin 3 than for bin 2. This is due to the contribution from shot noise to the error and is caused by the lower number density of sources in the highest redshift bin.For each redshift bin, we applied the ℓ max cut defined in Section 3.3.Dashed lines represent the theoretical predictions obtained, for each bin, by considering the best-fit bias parameters for the GMV case.The estimated power spectra are plotted with an offset along the x-axis for a better visualization.
Nevertheless, all power spectra are detected with a high signal-to-noise ratio in all bins.
The detection significance is estimated as χ 2 null − χ 2 bf , where χ 2 null is the value of χ 2 (see Equation 3.14) for a null prediction (t = 0), and χ 2 bf is the best-fit value in the CosmoMarg model (see Section 3.4).The detection significances of each individual power spectrum are summarised in Table 2.In most cases, we find that the cross-spectrum has a larger statistical significance than the auto-spectrum, since the latter is more strongly affected by shot noise.The most sensitive cross-correlation measurements correspond, unsurprisingly, to the GMV κ map.Since the sensitivity of the Planck CMB lensing data is dominated by temperature information, the TT-only and TT-noSZ cross-correlations achieve better sensitivity than the Pol-only and No-TT maps.Nevertheless, the No-TT map, which recovers some of  the temperature information, while remaining largely immune to unpolarized foregrounds, is significantly more sensitive than the Pol-only map employed in [25] to quantify a potential foreground contribution to the Quaia cross-correlation at high redshifts.Thus, our fiducial results will largely focus on the GMV and No-TT maps, as they represent the most sensitive CMB lensing reconstructions, achieving the best compromise between sensitivity and robustness to potential extragalactic foregrounds, respectively.However, we will also discuss the results obtained with the other CMB lensing reconstruction methods in Section 4.4 As hinted at in Section 2.2, and in agreement with the findings of [25], we can see that, particularly in the case of bin 2, the amplitude of the cross-correlation with the No-TT and Pol-only κ maps is consistently higher than all other κ maps.We will explore this effect further in Section 4.4.Appendix A presents the constraints on standard ΛCDM parameters found from these measurements, and quantifies the differences with respect to the results presented in [25] due to changes in the Quaia catalogue and selection function, the choice of redshift binning, and the impact of potential Galactic and extragalactic contamination (see Section 4.4).

Growth reconstruction
We use the power spectrum measurements presented in the previous section, to reconstruct the growth history following the procedure described in Section 3.2.As described there, we constrain the deviations with respect to the fiducial ΛCDM growth history ∆(z i ) and, from them, place constraints on σ 8 (z).The results are shown in Figure 5, where the green band shows the 1σ constraints found for the GMV CMB lens reconstruction, while the magenta hexagons and purple circles show results for the GMV-Gal.mask, and No-TT cases, respectively.The dashed line shows the constraints of Planck ( [33]), while the navy band shows the constraints found by [42] by combining a larger suite of galaxy clustering, cosmic shear, and CMB lensing datasets, covering mostly lower redshifts than those explored here.With Quaia, we are able to extend the constraints of [25] to higher redshifts, adding a new data point at z ∼ 2.7.We find that, throughout the entire redshift range covered by Quaia (z ≲ 3), our constraints are in reasonably good agreement with the Planck best-fit model.
The multi-dimensional constraints on the three ∆(z i ) parameters recovered by this analysis are shown in Figure 6.The marginalized constraints on these parameters for the three redshift bins are: ∆(z i ) GMV = {0.15± 0.17, 0.01 ± 0.16, −0.21 ± 0.21}, ∆(z i ) GMV,Gal.mask = {0.01 ± 0.17, 0.11 ± 0.17 Figure 5: Amplitude of the matter fluctuations σ 8 as a function of redshift.The dashed line is obtained with the fiducial Planck cosmological parameters, while the navy shadowed area represents the constraints found by [42].The points in the plot represent the measurements obtained with the Gaussian likelihood for three different CMB lensing reconstructions: GMV (green diamonds), GMV Gal.mask (magenta hexagons) and No-TT (purple dots).
The agreement with the growth history predicted by Planck can be further quantified by fitting the recovered values of σ 8 (z) at the three redshift nodes (i.e. at zi ) with the prediction of Planck with a free amplitude parameter A σ 8 , assuming that the measured values of ∆(z i ) are Gaussian distributed, with a covariance matrix recovered from the MCMC chains.The resulting values of A σ 8 are = 1.12 ± 0.09 compatible with A σ 8 = 1 within ∼ 1σ.More in detail, the upward shift in σ 8 when discarding temperature auto-correlations in the CMB lensing reconstruction, which was described in the previous Section, is only evident here in the intermediate redshift bin (z ∼ 1.6), in which the No-TT analysis recovers a value of ∆(z) that is ∼ 1σ higher than the GMV value 8 .This 0.5 0.0 0.5  weak evidence disappears at both lower and higher redshifts.In fact, the new measurement at z ∼ 2.7 is ∼ 1σ low compared to Planck for the GMV map, and ∼ 1.3σ low for the No-TT map.Using the more restrictive Galactic mask does not lead to significant changes in these results.We explore this in more detail in Section 4.4.

No-TT GMV
Overall, we find that the growth of structure over the redshift range covered by Quaia, which spans the last ∼ 11.8 billion years of cosmic history, is well described by the best-fit ΛCDM model found by Planck.

Bias evolution
Considering the fiducial N (z) for the Quaia bins to be the true description of the quasar redshift distribution, we can constrain the bias evolution of the sources by extracting information from the power spectra.We use the likelihood defined in Section 3 and the bias is evaluated in each redshift bin using Equation 3.6.In particular, we perform two different tests, depending on the data vector and the cosmological model considered.In the first one, the data vector d of Equation 3.14 is given by the combination of the cross-and auto-power spectra of the Quaia bins, and we assumed a CosmoMarg model.The results are shown with filled markers in Figure 7 for the GMV and No-TT cases.The evolution of the bias is generally in qualitative agreement with the fiducial eBOSS-based model (dashed black line), although some differences must be noted.First, the values recovered at low redshifts are consistently  These bias models display a milder evolution at high redshifts (see the text for more details).
higher than those in the eBOSS model.The fact that the Quaia bias is potentially higher than that of the eBOSS model is not entirely surprising, since the Quaia sample is slightly brighter than the eBOSS sample.Perhaps more interestingly, we find that the quasar bias grows more slowly than the eBOSS model at higher redshifts.This, again, is likely due to differences between the eBOSS and Quaia samples.
Given the previous discussion, it is interesting to consider the predicted value and evolution of the quasar bias under the assumption that the best-fit cosmological model found by Planck is the true underlying cosmology.In particular, fixing cosmology allows us to measure the quasar bias using only the quasar-CMB lensing cross-correlation, which is potentially less sensitive to systematic contamination in the quasar overdensity maps, as well as to uncertainties in the redshift distribution of the sample (which may be significant at high redshifts) [49].The outcome of this analysis is shown with hollow markers in Figure 7, squares for the combination of C gg ℓ and C κg ℓ , and circles for C κg ℓ alone.We find that the main effect of fixing the cosmological model is an overall downward shift in the recovered b(z) by about 1 to 1.5σ.This makes sense, since the value of σ 8 recovered by the Quaia data is ∼ 1σ lower than the Planck best fit (see [25] and Appendix A), which must be compensated for by lowering the amplitude of b(z).It is worth noting that this effect is present for both C gg ℓ + C κg ℓ and for C κg ℓ alone, showing that the amplitudes of both power spectra are compatible within the preferred ΛCDM model.In this case, we also recover a slower evolution for b(z) at high redshifts when compared with the eBOSS model.The results found using the No-TT κ map (shown in purple) are generally in good agreement with our fiducial constraints, although with a ∼ 0.5 − 1σ downwards shift, compatible with the higher value of σ 8 recovered by this map.Thus, we conclude that this result is not strongly driven by contamination from extragalactic foregrounds in the CMB lensing map.Repeating our analysis imposing the 40% Galactic mask (not shown in the figure) recovers results that are in good agreement with the fiducial case, reassuring us that our results are not significantly affected by Galactic systematics (neither in Quaia nor in the κ maps).
Finally, we provide a model for the bias of the Quaia sample based on a simpler parametrization of the form b(z) = b 0 /D fid (z), where D fid (z) is the linear growth factor for the best-fit Planck cosmology, and b 0 is a free parameter.This may be convenient for future studies using this sample.For the combination of gg and κg with free cosmological parameters, we find: The resulting best-fit models are shown in Figure 7 with dotted and dashed lines.Using only the galaxy-lensing cross-correlation, and fixing all cosmological parameters to the Planck best fit, we find

Foreground systematics in CMB lensing reconstruction
We now perform an analysis which aims to determine the possible presence of the foreground contamination in the Planck CMB lensing maps, as well as its potential origin.Ref. [25] reported a potentially significant difference in the amplitude of the Quaia-CMB lensing crosscorrelation using different κ maps at high redshifts, which we also saw hints of in the previous sections, and in Appendix A. Isolating the range of redshifts over which this contamination is more relevant would be interesting in order to obtain a more significant detection of it, and to establish its potential extragalactic origin.For this reason, in this case we divide the total redshift distribution of the Quaia sources into six redshift bins, ensuring that each bin contains a suitable number of sources to mitigate excessive shot noise.Specifically, the redshift ranges used to define these bins are: z edges = [0.0,0.5, 1.0, 1.5, 2.0, 2.5, 5].For each redshift bin, we also computed the associated selection function w to take into account the new redshift ranges.
In this case, we did not perform a full likelihood analysis for each map, but instead we simply used the galaxy-κ cross-correlation with fixed cosmological parameters to measure the galaxy bias as a way to quantify the amplitude of this cross-correlation.Since, in this case, the model is linear in b g , it can be computed analytically as: where t is the theoretical cross-power spectrum calculated for b g = 1, d is the measured cross-spectrum, and Cov is the associated covariance matrix.
The results for all eight different CMB κ maps are shown in Figure 8.To aid in visualizing these results, we present our measurements with a normalization factor determined as follows: Colored markers represent constraints from different CMB lensing reconstructions, as indicated in the legend.Each cloud of points is centered on the mean redshift value of the redshift bin, as indicated by the x-axis label.In the fourth redshift bin, centered at z = 1.7, we also display results for the Pol-only map imposing the 40% Galactic mask (dashed error bar), and a more conservative selection function with w > 0.6 (dotted error bar) (see text for more details).
we consider the No-TT measurements (purple points in Figure 8) and fit a second-order polynomial function in redshift to them.Subsequently, we use the best-fit b(z) polynomial obtained as the normalization factor for all the b g measurements.From this figure, it is visually evident that different CMB κ maps are largely consistent with one another throughout the whole redshift range except in the fourth bin, centered around z ≃ 1.7.As in the case of the constraints on σ 8 found by [25] (see also Appendix A), we find that the CMB lensing reconstructions involving T T correlations recover a lower value of b g than Pol-only and No-TT, seemingly at ∼ 2 − 3σ significance.
The actual significance of these differences is hard to assess just from the statistical uncertainties shown in this figure, since the different CMB lensing maps were constructed from the same set of CMB maps, and thus the different estimates of b g are correlated.To quantify the statistical significance of these differences consistently, we made use of correlated simulations as follows.The Planck PR4 lensing analysis delivered 480 independent realizations of simulated lensing reconstructions based on the signal and sky model of the Planck FFP109 simulations and the noise properties of the NPIPE maps.The 8 different types of lensing reconstructions explored here are run on the same simulated data sets.As such, they share the underlying signal and noise in the CMB maps, while the resulting lensing reconstruction noise differs.For each of the 480 available realizations of noisy lensing reconstructions, we retrieved the corresponding input κ signal realization.For each of the κ maps, we generated a correlated Gaussian realization of the quasar overdensity map, δ g .To do so, we assumed the theoretical auto and cross-correlation with κ, the same cosmology as in the FFP10 simulations, as well as the fiducial bias model and redshift distribution for each bin.To simulate the noise contribution for the quasars, we added to δ g a shot noise component modeled using random realizations from the Quaia catalogue matching the number of objects included in  each redshift bin.We then computed the cross-spectrum C κg ℓ for each realization and κ reconstruction method, and estimated the value of b g as in Equation 4.4.Finally, we quantify the level of disagreement between different CMB lensing maps by calculating the fraction of simulated realizations for which a difference in the value of b g with respect to the one found for the corresponding GMV map is larger than the value measured in the data.This corresponds to a simulation-based estimate of the consistency between the cross-correlations with different κ maps in terms of a probability to exceed (PTE).When computing the PTE, we take into account the sign of the difference in the value of b g and, therefore, we considered the 2-sided probability.The probability-to-exceed values are shown in Figure 9.We find that the level of tension in the 4-th bin between the No-TT map and the GMV map, and between the Pol-only map and the GMV map, is relatively high, with probabilities of 0.03 and 0.006 respectively, corresponding to a 2-3σ significance.
Interestingly, this tension is not reproduced by the source-hardened maps (MVh and TTh) and, although the TE map does recover a marginally larger value of b g (by about 0.5σ in the same direction as the Pol-only and No-TT maps), the level of disagreement with the GMV map is much smaller.Although this still does not rule out the CIB or other extragalactic foregrounds a potential cause of this disagreement (since their clustered component might be immune to source hardening), it is worth exploring Galactic foreground contamination in the polarization maps as an alternative possibility to explain this tension.To do so, we repeat our analysis using the 40% Galactic mask from Planck (see Figure 1).The resulting value of b g from the Pol-only map, shown as a blue triangle with dashed error bars, is in good agreement with the result found with our fiducial mask.We performed a similar test using a more restrictive mask, defined by nulling pixels in which the selection function was w ≤ 0.6, and removing a substantially larger sky fraction.Although the value of b g recovered in this case (shown in Figure 8 as a triangle with dotted error bars) is in better agreement with the GMV result, we find that this is mostly driven by statistical fluctuations at ℓ ≳ 300, where the cross-correlation is compatible with zero, whereas at low ℓ, where the cross-correlation is actually detected, the measurement is in better agreement with the No-TT and Pol-only maps on the fiducial mask than with the GMV map.Because of this, together with the significantly larger error bars of this measurement, due to the small sky fraction, we consider this to be compatible with a statistical fluctuation.
We thus find no conclusive evidence that the differences found between different redshift bins are caused by Galactic contamination in the CMB lensing maps.On the other hand, if extragalactic in origin, this potential systematic would evade source-hardened estimators, thus requiring it to have a significant clustered component.Although the CIB is an attractive candidate, it is not clear why its presence is not evident at higher redshifts, and our investigation is ultimately inconclusive.Ultimately, repeating this analysis with CMB lensing maps constructed from other CMB experiments, such as the Atacama Cosmology Telescope (ACT, [50]), or the South Pole Telescope (SPT, [51]), will be vital to address this issue.We leave this analysis for future work.

Conclusions
In this work, we have studied the growth of cosmic structure as a function of time during the last 11.8 billion years, using the auto-correlation of the Quaia quasar sample in combination with its cross-correlation with CMB lensing from Planck.In addition to this, we have studied the redshift dependence of the linear bias of this sample, and quantified the potential presence of Galactic and extragalactic foreground contamination on the Planck CMB lensing maps.
Taking advantage of the wide Quaia redshift coverage, we measured the growth of structure at z ≃ 2.7, achieving one of the highest-redshift constraints on σ 8 in the literature (σ 8 (z = 2.7) = 0.22 ± 0.06 -see Equation 4.1).Our results, together with other relevant growth constraints from [21,26,27,52], are shown in Figure 10.Our measurements of σ 8 (z) are in reasonable agreement (within 1σ) with the evolution predicted by Planck.This result is robust against potential Galactic and extragalactic contamination in the CMB lensing maps, recovering compatible constraints using lensing reconstruction algorithms with varying sensitivity to such contamination.Although our constraints are less precise at z ≲ 1 than others found in the literature, this result is non-trivial: it confirms the validity of the standard cosmological paradigm over a broad range of redshifts extending to a largely unexplored regime (z ≳ 2).As discussed in Appendix A, our constraints on ΛCDM parameters using the newer Quaia catalogue find a value of Ω M that is mildly in tension with the best-fit found by Planck (∼ 2σ higher).
We also measure the bias of the quasars and its redshift evolution.We find that, while the overall amplitude of b(z) is comparable to that of previous quasar samples (e.g.eBOSS), it seems to display a milder evolution at high redshifts, more compatible with a scaling b(z) ∝ 1/D(z), where D(z) is the linear growth factor.This qualitative behavior seems to be dependent of the cosmological model assumed (or simultaneously constrained), and of the subset of the data considered.The amplitude of the b(z) curve does show some sensitivity to the assumptions made about the cosmological model.Since the Quaia data favors a model with a lower value of σ 8 than Planck (although not in significant tension with it), when jointly constraining cosmology and b(z), a larger value of the latter is recovered (by ∼ 1σ) than that 0.0 0.5 found when fixing all cosmological parameters to the best-fit Planck values.This effect is less relevant when using the No-TT κ map, given its better agreement with Planck.Despite these shifts, our results are nevertheless largely compatible between different analysis choices.
Inspired by the slightly different amplitude of the CMB lensing cross-correlation between different reconstruction algorithms (first pointed out in [25]), we investigated the possible presence of foreground contamination on different Planck κ maps.If sourced by extragalactic foregrounds, this contamination could affect other cross-correlation analyses.In order to ascertain the redshift range over which this contamination might be most prominent, we divided the total Quaia redshift distribution into six bins and quantified the amplitude of the quasar-κ cross-correlation in each of them for 8 different κ maps.We find that the effect seems to be localized around redshifts z ∼ 1.7, where results from different κ maps (in particular those using or discarding T T correlations) differ at the 2 − 3σ level.We find no evidence that this tension is the result of Galactic contamination, and source hardening does not ameliorate it either.This leaves extragalactic contamination from a clustered component as perhaps the most likely explanation, although the nature of this contaminant is not entirely clear.Nevertheless, we find that our constraints on the growth history are not significantly affected by this potential systematic.
Our analysis highlights the potential of combining CMB lensing data with high-redshift galaxy samples to constrain and validate the ΛCDM model over the range of cosmic times not directly probed by the CMB or by the low-redshift (z ≲ 1) optical surveys that have so far dominated large-scale structure cosmology.In the future, further progress in this regard may be possible thanks to more extensive optical/IR quasar catalogues [54,55], samples of Lyman-break galaxies and similar dropout populations [21,56], radio continuum surveys, and, potentially, 21cm intensity mapping [57][58][59].In this respect, the cross-correlation with CMB lensing maps will be a vital element to produce reliable cosmological constraints and advances from existing and future experiments, such as the Atacama Cosmology Telescope [13], the Simons Observatory [60], and CMB Stage-4 [61].In particular, the ability to probe significantly smaller angular scales with higher sensitivity, and over a wide range of frequencies, will make it possible to improve the robustness of these constraints to both Galactic and extragalactic systematics.analysis choices, namely the use of 3 redshift bins instead of 2 (to better reconstruct the growth history and isolate the high-redshift contribution), and the study of potential residual Galactic contamination.To quantify the impact of these differences on our analysis, we have derived ΛCDM from our data and compared them with the results of [25]  Figure 11 shows the constraints on Ω M , σ 8 , and S 8 obtained from the old catalogue and the new catalogue under the same analysis choices (2 redshift bins), but using different selection functions (blue and black contours).We find identical constraints on σ 8 , although we observe a shift in 0.7σ upward in Ω M .Note that the changes in the posteriors are primarily determined by differences in the old and updated selection functions, rather than by the catalogue itself.The observed shift leads to a mild tension with the value of this parameter preferred by Planck (shown in orange in the figure) at the level of 2.3σ.To determine if this might be caused by potential Galactic contamination, we repeat the analysis imposing the 40% Galactic mask from Planck (red contours).This results in virtually the same shift in Ω M , albeit with ∼ 15% larger error bars, commensurate with the loss of sky area.We thus conclude that this mild tension is not caused by Galactic contamination.The right panel of Figure 11 then compares the results obtained with the new catalogue using 2 bins (blue) and 3 bins (black), as well as with the 3-bin case using a Galactic mask (red).The use of 3 bins leads to a 0.5σ upward shift in σ 8 , and a 0.3σ 8 downward shift in Ω M , reducing the tension in Ω M to 1.9σ.Unlike in the 2-bin case, applying the Galactic mask shifts both σ 8 and Ω M upward, although the tension with Planck remains at the same (mild) level.The use of the No-TT CMB lensing map (not shown in the figure) does not lead to significant changes in the Ω M tension, but leads to a ∼ 1σ upward shift in σ 8 .This is compatible with the results found by [25] with the Pol-only map and, as hinted at in Section 4.4, may be caused by residual extragalactic contamination in the κ maps using temperature correlations.

Figure 1 :
Figure 1: Mollweide projection in Equatorial coordinates of the masked area of the CMB lensing map (purple), the first Quaia redshift bin (blue), and of the 40% Galactic mask made available by Planck (green).

Figure 2 :
Figure 2: Normalized redshift distribution of the three Quaia redshift bins, defined in Section 2.1.The grey-shaded area represents the CMB lensing kernel W κ (z), with an arbitrary normalisation for visualisation purposes..

Figure 3 :
Figure3: Auto angular power spectrum of the quasars for three Quaia bins, as indicated in the legend, in the range of multipoles from ℓ > 30.For each redshift bin, we applied the ℓ max cut defined in Section 3.3.Dashed lines represent the theoretical predictions obtained, for each bin, by considering the best-fit bias parameters for the GMV case.The estimated power spectra are plotted with an offset along the x-axis for a better visualization.

Figure 4 :
Figure4: Quaia-CMB lensing cross angular power spectrum for different redshift bins (bin1: top panel, bin2: central panel and bin3: bottom panel).In each plot, different colours show the cross-correlation with the various CMB lensing reconstructions considered in the analysis while solid black lines represent the theoretical predictions.The estimated power spectra are plotted with an offset along the x-axis for a better visualisation.

Figure 6 :
Figure 6: Marginalized posterior distribution for ∆ i parameters obtained from the GMV (green) and No-TT (purple) CMB lensing maps.Results are obtained within the CosmoGrowth model.Dashed lines represent the value of ∆(z i ) = 0 predicted by ΛCDM cosmology.

Figure 7 :
Figure 7: Quasar bias evolution as a function of redshift.Different colours show results obtained with different CMB lensing maps (GMV in green, No-TT in purple).Filled markers represent results obtained with the quasars' auto-and cross-correlation with CMB lensing, while empty diamonds and circles are obtained with the CMB lensing cross-correlation alone.Lastly, empty square markers are obtained using the full data vector, namely κg + gg but within the CosmoFix cosmological model.The bias is evaluated at the mean redshift of each bin.The dashed black line represents the fiducial eBOSS bias model of Equation 3.6 for which b g = 1, while the coloured dashed lines represent the bias model proposed here b(z) = b 0 /D fid (z), where b 0 is the best-fit for the CosmoMarg case obtained considering κg+gg.These bias models display a milder evolution at high redshifts (see the text for more details).

6 Figure 8 :
Figure8: Relative amplitude of the CMB lensing cross-correlation as function of mean redshift.Colored markers represent constraints from different CMB lensing reconstructions, as indicated in the legend.Each cloud of points is centered on the mean redshift value of the redshift bin, as indicated by the x-axis label.In the fourth redshift bin, centered at z = 1.7, we also display results for the Pol-only map imposing the 40% Galactic mask (dashed error bar), and a more conservative selection function with w > 0.6 (dotted error bar) (see text for more details).

Figure 9 :
Figure 9: Probability to exceed of the difference between the amplitude of the κg crosscorrelation found with the GMV κ map and all other lensing reconstruction maps considered in Figure 8, as a function of redshift.

Figure 11 :
Figure 11: Constraints on ΛCDM parameters Ω M , σ 8 , and S 8 under different analysis setups, comparing the previous and latest Quaia catalogues, and quantifying the impact of redshift binning and Galactic mask.See main text for details.

Table 1 :
Prior distribution for the parameters used for the different analyses carried out in this work.U (a, b) represents a uniform distribution.thatat times we will study only subsets of this data vector (e.g.only cross-correlations), or different numbers of redshift bins (e.g.see Section 4.4).
3. The length of d is therefore N d = 72 in our fiducial analysis.Note, however,

Table 2 :
Significance of the detection S/N ≡ χ 2 null − χ 2 bf in unit of equivalent σ for a Gaussian distribution for the different bins.Second column: results for the auto-spectrum.Third to seventh column: results for the cross-power spectra between Quaia bins and the lensing reconstructions used in the analysis.Values in bold indicate the highest significance obtained with the different CMB lensing reconstruction.
Figure 10: σ 8 (z) measurements as a function of redshift.We compare our results (magenta stars) with other measurements from previous studies ([21, 42, 52, 53]), as indicated in the legend.The dashed black line shows the ΛCDM prediction when using the Planck best-fit cosmological parameters.
. The left panel of