Constraints on the origin of the radio synchrotron background via angular correlations

The origin of the radio synchrotron background (RSB) is currently unknown. Its understanding might have profound implications in fundamental physics or might reveal a new class of radio emitters. In this work, we consider the scenario in which the RSB is due to extragalactic radio sources and measure the angular cross-correlation of LOFAR images of the diffuse radio sky with matter tracers at different redshifts, provided by galaxy catalogs and CMB lensing. We compare these measured cross-correlations to those expected for models of RSB sources. We find that low-redshift populations of discrete sources are excluded by the data, while higher redshift explanations are compatible with available observations. We also conclude that at least 20\% of the RSB surface brightness level must originate from populations tracing the large-scale distribution of matter in the universe, indicating that at least this fraction of the RSB is of extragalactic origin. Future measurements of the correlation between the RSB and tracers of high-redshift sources will be crucial to constraining the source population of the RSB.


INTRODUCTION
The background level of radio sky brightness, dominated by steepspectrum synchrotron-like emission below 5 GHz and so here referred to as the radio synchrotron background (RSB), is due to some as of now unknown combination of the integrated emission from extragalactic radio sources and a possible large-scale Galactic halo,  RSB =  G +  E .Direct measurement of the level of the radio background, combining Absolute Radiometer for Cosmology, Astrophysics, and Diffuse Emission 2 (ARCADE 2) measurements from 3 to 90 GHz (Fixsen et al. (2011)) with several radio maps at lower frequencies from which an absolute zero level has been inferred (recently summarized in Singal et al. (2023)), reveals a brightness spectrum (Dowell & Taylor (2018)): RSB () = 30.4± 2.6K  310 MHz −2.66±0.04 . (1) ★ E-mail: elisamaria.todarello@unito.it† E-mail: marco.regis@unito.it At  = 140 MHz, this radio background amounts to  RSB ∼ 252 K.
As summarized in, e.g.Singal et al. (2023), it has been apparent for nearly 15 years that this level of surface brightness is several times higher than can be accounted for by known populations of extragalactic radio sources and accepted Galactic emission characteristics.
To estimate the Galactic component we consider Fornengo et al. (2014), and from their best-fits, we obtain  E ∼ 180 K.The contribution from known sources down to a threshold of ∼5 mJy (which is the one of the images used in this work) is ∼ 30 K, and can be obtained by integrating, e.g., the / curve provided in Intema et al. (2017).
Therefore, if the RSB is of extragalactic origin, there remains a contribution from below-threshold sources of  bts ∼ 150 K at 140 MHz, which is much higher than typical extrapolations for known radio populations.We will refer to this unexplained contribution to the surface brightness as the "RSB excess". 1 The observational sta-tistical error on this estimate is small, of order of 10%, whilst there are possible significant systematic errors related to the estimate of the Galactic contribution or to the zero-level calibration of the observations.Discussing those issues is not the goal of this work, and we will normalize the radio background contribution from extragalactic sources below 5 mJy to provide  bts = 150 K.
It is natural to expect that extragalactic sources providing the RSB would follow the matter distribution in the Universe.Their radio emission should therefore exhibit a certain level of correlation with cosmological matter tracers.Analyses of the anisotropy angular power spectrum of the radio background include Offringa et al. (2021) and Cowie et al. (2023) at the same background radio frequencies considered in this analysis, and Holder (2014) at much higher background frequencies (above 4 GHz).The latter work indicated that the observed anisotropy power of the radio background was lower than that which could result from sources following the large-scale distribution of matter in the universe, while the former works indicated, in contrast, an excess of anisotropy power above that which would be produced by known classes of extragalactic radio sources.
In this work, we constrain the properties of radio source populations explaining the RSB excess by computing their cross-correlation with the CMB lensing and with different galaxy catalogs.We measure the angular cross-correlation of LOFAR images at 140 MHz (Offringa et al. 2021;Cowie et al. 2023) with the Planck CMB lensing (Carron et al. 2022), and the 2MPZ (Bilicki et al. 2014), SDSS main (Data Release 12) (Alam et al. 2015) and SDSS quasar (Data Release 16) (Lyke et al. 2020) catalogs.By comparing those measurements to a simple but flexible phenomenological model of extragalactic radio emitters, we derive bounds on the redshift behavior required for a population that successfully provides the RSB excess.
The cross-correlation between the radio sky and low- galaxy catalogs has been discussed before in (Brown et al. 2010;Vernstrom et al. 2017), even though with different goals than here and with shallower images from the Bonn survey (Brown et al. (2010)) and Murchison Widefield Array (Vernstrom et al. 2017).
The paper is organized as follows.In Section 2, we describe the observational data that we are going to use to perform the crosscorrelation measurements.Section 3 introduces the phenomenological models adopted throughout the paper to describe extra-galactic radio sources.The analysis, results, and discussion related to the angular autocorrelation, cross-correlation with CMB lensing, and cross-correlation with galaxy catalogs are provided in Sections 4, 5 and 6, respectively.In Section 7, we summarize and discuss our findings.In Appendix A, we use the web tool Tomographer to further validate our results regarding the radio-catalog cross-correlation.

OBSERVATIONAL DATA
We consider radio maps from LOFAR observations obtained in highband antenna dual mode.Two of the maps were the result of a targeted observation at 142 MHz used for the analysis of Ref. Offringa et al. (2021).The maps are centered at 9ℎ 38 41, 30 • 49 ′ 12 ′′ (referred to as coldest patch) and 10ℎ 25 00, 30 • 00 ′ 00 ′′ (referred to as reference field).Each image has 25002 pixels and covers an area of approximately 48 deg 2 .In these maps, bright sources with a flux larger than 7 times the rms fluctuation (i.e., with a flux ≳ 5 mJy) were selected in a uniformly weighted image that was tapered to make the synthesised beam a Gaussian with a full width at half maximum (FWHM) of 30" resolution.These sources have been removed from the visibilities.The residuals were then imaged with natural weighting, resulting in a synthesised beam that has a FWHM of ∼ 3.5 arcmin.Further details on observation and data reduction can be found in Offringa et al. (2021).We also consider another map, coming from multi-beaming observations from legacy LOFAR data and being a mosaic composed by 7 adjacent fields, as described in Cowie et al. (2023).The mosaic area is a 16 × 16 deg 2 , centered on the North Celestial Pole (NCP).The synthesised beam has a size of approximately 3.3 arcmin.Sources have been subtracted off, as for the other images, again up to a threshold of around 5 mJy.The images are naturally weighted to enhance sensitivity to the angular power spectrum.
As tracers of the large-scale structure, we consider three galaxy catalogs.The Two Micron All Sky Photometric Redshift catalog (2MPZ) (Bilicki et al. 2014) is peaked at very low redshift,   ≃ 0.07, and contains ∼ 10 3 galaxies in the region of a LOFAR field.The Sloan Digital Sky Survey (SDSS) main sample (Data Release 12) (Alam et al. 2015) has a wider distribution peaked at intermediate redshift,   ≃ 0.35, and offers higher statistics since the number of galaxies in the region of a LOFAR field amounts to ∼ 10 5 .The content of the SDSS main sample in the coldest patch field is shown in Figure 1.Finally, we consider the SDSS quasar catalog (Data Release 16) (Lyke et al. 2020), so to cover higher redshifts.To this aim, namely, in order to test the correlation of the radio images with a population of sources at high redshift, we further remove sources with  < 1 from the quasar catalog.In the region of a single LOFAR field, it then contains about 8 × 10 3 objects.The redshift distribution of the sources for the different catalogs is shown on the right panel of Figure 3.To perform the angular cross-correlation measurement between a catalog and a radio image, we will actually consider sources from the catalog that are located within an area four times larger than the image, sharing the same center.
We also consider the publicly available cosmic microwave background (CMB) lensing convergence map  CMB , which is reconstructed from Planck PR4 data 2 (Carron et al. 2022).This lensing convergence map is extracted using an improved lensing quadratic estimator applied to the reprocessed Planck PR4 NPIPE foregroundcleaned CMB temperature and polarization maps.These NPIPE maps incorporate about 8% more data than the previous 2018 Planck PR3 release.For the CMB lensing reconstruction, CMB angular scales over the 100 ≤ ℓ ≤ 2048 range are optimally filtered and fed to the quadratic estimator.We show in the left panel of Fig. 2 the Wienerfiltered Planck CMB lensing convergence map over a patch of the sky centered on the North Celestial Pole.

PHENOMENOLOGICAL MODELS OF RADIO SOURCE POPULATIONS
We model the contribution of radio sources by assuming they populate halos in the Universe.The halo distribution is described adopting the halo mass function d/(dd) from Sheth & Tormen (1999).
Then for each halo, we define an average radio flux density which depends on the halo mass and on the redshift.For simplicity, we take the two dependencies to be separable, i.e., the average radio flux density at 140 MHz of a halo of mass  at redshift  is S(, ) = S 0  () (), with the mass dependence described by a power-law and the redshift behavior described by a Gaussian func-  tion: This completely defines our models.No additional ingredients, besides standard cosmological assumptions, will be introduced in the following.Then, in order to make predictions, we will only have to specify the values of the model parameters ,  0 , and   (while  0 will be set by the requirement of reproducing the RSB excess, as detailed below).
A different approach to predict the angular correlations of a specific population can be adopted if models of its d/d, luminosity function, and halo occupation are available.Our aim is instead to derive more general conclusions, without focusing on a given radio source population, and in this case, there would be three different functions and many more parameters to be arbitrarily defined.The idea behind the function used in Eq. ( 2) is that typically the luminosity grows with the hosting halo mass up to a certain cutoff (in our case large masses will be effectively suppressed by the d/d) and that radio emissions are typically associated to star formation which starts at a certain  max and possibly ends at a certain  min .De-spite the function used being an arbitrary choice, we will see that by varying its parameters, one can obtain very different number counts, luminosity functions, and bias, making it hard to imagine a physical scenario that is not represented, at least in an approximate way, in our sample.
For illustrative purposes, we will show plots for weak, medium and strong scalings with the mass by taking  equal to 0.1 (magenta), 0.5 (black) and 1 (cyan), and for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).We define the maximum redshift of the sources to be  max = 5.The latter value is an arbitrary choice.It can be seen also as a physical assumption, since the most distant tracers of the radio source population, such as quasars, fall off beyond  ∼ 4.However, results are insensitive to the choice.
The above models are chosen in order to cover a variety of physical scenarios.They are not intended to necessarily provide a good fit to anisotropy data (and most of them will not, as we will see later on), but to help understand the key physical requirements for a model to work.The values of the parameters that are preferred by our measurements will be derived in Section 6.  3) and corresponding to models arising from Eq. ( 2) with weak, medium and strong scalings with the mass by taking  equal to 0.1 (magenta), 0.5 (black) and 1 (cyan), and for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).The bold black line at high fluxes shows the number counts measured from detected sources at 150 MHz Intema et al. (2017).Right: Window functions of radio emitters (same color/style coding of the left panel, -axis labels on the left); of galaxy/quasar catalogs and of CMB lensing (black solid curves, -axis labels on the right).
The number counts are obtained through where we use the notation for differential counts customary in the literature, i.e., without reporting dΩ at the denominator.The cosmological volume factor is obtained from a ΛCDM model with parameters from Aghanim et al. (2020).We consider halos between  min = 10 6  ⊙ and  max = 10 16  ⊙ (the precise values of the boundaries of the mass range do not affect our results in a significant way).With the functional form of Eq. ( 2), the last factor is just d/dS = /( S).
The contribution to the radio background from extragalactic sources below the detection threshold is where we set S thr = 5 mJy for the LOFAR images considered in this work (Offringa et al. 2021;Cowie et al. 2023).For each model, i.e., each combination of ,  0 and   , we obtain S 0 in Eq. ( 2) by requiring  bts = 150 K, namely, by normalizing each model to reproduce the RSB excess discussed in the Introduction.
We show the source counts of the nine models mentioned above in Fig. 3 (left panel).They span a variety of scenarios, with excess counts just below the threshold of detected sources (cyan) or peaking at much lower fluxes, with/without (black/magenta) a tail at intermediate fluxes (i.e., in ranges possibly accessible by near future telescopes).
The expected auto-angular power spectrum of the unresolved radio sources and the cross-spectrum with extragalactic tracers of the matter distribution can both be written as where () is the comoving distance to redshift , obeying d/d =  () with  () the Hubble rate,  is the window function described below, and  radio, is the three-dimensional power spectrum.
The subscript  means that Eq. ( 5) is used both for computing the autocorrelation ( =radio) and cross-correlation with galaxy catalogs ( =galaxy) or CMB lensing ( =lensing).For radio sources, the window function can be written as: ) so again using the ingredients described at the beginning of this Section.The upper integration limit is  ′ max = min[ max ,  (S thr )].We show the window functions for our reference models in Fig. 3 (right panel), where one can also note that, since the mass dependence is integrated, see Eq. ( 6), there is very little dependence on .The galaxy window function is described in Section 6, and the lens-ing case in Section 5.The three-dimensional power spectrum  radio, can be computed in the framework of the halo model as the sum of two terms:  =  1ℎ +  2ℎ .More details on the power spectrum for the different cases are provided in the next Sections.

RADIO AUTOCORRELATION
In the case of the radio angular autocorrelation, we do not perform a statistical analysis, but nevertheless, we report an illustrative comparison of the predictions from the models considered in this work with the measurement in Offringa et al. (2021).
Assuming that the radio sources of our interest are point-like, the one-halo term of the power spectrum is given by The assumption of point-like sources makes it -independent.The associated angular power spectrum is therefore flat in harmonic space and sometimes called Poisson-noise term   .Combining the above ingredients one gets Eq. ( 8) provides the dominant angular autocorrelation at very small scales.
The two-halo term is given by  2ℎ (, ) = ⟨ radio ()⟩ 2  lin (, ), where  lin is the linear matter power spectrum.The radio effective bias is computed within the halo model approach, and averaging the halo bias  ℎ [, ] (taken from Sheth & Tormen (1999)) with the radio flux distribution In Fig. 4, we show the autocorrelation of the source models compared with the measurement presented in Offringa et al. (2021).Models with significant emission at low redshift are clearly disfavoured.The flat  case (extending to  max = 5) does not overshoot the data, even though some additional contribution might be needed to explain the rise of the measured power spectrum at ℓ ≳ 1000.Despite no error bars are provided in Fig. 4, we mention here that very high multipoles, ℓ ≳ 3000 (which corresponds to 3.5 ′ ), refer to scales below the radio beam size and thus can suffer from significant uncertainty.Note that since the 2-halo component of the 3D power spectra differ only for the bias term which is of order 1, and the radio window functions are nearly independent from , it is mostly the different redshift dependence that sets the amplitude of the predicted signal at large scales.
Let us also stress again that the quantitative analysis of this work will be focused only on the cross-correlation measurements presented below, which typically suffer from less systematic uncertainties than the autocorrelation.

CROSS CORRELATION BETWEEN RADIO MAPS AND CMB LENSING
Now we move to a quantitative analysis based on cross-correlations, starting from the cross-correlation of the Planck CMB lensing convergence map with the LOFAR 140 MHz NCP mosaic image.We use the NCP mosaic image instead of the individual LOFAR fields in order to maximize the area of overlap, so to enhance the crosscorrelation signal.The path of CMB photons is deflected by the intervening gravitational potentials associated with the large-scale structure (see, e.g., Lewis & Challinor (2006) for a review).This induces a remapping of the primary CMB anisotropies  ( n) according to X ( n) =  ( n + ∇( n)), where X ( n) are the observed CMB maps and ( n) is the CMB lensing potential.These deflections, typically of a few arcminutes, introduce correlations between modes of the CMB anisotropies that can be exploited to reconstruct the projected gravitational potential ( n).As a result, CMB lensing maps contain information on the integrated total matter distribution along the line-of-sight, with peak sensitivity around redshift  ∼ 2 (see Fig. 3).Conversely, radio galaxies are biased signposts of the same dark matter haloes that act as lenses for CMB photons.It is therefore natural to expect a degree of correlation between maps of CMB lensing and radio data since they respond to the underlying dark matter field in complementary ways.In fact, many statistically significant crosscorrelation measurements between CMB lensing maps and galaxy catalogs extracted from radio surveys have been reported to date (e.g., Smith et al. 2007;Allison et al. 2015;Alonso et al. 2021;Piccirilli et al. 2023).In this paper, we instead focus on the crosscorrelation between CMB lensing convergence  = − 1 2 ∇ 2  and the spatial anisotropies of the radio background at 140 MHz.
The CMB lensing window function is given by Bartelmann (2010): where  0 is the Hubble constant today, Ω m is the matter-density parameter, and  * is the comoving distance to the last-scattering surface, all taken from Aghanim et al. (2020).The shape of  lensing is shown in Fig. 3 (right).Note that it extends over a large range of redshifts.The lensing bias is obtained by averaging the halo bias over the matter distribution where ũ( |, ) is the Fourier transform of the matter halo density profile.
The cross-correlation between radio data from LOFAR and Planck CMB lensing is measured in harmonic space.We use the pseudo- ℓ estimator implemented in the NaMaster3 package (Alonso et al. 2019) to extract the angular cross-power spectrum  Δ ℓ .Given two masks to weight the observed maps, NaMaster calculates the modecoupling matrix  ℓℓ ′ and undoes the effects of masking-induced statistical coupling between different harmonic modes.To improve the condition number of the mode-coupling matrix, we bin all spectra into 15 linearly-spaced bandpowers between 50 ≤ ℓ ≤ 2000 (approximately corresponding to Δℓ ≃ 139) and return an unbiased estimate of the cross-power spectrum.We focus on 50 ≤ ℓ ≤ 2000 since the uncertainty becomes large outside this range (compared to model predictions) and there is no statistical gain in considering additional multipoles.When estimating the cross-power spectrum, we mask the lensing convergence map using the fiducial mask provided by the Planck team which removes regions dominated by Galactic foregrounds as well as a number of extragalactic sources (including infrared galaxies and Sunyaev-Zel'dovich clusters).The LOFAR radio data are similarly weighted by a map of their sensitivity over the NCP mosaic image.Due to the small size of the observed overlapping The color/style coding for the theoretical models is as in Fig. 3, namely arising from Eq. ( 2) with weak, medium, and strong scalings with the mass by taking  equal to 0.1 (magenta), 0.5 (black) and 1 (cyan), and for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).
patch (approximately 100 deg 2 ), we work in the flat-sky limit and reproject the original Planck data from the HEALPix (Górski et al. 2005) pixelation scheme to the same LOFAR orthographic projection using the reproject package. 4We additionally correct the crosspower spectrum for the pixel window function, reprojection effects, and a small (∼ 5%) multiplicative CMB lensing transfer function to account for the spatial-dependent normalization of the lensing map (see, e.g., Carron et al. (2022); Alonso et al. ( 2023)).We calculate the covariance matrix of the band powers by cross-correlating 480 realizations of the CMB lensing convergence maps as reconstructed by the Planck pipeline with the real LOFAR radio map.This implicitly assumes the two maps to be uncorrelated, which is a well-founded assumption given the noise levels of the respective fields (especially CMB lensing).We take the square root of the diagonal of the covariance matrix C ℓℓ ′ as the statistical uncertainties on the recovered cross-power spectrum.
We do not detect a statistically significant cross-power spectrum.Defining the null-hypothesis as the absence of correlation between the CMB lensing and the radio background,  Δ ℓ = 0, the chi-square value can be evaluated as responding to a probability-to-exceed of 12%.5 Therefore, we cannot rule out the null hypothesis.A comparison between our reference models and the measurement is provided in Fig. 5 (left).Note that, since the CMB lensing window function is rather flat in redshift, the predictions from the different models are rather similar (remember they have the same overall "normalization" given by matching Eq. ( 4) with  RSB ).This highlights the importance of this cross-correlation which has the potential to constrain different interpretations of the RSB excess with little dependence on the specific model.On the other hand, the data in Fig. 5 are poorly constraining.A measurement with a larger radio field, so to increase the statistics and reduce errors, is in order to fully exploit this signature for constraining the RSB excess.

CROSS CORRELATION BETWEEN RADIO MAPS AND GALAXY CATALOGS
In this Section, we analyze the cross-correlation of the RSB with galaxy catalogs.Contrary to the CMB lensing case, where we used the NCP mosaic image, here we consider the coldest patch and reference individual LOFAR fields.Despite having a smaller area, they are more suitable for the purposes of this Section, i.e., they provide a more significant signal, because of the lower noise of the images and since the NCP location of the mosaic image is not ideal for galaxy surveys.
In the case of galaxies, the window function entering Eq. ( 5) simply reads  galaxies () =  ()   /, where   / is the redshift distribution of galaxies (normalized to unity).The window functions used in this work are derived directly from the redshift of the sources reported in the catalogs.They are shown with black solid lines in Fig. 3 (right panel).
At the scale of interest, and considering sufficiently small systems (i.e., galaxy sizes), the one-halo term of the 3D power spectrum is -independent, leading to an ℓ-independent angular power spectrum.This might be not the case for cluster emission and very low redshifts (as the ones in the 2MPZ catalog).However, an explanation of the RSB excess in terms of clusters is viable only if we assume they are surrounded by a giant halo which would be at least partially resolved out by the LOFAR interferometer (and beyond extended emissions detected so far around clusters).On the contrary, if we assume their size is such that they are fully detected in LOFAR observations, clusters are too rare to offer a significant contribution to the RSB excess.Therefore, with the data considered in this work, we do not evaluate clusters and we can simplify our treatment assuming point-like sources, i.e.,  1ℎ ℓ =const.Since the estimate of  1ℎ ℓ heavily depends on details of the population under investigation, i.e., on the matching between the type of galaxies in the catalog under consideration and the source of radio emission, the prediction of the cross-correlation signal is highly model-dependent, i.e., it would require a precise estimate of the halo occupation distribution of both galaxies and radio emitters.Given the spirit of this work of deriving general conclusions, we decided to leave the effective shot noise level, i.e.,  1ℎ ℓ , as a free parameter, determined by the fit.To better understand this point, consider, e.g., two radio source populations providing the same RSB, but in one case given by a "few" bright sources and in the other by many dim sources.In the first case, if all the sources are hosted by the galaxies of the catalog involved in the cross-correlation, then the one-halo term would be large.On the contrary, if they are hosted by different galaxies (so do not belong to the same halo), there is no one-halo correlation.In the case of the numerous and faint population, they are distributed in so many halos that a single galaxy catalog includes only a fraction of them and in turn, the associated one-halo term would be small.Given the awkward link between the RSB contribution and the correlation, the one-halo term is not the main subject of this work.
The two-halo term can be instead computed in a much more reliable way and it is given by  2ℎ radio,gal (, ) = ⟨ gal ()⟩ ⟨ radio ()⟩ lin (, ), with the effective bias of galaxies taken from the literature: for 2MPZ we follow Balaguera-Antolínez et al. ( 2018), for SDSS main Repp & Szapudi (2020), and for SDSS quasars Laurent et al. (2017).
We express the results of the cross-correlation in terms of the correlation function , instead of the angular power spectrum  ℓ , the two being related by the transformation: where  ℓ are the Legendre polynomials, and Cℓ =  ℓ Wradio ℓ with  ℓ from Eq. ( 5) and Wradio ℓ being the radio beam window function (we neglect the galaxy beam since it is much smaller than the radio one).
We already saw that the dependence on  was mild in the case of autocorrelation.It is even milder for the cross-correlation, where the bias enters just linearly (instead of quadratically as in the autocorrelation).The illustrative plots of this Section will refer only to the case with  = 0.5.
We perform the measurement of the cross-correlation between LO-FAR images at 140 MHz and galaxy catalogs described in Section 2 using the "NKCorrelation" class from the TreeCorr package (Jarvis et al. 2004).This class handles two-point correlation functions between a catalog of source positions and a scalar quantity, in our case the radio flux in each pixel of the image under consideration.The cross-correlation is computed in 20 logarithmically spaced bins, with centers ranging from 0.12 to 420 arcmin.The correlation estimator  () is defined such that ∫   () = 0.6 The covariance matrix is then estimated using the jackknife method with 300 patches.
The measured cross-correlation of the LOFAR images with SDSS main catalog is shown in Fig. 5 (right) and with the 2MPZ and SDSS quasar catalogs in Fig. 6.We show the predicted cross-correlation for the different models with black lines.We also superimpose a Poisson-noise term, i.e., a contribution with  ℓ =const, (red and blue solid lines), arbitrarily normalized, since, as already mentioned, our modeling aims to be very general but cannot precisely predict the one-halo component (i.e., the correlation at angular size ≲ radio beam size).Note that a theoretical term  ℓ = const has to be convolved with the telescope beam to be compared with observations, and therefore when transformed to real space it has exactly the beam shape.
It is clear that low redshift scenarios, such as the case with ( 0 , ) = (0.3, 0.15), are excluded by SDSS main data.In Fig. 6 (left), we add a case with ( 0 , ) = (0.1, 0.05), shown with dasheddotted line, to illustrate that 2MPZ is very effective in excluding very low- models.Therefore, as for the case of autocorrelation, we find that the low- cases are challenged by data.High redshift scenarios are instead compatible with our cross-correlation measurements.
To be more quantitative about this last statement and about whether a scenario involving spatial correlation of the RSB origin population with that of galaxy catalogs is preferred to a scenario without such spatial correlation, we perform a global fit including all four cross-correlation measurements.Six free parameters ( 0 ,   ,  1ℎ 2  ,  1ℎ  ,  1ℎ   1ℎ  ) are present in the fit.We set  to three different cases  = 0.1, 0.5, 1.The  2 is computed using the covariance matrix discussed above. 7In Fig. 7 (left), we show as a red-filled area the 99% C.L. region of ( 0 ,   ) obtained after profiling out the four one-halo parameters  1ℎ  .The region at low ( 0 ,   ), i.e., below the preferred region, overshoots low- measurements and is excluded.The region at large ( 0 ,   ), i.e., above/right of the preferred region, undershoots the measurements.So, a source population peaked at high-z would not explain the SDSS Black lines are the predictions for the models described in the text (including only the clustering, i.e, 2-halo, term), with style coding as in Fig. 3, namely for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).Green points are the measurements at different multipoles with dark green for positive values and light green for negative ones.For illustrative purposes, we show the measurement also in linear scale in the inset plot.Right: Measured cross-correlation angular power of the LOFAR images of the reference field (red) and coldest patch (blue) with SDSS main catalog.Solid blue and red lines show the LOFAR radio beams in the two fields.Black lines are the predictions for the same models as in the left panel.
main measurement, but, at the same time, it is also not constrained by our analysis and could contribute to the RSB.We note here that cases with  ≲ 0.5 provide nearly identical contour regions.For larger values,  ≃ 1, the contribution from massive halos increases and since they are more numerous at low redshift than at high redshift, the bounds become somewhat more constraining (see dotted line with respect to dashed line and red band).The two LOFAR fields, reference (red) and coldest (blue) provide compatible results, with the latter being slightly more constraining, as expected.
To better understand whether or not we need an extragalactic component from large-scale clustering to fit the individual crosscorrelation measurements of galaxies with the LOFAR fields, we compare a model with only a (beam-like) Poisson-noise term with a model where we also add the two-halo component coming from clustering.For the latter term, since the angular dependence of the correlation function is not dramatically different for different models, see Figs. 5 and 6, we consider the case with  = 0.5 and ( 0 , ) = (1.0,0.5), with a parameter providing the overall normalization and determined by the fit.Thus we have one normalization parameter for the beam-like model (which can be seen as parametrizing our ignorance about the shot noise contribution) and two parameters when adding the 2-halo term (the second one essentially parametrizing our ignorance about the clustering bias).The Δ 2 can be thus taken with one degree of freedom.Values are reported in Table 1.
The cross-correlation with the SDSS main catalog shows clear evidence for an extragalactic clustered component.The fact that we obtain closed contours in Fig. 7 (left), is mostly due to the correlation with the SDSS main catalog.The presence of a two-halo term in the cases of 2MPZ and SDSS quasar is more uncertain.shows that the evidence of the signal is at a statistical significance ≳ 7.This implies that at least a fraction of the RSB must be of extragalactic origin in order to provide such a correlation.We estimate the minimum RSB fraction required, by considering a radio population with a window function matching the one of the SDSS main catalog, shown in Fig. 3, imposing its cross-correlation to be compatible at 99% C.L. with our measurement, and exploring different bias obtained from  ≤1.Then we compute  bts from Eq. ( 4) and find values always larger than 20% of the RSB excess.In other words, we conclude that, in order to explain the cross-correlation of the LOFAR images with the SDSS main catalog, at least 20% of the RSB at 140 MHz has to come from an extragalactic and clustered source population.
In contrast, we see that the RSB excess cannot predominantly come from sources at low redshift.On the left panel of Fig. 7 we draw the 8 contour (solid red line) which is indeed an open contour, excluding only models at low redshift but not the high redshift ones.
In Fig. 7 (right), we summarize the fraction of the RSB excess that can come from a certain redshift, by plotting d/d as a function of  for all the models within the 99% C.L. region of the profiled ( 0 , ) distribution (i.e., within the regions like the ones reported in the left panel for  = 0.1, 0.5, 1 in the case of the reference field and for  = 0.5 for the coldest patch).(13) By taking at the upper boundary of the red band, we find that the contribution for  ≲ 0.5 amounts to ≲ 15% of the RSB.One might wonder if the detected correlation can be explained by known source populations, i.e., by their components below the detection threshold.Using the model of counts in Intema et al. (2017), we found that the RSB contribution from below threshold sources amounts to 2 K.The counts essentially refer to AGNs.At 140 MHz the number counts of star-forming galaxies (SFGs) has been not determined.We can however use information from deeper surveys at higher frequencies to estimate it.From the model in Gervasi et al. (2008), we found an RSB contribution from SFGs of about 10 K at 140 MHz.Therefore the overall contribution from the unresolved components of known sources is ∼ 10% of the total RSB, not far from the 20% mentioned above.The redshift behavior plotted in Fig. 7 is compatible with an SFG population.Therefore, considering that the extrapolation of the counts is an oversimplified estimate, it is not inconceivable that sub-threshold sources of known type could account for the measured cross-correlation.

CONCLUSION
We have investigated the angular cross-correlation power spectra of LOFAR images of the diffuse radio emission with galaxy and quasar catalogs from SDSS and with the map of CMB lensing from Planck.We have considered the cross-correlation with source-subtracted radio maps, as opposed to radio sources extracted from maps, with the goal of constraining extragalactic origin scenarios for the RSB.
As discussed in §5, we do not detect any significant correlation between the radio images and CMB lensing, likely in part because of the limited region of effective overlap between the two maps.Still, since the CMB lensing extends to a wide redshift range and can overlap with very different RSB interpretations, we believe that crosscorrelations of CMB lensing with radio maps are a useful avenue to pursue further.
As discussed in §6, we find that models of the RSB resulting from very low redshift sources ( ∼ 0.1-0.3)greatly overproduce the measured cross-correlation power spectrum between the LOFAR radio maps and galaxy catalogs, excluding those redshifts of origin for the predominant source classes of the RSB.However, models of the RSB resulting from intermediate or higher redshift sources ( ≳0.5) are not incompatible with the measured cross-correlation power.This is true even for a wide range of faint radio source count (/) distributions.Generally, as discussed in §6, a scenario involving spatial correlation of the RSB origin population with that of the SDSS main galaxy catalog is preferred to a scenario without such spatial correlation, favoring hypothetical source populations that trace the matter distribution in the universe.In §6, we calculate that at a minimum 20% of the RSB surface brightness level must originate from populations tracing the large-scale distribution of matter in the universe, indicating that at least this fraction of the RSB is of extragalactic origin.
These results are compatible with the analyses in Offringa et al. ( 2021) and Cowie et al. (2023), which both showed an excess of radio anisotropy power above that which would result from known extragalactic radio source classes.The latter work did raise the possibility of significantly extended sources as an important source class which this analysis does not constrain.Here we see that if 100% of the RSB is produced by a new class of numerous, low-flux radio sources, these sources should be at intermediate or higher redshifts.This class of sources has to provide a strong evolution of the far infrared-radio correlation (Ysard & Lagache 2012) so as not to overproduce the measured level of the far infrared background.Whether this can be compatible with constraints derived for known source classes (see, e.g., (Magliocchetti 2022)) depends on the redshift distribution of their radio emission.Looking forward, in order to ascribe the RSB surface brightness excess to a specific particular extragalactic population with angular cross-correlation analyses, it would be crucial to measure a cross-correlation at  ∼ 1 and angular scales much greater than radio beam size. of redshift distributions from measurements of the spatial clustering of arbitrary sources (or diffuse intensity maps), using a set of reference objects for which redshifts are known (see, e.g., Newman (2008); Ménard et al. (2013); Rahman et al. (2015); Chiang et al. (2019) for a thorough discussion of the method's principles and validation).The key idea hinges on the spatial clustering observed on the sky when matter tracers share overlapping redshift distributions (neglecting gravitational lensing effects).In the case of Tomographer, the spectroscopic reference sample is constructed from a compilation of about two million Sloan Digitial Sky Survey sources (including galaxies and quasars) distributed over approximately 10,000 deg 2 and extending out to redshift  ≲ 4.
Given a sky map (or a source catalog), Tomographer returns a measurement of the bias-weighted redshift distribution of the photons contributing to the radio background as well as an estimate of its uncertainties (estimated through a resampling approach).More specifically, Tomographer computes the spatial correlation function () between the pixelized density maps (at an HEALPix resolution of  side = 2048) of the reference spectroscopic sample and of the input diffuse map (or source catalog).To enable a fast calculation of the results, the redshift and angular binnings are pre-set to allow the algorithm to use a number of pre-computed quantities.Schematically, the cross-correlation amplitude w between the input data  and the reference data  calculated in the reference redshift bin   is given by w (  ) =   (  )  (  )  (  ) wDM (  ), where   and   are the effective linear bias parameters for the input/reference data, wDM is the dark matter clustering, and   (  ) is the quantity of interest.Note that the correlation function amplitudes w are obtained by integrating the correlation functions () over the angles that correspond to physical scales 2 ≤  ≤ 8 Mpc and that the scales larger than 2 deg are filtered out to mitigate potential large-scale systematics (including stellar contamination and dust extinction).By looking at Eq. A1, we can see how the ratio between the cross-correlation and reference sample auto-correlation is sensitive to the bias-weighted redshift distribution of the unknown sample, w / w ∝   ×   .In Fig. A1, we show the recovered   × () estimated from the LOFAR coldest patch.As can be seen, the inferred distribution qualitatively agrees with that shown in Fig. 7, exhibiting a bulk contribution from structures at  around 1, once a (moderate) bias evolution with redshift is foreseen.

Figure 1 .
Figure 1.Scatter plot of the SDSS main sample (left) in the field of the LOFAR coldest patch map (right).

Figure 2 .
Figure 2. A 16 • × 16 • patch of sky centered on the North Celestial Pole as seen by the Widefield LOFAR High-Band Antenna at a radio frequency of 140 MHz (right panel), juxtaposed with a Wiener-filtered CMB lensing convergence map derived from the Planck PR4 dataset (left panel, dimensionless).Gridlines of constant declination are shown every 3 degrees.

Figure 3 .
Figure 3. Left: Illustrative examples of radio source count models (  /) considered in this work, computed through Eq. (3) and corresponding to models arising from Eq. (2) with weak, medium and strong scalings with the mass by taking  equal to 0.1 (magenta), 0.5 (black) and 1 (cyan), and for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).The bold black line at high fluxes shows the number counts measured from detected sources at 150 MHzIntema et al. (2017).Right: Window functions of radio emitters (same color/style coding of the left panel, -axis labels on the left); of galaxy/quasar catalogs and of CMB lensing (black solid curves, -axis labels on the right).

Figure 4 .
Figure 4.The measured radio autocorrelation angular power spectrum obtained from the LOFAR images of the reference field (red) and coldest patch (blue) fromOffringa et al. (2021), along with autocorrelation angular power spectra that would result from the theoretical models of radio sources discussed in Sec. 3. The color/style coding for the theoretical models is as in Fig.3, namely arising from Eq. (2) with weak, medium, and strong scalings with the mass by taking  equal to 0.1 (magenta), 0.5 (black) and 1 (cyan), and for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).

Figure 5 .
Figure5.Left: Measured cross-correlation angular power spectrum between the LOFAR NCP mosaic image and CMB lensing measured by Planck.Black lines are the predictions for the models described in the text (including only the clustering, i.e, 2-halo, term), with style coding as in Fig.3, namely for populations with distribution in redshift peaked at low or intermediate , or flat in , by taking ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted).Green points are the measurements at different multipoles with dark green for positive values and light green for negative ones.For illustrative purposes, we show the measurement also in linear scale in the inset plot.Right: Measured cross-correlation angular power of the LOFAR images of the reference field (red) and coldest patch (blue) with SDSS main catalog.Solid blue and red lines show the LOFAR radio beams in the two fields.Black lines are the predictions for the same models as in the left panel.

Figure 6 .
Figure6.Measured cross-correlation angular power of the LOFAR images of the reference field (red) and coldest patch (blue) with the 2MPZ (left) and SDSS quasar (right) catalogs, with error bars.Color/style coding as in Fig.5with ( 0 , ) equal to (0.3, 0.15) (dashed), (1, 0.5) (solid), and (0,+∞) (dotted), with the addition of a model with ( 0 , ) = (0.1, 0.05) (dashed-dotted line) in the left panel to emphasize that this measurement excludes very low- explanations of the RSB excess.Solid blue and red lines show the LOFAR radio beams in the two fields.

Figure 7 .
Figure 7. Left: Closed contours are the 99% C.L. regions of the two parameters ( 0 , ) describing the redshift behaviour of radio sources.The regions are obtained by fitting the cross-correlation measurements shown in Figs. 5 and 6 for the reference LOFAR field (red filled region for  = 0.5 and black line for  = 0.1 (dashed) and  = 1 (dotted)) and coldest patch (blue line).The red solid line shows the 8 (open) contour for  = 0.5 and the reference LOFAR field.Right: Contribution to the extragalactic radio background as a function of the redshift, taking the envelope of all models belonging to the preferred regions (see left panel), with  ∈ [0.1, 1] and for reference LOFAR field (red) and coldest patch (blue).

Figure A1 .
Figure A1.Inferred bias-weighted redshift distribution of the photons contributing to the radio background in the LOFAR coldest patch.The units are the same ones of the original LOFAR map, i.e.Jy/beam.

Table 1 .
Δ 2 values obtained from  2   −  2   2ℎ, where in the model providing  2   only a Poisson-noise term scaling as the radio beam is included in the fit, whilst for  2   2ℎ we add a clustering term, computed from the 2-halo correlation described in the text.