11institutetext:  Division of Physics, Mathematics and Astronomy, California
Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125
22institutetext: Department of Astrophysics, National Centre for Nuclear
Research, Pasteura 7, 02-093 Warsaw, Poland 33institutetext: Sub-department of
Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road,
Oxford, OX1 2DL, UK 44institutetext: Center for Theoretical Physics, Polish
Academy of Sciences, al. Lotników 32/46, 02-668 Warsaw, Poland
55institutetext: Fakultät für Physik, Universität Bielefeld, Postfach 100131,
33501 Bielefeld, Germany 66institutetext: Institute for Astronomy, School of
Physics and Astronomy, University of Edinburgh, Royal Observatory Edinburgh,
Blackford Hill, Edinburgh, EH9 3HJ, UK 77institutetext: Astronomical
Observatory of the Jagiellonian University, 31-007 Cracow, Poland
88institutetext: Institute for Theoretical Physics, Heidelberg University,
Philosophenweg 16, 69120 Heidelberg, Germany 99institutetext: National
Astronomical Observatories, Chinese Academy of Science, Beijing, 100101,
P.R.China 1010institutetext: Hamburg Observatory, University of Hamburg,
Gojenbersgweg 112, 21029 Hamburg, Germany 1111institutetext: Department of
Physics, University of the Western Cape, Bellville 7535, South Africa
1212institutetext: ASTRON, Netherlands Institute for Radio Astronomy, Oude
Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands 1313institutetext:
Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, The
Netherlands

# Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-correlation
with Cosmic Microwave Background

S.J. Nakoneczny Cosmology from LOFAR Two-metre Sky Survey Data Release 2:
Cross-correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-
metre Sky Survey Data Release 2: Cross-correlation with Cosmic Microwave
BackgroundCosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background D.
Alonso Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background M.
Bilicki Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background D.J.
Schwarz Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background C.L.
Hale Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background A.
Pollo Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave
BackgroundCosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background C.
Heneka Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background P.
Tiwari Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background J.
Zheng Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background M.
Brüggen Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background M.J.
Jarvis Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave
BackgroundCosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background T.W.
Shimwell Cosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave
BackgroundCosmology from LOFAR Two-metre Sky Survey Data Release 2: Cross-
correlation with Cosmic Microwave BackgroundCosmology from LOFAR Two-metre Sky
Survey Data Release 2: Cross-correlation with Cosmic Microwave Background
<EMAIL_ADDRESS>
###### Abstract

Aims. We combine the LOw-Frequency ARray (LOFAR) Two-metre Sky Survey (LoTSS)
second data release (DR2) catalogue with gravitational lensing maps from the
Cosmic Microwave Background (CMB) to place constraints on the bias evolution
of LoTSS detected radio galaxies, and on the amplitude of matter
perturbations.

Methods. We construct a flux-limited catalogue from LoTSS DR2, and analyse its
harmonic-space cross-correlation with CMB lensing maps from Planck,
$C_{\ell}^{g\kappa}$, as well as its auto-correlation, $C_{\ell}^{gg}$. We
explore the models describing the redshift evolution of the large-scale radio
galaxy bias, discriminating between them through the combination of both
$C_{\ell}^{g\kappa}$ and $C_{\ell}^{gg}$. Fixing the bias evolution, we then
use these data to place constraints on the amplitude of large scale density
fluctuations, parametrised by $\sigma_{8}$.

Results. We report the significance of the $C_{\ell}^{g\kappa}$ signal at a
level of 26.6$\sigma$. We determine that a linear bias evolution of the form
$b_{g}(z)=b_{g,D}/D(z)$, where $D(z)$ is the growth rate, is able to provide a
good description of the data, and measure ${b_{g,D}=1.41\pm{0.06}}$ for a
sample flux-limited at $1.5\,{\rm mJy}$, for scales $\ell<250$ for
$C_{\ell}^{gg}$, and $\ell<500$ for $C_{\ell}^{g\kappa}$. At the sample’s
median redshift, we obtain ${b(z=0.82)=2.34\pm{0.10}}$. Using $\sigma_{8}$ as
a free parameter, while keeping other cosmological parameters fixed to the
Planck values, we find fluctuations of ${\sigma_{8}=0.75^{+0.05}_{-0.04}}$.
The result is in agreement with weak lensing surveys, and at $1\sigma$
difference with Planck CMB constraints. We also attempt to detect the late-
time integrated Sachs-Wolfe effect with LOFAR data, but with the current sky
coverage, the cross-correlation with CMB temperature maps is consistent with
zero. Our results are an important step towards constraining cosmology with
radio continuum surveys from LOFAR and other future large radio surveys.

###### Key Words.:

cosmological parameters – large-scale structure of Universe – cosmic
background radiation – dark matter – Cosmology: observations – Radio
continuum: galaxies – Methods: statistical – Galaxy: evolution – Galaxies:
active – Galaxies: distances and redshifts – Galaxies: photometry –
Techniques: photometric

††offprints: S.J. Nakoneczny, .

##  1 Introduction

One of the main current goals of observational cosmology is constraining the
history of structure growth as a way to pin down the different components that
dominate the background expansion of the Universe at late times (Huterer,
2022). To do so, one must investigate probes of structure that, ideally,
satisfy a number of criteria. Among these:

  1. 1.

They should cover a large enough patch of the Universe to reconstruct the
density field on cosmological scales.

  2. 2.

They should also be easy to connect with fundamental quantities, such as the
matter overdensity.

  3. 3.

They should contain redshift information, allowing for an accurate
reconstruction of the structure growth history.

Unfortunately, virtually no cosmological probe, taken alone, is able to
fulfill these requirements. Although weak gravitational lensing, measured from
its effect on the shapes of background galaxies, or on the cosmic microwave
background (CMB) fluctuations (Bartelmann & Schneider, 2001) is an unbiased
tracer of the total matter fluctuations, it has significantly lower raw
statistical power and poorer ability to trace redshift evolution than
measurements of galaxy clustering. The latter, in turn, is cursed by the
problem of galaxy bias: the complicated relation between galaxy and matter
overdensities. However, if reasonably accurate redshifts are available (be
them spectroscopic or photometric), the growth of structure can be
reconstructed through redshift-space distortions (Guzzo et al., 2008; Blake et
al., 2011; de la Torre et al., 2013; Blake et al., 2013; Howlett et al., 2015;
Okumura et al., 2016; Pezzotta et al., 2017; Alam et al., 2021), or by
combining galaxy clustering and weak lensing (Hu, 2002; de la Torre et al.,
2017; Peacock & Bilicki, 2018; Wilson & White, 2019; Krolewski et al., 2020;
Heymans et al., 2021; White et al., 2022; García-García et al., 2021; Alonso
et al., 2023). The combination of different tracers of the large-scale
structure is thus able to overcome their individual shortcomings and fulfill
the requirements listed above. It is for this reason that multi-tracer large-
scale structure analyses have now become one of the staples of late-Universe
cosmology.

In this context, radio continuum surveys are an interesting and promising
probe. Due to the large instantaneous field of view of modern low frequency
radio interferometers, such surveys cover wide areas of the sky. Unencumbered
by dust extinction, they are able to cover large swathes of the Universe.
Their ability to recover clustering information on gigaparsec scales thus
makes them potentially valuable for specific cosmological science cases, such
as the search for primordial non-Gaussianity (Ferramacho et al., 2014; Alonso
& Ferreira, 2015; Gomes et al., 2020). Furthermore, radio continuum samples
are dominated by Active Galactic Nuclei (AGNs) and Star-Forming Galaxies
(SFGs), and thus their study can shed light on key processes in the formation
and evolution of galaxies. In this sense, the clustering of radio sources has
been used to place constraints on the properties of the different radio
populations, making use of a variety of datasets, such as SUMSS (Blake et al.,
2004b), NVSS (Blake & Wall, 2002; Overzier et al., 2003; Negrello et al.,
2006; Nusser & Tiwari, 2015; Chen & Schwarz, 2016), FIRST (Lindsay et al.,
2014), COSMOS 3GHz (Hale et al., 2018), and TGSS (Dolfi et al., 2019; Rana &
Bagla, 2019). In this work, we will study the clustering of radio galaxies in
the second data release of the LOw-Frequency ARray Two-metre Sky Survey (LoTSS
DR2 Shimwell et al., 2022), extending the analyses carried out making use of
the first data release (DR1, Siewert et al., 2020; Alonso et al., 2021; Tiwari
et al., 2022).

The dominant mechanism for radio emission is synchrotron radiation, which is
characterised by a featureless almost power-law spectrum (Condon, 1992). The
absence of bright emission lines or other sharp features in the radio
spectrum, therefore, precludes redshift measurement. This gives rise to two
key sources of uncertainty in the cosmological analysis of radio continuum
surveys: the evolution of the main properties of the sample (e.g. galaxy bias,
relative fractions of different source types) over the large range of
redshifts covered by the sample, and the detailed description of the sample’s
redshift distribution.

As in the case of optical surveys, some of these uncertainties can be overcome
or mitigated through the use of cross-correlations. For example, their cross-
correlation with optical catalogs can be used to infer the clustering
properties of the radio sample, and to constrain its redshift distribution via
tomography (e.g. Ménard et al., 2013). Our focus in this work will be the
cross-correlation with maps of the CMB lensing convergence (Planck
Collaboration et al., 2020b). CMB lensing, sourced at $z\sim 1100$, receives
contributions from density inhomogeneities covering a wide range of redshifts,
peaking at $z\sim 2$. As such, it is an interesting tracer to cross-correlate
with radio data, one of the few probes able to cover comparable volumes. In
fact, the cross-correlation with radio data from the NRAO VLA Sky Survey
(NVSS, Condon et al., 1998) was used to make the first detection of the CMB
lensing signal (Smith et al., 2007). Since then, this cross-correlation has
been used for the benefit of both probes, e.g. as a way to measure the radio
galaxy bias (Allison et al., 2015; Piccirilli et al., 2023), and in the
context of delensing (Namikawa et al., 2016). As the sensitivity of radio and
CMB experiments increases, and the statistical uncertainties of this cross-
correlation decrease, the joint analysis of radio continuum and CMB lensing
data becomes a more powerful tool for cosmological studies, able to not only
constrain amplitude-like parameters, but also to discriminate between more
nuanced details of the underlying astrophysical model. For example, the
analysis of the first LoTSS data release in combination with CMB lensing data
(Alonso et al., 2021, A21 hereafter) showed that the galaxy auto-correlation
and its cross-correlation with CMB lensing respond differently to changes in
the galaxy bias and to the width of the redshift distribution, two effects
that would otherwise be highly degenerate. The inclusion of CMB lensing data
can therefore shed light on the main systematic uncertainties affecting
continuum surveys as described above. In addition, it provides a way to
measure the global amplitude of matter fluctuations.

In order to constrain the redshift distribution of the LoTSS DR2 radio
sources, we make use of the cross-identifications of radio sources with
multiwavelength data in three LoTSS deep fields (Tasse et al., 2021; Sabater
et al., 2021; Kondapally et al., 2021). These allow for the measurement of
photometric redshifts for over 90 per cent of all deep field sources at flux
densities above $1.5\,{\rm mJy}$, of which about a quarter have spectroscopic
redshift as well (Duncan et al., 2021; Bhardwaj et al., 2023).

In this paper, we will address two main science questions: what is the galaxy
bias within flux-limited radio samples and its redshift evolution? And what is
the amplitude of density fluctuations as probed by radio data? To do so, we
will carry out a joint harmonic-space analysis of two radio samples extracted
from the LoTSS DR2, defined by different flux and signal-to-noise cuts,
together with maps of the CMB lensing convergence provided by the Planck
collaboration. The analysis of the LoTSS sample will follow closely the
treatment described in the companion paper (Hale et al., 2023, H23 hereafter),
focused on constraints from the real-space galaxy auto-correlation. Taking
advantage of the tools developed for that analysis, we will also apply them to
the cross-correlation between LoTSS and CMB primary anisotropies, in order to
place constraints on the integrated Sachs-Wolfe effect (ISW, Sachs & Wolfe,
1967).

This paper is structured as follows: Section 2 lays out the theoretical
background behind the auto- and cross-correlations we will use here. Sections
3 and 4 present the datasets used in this work, as well as the methods used to
analyse them. The measurements of the radio galaxy bias and the tentative
constraints on $\sigma_{8}$ are presented in Section 5. We summarise our
results and conclude in Section 7.

##  2 Theory

Our main observable is a sky map of the galaxy surface density

$\delta_{g}(\hat{\bf n})=\frac{N_{g}(\hat{\bf n})-\bar{N_{g}}}{\bar{N}_{g}},$
(1)

where $\hat{\bf n}$ is a unit vector pointing along a line of sight,
$N_{g}(\hat{\bf n})$ is the number of galaxies along $\hat{\bf n}$ per unit
solid angle, and $\bar{N}_{g}$ is the mean number of galaxies per unit solid
angle. The projected overdensity is related to the three-dimensional galaxy
overdensity $\Delta_{g}$ through (Peebles, 1980)

$\delta_{g}(\hat{\bf n})=\int_{0}^{\infty}dz\,p(z)\,\Delta_{g}(\chi(z)\hat{\bf
n},z),$ (2)

where $\chi(z)$ is the comoving radial distance, and $p(z)$ is the redshift
distribution of the galaxy sample, normalised to 1 when integrated over $z$.

In addition to $\delta_{g}$, we will study maps of the CMB lensing convergence
$\kappa(\hat{\bf n})$, which quantifies the distortion in the trajectories of
the CMB photons caused by the gravitational potential of the intervening
matter structures (Lewis & Challinor, 2006), and is proportional to the
divergence of the deflection in the photon arrival angle
$\boldsymbol{\alpha}$: $\kappa\equiv-\nabla_{\hat{\bf
n}}\cdot\boldsymbol{\alpha}/2$. As such, $\kappa$ is an unbiased tracer of the
matter density fluctuations $\Delta_{m}(\boldsymbol{x},z)$, and is related to
them through:

$\kappa(\hat{\bf
n})=\int_{0}^{\chi_{LSS}}d\chi\frac{3H_{0}^{2}\Omega_{m}}{2a(\chi)}\chi\frac{\chi_{LSS}-\chi}{\chi_{LSS}}\Delta_{m}(\chi\hat{\bf
n},z(\chi)),$ (3)

where $\Omega_{m}$ is the fractional matter density, $a=1/(1+z)$ is the scale
factor, $H_{0}$ is the Hubble constant today, and $\chi_{LSS}$ is the comoving
distance to the surface of last scattering.

Finally, we will also consider the cross-correlation with CMB temperature
anisotropies, which receives a contribution from the so-called Integrated
Sachs-Wolfe effect (Sachs & Wolfe, 1967). Caused by time-varying gravitational
potentials at late times, the ISW leads to an additional temperature
fluctuation of the form

$\left.\frac{\Delta T}{T}\right|_{\rm ISW}(\hat{\bf n})=2\int_{0}^{\chi_{\rm
LSS}}d\chi\,a\,\dot{\phi},$ (4)

where $\dot{\phi}$ is the derivative of the Newtonian potential with respect
to cosmic time. In Fourier space, $\dot{\phi}$ can be related to the matter
overdensity, assuming linear growth, via

$\dot{\phi}({\bf
k},t)=-\frac{3H_{0}^{2}\Omega_{m}}{2a}\frac{H}{k^{2}}(f-1)\,\Delta_{m}({\bf
k},t),$ (5)

where $f(a)\equiv d\log\Delta_{m}/d\log a$ is the “growth rate”, which is
scale-independent in the linear regime.

Consider a generic three dimensional field ($U$) projected onto a sphere with
a kernel $W_{u}$

$u(\hat{\bf n})=\int d\chi W_{u}(\chi)U(\chi\hat{\bf n},z(\chi)).$ (6)

Any such projected quantity can be decomposed in terms of its spherical
harmonic coefficients $u_{\ell m}$, the covariance of which with another field
$V$ is the so-called angular power spectrum ($C_{\ell}^{uv}$). The angular
power spectrum can be related to the power spectrum of the 3D fields
$P_{UV}(k,z)$ through

$C_{\ell}^{uv}=\int\frac{d\chi}{\chi^{2}}W_{u}(\chi)W_{v}(\chi)P_{UV}\left(k_{\ell}(\chi),z(\chi)\right),$
(7)

where $P_{UV}(k,z)$ is the covariance of the Fourier coefficients of $U$ and
$V$, and $k_{\ell}(\chi)\equiv(\ell+1/2)/\chi$. In this formalism, for the
three fields under consideration (galaxy overdensity, CMB lensing convergence,
and ISW), the radial kernels are given by

$\displaystyle W_{g}(\chi)=\frac{H(z)}{c}p(z),$ (8) $\displaystyle
W_{\kappa}(\chi)=f_{\ell}\frac{3H_{0}^{2}\Omega_{m}}{2a}\chi\frac{\chi_{LSS}-\chi}{\chi_{LSS}}\Theta(\chi_{LSS}-\chi),$
(9) $\displaystyle W_{\rm
ISW}(\chi)=\frac{3\,H_{0}^{2}\Omega_{m}}{k_{\ell}^{2}}H(z)(1-f),$ (10)

where $H(z)$ is the expansion rate, $\Theta(x)$ is the Heaviside function, $c$
is the speed of light, and $f_{\ell}$ is the scale-dependent prefactor given
by

$f_{\ell}=\frac{\ell(\ell+1)}{(\ell+1/2)^{2}}=1-\frac{1}{(2\ell+1)^{2}},$ (11)

which significantly differs from unity only for $\ell\lesssim 10$, and
accounts for the fact that $\kappa$ is related to $\Delta_{m}$ through the
angular Laplacian of the gravitational potential $\phi$. Eq. (7) is only valid
in the Limber approximation (Limber, 1953), which holds when the extent of the
radial kernels is much broader than the correlation scale of the matter
fluctuations (which is the case for $W_{g}$, $W_{\kappa}$ and $W_{\rm ISW}$ in
this work).

In order to use the galaxy distribution as a probe of structure, we need to
define its relation to the matter overdensities. This implies developing
models for the 3D power spectra $P_{gg}(k,z)$, and $P_{gm}(k,z)$, as well as
the matter power spectrum $P_{mm}(k,z)$. In this work, we will assume a simple
parametrisation, for which

$P_{gg}(k,z)=b_{g}^{2}(z)P_{mm}(k,z),\hskip
12.0ptP_{gm}(k,z)=b_{g}(z)P_{mm}(k,z),$ (12)

where $b_{g}(z)$ is the linear bias function. The models used to describe the
redshift distribution of our sample, and the redshift evolution of the bias
are described in Section 4.

To compute the linear matter power spectrum, we use the CAMB Boltzmann solver
(Lewis et al., 2000), and estimate the non-linear power spectrum from it using
HALOFIT (Smith et al., 2003; Takahashi et al., 2012). All other theoretical
calculations (e.g. Limber integrals) were carried out using the Core
Cosmological Library (CCL, Chisari et al., 2019). Unless stated otherwise, we
fix all cosmological parameters to the best-fit $\Lambda$CDM values found by
Planck Collaboration et al. (2020a): $\Omega_{c}=0.26503$,
$\Omega_{b}=0.04939$, $h=0.6732$, $\sigma_{8}=0.8111$, $n_{s}=0.96605$.

##  3 Data

###  3.1 LoTSS DR2

The LOw-Frequency ARray (LOFAR) Two-metre Sky Survey (LoTSS) second data
release (DR2) (Shimwell et al., 2022) covers 27% of the northern sky at
120-168 MHz. It consists of 841 pointings, split into two regions separated by
the Galactic plane, spanning 4178 and 1457 square degrees, respectively, and
shown in Figure 1 (top-left). Data reduction was performed using both
direction-dependent and independent calibration pipelines, and the source
catalogue was created with the source finder PyBDSF (Mohan & Rafferty, 2015).
The catalogue derived from the total intensity (Stokes I) maps contains
4,396,228 radio sources.

The completeness and spatial homogeneity of the sample have a complex
dependence on flux morphology and signal-to-noise ratio (SNR), defined as the
ratio of peak flux density per beam and root mean square noise per beam. The
fiducial sample used in this work comprises 1,136,219 galaxies with a flux
density brighter than $1.5\,{\rm mJy}$, and detected with ${\rm SNR}\geq 7.5$,
which H23 finds as the best balance between the number of sources and
variation in data compared to random samples. To test the robustness of the
cosmological constraints depending on the choice of the sample, we will also
present results for an alternative selection corresponding to galaxies above
$2.0\,{\rm mJy}$ with ${\rm SNR}\geq 5$ (958,438 objects). These cuts are
similar to those used in the analysis of the two-point correlation function
(Siewert et al., 2020) and the cross correlation with the CMB lensing of LoTSS
DR1 (Alonso et al., 2021).

As described in H23, a spatial mask was created by removing regions with tiles
that have not been mosaiced together, or have a large number of gaps due to
problems in the reduction process, which would lead to strong spatial
variations in the flux scale. These regions are mostly located by the outer
edges. The resulting unmasked footprint covers 4357 ${\rm deg}^{2}$,
corresponding to a sky fraction $f_{\rm sky}=0.11$. After data cuts and
masking, 896,637 objects remain in our fiducial catalogue ($1.5\,{\rm mJy}$
flux density cut, ${\rm SNR}\geq 7.5$), and 742,692 in the $2.0\,{\rm mJy}$,
${\rm SNR}\geq 5.0$ sample. We note that because we use pixel coordinates to
mask the resulting maps, the final number of objects is slightly different
than in H23 which uses the object coordinates to mask the catalogue.

Figure 1: Maps used in this work, smoothed over a scale of $1^{\circ}$,
strictly for a visualization purpose. Top left: LoTSS DR2 overdensity, top
right: LoTSS DR2 completeness based on randoms, bottom left: Planck CMB
lensing convergence, bottom right: Planck CMB temperature. The overdensity and
completeness maps include the LoTSS DR2 mask (Hale et al., 2023), while CMB
maps include this and an appropriate mask from the Planck survey.

###  3.2 Planck

We use CMB data released as part of the 2018 Planck analysis (Planck
Collaboration et al., 2020a). We use the “minimum variance” (MV) CMB lensing
convergence harmonic coefficients released in Planck Collaboration et al.
(2020b), together with the associated sky mask. The harmonic coefficients are
transformed into a HEALPix (Górski et al., 2005) map with resolution parameter
$N_{\rm side}$, after truncating them to $(\ell,m)<3\,N_{\rm side}$. In our
analysis, we will use a common resolution $N_{\rm side}=512$, corresponding to
a pixel size of $\sim 6.9\,{\rm arcmin}$. We repeated our analysis using
$N_{\rm side}=256$, finding compatible results. The lensing map covers a sky
fraction $f_{\rm sky}\simeq 0.67$, and overlaps with LoTSS over the whole
footprint of the latter.

For the ISW analysis described in Section 5.5, we use the foreground-cleaned
temperature fluctuation map produced through the SMICA component separation
method, described in Planck Collaboration et al. (2020c), as well as its
associated sky mask. Both mask and map were downgraded to the common
resolution $N_{\rm side}=512$.

##  4 Methodology

###  4.1 Maps

Figure 1 shows the maps used in this analysis. In the figure, we smooth them
for visualisation purposes with a Gaussian filter with a full width at half-
maximum of $1^{\circ}$. The LOFAR maps are limited to the region corresponding
to the LOFAR mask, while Planck maps are shown at the intersection between the
LOFAR mask and the corresponding CMB mask. Lower and higher resolution falls
outside of the scales of interest for this analysis, $50<\ell<800$ (see
Section 5.2).

We use catalogues of randomly-generated sources to account for the spatially-
varying survey depth. The object positions in the random catalogue should be
uncorrelated, while tracing the detection rate, which is not uniform across
the footprint. The simulated sources are based on a modified SKA Design
Studies Simulated Skies (SKADS, Wilman et al., 2008, 2010), to account for an
underestimated number of SFGs at the faintest flux densities (e.g. Hale et
al., 2023). The catalogue provides multiple observable properties for
simulated sources, which were used in combination with simulations from
Shimwell et al. (2022) to account for the effects of smearing, variation of
sensitivity due to elevation or declination, location within the mosaic,
proximity to bright sources or edge of the observed field, where the number of
mosaiced pointings is smaller. As we do not split the sources between their
type (AGN/SFG), redshift, or luminosity, it is the input flux density
distribution of the randoms which is the most important, and the modified
SKADS represents it well for 144 MHz sources, to below the source detection
limit of the survey. The process of generating the random catalogues is
described in detail in \citetalias{Hale:2023}. It provides us with simulated
“output” number count maps, where “output” represents galaxies as if they were
detected by LOFAR in an idealized case. We use this “output” map to correct
for depth fluctuations, and as a weight for the resulting galaxy overdensity
map when computing power spectra (i.e., $w_{g}(\hat{\bf n})$ constitutes the
mask of the galaxy overdensity map). The overdensity $\delta_{g}$ is computed
as

$\delta_{g}(\hat{\bf n})=\frac{N_{g}(\hat{\bf n})}{\bar{N}_{g}w_{g}(\hat{\bf
n})}-1,$ (13)

where $N_{g}(\hat{\bf n})$ is the number of galaxies in the pixel lying in the
direction of $\hat{\bf n}$. $\bar{N}_{g}$ is the mean number of objects per
pixel, and is estimated as $\bar{N}_{g}=\langle N_{g}(\hat{\bf
n})\rangle_{\hat{\bf n}}/\langle w_{g}(\hat{\bf n})\rangle_{\hat{\bf n}}$,
where $\langle\cdots\rangle_{\hat{\bf n}}$ represents a mean over all pixels
within the mask.

###  4.2 Redshift distribution

An important ingredient of the analysis is the redshift distribution, $p(z)$,
of the LoTSS sources, necessary to recover the three-dimensional clustering
parameters, which can then be compared with theoretical predictions (eq. 8).
For radio continuum objects, the individual redshifts are not known and cannot
be estimated from radio fluxes. At present, we do not have optical
identifications and photometric redshift estimates for most of the LoTSS DR2
sources. Therefore, we need to model the underlying $p(z)$ in a more indirect
way.

Figure 2: Redshift distribution based on the three deep fields located within
the LoTSS DR2 footprint, for $2\,{\rm mJy}$ and $1.5\,{\rm mJy}$ flux cuts.
The thick lines show the models fitted with Eq. 14, and the shaded areas are a
$1\sigma$ region from the deep fields measurements. The redshift distribution
is limited to $z<6$.

Extragalactic radio sources consist mostly of SFGs and AGNs, although their
fractions vary with both redshift and flux density (see Best et al., 2023).
However, limitations in the multi-wavelength coverage of the sample may lead
to some uncertainty in the redshifts and classification of sources. In order
to calibrate the redshift distribution of our sample, we make use of the LOFAR
deep fields observations (Tasse et al., 2021; Sabater et al., 2021). The Deep
Field data consist of three fields: Boötes, ELAIS and Lockman Hole. For each
field, a smaller region was defined for which there exists deep multi-
wavelength information, of an area equal to 8.6 deg2 in the Boötes field, 6.7
deg2 in ELAIS and 10.3 deg2 in the Lockman Hole field (Kondapally et al.,
2021; Duncan et al., 2021). A redshift and its probability density function
were associated with each source using a hybrid method that combined template
fitting and machine learning (further details can be found in Duncan et al.
2021). The photometric redshift quality is characterised by normalised median
absolute deviation ($\sigma_{\rm{NMAD}}$) ranging from 1.6 to 2% for galaxies
and 6.4 to 7% for AGNs, while the outlier fraction
($|z_{\rm{phot}}-z_{\rm{spec}}|/(1+z_{\rm{spec}})>0.15$) equals around 2% for
galaxies and 20% for AGNs. It is worth noting that $\sim 5\%$ of the sources
satisfying our sample cuts in the deep fields do not have an optical cross-
match.

We estimate the redshift distribution for each flux density cut catalogue
using a technique based on sampling redshift values from the probability
distributions of photometric redshifts, using spectroscopic redshifts where
available. Given the full probability distribution over a redshift range for
each photometric redshift measurement, we sample a single redshift value over
this probability for each object, and build a histogram of such a
distribution, binning in $\Delta z=0.05$. For objects with spectroscopic
redshifts available, we always take the reported value (i.e., equivalent to
zero photo-$z$ uncertainty). We repeat this procedure of histogram creation
for each deep field separately, and the number of histograms created for each
field is proportional to the number of objects in each field, which makes
fields with more observations more significant in the final estimate. We find
that the final results do not change after sampling at least 200 histograms in
total. The final distribution and its statistical uncertainty is given by the
mean and standard deviation calculated over all histogram realisations, and
then normalised to a unit integral over the redshift range $0<z<6$. This
approach to redshift distribution is also described in H23. The method is able
to combine both photometric and spectroscopic redshifts, and ensures a
reasonable estimate of the final uncertainty in the redshift distribution. The
uncertainty estimated with this method accounts both for errors in every
single measurement of the photometric redshift, and for differences in
redshift distributions between the three deep fields. We found that the errors
estimated with this method are significantly larger in comparison to bootstrap
sampling over the probability distributions of photometric redshifts, where
single redshift distributions are calculated as a sum of probability
distributions within the bootstrap samples, and uncertainty is taken as a
standard deviation within those.

Table 1: Constraints on the parameters of the redshift distribution, as given in the equation 14, for the $1.5\,{\rm mJy}$ and $2.0\,{\rm mJy}$ samples. Sample | $z_{0}$ | $r$ | $a$
---|---|---|---
$1.5\,{\rm mJy}$ | ${0.05\pm 0.01}$ | ${0.20\pm 0.03}$ | ${4.9\pm 0.1}$
$2.0\,{\rm mJy}$ | ${0.04\pm 0.01}$ | ${0.17\pm 0.03}$ | ${5.0\pm 0.1}$

We model the resulting redshift distribution using a functional form

$p(z)\propto\frac{z^{2}}{1+z}\left(\exp\left({\frac{-z}{z_{0}}}\right)+\frac{r^{2}}{(1+z)^{a}}\right),$
(14)

normalised to a unit integral over the redshift range $0<z<6$, with
$\\{z_{0},r,a\\}$ being free parameters. This form is motivated by the fact
that the LoTSS radio sources contain two main populations of objects, AGNs,
and SFGs. At low redshifts, we expect their numbers to grow proportionally to
the volume for both populations, which motivates the factor of $z^{2}$, which
would be exact for any redshift in a de Sitter model. The factor $1/(1+z)$
provides a simple correction for a $\Lambda$CDM model, and gives a good
approximation up to the redshift of $\sim 0.2$. For higher redshifts, the flux
density limitation of the sample becomes the dominant aspect, and the form of
the luminosity function for each population starts to be important. The AGN
radio luminosity function is typically approximated by a double power law,
which motivates the power law term, while the SFG radio luminosity function is
typically modelled as a Schechter function (Bonato et al., 2017), which
exhibits an exponential cut-off, and motivates the first term. The relative
fraction of both contributions is controlled by the parameter $r$. We verified
that this three-parameter model provides a good semi-empirical fit and is
superior to other simple parameterisations that have been tested. Table 1
shows the constrained parameters, based on the uncertainties mentioned above,
for the $1.5\,{\rm mJy}$ and $2.0\,{\rm mJy}$ samples, while figure 2 shows
the resulting redshift distributions. The blue and orange bands show the
$1\sigma$ constraints measured from the deep fields for the $1.5\,{\rm mJy}$
and $2\,{\rm mJy}$ cuts, respectively, with the corresponding solid lines
showing the best-fit model of Eq. (14) in each case.

###  4.3 Bias models

Given the wide range of redshifts covered by the samples studied, the
evolution of the linear galaxy bias over that range must be taken into
account. This is non-trivial, as the sample includes several types of
extragalactic sources (SFGs and AGNs, in the simplest description), and their
relative abundances and intrinsic galaxy biases evolve with $z$. To assess the
impact of our assumptions regarding the evolution of the effective bias of the
sample, we will consider three different models (Nusser & Tiwari, 2015; Alonso
et al., 2021):

  * •

A _constant bias_ model $b_{g}(z)=b_{g}$ represents the simplest case.
Although likely an unrealistic model, the corresponding value of $b_{g}$ can
be interpreted as the effective bias of the sample accounting for redshift
evolution.

  * •

A _constant amplitude_ model, in which the bias evolves inversely with the
linear growth factor $D(z)$

$b_{g}(z)=b_{g,D}/D(z).$ (15)

This model has the advantage of reproducing the expected rise in $b_{g}(z)$ at
high $z$ for a flux-limited sample (assuming a monotonic mass-luminosity
relation), while preserving the simplicity of the constant-bias model, with
only a single free parameter. In this model, the amplitude of $\Delta_{g}$
does not change over time at linear order (since $\Delta_{m}\propto D(z)$).
This would correspond to a galaxy distribution that is fixed at some early
time and preserves its large-scale properties unchanged (Bardeen et al., 1986;
Mo & White, 1996; Tegmark & Peebles, 1998; Coil et al., 2004).

  * •

The two previous models fix the redshift evolution of $b_{g}(z)$, allowing
only its overall amplitude to vary. As a more flexible alternative, we will
also use a _quadratic bias_ model, in which

$b_{g}(z)=b_{0}+b_{1}z+b_{2}z^{2},$ (16)

with $\\{b_{0},b_{1},b_{2}\\}$ free parameters.

###  4.4 Power spectra

We use NaMaster (Alonso et al., 2019) to compute the angular power spectra of
fields defined on a limited region of the sphere using the pseudo-$C_{\ell}$
estimator (Peebles, 1973; Hivon et al., 2002).

We calculate the shot-noise contribution to the galaxy auto-correlation before
inverting the pseudo-$C_{\ell}$ mode-coupling matrix, as (Nicola et al., 2020)

$\tilde{N}^{gg}_{\ell}=\frac{\langle w_{g}\rangle}{\bar{N}_{\Omega}},$ (17)

where $\bar{N}_{\Omega}$ is the mean angular number density of galaxies (in
units of ${\rm sr}^{-1}$), and $\langle w_{g}\rangle$ is the value of the mask
averaged across the sky.

There are reasons to expect departures from a purely Poisson shot noise
contribution to $N^{gg}_{\ell}$. Prominently, in the case of radio surveys, a
fraction of sources may have multi-component detections, which effectively
leads to a higher shot-noise amplitude than predicted by Poisson statistics
(Blake et al., 2004a; Tiwari et al., 2022). Additionally, stochastic and non-
local effects in galaxy formation, as well as effects such as halo exclusion,
can lead to similar departures from Poissonian shot noise (Baldauf et al.,
2016; Kokron et al., 2022). To account for these effects, we marginalise over
a free shot noise amplitude $A_{sn}$, which in practice makes the pipeline
sensitive only to non-flat contributions to the galaxy auto-correlation.

To calculate the statistical uncertainties of our measurements of
$C_{\ell}^{xy}$, we use a jackknife resampling procedure (Norberg et al.,
2009). We divide the LoTSS DR2 footprint into 54 similarly-sized rectangular
areas. We find this number of regions to provide a good balance between the
small scale and large scale errors. Then, removing one of these areas at a
time, we calculate the power spectra in the resulting footprint. The power
spectrum covariance is then calculated as

${\sf Cov}(C_{\ell}^{x},C_{\ell}^{y})=\frac{N_{\rm JK}-1}{N_{\rm
JK}}\sum_{i=1}^{N_{\rm
JK}}(C_{\ell}^{x,i}-\bar{C}_{\ell}^{x})(C_{\ell}^{y,i}-\bar{C}_{\ell}^{y}),$
(18)

where $x$ and $y$ stand for $(gg,g\kappa,gT)$, $N_{\rm JK}$ is the number of
jackknife samples, $C_{\ell}^{x,i}$ is the power spectrum measured in the
$i$-th sample, and $\bar{C}_{\ell}^{x}\equiv\sum_{i}C_{\ell}^{x,i}/N_{\rm JK}$
is the average over jackknife samples. To validate this estimate of the
covariance matrix, we compared it with the analytical prediction assuming all
fields studied can be described by Gaussian statistics (e.g. García-García et
al., 2019). Both estimates were found to be in good agreement. We also report
a correlation matrix

${\sf r}_{ij}={\sf Cov}_{ij}\,/\sqrt{{\sf Cov}_{ii}{\sf Cov}_{jj}},$ (19)

where ${\sf r}$ is the correlation coefficient, while $i$ and $j$ are
corresponding indices of the covariance matrix.

###  4.5 Likelihood inference

We assume the power spectrum data follow a Gaussian distribution, and we
estimate the log-likelihood as

$\chi^{2}\equiv-2\log p({\bf d}|{\bf q})=[{\bf d}-{\bf t}({\bf q})]^{T}{\sf
Cov}^{-1}[({\bf d}-{\bf t}({\bf q})],$ (20)

where ${\bf d}$ denotes the data vector, consisting of combinations of
$C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$, as well as the deep field
measurements of the redshift distribution described in Section 4.2. We do not
use $C_{\ell}^{gT}$ as part of the inference scheme, as we only report its
significance. ${\bf t}({\bf q})$ is the theoretical prediction for ${\bf d}$
given a set of parameters ${\bf q}$, describing both the power spectra and the
redshift distribution (as parameterised in Eq. 14). The covariance matrix
${\sf Cov}$ incorporates the correlated uncertainties of the different
elements of ${\bf d}$. We assume that the power spectrum and redshift
distribution measurements are uncorrelated, while retaining all potential
correlations between different power spectra (at different scales and for
different fields). In other words, the covariance matrix elements for
$(p(z),C_{\ell}^{xy})$ are set to zero. The covariance of the measured
redshift distribution was assumed to be diagonal, with errors estimated via
sampling as described in Section 4.2.

We report the significance of the $C^{g\kappa}_{\ell}$ power spectrum and the
ISW signal as the square root of the difference in $\chi^{2}$ between a null
hypothesis, defined as $C_{\ell}=0$, and the best-fit model
($\sqrt{\Delta\chi^{2}}$). We calculate the reduced chi-squared as:

$\chi^{2}_{\nu}=\frac{\chi^{2}}{\nu},$ (21)

where $\nu$ stands for the number of degrees of freedom equal to the number of
observations minus the number of fitted parameters. The number of observations
includes the data points from the $C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$,
while the number of fitted parameters includes the bias parameters, amplitude
of shot noise, and $\sigma_{8}$. Data points from deep fields and parameters
from the redshift distribution modelling are not included while reporting
those statistics. Finally, we report the “probability to exceed” (PTE),
calculated in terms of the $\chi^{2}$ as:

$\mathrm{PTE}(\chi^{2},\nu)=1-F(\chi^{2},\nu),$ (22)

where $F$ denotes the $\chi^{2}$ cumulative distribution function.

Table 2: List of all parameters used in the inference with the corresponding priors and initial values. We require $b_{0}$ and $b_{1}$ to be positive, in order to obtain positive $a$, and to be increasing at low redshifts. We allow $b_{2}$ to be negative, as it can both increase and decrease the bias evolution. The initial values are drawn uniformly from a range centred at the ‘middle point’ value, plus/minus 20%. |  | prior | middle point
---|---|---|---
constant bias | $b_{g}$ | positive | 2.0
const. amplitude bias | $b_{g,D}$ | positive | 1.5
quadratic bias | $b_{0}$ | positive | 1.5
| $b_{1}$ | positive | 1.0
| $b_{2}$ | none | 0.1
redshift distribution | $z_{0}$ | positive | 0.05
| $a$ | none | 5.0
| $r$ | positive | 0.2
shot noise amplitude | $A_{sn}$ | [0.8, 1.4] | 1.1
matter fluctuations | $\sigma_{8}$ | positive | 0.81

To explore the posterior distribution function we make use of rejection
sampling via Monte-Carlo Markov Chains (MCMC) as implemented in the public
emcee code (Foreman-Mackey et al., 2013). Table 2 lists the parameters of
interest explored in Section 5 together with their priors. We free up the
value of $\sigma_{8}$ only in Section 5.4. The MCMC chains were generated
using 32 walkers and a convergence condition ensuring that the number of
samples is equal to or higher than 40 times the mean of the auto-correlation
scale for all the inferred parameters.

##  5 Results

###  5.1 Power spectra

(a)  (b)

Figure 3: Comparison of the linear and HALOFIT matter power spectrum for the
auto- and cross-correlation (left and right, respectively). We note that the
shot noise is reported for each multipole separately, while the correlation
signal is calculated in the bins of 50 multipoles. The solid lines show the
best-fit results with different 3D power spectrum models, while dashed lines
show the models with the same resulting best-fit parameters as obtained for
solid lines, but with only matter power spectrum changed to the other model.
Hence, the difference between the corresponding solid and dashed lines stems
only from a difference between linear and HALOFIT models. The vertical dashed
lines mark the multipole ranges used in this analysis: the fiducial
$50\leq\ell\leq 250$ and $50\leq\ell\leq 500$, as well as larger
$50\leq\ell\leq 500$ and $50\leq\ell\leq 800$, for $C_{\ell}^{gg}$ and
$C_{\ell}^{g\kappa}$ respectively. The fits shown here were made on the
fiducial multipole range, where differences between the linear and HALOFIT
models are between $1\sigma-2\sigma$ of errors on data measurements.

Figure 4: Correlation matrix (${\sf r}_{ij}$, see eq. 19) for $C_{\ell}^{gg}$,
$C_{\ell}^{g\kappa}$ power spectra. Multipole ranges $50\leq\ell\leq 500$ and
$50\leq\ell\leq 800$ are shown, which are the largest ranges used in this
work, in bins of $\Delta\ell=50$.

Figure 3 shows the measurements of the LoTSS DR2 auto-spectrum, and its cross-
correlation with the Planck lensing map (left and right panels respectively).
The solid gray line in the left panel shows the expected contribution from
shot noise. As expected, given the broad redshift range covered by the sample,
the auto-correlation has a featureless, roughly power-law-like behaviour,
which is detected at relatively high significance over all the scales
explored. The cross-correlation is also clearly detected to scales $\ell\sim
800$. Quantifying the significance of this detection as described in Section
4.5 (from the $\chi^{2}$ difference between best-fit model and null
hypothesis), including these scales, we obtain a signal-to-noise of:

$\left(\frac{S}{N}\right)_{\ell\leq 800}={26.6}.$ (23)

This is one of the most significant detections of the cross-correlation
between radio galaxies and CMB lensing so far, comparable to the significance
of the correlation with the NVSS sample over a much larger area (Ade et al.,
2014). Considering only the fiducial scales $\ell\leq 500$ that we will
include in the analysis, the significance is

$\left(\frac{S}{N}\right)_{\ell\leq 500}={23.1}.$ (24)

This is a factor $\sim 3.6$ higher than the detection in A21 using LoTSS DR1,
in good agreement with the expectation given the relative increase in area
between both releases (assuming that $S/N$ scales as $\sqrt{f_{\rm sky}}$).
For the auto-correlation, we obtain the signal-to-noise ratio of
${34.6}\sigma$ at $\ell\leq 500$, and ${17.9}\sigma$ at $\ell\leq 250$.

Figure 3 also shows the best-fit predictions for both power spectra using the
constant amplitude bias model and linear (green line) and HALOFIT (orange
line) predictions, as well as corresponding models resulting from using the
same best-fit parameters as found earlier, but changing only the matter power
spectrum (both dashed lines). Comparing the same colour solid and dashed
lines, linear and HALOFIT predictions begin to differ from one another by more
than $2\sigma$ of the statistical uncertainties in our measurements of
$C_{\ell}^{gg}$ at $\ell=250$, and about $1\sigma$ of the error of
$C_{\ell}^{g\kappa}$ at $\ell=500$ in case of the cross-correlation. Using the
approximation $\theta\simeq 180^{\circ}/\ell$, those scales translate to $\sim
0.72^{\circ}$ and $\sim 0.36^{\circ}$, respectively. On the other hand, at the
median redshift $z_{\rm med}\simeq 0.82$, based on the fitted distribution,
these angular scales correspond to wave numbers $k=0.13\,h\,{\rm Mpc}^{-1}$
and $0.26\,h\,{\rm Mpc}^{-1}$ respectively. We can thus use this to define
conservative scale cuts for which the linear bias model can be considered. Our
fiducial scale cuts will therefore be reported here as $\ell<(250,500)$. To
eliminate any residual systematics associated with large-scale survey depth
variations, we will also remove the first bandpower in the galaxy auto-
correlation.

Table 3: Comparison of bias estimates for constant bias and constant amplitude models, using different power spectrum measurements at the fiducial $\ell<(250,500)$ scale cuts, together with the HALOFIT matter power spectrum. | $b_{g}(z)=b_{g}$ | $b_{g}(z)=b_{g,D}/D(z)$
---|---|---
| $b_{g}$ | $A_{sn}$ | $\chi^{2}_{\nu}$ | PTE | $b_{g,D}$ | $A_{sn}$ | $\chi^{2}_{\nu}$ | PTE
$C_{\ell}^{gg}$ | $1.86^{+0.14}_{-0.14}$ | $0.93^{+0.10}_{-0.08}$ | 1.7 | 19% | $1.53^{+0.09}_{-0.11}$ | $0.93^{+0.11}_{-0.08}$ | 1.7 | 19%
$C_{\ell}^{g\kappa}$ | $2.16^{+0.10}_{-0.09}$ |  | 1.2 | 30% | $1.39^{+0.06}_{-0.06}$ |  | 1.2 | 32%
$C_{\ell}^{gg}$ & $C_{\ell}^{g\kappa}$ | $2.08^{+0.09}_{-0.09}$ | $0.89^{+0.08}_{-0.06}$ | 1.4 | 18% | $1.41^{+0.06}_{-0.05}$ | $1.01^{+0.08}_{-0.09}$ | 1.2 | 25%

It is worth noting, however, that using a linear bias model applied to the
non-linear matter power spectrum has been empirically found to extend the
validity of the model to mildly non-linear scales (Pandey et al., 2020;
Sugiyama et al., 2022; Porredon et al., 2022). In some cases, we will
therefore also report results for less conservative scale cuts
$\ell<(500,800)$ (corresponding to $k_{\rm max}=0.26\,h\,{\rm Mpc}^{-1}$ and
$0.42\,h\,{\rm Mpc}^{-1}$ respectively). We stress, however, that these
results should be interpreted with care, since they rely on the validity of
the linear bias model over mildly non-linear scales. This could be quantified
via numerical simulations including a physics-based model for the galaxy-halo
relation for SFGs and AGNs, but this lies beyond the scope of this work. The
correlation matrix of the joint $(C_{\ell}^{gg},C^{g\kappa}_{\ell})$ data
vector after imposing these scale cuts is shown in Fig. 4. As evidenced by
this plot, the uncertainties between different bandpowers are largely
uncorrelated.

###  5.2 Constraining bias

We now use the measurements presented in the previous section to constrain the
bias of radio sources in the LoTSS DR2 sample. For now, we will fix all
cosmological parameters to the best-fit Planck cosmology (Planck Collaboration
et al., 2020a), and will only vary bias and $p(z)$ parameters.

We begin by comparing our two 1-parameter bias models, the constant-bias and
constant-amplitude (or $1/D(z)$) parametrisations, when constrained by
different combinations of correlation functions, but in all cases using our
fiducial scale cuts $\ell<(250,500)$, and the HALOFIT matter power spectrum.
Table 3 shows the constraints on the bias and the shot-noise amplitude
$A_{sn}$ obtained using our fiducial scale cuts when including only the galaxy
auto-correlation (first row), the cross-correlation (second row), and both
(third row). We find that, while both models are able to provide a good fit to
the auto- and cross-correlations separately, the $1/D(z)$ model provides a
better representation of both signals simultaneously. The combined constraint
${b_{g,D}=1.41\pm{0.06}}$ is compatible with the individual constraints from
$C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$, which are also in agreement with
each other. The constant-bias model, in turn, finds broadly incompatible best-
fit values for $b_{g}$, and the model is a worse fit for the combined data
vector than the $1/D(z)$ model. It is also interesting to note that, with the
conservative scale cuts applied here, the bias is better constrained with
$C_{\ell}^{g\kappa}$ than with the galaxy auto-spectrum. When including the
CMB lensing cross-correlation, there is then significant evidence that the
effective bias of the sample grows with redshift, as would be expected for
most flux-limited samples.

Table 4: Comparison of bias estimates for constant amplitude and quadratic models, using different multipole ranges and modelling of matter power spectrum, obtained with both $C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$. The last column shows the number of data points. |  | $b_{g}(z)=b_{g,D}/D(z)$ | $b_{g}(z)=b_{0}+b_{1}z+b_{2}z^{2}$
---|---|---|---
|  | $b_{g,D}$ | $A_{sn}$ | $\chi^{2}_{\nu}$ | PTE | $b_{0}$ | $b_{1}$ | $b_{2}$ | $A_{sn}$ | $\chi^{2}_{\nu}$ | PTE
$\ell<(250,500)$ | linear | $1.54^{+0.06}_{-0.06}$ | $1.19^{+0.06}_{-0.06}$ | 1.4 | 16% | $1.54^{+0.21}_{-0.24}$ | $0.67^{+0.67}_{-0.46}$ | $0.19^{+0.25}_{-0.28}$ | $1.21^{+0.06}_{-0.07}$ | 1.8 | 6.8%
| HALOFIT | $1.41^{+0.06}_{-0.05}$ | $1.01^{+0.08}_{-0.09}$ | 1.2 | 25% | $1.56^{+0.19}_{-0.21}$ | $0.57^{+0.50}_{-0.39}$ | $0.06^{+0.20}_{-0.17}$ | $0.98^{+0.10}_{-0.09}$ | 1.4 | 17%
$\ell<(500,800)$ | linear | $1.65^{+0.04}_{-0.04}$ | $1.14^{+0.02}_{-0.02}$ | 1.8 | 1.6% | $1.60^{+0.19}_{-0.26}$ | $0.83^{+0.77}_{-0.59}$ | $0.23^{+0.28}_{-0.34}$ | $1.15^{+0.02}_{-0.02}$ | 1.9 | 0.9%
| HALOFIT | $1.44^{+0.04}_{-0.04}$ | $1.04^{+0.02}_{-0.02}$ | 1.5 | 7.8% | $1.51^{+0.16}_{-0.18}$ | $0.65^{+0.53}_{-0.45}$ | $0.11^{+0.23}_{-0.22}$ | $1.04^{+0.02}_{-0.02}$ | 1.6 | 3.9%

Having established that the constant bias model is mildly disfavoured, we use
constant-amplitude and quadratic bias models to further compare the bias
estimates between the linear and HALOFIT models of the power spectrum at
different scales. The results are shown in Table 4. The linear and HALOFIT
models are in broad agreement within $\sim 2\sigma$ at the fiducial scale cuts
$\ell<(250,500)$. At the less conservative scale cuts $\ell<(500,800)$, the
agreement is significantly poorer, and neither model is able to provide a good
fit to the data (with PTEs at around 3% and 13% for linear and HALOFIT models
respectively). This shows that the linear bias assumption employed here is not
a reliable representation of the data on these mildly non-linear scales. It is
worth noting that, for either choice of scale cuts, the linear power spectrum
model achieves a consistently poorer $\chi^{2}$ than HALOFIT
($\Delta\chi^{2}\simeq 6$ for $\ell<(500,800)$, with a significantly worse
PTE). The linear power spectrum model also generally prefers a $10\%$ to
$20\%$ higher value of the shot-noise amplitude $A_{sn}$, to compensate for
the lower small-scale power in comparison with HALOFIT.

Figure 5: Bias constraints for three different models based on $C_{\ell}^{gg}$
and $C_{\ell}^{g\kappa}$ at the fiducial $\ell<(250,500)$, and for the HALOFIT
matter power spectrum. The blue vertical line shows the LoTSS DR2 median
redshift and 68 percentiles based on the deep fields $N(z)$. The markers show
bias estimates for radio galaxies known from the literature (Nusser & Tiwari,
2015; Hale et al., 2018; Alonso et al., 2021; Chakraborty et al., 2020;
Mazumder et al., 2022; Hale et al., 2023). Markers with white fillings and
grey/black borders stand for AGNs/SFGs respectively, while markers with grey
fillings and black borders denote mixed populations.

The combination of $C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$ allows us to
successfully constrain the quadratic bias model. The constraints on the bias
evolution obtained from the joint data vector, analysed under the HALOFIT
model, for the three bias evolution models explored here, are shown in Fig. 5.
The figure also shows other existing estimates of the bias of various radio
galaxy samples (Nusser & Tiwari, 2015; Hale et al., 2018; Chakraborty et al.,
2020; Mazumder et al., 2022), which are of different depth and different
ratios of AGNs and SFGs, as well as the results obtained in H23 using the
galaxy correlation function. The figure also shows, as a vertical shaded band,
the median and 68 percentiles of the redshift distribution estimated from the
deep fields. Using the median survey redshift, the $1/D(z)$ model predicts a
value of the bias ${b(z=0.82)=2.34\pm{0.10}}$. This is in good agreement with
the prediction from the other two bias models at the same redshift. Future
releases of LoTSS data may allow us to better constrain the values of the
quadratic bias model, which imposes fewer assumptions on bias evolution. With
the current data, the predictions of the quadratic bias model are in good
agreement with those of the $1/D(z)$ model. From these results, we are led to
conclude that the more reliable setup to carry out cosmological analyses with
the LoTSS data is to adopt the conservative scale cuts $\ell<(250,500)$ paired
with the HALOFIT matter power spectrum and the $1/D(z)$ bias model, as the
simplest parametrisation able to describe all the data used in the analysis.

The different results from the literature shown in Fig. 5 are in rough
agreement with our measurements. It is worth noting, however, that each of
these works was carried out on samples of radio galaxies with different depths
and ratios of AGNs and SFGs, also using different bias parametrisations, and
hence a direct comparison is not possible. The most direct comparisons can be
made with the measurements of A21, using a sample at flux higher than $2\,{\rm
mJy}$ and an SNR higher than $5$ defined on the LoTSS DR1 catalog, and with
the estimate of H23, based on the angular auto-correlation function of LoTSS
DR2 for a flux limit of $1.5\,{\rm mJy}$ and an SNR cut of $7.5$. A21 found
the joint constraints on galaxy bias to depend strongly on the assumed
redshift distribution of the sample, although the cross-correlation alone was
extremely robust against this systematic. In this case, for the $1/D(z)$
model, A21 finds $b_{g,D}=1.46\pm 0.28$, in agreement with our findings
(${b_{g,D}=1.40\pm 0.07}$ at $2.0\,{\rm mJy}$, at $5\,\rm{SNR}$). H23 used a
setup that is much more similar to ours, although only studying the LoTSS
auto-correlation. Using the linear matter power spectrum and a scale range
$36<\ell<360$ (assuming a conversion $\ell=180^{\circ}/\theta$ between real
and harmonic scales), H23 found $b_{g,D}=1.79^{+0.15}_{-0.14}$, which is at
$1.3\sigma$ difference from our corresponding setup which yields
$b_{g,D}=1.61\pm 0.11$, while using only the auto-correlation at $50<\ell<250$
and the linear matter power spectrum. Our bias estimate becomes lower,
$b_{g,D}=1.54\pm 0.06$, and the difference becomes larger, at a level of
$1.8\sigma$, if we add the cross-correlation, as shown in the first row of
table 4. The difference becomes even larger, at a level of $2.3\sigma$, if we
assume the HALOFIT matter power spectrum, which gives better fits in our case
and evaluates to ${b_{g,D}=1.41\pm{0.06}}$, in comparison to
$b_{g,D}=1.75^{+0.16}_{-0.15}$ from H23, while also using the HALOFIT. We
follow upon this in Sect. 6.

###  5.3 Bias evolution and clustering redshifts

Figure 6: Comparison of our fit to redshift distribution and bias from our
fiducial approach and $b_{g,D}/D(z)$ bias modelling, against the results from
Tomographer.

Since all the bias models explored here assume some form of bias evolution, we
carry out one additional test of their validity, making use of the clustering
redshifts technique (Newman, 2008; Ménard et al., 2013; Scottez et al., 2016).
Clustering-based redshift estimation uses a set of angular cross-correlations
between a sample for which the redshift distribution is unknown and a
reference spectroscopic sample with known redshifts, to infer the unknown
redshift distribution. Due to the degeneracy between the redshift distribution
and the galaxy bias of the target sample in setting the amplitude of the
cross-correlations, the technique is in fact only able to constrain the
combination $b_{g}(z)\,p(z)$. Although this degeneracy with bias evolution is
one of the drawbacks of the clustering redshifts technique, we can use it to
our advantage in order to validate our assumptions regarding bias evolution.
We estimate $b_{g}(z)\,p(z)$ for our sample using the public tool
Tomographer111http://tomographer.org.(Chiang & Ménard, 2019), which uses about
2 million spectroscopic objects covering about 10,000 square degrees, based on
samples of galaxies and quasars from the Sloan Digital Sky Survey (SDSS,
Strauss et al., 2002; Blanton et al., 2005; Schneider et al., 2010; Reid et
al., 2016; Pâris et al., 2017; Ata et al., 2018; Bautista et al., 2018). We
compare the result from Tomographer with the product of the redshift
distribution obtained from the deep fields and the best-fit $1/D(z)$ bias
model. Figure 6 shows the results, with Tomographer measurements shown as
points with error bars, and the $68\%$ confidence interval of our estimated
$b_{g}(z)\,p(z)$ shown as an orange band. Both results are in broad agreement,
but there is some potential evidence of a higher bias at $z\gtrsim 1.5$, which
could be confirmed with future LoTSS data releases, or with a dedicated cross-
correlation analysis involving a dense optical galaxy sample at those
redshifts (e.g. Meisner et al., 2018; Storey-Fisher et al., 2023).

###  5.4 Constraining $\sigma_{8}$

Figure 7: Constraints on $\sigma_{8}$ using $C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$ at the fiducial $\ell<(250,500)$, the $b_{g,D}/D(z)$ bias modelling, HALOFIT matter power spectrum, and Planck cosmology assumed for parameters other than $\sigma_{8}$, The top bar shows constraints form the LoTSS DR1 (Alonso et al., 2021), and the three bottom bars present Planck (Planck Collaboration et al., 2020a), KiDS (Heymans et al., 2021), and DES (Abbott et al., 2022). Table 5: Comparison of $\sigma_{8}$ estimates at two different choices of data cuts, using both $C_{\ell}^{gg}$ and $C_{\ell}^{g\kappa}$ and the fiducial scale cut $\ell<(250,500)$ and HALOFIT matter power spectrum. The number of data points in the correlation functions is 13. Sample | Matter power spectrum | $\sigma_{8}$ | $b_{g,D}$ | $A_{sn}$ | $\chi^{2}$ | PTE
---|---|---|---|---|---|---
($1.5\,{\rm mJy}$, ${\rm SNR}>7.5$) | HALOFIT | $0.75^{+0.05}_{-0.04}$ | $1.62^{+0.21}_{-0.19}$ | $0.97^{+0.09}_{-0.09}$ | 1.2 | 26%
| linear | $0.79^{+0.05}_{-0.05}$ | $1.62^{+0.20}_{-0.18}$ | $1.17^{+0.08}_{-0.08}$ | 1.5 | 12%
($2.0\,{\rm mJy}$, ${\rm SNR}>5.0$) | HALOFIT | $0.82^{+0.08}_{-0.07}$ | $1.38^{+0.25}_{-0.22}$ | $1.17^{+0.10}_{-0.09}$ | 1.4 | 20%
| linear | $0.82^{+0.06}_{-0.06}$ | $1.49^{+0.20}_{-0.18}$ | $1.30^{+0.06}_{-0.08}$ | 1.8 | 5.7%

We put constraints on the $\sigma_{8}$ parameter by using $C_{\ell}^{gg}$ and
$C_{\ell}^{g\kappa}$ at the fiducial $\ell<(250,500)$ scale cuts, together
with the deep fields $p(z)$, the HALOFIT matter power spectrum, and the
$1/D(z)$ bias model, as justified in previous sections. Our data is not yet
powerful enough to break the degeneracy between different cosmological
parameters, and therefore we only vary the amplitude of matter fluctuations,
parametrised by $\sigma_{8}$. Our measurement thus corresponds to an
independent constraint on the growth of structure at low redshifts, assuming
that CMB data can reliably constrain all background evolution parameters
($\Omega_{c}$, $\Omega_{b}$, $H_{0}$, etc.). Combinations with other datasets
(e.g. BAO measurements) may allow us to break these degeneracies independently
from the CMB, but we leave this analysis for future work.

The resulting $68\%$ CL constraints on $\sigma_{8}$ are

${\sigma_{8}=0.75^{+0.05}_{-0.04}}.$ (25)

The full marginalised distribution is shown in Figure 7, together with the
constraints from Planck (Planck Collaboration et al., 2020a), as well as the
Kilo-Degree Survey (KiDS, Heymans et al., 2021), and the Dark Energy Survey
(DES, Abbott et al., 2022). The fiducial measurement (first from the top in
Fig. 7) is in agreement with the weak lensing surveys, and at $1.2\sigma$
difference from the CMB constraints by Planck. The measurement using linear
matter power spectrum (second from the top in Figure 7) is in agreement with
the fiducial setup, which uses the HALOFIT modelling, but it is closer than
the HALOFIT to the measurements from Planck. To test the robustness of this
result to the choice of galaxy sample, we repeat the analysis for the fiducial
flux and SNR cuts of A21 ($2.0\,{\rm mJy}$, 5.0 respectively, third and fourth
from the top in Figure 7). We obtain ${\sigma_{8}=0.82^{+0.08}_{-0.07}}$, in
agreement with the result found for the fiducial sample, but we note higher
uncertainty of the estimations using this additional sample. These results are
summarised in Table 5. The full posterior distribution of all model parameters
is shown in Figure 9.

###  5.5 Integrated Sachs-Wolfe effect

Figure 8: Cross-correlation with the CMB temperature, based on fiducial bias
constraints from the section 5.2, using scales $2<\ell<50$. The signal-to-
noise ratio is consistent with zero.

The ISW signal is especially useful for cosmology, as it is sensitive to the
dark energy equation of state. However, sufficiently large sky coverage of the
galaxy sample is needed for it to be detected, which is not yet the case for
LoTSS. Figure 8 shows the measured cross-correlation between our LoTSS DR2
sample limited at $1.5\,{\rm mJy}$ and the CMB temperature anisotropies
measured by Planck. This measurement was carried out using thinner $\ell$ bins
($\Delta\ell=16$) to concentrate on the largest scales, where the ISW signal
is the most significant. The orange lines in the same plot show the
theoretical prediction for values of the galaxy bias selected by the MCMC
chains run in Section 5.2 for the $1/D(z)$ model using HALOFIT. Fixing the
galaxy bias to the best-fit value found in Section 5.2, and comparing the
$\chi^{2}$ value of the measured $C_{\ell}^{gT}$ with respect to the null-
hypothesis and the best-fit model, we determine that the ISW signal is not
significantly detected. A higher-significance measurement of this signal can
be expected with future releases of LoTSS covering the full northern sky.

##  6 Discussion

We have shown that our results are in reasonable agreement with previous
measurements of the galaxy bias for various radio samples, including previous
CMB lensing cross-correlation analyses. The comparison with the real-space
analysis of LoTSS DR2 carried out in H23 shows $1.3\sigma$ difference if we
use only the auto-correlation and linear matter power spectrum, $1.8\sigma$
difference if we add the cross-correlation, and $2.3\sigma$ difference if we
assume the HALOFIT modelling for both approaches, which provides better fits
in our case. As shown in Hamana et al. (2022), the difference in
$S_{8}=\Omega_{m}^{1/2}\sigma_{8}$ estimates between the real and harmonic
space can be even larger than $1\sigma$. In our case, there are several
possible reasons for the resulting difference.

  * •

The pure sample variance is due to the fact that both analyses actually do not
use the same modes. We only expect harmonic space and real space methods to
agree for full sky coverage or isotropic sampling of a statistically isotropic
universe. The second aspect, the isotropic sampling, is indeed violated, and
we tried to correct it by means of the weight mask and data cuts that we
applied. However, the difference between the real and harmonic space analyses
can point out that we did not correct for all the large-scale systematics.

  * •

Both approaches treat the multi-component sources in different ways. In our
case, it is a marginalisation over the amplitude of shot noise, whereas H23
selects the scales which allow avoiding the effects resulting from the multi-
component sources.

  * •

The angular two-point correlation function can be affected by contamination
from a dipole. Chen & Schwarz (2016) show that an excess of two-point
correlation at the degree scale in the NVSS data set can be removed by
properly removing the NVSS dipole before analysing the two-point correlation.
A study of that issue in the scope of the LoTSS survey will be published in
Böhme et al. (2023).

##  7 Conclusion

We combined LoTSS DR2 wide field and LoTSS DR1 deep fields, supplemented by
multiwavelength data, with gravitational lensing from the Planck Cosmic
Microwave Background (CMB) to place constraints on the bias and its evolution
for radio galaxies, and on the amplitude of matter perturbations. Our main
results can be summarised as follows:

  * •

We obtain one of the most significant detections of the cross-correlation
between radio and CMB lensing data, resulting in the signal-to-noise ratio at
a level of ${26.6}\sigma$.

  * •

We show that the inclusion of CMB lensing information leads to a clear
preference for an evolving galaxy bias, growing towards higher redshifts, as
expected from linear theory. We determine that a linear bias evolution of the
form $b_{g}(z)=b_{g,D}/D(z)$, where $D(z)$ is the linear growth factor, is
able to consistently provide a good description of different sectors of the
data. This allows us to measure ${b_{g,D}=1.41\pm{0.06}}$, which evaluates to
${b(z=0.82)=2.34\pm{0.10}}$ at the median survey redshift, for a sample flux-
limited at $1.5\,{\rm mJy}$. These results are also in good agreement with
more flexible bias parametrisations (e.g. a quadratic polynomial in
redshifts), which lead to similar constraints.

  * •

Freeing up the value of $\sigma_{8}$, we are able to constrain it to
${\sigma_{8}=0.75^{+0.05}_{-0.04}}$ using our fiducial sample. The result is
in good agreement with weak lensing surveys – KiDS (Asgari et al., 2021;
Heymans et al., 2021) and DES (Abbott et al., 2022), as well as CMB data from
Planck.

  * •

We attempt a first measurement of the ISW signal with LOFAR data, but find
that the signal is compatible with zero.

Throughout this analysis, we used conservative scale cuts, $\ell<250$ for
$C_{\ell}^{gg}$, and $\ell<500$ for $C_{\ell}^{g\kappa}$, and showed that for
more permissive cuts ($\ell<(500,800)$), including mildly non-linear scales,
the simple linear bias models used here are not able to fit the data
adequately. More work is needed in order to provide a robust model for the
bias of radio galaxies that extends to non-linear scales. This could be done
by making use of perturbative bias expansions (Matsubara, 2008; Desjacques et
al., 2018), or phenomenological halo-based models (Peacock & Smith, 2000;
Berlind & Weinberg, 2002). However, additional information, in the form of
cross correlations with other tracers (e.g. optical redshift surveys,
tomographic cosmic shear data) will be necessary in order to disentangle non-
linear bias from evolutionary effects.

The strength of our approach comes from the ability of high-resolution radio
surveys to detect galaxies at high redshifts, and from the combination of both
auto- and cross-correlations with CMB lensing, which allows us to break
degeneracies between the amplitude of matter perturbations and galaxy bias,
and to potentially constrain the redshift evolution of the latter. However, by
far the largest source of systematic uncertainty in our results is the lack of
redshift information for radio continuum samples. We have modelled and
calibrated the redshift distribution of our sample using three deep fields
within the LoTSS DR2 footprint. The resulting redshift distribution is still
subject to caveats, due to the use of photometric redshifts, the small area
covered by the deep fields, and the uncertainty due to radio sources with no
optical cross-matches (around $5\%$ of our sample). In our analysis, however,
we have propagated the $p(z)$ calibration uncertainties, by making the
redshift distribution measurements part of the data vector, modelled together
with the galaxy and CMB lensing power spectra. Future LoTSS data releases will
include one additional deep field and even deeper observations of those
fields, which will likely help reduce this source of uncertainty.

The current LoTSS catalog does not allow us to make a significant detection of
the ISW effect. However, future data releases covering the majority of
northern sky, should allow us to improve this result. If included in the
cosmological analysis, this measurement could help improve cosmological
constraints, particularly in the context of dark energy, which is an important
source of the ISW signal.

In this work, we have demonstrated the significant improvement on cosmological
and astrophysical constraints from radio continuum data enabled by the
inclusion of CMB lensing cross-correlations. This will allow future LOFAR data
releases to start providing meaningful constraints on cosmological parameters,
on par with (and in combination with) other probes of the large-scale
structure.

###### Acknowledgements.

SJN is supported by the US National Science Foundation (NSF) through grant
AST-2108402, and the Polish National Science Centre through grant
UMO-2018/31/N/ST9/03975. DA acknowledges support from the Beecroft Trust, and
from the Science and Technology Facilities Council through an Ernest
Rutherford Fellowship, grant reference ST/P004474. MBi is supported by the
Polish National Science Center through grants no. 2020/38/E/ST9/00395,
2018/30/E/ST9/00698, 2018/31/G/ST9/03388 and 2020/39/B/ST9/03494. DJS
acknowledges support from the Bundesministerium für Bildung und Forschung
(BMBF) ErUM-Pro grant 05A20PB1 and Ministerium für Kultur und Wissenschaft des
Landes Nordrhein-Westfahlen Profilbildung 2020 grant B3D. CLH acknowledges
support from the Leverhulme Trust through an Early Career Research Fellowship.
AP is supported by the Polish National Science Centre grant
UMO-2018/30/M/ST9/00757, and by COST (European Cooperation in Science and
Technology) through grant COST Action CA21136 – “Addressing observational
tensions in cosmology with systematics and fundamental physics (CosmoVerse)”.
AP and MBi acknowledge support from the Polish Ministry for Science and Higher
Education through grant DIR/WK/2018/12. CSH’s work is funded by the Volkswagen
Foundation. CSH acknowledges additional support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306 and EXC 2181/1 -
390900948 (the Heidelberg STRUCTURES Excellence Cluster). CLH acknowledges
support from the Leverhulme Trust through an Early Career Research Fellowship.
PT acknowledges the support of the RFIS grant (No. 12150410322) by the
National Natural Science Foundation of China (NSFC). JZ is supported by the
project “NRW-Cluster for data intensive radio astronomy: Big Bang to Big Data
(B3D)“funded through the programme “Profilbildung 2020”, an initiative of the
Ministry of Culture and Science of the State of North Rhine-Westphalia. MBr
acknowledges support from the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy - EXC 2121 “Quantum Universe” - 390833306 and from the
BMBF ErUM-Pro grant 05A2023. MJJ acknowledges support of the STFC consolidated
grant [ST/S000488/1] and [ST/W000903/1] and from a UKRI Frontiers Research
Grant [EP/X026639/1]. MJJ also acknowledges support from the Oxford Hintze
Centre for Astrophysical Surveys which is funded through generous support from
the Hintze Family Charitable Foundation. LOFAR is the Low Frequency Array
designed and constructed by ASTRON. It has observing, data processing, and
data storage facilities in several countries, which are owned by various par-
ties (each with their own funding sources), and which are collectively
operated by the ILT foundation under a joint scientific policy. The ILT
resources have benefited from the following recent major funding sources:
CNRS-INSU, Observatoire de Paris and Univer- sité d’Orléans, France; BMBF,
MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of
Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The
Science and Technology Facilities Council, UK; Ministry of Science and Higher
Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. This
research was carried out using Python 3 (Van Rossum & Drake, 2009) and a
number of software packages from which we enumerate the most significant ones
for our analysis: healpy (Zonca et al., 2019), HEALPix (Górski et al., 2005),
Astropy (Astropy Collaboration et al., 2013, 2018, 2022), pymaster (Alonso et
al., 2019), pyccl (Chisari et al., 2019), emcee (Foreman-Mackey et al., 2013),
getdist (Lewis, 2019), NumPy (Harris et al., 2020), SciPy (Virtanen et al.,
2020), IPython (Pérez & Granger, 2007), Pandas (Wes McKinney, 2010),
Matplotlib (Hunter, 2007), seaborn (Waskom, 2021), tqdm (da Costa-Luis, 2019).

## References

  * Abbott et al. (2022) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520
  * Ade et al. (2014) Ade, P. A. R. et al. 2014, Astron. Astrophys., 571, A17
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533
  * Allison et al. (2015) Allison, R. et al. 2015, Mon. Not. Roy. Astron. Soc., 451, 849
  * Alonso et al. (2021) Alonso, D., Bellini, E., Hale, C., Jarvis, M. J., & Schwarz, D. J. 2021, Mon. Not. Roy. Astron. Soc., 502, 876
  * Alonso et al. (2023) Alonso, D., Fabbian, G., Storey-Fisher, K., et al. 2023, arXiv e-prints, arXiv:2306.17748
  * Alonso & Ferreira (2015) Alonso, D. & Ferreira, P. G. 2015, Phys. Rev. D, 92, 063525
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127
  * Asgari et al. (2021) Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Ata et al. (2018) Ata, M., Baumgarten, F., Bautista, J., et al. 2018, MNRAS, 473, 4773
  * Baldauf et al. (2016) Baldauf, T., Codis, S., Desjacques, V., & Pichon, C. 2016, MNRAS, 456, 3985
  * Bardeen et al. (1986) Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15
  * Bartelmann & Schneider (2001) Bartelmann, M. & Schneider, P. 2001, Phys. Rep, 340, 291
  * Bautista et al. (2018) Bautista, J. E., Vargas-Magaña, M., Dawson, K. S., et al. 2018, ApJ, 863, 110
  * Berlind & Weinberg (2002) Berlind, A. A. & Weinberg, D. H. 2002, Astrophys. J., 575, 587
  * Best et al. (2023) Best, P. N., Kondapally, R., Williams, W. L., et al. 2023, MNRAS, 523, 1729
  * Bhardwaj et al. (2023) Bhardwaj, N. et al. 2023, In preparation
  * Blake et al. (2013) Blake, C., Baldry, I. K., Bland-Hawthorn, J., et al. 2013, MNRAS, 436, 3089
  * Blake et al. (2011) Blake, C., Brough, S., Colless, M., et al. 2011, MNRAS, 415, 2876
  * Blake et al. (2004a) Blake, C., Ferreira, P. G., & Borrill, J. 2004a, MNRAS, 351, 923
  * Blake et al. (2004b) Blake, C., Mauch, T., & Sadler, E. M. 2004b, MNRAS, 347, 787
  * Blake & Wall (2002) Blake, C. & Wall, J. 2002, MNRAS, 337, 993
  * Blanton et al. (2005) Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208
  * Bonato et al. (2017) Bonato, M., Negrello, M., Mancuso, C., et al. 2017, MNRAS, 469, 1912
  * Böhme et al. (2023) Böhme, L. et al. 2023, in preparation
  * Chakraborty et al. (2020) Chakraborty, A., Dutta, P., Datta, A., & Roy, N. 2020, MNRAS, 494, 3392
  * Chen & Schwarz (2016) Chen, S. & Schwarz, D. J. 2016, A&A, 591, A135
  * Chiang & Ménard (2019) Chiang, Y.-K. & Ménard, B. 2019, ApJ, 870, 120
  * Chisari et al. (2019) Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2
  * Coil et al. (2004) Coil, A. L., Davis, M., Madgwick, D. S., et al. 2004, ApJ, 609, 525
  * Condon (1992) Condon, J. J. 1992, ARA&A, 30, 575
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693
  * da Costa-Luis (2019) da Costa-Luis, C. 2019, The Journal of Open Source Software, 4, 1277
  * de la Torre et al. (2013) de la Torre, S., Guzzo, L., Peacock, J. A., et al. 2013, A&A, 557, A54
  * de la Torre et al. (2017) de la Torre, S., Jullo, E., Giocoli, C., et al. 2017, A&A, 608, A44
  * Desjacques et al. (2018) Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rept., 733, 1
  * Dolfi et al. (2019) Dolfi, A., Branchini, E., Bilicki, M., et al. 2019, A&A, 623, A148
  * Duncan et al. (2021) Duncan, K. J., Kondapally, R., Brown, M. J. I., et al. 2021, A&A, 648, A4
  * Ferramacho et al. (2014) Ferramacho, L. D., Santos, M. G., Jarvis, M. J., & Camera, S. 2014, MNRAS, 442, 2511
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * García-García et al. (2019) García-García, C., Alonso, D., & Bellini, E. 2019, JCAP, 11, 043
  * García-García et al. (2021) García-García, C., Ruiz-Zapatero, J., Alonso, D., et al. 2021, J. Cosmology Astropart. Phys., 2021, 030
  * Gomes et al. (2020) Gomes, Z., Camera, S., Jarvis, M. J., Hale, C., & Fonseca, J. 2020, MNRAS, 492, 1513
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Guzzo et al. (2008) Guzzo, L., Pierleoni, M., Meneux, B., et al. 2008, Nature, 451, 541
  * Hale et al. (2018) Hale, C. L., Jarvis, M. J., Delvecchio, I., et al. 2018, MNRAS, 474, 4133
  * Hale et al. (2023) Hale, C. L., Whittam, I. H., Jarvis, M. J., et al. 2023, MNRAS, 520, 2668
  * Hale et al. (2023) Hale, C. L. et al. 2023, In preparation
  * Hamana et al. (2022) Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2022, PASJ, 74, 488
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Heymans et al. (2021) Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140
  * Hivon et al. (2002) Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2
  * Howlett et al. (2015) Howlett, C., Ross, A. J., Samushia, L., Percival, W. J., & Manera, M. 2015, MNRAS, 449, 848
  * Hu (2002) Hu, W. 2002, Phys. Rev. D, 66, 083515
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Huterer (2022) Huterer, D. 2022, arXiv e-prints, arXiv:2212.05003
  * Kokron et al. (2022) Kokron, N., DeRose, J., Chen, S.-F., White, M., & Wechsler, R. H. 2022, MNRAS, 514, 2198
  * Kondapally et al. (2021) Kondapally, R., Best, P. N., Hardcastle, M. J., et al. 2021, A&A, 648, A3
  * Krolewski et al. (2020) Krolewski, A., Ferraro, S., Schlafly, E. F., & White, M. 2020, J. Cosmology Astropart. Phys., 2020, 047
  * Lewis (2019) Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970
  * Lewis & Challinor (2006) Lewis, A. & Challinor, A. 2006, Phys. Rep, 429, 1
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
  * Limber (1953) Limber, D. N. 1953, ApJ, 117, 134
  * Lindsay et al. (2014) Lindsay, S. N., Jarvis, M. J., Santos, M. G., et al. 2014, MNRAS, 440, 1527
  * Matsubara (2008) Matsubara, T. 2008, Phys. Rev. D, 78, 083519, [Erratum: Phys.Rev.D 78, 109901 (2008)]
  * Mazumder et al. (2022) Mazumder, A., Chakraborty, A., & Datta, A. 2022, MNRAS, 517, 3407
  * Meisner et al. (2018) Meisner, A. M., Lang, D., & Schlegel, D. J. 2018, Research Notes of the American Astronomical Society, 2, 1
  * Ménard et al. (2013) Ménard, B., Scranton, R., Schmidt, S., et al. 2013, arXiv e-prints, arXiv:1303.4722
  * Mo & White (1996) Mo, H. J. & White, S. D. M. 1996, MNRAS, 282, 347
  * Mohan & Rafferty (2015) Mohan, N. & Rafferty, D. 2015, PyBDSF: Python Blob Detection and Source Finder, Astrophysics Source Code Library, record ascl:1502.007
  * Namikawa et al. (2016) Namikawa, T., Yamauchi, D., Sherwin, B., & Nagata, R. 2016, Phys. Rev. D, 93, 043527
  * Negrello et al. (2006) Negrello, M., Magliocchetti, M., & De Zotti, G. 2006, MNRAS, 368, 935
  * Newman (2008) Newman, J. A. 2008, ApJ, 684, 88
  * Nicola et al. (2020) Nicola, A., Alonso, D., Sánchez, J., et al. 2020, J. Cosmology Astropart. Phys., 2020, 044
  * Norberg et al. (2009) Norberg, P., Baugh, C. M., Gaztañaga, E., & Croton, D. J. 2009, MNRAS, 396, 19
  * Nusser & Tiwari (2015) Nusser, A. & Tiwari, P. 2015, ApJ, 812, 85
  * Okumura et al. (2016) Okumura, T., Hikage, C., Totani, T., et al. 2016, PASJ, 68, 38
  * Overzier et al. (2003) Overzier, R. A., Röttgering, H. J. A., Rengelink, R. B., & Wilman, R. J. 2003, A&A, 405, 53
  * Pandey et al. (2020) Pandey, S., Krause, E., Jain, B., et al. 2020, Phys. Rev. D, 102, 123522
  * Pâris et al. (2017) Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79
  * Peacock & Bilicki (2018) Peacock, J. A. & Bilicki, M. 2018, MNRAS, 481, 1133
  * Peacock & Smith (2000) Peacock, J. A. & Smith, R. E. 2000, Mon. Not. Roy. Astron. Soc., 318, 1144
  * Peebles (1973) Peebles, P. J. E. 1973, ApJ, 185, 413
  * Peebles (1980) Peebles, P. J. E. 1980, The large-scale structure of the universe
  * Pérez & Granger (2007) Pérez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * Pezzotta et al. (2017) Pezzotta, A., de la Torre, S., Bel, J., et al. 2017, A&A, 604, A33
  * Piccirilli et al. (2023) Piccirilli, G., Migliaccio, M., Branchini, E., & Dolfi, A. 2023, A&A, 671, A42
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A1
  * Planck Collaboration et al. (2020b) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A8
  * Planck Collaboration et al. (2020c) Planck Collaboration, Akrami, Y., Ashdown, M., et al. 2020c, A&A, 641, A4
  * Porredon et al. (2022) Porredon, A., Crocce, M., Elvin-Poole, J., et al. 2022, Phys. Rev. D, 106, 103530
  * Rana & Bagla (2019) Rana, S. & Bagla, J. S. 2019, MNRAS, 485, 5891
  * Reid et al. (2016) Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553
  * Sabater et al. (2021) Sabater, J., Best, P. N., Tasse, C., et al. 2021, A&A, 648, A2
  * Sachs & Wolfe (1967) Sachs, R. K. & Wolfe, A. M. 1967, ApJ, 147, 73
  * Schneider et al. (2010) Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, AJ, 139, 2360
  * Scottez et al. (2016) Scottez, V., Mellier, Y., Granett, B. R., et al. 2016, MNRAS, 462, 1683
  * Shimwell et al. (2022) Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1
  * Siewert et al. (2020) Siewert, T. M., Hale, C., Bhardwaj, N., et al. 2020, A&A, 643, A100
  * Smith et al. (2007) Smith, K. M., Zahn, O., & Dore, O. 2007, Phys. Rev. D, 76, 043510
  * Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311
  * Storey-Fisher et al. (2023) Storey-Fisher, K., Hogg, D. W., Rix, H.-W., et al. 2023, arXiv e-prints, arXiv:2306.17749
  * Strauss et al. (2002) Strauss, M. A., Weinberg, D. H., Lupton, R. H., et al. 2002, AJ, 124, 1810
  * Sugiyama et al. (2022) Sugiyama, S., Takada, M., Miyatake, H., et al. 2022, Phys. Rev. D, 105, 123537
  * Takahashi et al. (2012) Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152
  * Tasse et al. (2021) Tasse, C., Shimwell, T., Hardcastle, M. J., et al. 2021, A&A, 648, A1
  * Tegmark & Peebles (1998) Tegmark, M. & Peebles, P. J. E. 1998, ApJ, 500, L79
  * Tiwari et al. (2022) Tiwari, P., Zhao, R., Zheng, J., et al. 2022, ApJ, 928, 38
  * Van Rossum & Drake (2009) Van Rossum, G. & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61
  * White et al. (2022) White, M., Zhou, R., DeRose, J., et al. 2022, J. Cosmology Astropart. Phys., 2022, 007
  * Wilman et al. (2010) Wilman, R. J., Jarvis, M. J., Mauch, T., Rawlings, S., & Hickey, S. 2010, MNRAS, 405, 447
  * Wilman et al. (2008) Wilman, R. J., Miller, L., Jarvis, M. J., et al. 2008, MNRAS, 388, 1335
  * Wilson & White (2019) Wilson, M. J. & White, M. 2019, J. Cosmology Astropart. Phys., 2019, 015
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298

##  Appendix A Posterior probability distributions

Figs. 9 and 10 show the constraints on all model parameters for the 1.5 mJy
and 2.0 mJy samples respectively.

Figure 9: Probability density for constrained parameters, using the fiducial
approach of flux density brighter than $1.5\,\rm{mJy}$, SNR higher than $7.5$,
$\ell<(250,500)$, and HALOFIT matter power spectrum. The vertical line marks
the Planck constraints on $\sigma_{8}=0.81$. Figure 10: Probability density
for constrained parameters, using flux density brighter than $2.0\,\rm{mJy}$,
SNR higher than $5.0$, and the fiducial approach of $\ell<(250,500)$, and
HALOFIT matter power spectrum. The vertical line marks the Planck constraints
on $\sigma_{8}=0.81$.