Cross-correlating radio continuum surveys and CMB lensing: constraining redshift distributions, galaxy bias and cosmology

We measure the harmonic-space auto-power spectrum of the galaxy overdensity in the LOFAR Two-metre Sky Survey (LoTSS) First Data Release and its cross correlation with the map of the lensing convergence of the cosmic microwave background (CMB) from the Planck collaboration. We report a $\sim5\sigma$ detection of the cross-correlation. We show that the combination of the clustering power spectrum and CMB lensing cross-correlation allows us to place constraints on the high-redshift tail of the redshift distribution, one of the largest sources of uncertainty in the use of continuum surveys for cosmology. Our analysis shows a preference for a broader redshift tail than that predicted by the photometric redshifts contained in the LoTSS value added catalog, as expected, and more compatible with predictions from simulations and spectroscopic data. Although the ability of CMB lensing to constrain the width and tail of the redshift distribution could also be valuable for the analysis of current and future photometric weak lensing surveys, we show that its performance relies strongly on the redshift evolution of the galaxy bias. Assuming the redshift distribution predicted by the Square Kilometre Array Design simulations, we use our measurements to place constraints on the linear bias of radio galaxies and the amplitude of matter inhomogeneities $\sigma_8$, finding $\sigma_8=0.69^{+0.14}_{-0.21}$ assuming the galaxy bias scales with the inverse of the linear growth factor, and $\sigma_8=0.79^{+0.17}_{-0.32}$ assuming a constant bias.


INTRODUCTION
Radio continuum surveys have long attracted the attention of the cosmological community as deep tracers of the large-scale structure, able to probe enormous swathes of the Universe. Several forecasting studies have been carried out to explore the cosmological potential of these probes with the Square Kilometre Array (SKA) (Raccanelli et al. 2012;Jarvis et al. 2015). In particular, it has been argued that, using different radio galaxy populations as individual tracers, continuum surveys could be useful in detecting ultra-large scale effects, such as relativistic lightcone effects and the impact of primordial non-Gaussianity (Ferramacho et  Extragalactic radio continuum surveys produce a very different view of the Universe compared to those at shorter wavelengths (e.g. in the optical and IR). At ∼0.1-10 GHz, the dominant mechanism for radio continuum emission is from synchrotron radiation (Condon 1992), which originates from relativistic electrons spiralling in the magnetic fields. Due to this, radio observations are typically dominated by two main galaxy populations: Star Forming Galaxies (SFGs) and Active Galactic Nuclei (AGN). SFGs are often classed into normal starforming galaxies (with star formation rates, SFR 100 M yr −1 ) as well as starburst galaxies, where intense periods of star formation are present (SFR 100 M yr −1 ). For AGN, there is a diverse polychotomy of sources as well as classification schemes. This includes classifying sources based on orientation to the observer (Antonucci 1993;Urry & Padovani 1995), morphol-ogy (e.g. Faranoff-Riley Type I and II sources Fanaroff & Riley 1974), or the accretion mechanism of material onto the central supermassive black hole (High/Low excitation radio galaxies, e.g. Best & Heckman 2012;Heckman & Best 2014). A key advantage in observing both AGN and SFGs at these frequencies is that attenuation of emission by dust, in the inter-galactic medium and along the line of sight, is negligible. This is especially useful for the study of SFGs, where the radio emission can be used as an unbiased measurement of SFR (see e.g. Bell 2003;Davies et al. 2017;Gürkan et al. 2018). Moreover, at these low frequencies, the instantaneous field of view from radio surveys is often large. For example, the LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) has a field of view of ∼30 sq. deg at 150 MHz, the Meer Karoo Array Telescope (MeerKAT; Jonas 2009) has a field of view of ∼1 sq. deg at 1.2 GHz and the Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al. 2007) has an instantaneous field of view of ∼30 sq deg, due to its use of phased array feeds (PAFs). These large fields of view for radio interferometers are crucial for producing large, contiguous observations of the celestial sphere, which is especially advantageous for studying large-scale structure at large angular separations. As such, there exist many large sky surveys at radio frequencies, including: the NRAO VLA Sky Survey (NVSS; Condon et al. 1998, at 1.4 GHz), TIFR GMRT Sky Survey (TGSS-ADR; Intema et al. 2017, at 150 MHz), Sydney University Molongolo Sky Survey (SUMSS; Mauch et al. 2003). Further surveys such as the Evolutionary Map of the Universe (EMU; Norris et al. 2011) and the LO-FAR Two-metre Sky Survey (LoTSS Shimwell et al. 2019), will soon provide a huge leap by combining both depth and sky coverage (e.g. see Fig. 1 of Shimwell et al. 2019).
Already though, the first data release of LOFAR Twometre Sky Survey (Shimwell et al. 2019) has generated a low frequency radio continuum survey, covering 424 sq. deg at 144 MHz to a depth of ∼70 µJy/beam. This is one of the deepest surveys available at this frequency that covers a large area of sky, ideal for studies of large scale structure (see Siewert et al. 2019).
However, the main complication in using continuum data for cosmology is the absence of reliable redshifts for a large fraction of sources. This leads to large uncertainties in the redshift distribution of different samples, and in general on the redshift evolution of any of their properties, such as the galaxy bias or host halo masses. The large-scale structure analysis of continuum samples with optical matches, has been a useful method to address these issues and improve our understanding of the clustering properties of radio galaxies Hale et al. 2018;Siewert et al. 2019).
Further information can be gained through crosscorrelations with other tracers of the large-scale structure (e.g. Lindsay et al. 2014). In particular, since the gravitational lensing of the Cosmic Microwave Background (CMB) is sensitive to the distribution of matter inhomogeneities to very high redshifts, it constitutes an ideal probe to crosscorrelate with deep radio continuum data (e.g. Planck Collaboration et al. 2014; Allison et al. 2015), including in the context of de-lensing (Namikawa et al. 2016). In this paper, we will explore the combination of the auto-correlation of the LoTSS survey and its cross-correlation with CMB lens-ing data from the Planck satellite. In particular, we will show how the different dependence of both correlations on the galaxy bias and redshift distribution of the sample allows us to use this combination to simultaneously constrain the bias and the high-redshift tail of the redshift distribution. This possibility is also relevant in the context of weak lensing analyses with current and future deep imaging surveys (Wright et al. 2020;Hildebrandt et al. 2020), for which the robustness of their cosmological constraints will rely heavily on their ability to calibrate the redshift distributions of their samples.
This paper is structured as follows: in Section 2 we present the theoretical background describing the auto-and cross-correlation signals. Section 3 presents the different datasets used in our analysis. The methods used to extract the auto-and cross-correlations, and to analyse them, are described in Section 4. The main results are then presented and discussed in Section 5, and we conclude in Section 6.

Galaxy overdensity and CMB lensing
The two probes studied here are the projected overdensity of galaxies in the LoTSS sample, δg, and the gravitational lensing convergence of the CMB, κ.
The projected overdensity quantifies the overabundance of galaxies in a given sky positionn, Ng(n) with respect to the sky averageNg and is related to the three-dimensional galaxy overdensity, ∆g(x, z) (defined in a similar way as the fluctuation in the comoving number density of galaxies at redshift z) through δg(n) = dz dp dz ∆g(χ(z)n, z), where dp/dz is the redshift distribution of the galaxy sample normalized to 1 when integrated over z, and χ is the comoving radial distance. The lensing convergence κ(n) quantifies the distortion in the trajectories of the CMB photons caused by the gravitational potential of the intervening matter structures, and is defined to be proportional to the divergence of the deflection in the photon arrival angle α α α: κ ≡ −∇ · α α α/2. As such, κ is an unbiased tracer of the matter density fluctuations ∆m(x, z), and is related to them through: where H0 is the Hubble constant, Ωm is the fractional matter density, a = 1/(1 + z) is the scale factor, and χLSS is the comoving distance to the surface of last scattering. We will use the harmonic-space correlation between κ and δg to study the connection between ∆g and ∆m.

Angular power spectra
Both κ and δg can be generically described as a 2-dimensional field u(n) related to an underlying 3dimensional quantity U (x, z) through a projection onto the two-dimensional sphere with a given radial kernel Wu(χ): u(n) = dχ Wu(χ) U (χn, z(χ)). (4) Any such projected quantity can be decomposed in terms of its spherical harmonic coefficients u m , the covariance of which is the so-called angular power spectrum: The angular power spectrum can be related to the power spectrum of the 3D fields PUV through where PUV (k, z) is the variance of the Fourier coefficients of U and V : We use the following definition for the Fourier transform: For the two fields under consideration, the radial kernels are given by where H(z) is the expansion rate, and Θ(x) is the Heaviside function. Note that we have included a scale-dependent prefactor f in the lensing kernel, given by The presence of this factor is only relevant for 10, and accounts for the fact that κ is related to ∆m through the angular Laplacian of the gravitational potential Φ. We have also neglected the effects of lensing magnification, which are much smaller than our statistical uncertainties (Allison et al. 2015) 1 .
It is worth noting that Eq. 6 is stictly valid only in the Limber approximation, where the spherical Bessel functions are approximated by This is accurate when the radial kernels Wu are broader than the typical correlation length of the density inhomogeneities, as is the case for both δg and κ.
The final ingredient of the model needed to describe the measured angular power spectra is the auto-and crossspectra between the matter and galaxy overdensity fields, Pgg(k, z), Pgm(k, z) and Pmm(k, z). On the scales k 0.2 h Mpc −1 we probe here, we assume a linear, scaleindependent bias model between both fields: ∆g(k, z) = bg(z)∆κ(k, z), such that We explore two different models for the redshift dependence of the bias: • A constant bias model, in which bg does not evolve, and the clustering of galaxies follows the same growth as a function of time as the matter inhomogeneities. • A constant amplitude model, in which bg evolves inversely with the linear growth factor D(z) where the growth factor D(z) ≡ ∆m(k, z)/∆m(k, z = 0) describes the evolution of the matter fluctuations on linear scales (k → 0). In this scenario, the amplitude of ∆g does not vary as a function of time on linear scales. This would correspond to a galaxy distribution that is fixed at some early time and preserves its large-scale properties unchanged.
These two models represent the two extremes of the expected redshift evolution of bg(z). Flux-density limited samples such as the one studied here are expected to have bg grow with z, since more distant sources are intrinsically brighter and are therefore likely to reside in more massive haloes. At the same time, the LoTSS sample contains a wide variety of galaxy types, with different biases and different relative abundances as a function of redshift, and therefore the growth of bg(z) could be slower than the 1/D(z) model. Finally, we model the matter power spectrum Pmm(k, z) using the HALOFIT parameterization of Takahashi et al. (2012) as implemented in the CAMB Boltzmann code (Lewis et al. 2000). We use the Core Cosmology Library (CCL; Chisari et al. 2019) for all theoretical calculations.

DATA
3.1 The LOFAR Two-metre Sky Survey data release 1 The radio continuum sample used here is contained in the first data release (DR1) of the LOFAR Two-metre Sky Survey (LoTSS). Whilst the full LoTSS survey will eventually cover the entire northern sky, this first data release covers a contiguous patch of 424 square degrees in the northern hemisphere over the HETDEX spring field, and contains a total of 325,694 sources. The DR1 is described in detail in Shimwell et al. (2019), and the clustering properties of the sample were studied for the first time in Siewert et al. (2019). DR1 consists of 58 individual pointings, each covering an area approximately the shape of a circle of radius 1.7 • at 6 arcsecond resolution. Each of the 58 pointings was observed for a duration of 8 hours, resulting in an rms noise of ∼70 µJy/beam. Siewert et al. (2019) found five of these pointings to result in a poor completeness (at the level of ∼ 90% at flux S 1 mJy), we discard these pointings in our analysis as well. The average point-source completeness of the remaining area is approximately 99% level for S > 0.8 mJy.
Through a combination of likelihood-ratio matching and visual identification, a value added catalogue (VAC) was provided as part of the DR1, which we use for our analysis (Williams et al. 2019). The VAC contains 318,520 sources, of which 73% have additional optical and or infrared information, primarily from cross-matches with Pan-STARRS and/or WISE. For approximately half of the sample, redshifts were also assigned to sources within the VAC (Duncan The green histogram shows the estimate from the spectroscopic VLA-COSMOS sample at 3 GHz extrapolated to the LOFAR band. The blue histogram shows the distribution obtained from the SKADS simulation, based on existing measurements of the radio luminosity function. The dashed lines are empirical fits to the three distributions using the expression in Eq. 16 with a single free parameter z tail . The black solid line shows the CMB lensing kernel normalized to make it viewable on the same scale. The distributions shown here are plotted as probability distributions dp/dz. et al. 2019). The majority of these are photometric redshifts, whilst a minority of sources have associated spectroscopic redshifts. Of those sources included in the VAC, we use a fiducial flux limited sample defined by a cut in total flux of S > 2 mJy. We additionally discard all sources with a signal-to-noise ratio S/N < 5. The remaining sample contains 57,928 sources.

The Planck lensing maps
We use the publicly available CMB lensing convergence map provided by the Planck Collaboration (Planck Collaboration et al. 2018b). The data are provided in the form of spherical harmonic coefficients κ m , which we transform into a HEALPix map (Górski et al. 2005) with resolution parameter N side = 2048, corresponding to an angular size of ∼ 2 . We verified that this choice of resolution has a negligible impact on our results by re-calculating the δg-κ crosscorrelation presented in Section 5 for N side = 256. When estimating power spectra, we mask this convergence map using the mask provided by Planck, covering around 67% of the sky. This mask was designed to remove regions dominated by galactic foregrounds as well as strong Sunyaev-Zel'dovich sources. This mask removes around ∼ 2% of the LoTSS footprint.

Redshift distributions
The redshift distribution of the LoTSS sample used here is a key ingredient in order to obtain a theoretical prediction for the clustering auto-spectrum and the CMB lensing crossspectrum (dp/dz in Eq 2). The LoTSS value-added catalog provides photometric redshift (photo-z) information for about 51% of the total sample (and ∼ 48% of our flux limited sample), which can be used to estimate this redshift distribution. There are however significant uncertainties associated with these photometric redshifts. Of particular concern is the level to which the subset of the catalog with measured redshifts is representative of the full sample, especially at high redshifts. Although the distribution of radio fluxes in the two samples (the full sample and the subset with measured redshifts), shows no significant differences, this is no guarantee that other selection effects (e.g. cuts in the optical photometry) do not induce significant differences in their redshift distributions. This may be especially important in the tail of high redshift AGN which is typically found in radio surveys. It is likely that the N (z) distribution from Duncan et al. (2019) may underestimate the number of high-redshift sources, especially given radio emission does not suffer dust attenuation whilst the cross-matched catalogues is affected by this, especially for distant objects. Another potential source of error would be mis-identification of high-redshift AGNs as low-redshift SFGs. Also, as noted in Siewert et al. (2019), the low-and high-redshift differential source count distributions show markedly different shapes. Therefore other estimates of dp/dz are needed in order to assess the reliability of our results. The redshift distribution estimated from the LoTSS VAC is shown in red in Fig. 1 as a solid line.
We obtain a second estimate of the redshift distribution from the Very Large Array Cosmic Evolution Survey (VLA-COSMOS) 3 GHz catalog (Smolčić et al. 2017), which comprises ∼ 10, 000 radio sources in a 2 square-degree patch covering the COSMOS field in combination with optical and infrared data. Due to the vast wealth of deep multi-wavelength data in the COSMOS field, 8995 of these sources have measured redshifts. This is likely to therefore provide a much better estimate of the high redshift tail of radio sources. To obtain an estimate of the redshift distribution of the LoTSS sample, we extrapolate the radio fluxes of these objects to the LOFAR 144 MHz band assuming a spectral index α = −0.7 2 and impose the same flux cut of used for our sample (see Section 3.1). Due to this relatively bright flux cut, the resulting sample contains only 378 galaxies. Although the small size of this sample prevents us from obtaining a well-resolved measurement of the redshift distribution, it is enough to quantify the amplitude of the high-redshift tail in comparison with the LoTSS VAC measurement. The resulting redshift distribution, estimated by binning these galaxies into 20 linear redshift bins in the range 0 < z < 5, is shown in green in Fig. 1. Although the VLA-COSMOS catalog potentially underestimates the nearest extended objects due to its baseline sensitivity to extended emission, it displays a significantly larger high-redshift tail than the LoTSS VAC. As we will see, this high-redshift tail plays a significant role in determining the relative amplitudes of C gg and C gκ .
Finally, we use the Square Kilometre Array Design Study Simulated Skies (SKADS), in particular the semi-empirical simulation (Wilman et al. 2008), to obtain a third estimate of the redshift distribution of our sample. This simulation contains a realistic distribution of radio galaxies, including AGN and star-forming galaxies, down to flux densities of 10 nJy based on estimates of the radio luminosity function for different source types, on a square 100 deg 2 patch. Each source has measured flux densities at several frequencies in the range 151 MHz to 18 GHz, which we use to infer their flux in the LOFAR band. The redshift distribution is then computed by binning all sources in the simulation with fluxes above our cut of 2 mJy. The resulting distribution is shown in blue in Fig. 1. Much like the distribution inferred from the VLA-COSMOS catalog, the SKADS estimate displays a long high-redshift tail when compared with the LoTSS dp/dz. There is evidence though, that SKADS may underestimate the number of SFGs below flux densities of S1.4GHz < 0.1 mJy (Smolčić et al. 2017), however as this is below the survey limits considered in this work, the effect of this is likely to be minimal. Estimates of dp/dz will be improved in the future with deeper observations over well-studied multi-wavelength fields, such as through the MIGHTEE (Jarvis et al. 2016) or deep-field LOFAR observations. This large uncertainty in the true redshift distribution for radio continuum surveys motivates the use of the crosscorrelation with CMB lensing as a means to mitigate it. To aid in this study, we fit these three redshift distributions to the same analytic expression dp dz ∝ (z/z0) 2 1 + (z/z0) 2 The normalising prefactor is determined by the requirement that the integral over dp/dz be unity. We find that z0 = 0.1 and γ = 3.5 provide a good visual fit to the three redshift distributions. The free variable z tail parametrizes the extent of the high-redshift tail of the distribution, and therefore its width. We find that this expression is able to provide a rough fit to the three redshift distributions with z tail = 0.8, 1.5 and 2.0 for the LoTSS, VLA-COSMOS and SKADS N (z)s respectively. These analytic fits are shown as dashed curves in Fig. 1. The figure also shows, as a black solid line, the shape of the CMB lensing kernel, which peaks at z ∼ 2 but extends to significantly higher redshifts.

Mean density, sky mask and systematics
A crucial aspect of any analysis involving the fluctuations in galaxy number counts δg, is estimating the expected mean number of objects in each sky pixel. This involves accurately characterizing the geometry of the observed footprint, as well as the expected fluctuations in the number of observed sources due to variations in image noise. We will label these two the mask and mean density map respectively.

Sky mask
We produce masks for each of the 58 pointings from their associated low-resolution images at resolution N side = 2048. Each mask is a binary sky map containing 0s or 1s depending on whether the center of a given pixel lies inside the pointing. The final mask is then the logical addition of the 58 pointing masks.
When computing power spectra with the lower resolution maps, the associated mask is computed by averaging down the N side = 2048 masks. The resulting low-resolution mask contains, in each pixel, the fraction of its area that was observed. This mask is shown in the top left panel of Fig. 2.

Mean density
The expected mean number of observed objects in a given pixel varies across the sky as a result of spatial variations in observing conditions, resulting on fluctuations in flux noise. If not accounted for (i.e. if interpreted as intrinsic fluctuations in the galaxy number density), they will bias the estimated two-point functions and, through them, the final parameter constraints. We used a semi-analytic method to reconstruct the expected fluctuations in the mean density for a given sample that uses the same noise model used by Siewert et al. (2019).
Let S obs be the observed flux of a given object with true flux Strue. The probability to detect objects with observed fluxes above a given threshold S thr is: where p(Strue) is the distribution of true fluxes, and Π(Strue, S thr ) is the probability for a source with true flux Strue to have an observed flux above S thr , given in terms of the conditional distribution of observed fluxes p(S obs |Strue) as Assuming the noise in the measured fluxes has a Gaussian distribution with standard deviation σS (as done in Siewert et al. 2019), this is simply where erfc is the complementary error function. We therefore need two ingredients in order to compute P (> S thr ): • An estimate of the noise rms in each pixel σN (n). As done in Siewert et al. (2019), we make histograms of the rms noise for all sources in a given pixel, and use them to produce maps of the mean and median rms noise per pixel. The rms noise of each source is given as the averaged background rms of the corresponding island. Note that, in order to have enough galaxies in each pixel on average, we produce these maps with resolution N side = 256. • An estimate of the distributon of true fluxes p(Strue). We obtain this from the SKADS semi-empirical simulation (Wilman et al. 2008). We queried the SKADS database to retireve the flux distribution at 151 MHz, which we scaled to 144 MHz assuming a spectral index α = −0.7.
We use this method to produce a mean density map for our sample. Since we only consider sources detected above 5σ, the threshold flux is S thr = max(5σN , Scut), where Scut is the flux cut defining the sample. Figure 2 shows, in the top right panel, the resulting map for our fiducial Scut = 2 mJy. The sample is highly homogeneous, with a mean densitȳ n = 0.12 arcmin −2 that varies by ∼ 1.4% across the footprint (1 standard deviation).

Noise variations
The number of detected sources is likely to depend strongly on the depth of a given sky region. As a tracer for those depth variations, we map the coadded inverse background noise variance. For each pointing, we estimate its noise rms σN as the standard deviation of the low-resolution residual map made available with the DR1. The final inverse noise variance map is generated by adding the values of σ −2 N for all pointings overlapping on each pixel. It is worth noting that the procedure used to generate this map does not reflect exactly the mosaicing process used in LoTSS. Additionally, there are probably fluctuations of the flux density scale between pointings that have not been quantified. It is expected that these shortcomings will improve with future data releases. This map is shown in middle left panel of Figure 2, and will be used in the next section to study the impact of this systematic on the power spectrum.

Maps and power spectra
Before computing the power spectra, we first generate a map of the overdensity field δg(n) as: where N (n) is the number of galaxies observed in the pixel with positionn, and wg(n) is a weight map given by the product of the sky mask described in Section 4.1.1, and the mean density map described in Section 4.1.2. wg(n) therefore quantifies the fraction of galaxies passing our cuts at a given pixel that should have been observed. The quantitȳ N is the mean number of objects per pixel in the sample, estimated asN = N (n) θ / wg(n) θ , where · θ denotes an average over all pixels in the map. Finally, to avoid using heavily masked pixels, we set to zero all pixels in δg(n) and wg(n) where wg(n) < 0.5. Once we have the two maps δg(n) and κ(n) (both shown in Fig. 2, and their associated masks wg(n) and wκ(n) (the latter described in Section 3.2), we compute the autocorrelation of δg (C gg ), and its cross-correlation with κ (C gκ ) using the so-called pseudo-C estimator (Peebles 1973;Hivon et al. 2002) using the NaMaster code 3 (Alonso et al. 2019).
The pseudo-C algorithm uses analytical methods to calculate the coupling of different multipoles due to the incomplete sky coverage. Although the method is in principle less optimal than minimum variance quadratic estimators, it is able to achieve almost equivalent uncertainties for sufficiently flat power spectra (as is the case here) with a much higher computational speed. We bin both power spectra in bands of ∆ = 50 starting at = 2.
The estimated overdensity map contains shot noise due to the discrete nature of galaxy number counts. The associated contribution to the angular power spectrum (commonly called the "noise bias"), can be estimated analytically as described in Alonso et al. (2019). We will assess the validity of this calculation in Section 5.1.
We also use NaMaster to estimate the covariance matrix of the power spectra following García-García et al. (2019). This is calculated as the so-called Gaussian covariance, which approximates both κ and δg as Gaussian random fields. The covariance matrix estimator implemented in NaMaster accurately accounts for the correlation between different bins induced by the incomplete sky coverage. Since the estimator requires a best-guess estimate of the underlying power spectra (C gg , C gκ and C κκ ), the covariance is estimated in two steps: first, we compute theoretical power spectra for cosmological parameters fixed to the best-fit values found by Planck (Planck Collaboration et al. 2018a) and a constant galaxy bias bg = 1.3 assuming the LoTSS redshift distribution, which provides a good visual fit to the data. The resulting covariance is then used in the likelihood described in Section 4.3 to find the best-fit parameters. These are used to estimate new theory power spectra that are then used by NaMaster to estimate our final covariance matrix. Note that the auto-spectra C gg and C κκ should contain both the signal and noise contributions. For C gg we use the shot noise component described above, while for C κκ we use the noise curves provided in the Planck data release.
It is worth noting that, although the Gaussian covariance described above assumes Gaussian statistics for both δg and κ, which is known to be inaccurate due to the nonlinear growth of structure, the covariance estimated using this method accounts for the largest fraction of the statistical uncertainty (Barreira et al. 2018;Nicola et al. 2020), and is an excellent approximation in the range of scales studied here.

Likelihood
In order to extract information from the measured C gg and C gκ we use a Gaussian likelihood of the form: where the data vector d denotes all measured power spectra and t(q) is the theoretical prediction for d given a set of parameters q.
We will present results for different choices of data vector and parameters. Specifically, we will present constraints based on C gg and C gκ alone, as well as from the combination of both. We will also explore different combinations of three free parameters: • The galaxy bias bg (within the two redshift evolution models described in Section 2.2). • z tail , which quantifies the extent of the redshift distribution tail when parametrized according to Eq. 16.  Figure 3. Galaxy auto-correlation C gg (top panel) and its cross-correlation with the CMB lensing convergence C gκ (bottom panel). The measured power spectra are shown as black points with error bars. The solid red curve shows the best-fit theory prediction assuming the SKADS redshift distribution, the Planck cosmological parameters, and a galaxy bias that grows inversely with the linear growth factor (Eq. 15), with an amplitude bg = 2.1. The noise bias due to shot noise in the auto-correlation is shown as a gray line in the top panel. The gray points and error bars in the bottom panel show a null cross-correlation calculated by correlating the galaxy overdensity map with a randomly rotated CMB convergence map.
• The amplitude of matter fluctuations as parametrized by σ8.
We fix all cosmological parameters (including σ8 when not used as a free parameter) to the best-fit values found by Planck (Planck Collaboration et al. 2018a). We will only use the multipoles smaller than max = 500, corresponding to a wavenumber kmax 0.15 Mpc −1 at z 1. Thus we make sure that we only use modes where a linear, scale-dependent bias relation is a good approximation, and where the non-Gaussian contributions to the covariance matrix can be neglected.

Power spectra and covariances
The power spectra C gg and C gκ , measured using the methods described in Section 4.2, are shown as black points with error bars in the upper and lower panels of Figure 3. The correlation matrix rij = Covij/ CoviiCovjj associated with the covariance matrix of the full data vector is shown in Fig. 4. The lowest multipoles of C gg and C gκ are ∼ 40% correlated.
Before using these measurements to extract parameter constraints, we perform a number of sanity checks and null tests. The main concern regarding C gg is the presence of residual systematics in the galaxy overdensity producing extra power on large scales (Ross et al. 2012;Leistedt et al. 2013Leistedt et al. , 2016. In our case, the most likely source of such systematics is fluctuations in survey completeness caused by variations in survey properties. As a proxy for these residual systematics, we use the maps of pointing noise variation and the fluctuations in the mean density map itself. We then compare the power spectrumĈ gg estimated from the original overdensity map and from a "systematics-deprojected" map, in which these systematic templates are projected out at the map level (see Alonso et al. 2019, for details). We also verify that the power spectra are stable against the choice of pixelization by recomputing them with a resolution parameter N side = 256. The result of these two tests is displayed in Fig. 5. The figure shows the relative difference between the fiducial and alternative power spectra as a fraction of the 1σ uncertainties. The power spectra are robust to the choice of pixelization and to the potential contamination from variations in survey completeness within the range of scales explored here. In order to validate the cross-correlationĈ gκ , and to test for a potential mis-estimate of the statistical uncertainties, we perform one further null test by cross correlating the overdensity map δg with a random convergence map generated by arbitrarily rotating the original Planck map. The measured cross correlation is shown by the grey points in the lower panel of Fig. 3. Fitting a constant amplitude to these measurements we find the best-fit value A null = (−2.9 ± 4.1) × 10 −8 , consistent with zero.
The galaxy auto-correlationĈ gg could still be affected by unknown systematics that are not well described by the two template maps used above. These systematics would typically affect the largest scales, comparable with the size of a pointing ( ∼ π/δθ ∼ 50), and therefore we can test for the relevance of unknown contaminants by removing those scales from the analysis. Another source of systematic uncertainty that could bias our results is a mis-estimate of the noise bias inĈ gg . This noise bias is shown as a gray solid line in the top panel of Fig. 3, and was calculated analytically. An error in this estimate would mostly affect the smallest scales, which are dominated by shot noise, potentially biasing the inferred parameters. To study the impact of these two systematics, we constrain the galaxy bias parameter bg from the auto-correlationĈ gg alone in three cases: (i) A fiducial case, using all the data and the fiducial Gaussian covariance matrix. We find a bias value bg = 2.10 ± 0.10. (ii) We remove the first bandpower, covering scales < 50. The inferred bias is bg = 1.98 ± 0.12. (iii) Instead of analytically computing the noise bias and subtracting it from the data, we marginalize over a constant offset toĈ gg with a free amplitude parameter. This is done analytically by simply modifying the inverse covariance as (Rybicki & Press 1992): where n is our fiducial estimate of the noise bias. The corresponding constraint is bg = 2.23 ± 0.16.
In all cases, we parametrize the redshift dependence of bg(z) using the "constant amplitude" model, and we assume the SKADS redshift distribution. Given the restricted range of scales used in this analysis, marginalizing over the noise bias leads to a 60% increase in the final uncertainties. Nevertheless, the inferred values of bg are compatible with each other in all cases, and therefore we conclude that the measured spectra are robust to unknown sky contaminants and uncertainties in the shot noise level. . Measured auto-and cross-correlation (black dots with error bars int the top and bottom panels respectively), together with the theory prediction for different values of the galaxy bias bg and the high redshift tail z tail (left and right panels respectively), both in the range (0.5, 2.0). A "constant amplitude" model is assumed for the redshift evolution of the galaxy bias. bg is fixed to 1.3 in the right panel, while z tail = 1.1 in the left one. While both bg and z tail affect the amplitude of the auto-correlation, the cross-correlation depends only mildly on the high-redshift tail, making it possible to break the degeneracy between both parameters by combining C gg and C gκ .

Cross-correlation significance
The measured cross-correlationĈ gκ is shown in the lower panel of Figure 3. We can quantify the significance of this detection in two different ways. First, the probability to exceed associated with the χ 2 of the measurement ofĈ gκ with respect to the null hypothesis (t = 0) is p(> χ 2 ) = 8 × 10 −5 , corresponding to a ∼ 4σ detection. Secondly, we fit the measured cross-correlation to a single parameter model of the form where C gκ | b=1 is the predicted signal assuming the Planck best-fit cosmological parameters, the SKADS redshift distribution and a unit galaxy bias parameter bg = 1 in the constant-amplitude redshift dependent model. The constraint on the free amplitude bg is bg = 1.59 ± 0.31, corresponding to a 5.2σ measurement of the cross-correlation. Note that this is mathematically equivalent to quantifying the significance as the square root of the difference in χ 2 between the null hypothesis and this best-fit model ∆χ 2 .

Galaxy bias and high-redshift tails
As described in Section 3.3, we have at our disposal at least three different estimates of the redshift distribution for our sample, which differ mostly in terms of the size of their high-redshift tails. It is therefore important to understand the level to which the galaxy bias constraints depend on the redshift distribution uncertainties. Assuming the Planck cosmological parameters and a constant-amplitude redshift dependent model for the galaxy bias, the constraints on bg  . Constraints on galaxy bias bg and the high-redshift tail parameter z tail from the combination of C gg and C gκ . Results are shown for a constant galaxy bias as a function of redshift, the "constant bias" model (black), and for a bias that scales with the inverse of the linear growth factor, the "constant amplitude" model (orange).
assuming these three different redshift distributions, and using bothĈ gg andĈ gκ are: bg = 2.10 ± 0.10, (SKADS).  , and for a bias that grows with the inverse of the linear growth factor (orange region). The redshift distributions inferred from the photometric redshifts in the LoTSS value-added catalog, from the 3 GHz VLA-COSMOS catalog, and from the SKADS simulations are shown in red, green and blue respecitvely for comparison. All redshift distribution are normalized to have the same maximum amplitude. The constraints show a preference for longer redshift tails than the one predicted by the LoTSS photo-zs.
While the VLA-COSMOS and SKADS estimates are roughly compatible with each other, they differ from the LoTSS-VAC result by more than 7σ. The differences between the different redshift distribution estimates are therefore highly significant in this context. In order to investigate this further, we repeat the exercise using only the CMB lensing cross-correlation. In this case the constraints are fully compatible with each other (albeit with larger error bars): bg = 1.59 ± 0.31, (SKADS).
The reason for the disagreement between the different inferred bias values with the full data vector is the wellknown fact that both the galaxy bias and the width of the redshift distribution have a degenerate effect on the amplitude of the galaxy auto-correlation. Structure is washed out in samples with broader redshift distributions, which can be compensated with a larger bg.
To understand the role played by both parameters in the auto-and cross-correlation, Figure 6 shows the measured power spectra together with the theoretical predictions for varying values of the galaxy bias bg (left panel) and the high-redshift tail z tail (right panel). While the galaxy bias affects the amplitude of the auto-and cross-correlations (quadratically and linearly respectively), an increasing highredshift tail lowers the auto-correlation while leaving the cross-correlation almost constant. We can understand the latter result as follows: since the CMB lensing kernel Wκ extends to very high redshifts, a variation in the width of the galaxy redshift distribution leaves the overlap between Wg H18, All H18, SFG H18, AGN Figure 9. 1σ constraints on the galaxy bias bg(z) from the combination of C gg and C gκ . Results are shown for a constant bias model (gray region), and for a bias that grows with the inverse of the linear growth factor (orange region). The figure also shows measurements of the bias for different radio galaxy samples in Hale et al. (2018) (black symbols with error bars), and for the NVSS sample of Nusser & Tiwari (2015).
and Wκ almost unchanged. On the one hand, this makes the cross-correlation between CMB lensing and galaxy clustering robust to uncertainties in the redshift distribution width.
On the other hand, the different response of the auto-and cross-correlation to z tail should allow us to break its degeneracy with bg, allowing us to simultaneously constrain both parameters by combiningĈ gg andĈ gκ .
To explore this, we parametrize the redshift distribution according to Eq. 16 and derive constraints on two free parameters, bg and z tail , with flat priors bg ∈ (0.6, 6) and z tail ∈ (0.1, 5). The resulting constraints are shown in Fig. 7 for the constant-amplitude (orange contours) and constant-bias (black-gray contours) redshift dependent models of bg(z). Looking at the constant-amplitude constraints, the uncertainty on bg (σ(bg) = 0.28) degrades significantly with respect to the previous results where the redshift distribution was fixed to the SKADS estimate. The high-redshift tail, on the other hand, is constrained to be z tail = 1.30 +0.27 −0.40 . In the constant-bias case, a larger value of bg is preferred, to compensate for the lack of growth as a function of redshift, which is accompanied by a preference for larger values of z tail and overall larger uncertainties on both parameters. In both cases, however, the data show a preference for larger redshift tails than that implied by the photometric redshifts included in the LoTSS value-added catalog. This can be seen explicitly in Figure 8, which shows the 1σ bounds on the recovered redshift distribution in comparison with the estimates from the LoTSS VAC, VLA-COSMOS and SKADS.
The corresponding constraints on bg(z) are shown in Fig. 9 in both cases. The true redshift evolution of the effective galaxy bias for the LoTSS sample analysed here is most likely not constant but potentially less steep than the constant-amplitude model, therefore lying somewhere between these two extremes. To illustrate this, the figure also shows the bias values measured for different radio popula-tions by Hale et al. (2018), and the NVSS sample of Nusser & Tiwari (2015). The results above therefore imply that the true underlying distribution is probably more compatible with the SKADS or VLA-COSMOS estimates than the distribution of photometric redshifts in the LoTSS VAC.

Constraints on bias and σ8
Under the assumption that the true underlying redshift distribution of the LoTSS sample is well described by the SKADS estimate, we can use the combination ofĈ gg and C gκ to break the degeneracy between bg and the amplitude of matter fluctuations parametrized by σ8. Given the existing redshift distribution uncertainties, and the degeneracy with other cosmological parameters, the resulting constraints would be neither robust nor competitive. The aim of this exercise is therefore twofold: a sanity check to verify that our data is broadly compatible with the standard cosmological model, and a demonstration of the use of continuum surveys for cosmology.
Both values are in agreement with each other as well as with those those found by Planck (σ8 = 0.811 ± 0.006, shown in blue in the figure), with a preference for lower σ8 values along the bg-σ8 degeneracy direction.

CONCLUSIONS
Radio continuum surveys are a valuable probe to study the properties and evolution of star-forming galaxies and AGNs. Unimpeded by dust attenuation, radio surveys cover a much broader range of redshifts, and correspondingly larger volumes, than their optical counterparts. Because of this, continuum surveys also constitute a potential avenue to reconstruct the growth of structure over the largest scales, and therefore have drawn the interest of the cosmology community. However, the lack of redshift information from radio continuum data makes these surveys reliant on matching to deep optical catalogs to constrain the radial distribution of the sources, and to reconstruct the redshift evolution of their properties. The potential incompleteness of the optical crossmatches makes redshift calibration one of the largest sources of systematic uncertainty in the potential use of radio continuum surveys for cosmology. This is akin to the challenges faced by ongoing and future photometric weak lensing surveys when characterising and propagating the uncertainties in the redshift distribution of galaxies (Sánchez et al. 2020;Hildebrandt et al. 2020;Schaan et al. 2020).
In this paper we have studied the clustering of galaxies on the flux-density limited LoTSS first data release, both through the harmonic-space auto-correlation, and through their cross-correlation with the lensing convergence of the CMB, as measured by Planck. The cross-correlation is detected at the 5σ level.  To illustrate the challenge posed by the uncertainties in the redshift distribution dp/dz, we have considered three estimates of this quantity: from the LoTSS value-added catalog, from the VLA-COSMOS cross-matched sample (which includes both spectroscopic and photometric redshifts), and from the SKADS simulations. We have shown that the truncated high-redshift tail from the LoTSS-VAC photometric redshifts leads to radically different interpretations of the clustering amplitude when compared with the results from the VLA-COSMOS or SKADS redshift distributions. On the other hand, the CMB lensing cross-correlation alone is fairly insensitive to variations in dp/dz, and leads to consistent measurements of the galaxy bias.
The robustness of the CMB-lensing cross-correlation to redshift distribution uncertainties, allows us to break the degeneracy between the galaxy bias and the width of the dp/dz by combining it with the clustering auto-correlation. Through this joint analysis, we are able place constraints on the high-redshift tail of the distribution, showing that it is underestimated by the LoTSS-VAC photo-zs, and better represented by the VLA-COSMOS and SKADS estimates. To our knowledge, this is the first attempt in the literature at calibrating redshift distributions through cross-correlations with CMB lensing. This is akin to the "cross-correlation redshifts" approach, based on using cross-correlations against samples with a known redshift distribution (Newman 2008;Alonso et al. 2017;Gatti et al. 2018) (in this case, the CMB lensing kernel). Given the broad range of redshifts covered by the CMB lensing kernel, it is unlikely that this approach will be able to constrain a shift in the mean redshift of the target distribution (e.g. see Alonso et al. 2017), as is usually done in weak lensing analyses ), but it should be useful to calibrate the width of the distribution or the presence of photo-z outliers. This could provide a robust way to extract cosmological information from samples with poor spectroscopic coverage in future surveys, such as the Legacy Survey of Space and Time (LSST; The LSST Dark Energy Science Collaboration et al. 2018).
Finally, assuming that the redshift distribution is well described by the SKADS estimate (which is broadly consistent with the VLA-COSMOS redshift distribution), we use the combination of clustering and CMB lensing to measure the amplitude of matter inhomogeneities σ8. Although the resulting constraints are neither competitive nor robust, given the existing uncertainties on dp/dz and on the redshift dependence of the galaxy bias, this exercise allows us to demonstrate the use of continuum surveys for cosmology. The measured value of σ8 = 0.79 +0.17 −0.32 , assuming a constant bias, agrees well with the constraints found by CMB and large-scale structure experiments (Planck Collaboration et al. 2018a;Abbott et al. 2018;Hikage et al. 2019;Asgari et al. 2020).
This paper highlights the novel science that can be carried out by combining CMB and radio-continuum survey using the first data release from LoTSS, which covers 424 square degrees. LoTSS itself should cover the vast majority of the northern sky when completed (∼ 20, 000 square degrees, allowing much stronger constraints on the various results presented here. Furthermore, the EMU Survey (Norris et al. 2011) will provide a complementary southern hemisphere survey, down to a similar equivalent flux-density limit at GHz frequencies. However, as we have shown in this work, one of the key uncertainties in using radio continuum surveys for cosmology arises from the lack of a well constrained redshift distribution. Fortunately, there are a wealth of surveys that will aid in pinning this down on the timescale of these all-sky surveys. Future LoTSS data releases will be able to make use of deep fields with ∼ 90% photometric redshift coverage. Furthermore, the WEAVE-LOFAR Survey (Smith et al. 2016) will provide a huge number of spectroscopic redshifts for the LoTSS sources, combining deep spectroscopy in the LOFAR deep survey fields, with spectroscopy of the rarer brighter sources, through a tiered radio-continuum selected survey. However, spectroscopic surveys are rarely complete, and tend to be limited by the optical magnitude depth or the ability to detect emission lines (an advantage for radio continuum selected surveys). As such, further important information will be gleaned from measuring the redshift distribution of radio sources in the deep extragalactic survey fields, similar to how we have used the VLA-COSMOS survey here. However, a single deep field is inherently limited by sample variance, as such the MeerKAT International Giga-Hertz Tiered Extragalactic Exploration (MIGHTEE; Jarvis et al. 2016) Survey, which will cover ∼ 20 square degrees over four extragalactic deep fields, accessible from the southern hemisphere, will help pin down the redshift distribution of radio continuum sources to very deep flux-density limits using both photometric (e.g. Jarvis  Taken together these advances will allow us to begin to divide the radio continuum source population into distinct sub-samples, e.g. SFGs and AGN. This would then allow us to introduce more realistic bias evolution models for these different populations, with different redshift distributions, which is a limitation on our current work due to the relatively small survey area of the LoTSS DR1. Looking further in the future, the sky surveys planned for the Square Kilometre Array will lead to higher number density of objects across huge swathes of sky. Although this may not substantially increase the number of AGN at high redshift (z > 2), these surveys will increase the number of SFGs across all redshifts, and allow the possibility of a more robust decoupling of the SFG and AGN populations (e.g. Makhathini et al. 2015;Muxlow et al. 2020), given the much higher angular resolution (< 1 arcsec). This could also be achieved by LOFAR using its international baselines.

DATA AVAILABILITY STATEMENT
The data underlying this article are available in the Planck legacy archive https://www.cosmos.esa.int/web/planck/ pla, the LOFAR surveys website https://lofar-surveys. org/releases.html, and the VLA-COSMOS 3 GHz repository at the Strasbourg astronomical Data Center http: //cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/602/A2. data based on observations obtained with Planck (http: //www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. Contour plots were generated using the GetDist package (Lewis 2019). Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. We made extensive use of the scipy (Virtanen et al. 2020), healpy (Zonca et al. 2019), and matplotlib (Hunter 2007) python packages.