© The Authors 2024

1. Introduction

The excellent multi-frequency coverage of the LOFAR Two-metre Sky Survey (LoTSS) deep fields provides an opportunity to learn about the redshift distribution of low-frequency radio sources (120–168 MHz). In turn the redshift distribution of radio sources is an essential ingredient in the study of the spatial clustering of radio sources and their evolution. The LoTSS deep fields first data release (LoTSS-DF-DR1; Tasse et al. 2021; Sabater et al. 2021) offers information such as the type of sources (Best et al. 2023), cross-matching with multi-frequency observations (Kondapally et al. 2021), and an improved approach on estimating the probability distribution functions (pdfs) of photometric redshifts (photo-z) of radio sources (Duncan et al. 2021), resulting in a high level of photo-z completeness. Photo-z information derived from the full posterior pdf and spectroscopic redshift information if available have been included in released catalogues¹; the full posterior redshift distributions for the individual sources are not publicly available.

Continuum radio surveys enable us to study the angular distribution of the large-scale structure of the Universe, see Longair (1978), Becker et al. (1995), Condon et al. (1998), Miley & De Breuck (2008) for early works. Obtaining redshift information is essential for the study of the corresponding spatial large-scale distribution and to extract cosmological parameters (Camera et al. 2012), especially as radio galaxies exist over a wide range in redshift (Drinkwater & Schmidt 1996; Smolčić et al. 2017; Best et al. 2023). Thus, one needs to consider large radio surveys and complement them with redshift measurements for – if possible – all radio sources. This is a challenging task and one is limited by the number of sources in such surveys which have an estimate of their photometric or spectroscopic redshift. Usually, redshift estimates cannot be obtained for each source in wide area surveys, but are available for smaller deep fields with good multi-frequency or good spectroscopic coverage (see e.g. Smolčić et al. 2017; Kondapally et al. 2021). While flux density-limited continuum radio surveys can provide angular positions for a wide area of sky, the full three dimensional analysis of the cosmic structure requires to measure at least the statistical distribution of radio sources as a function of redshift or distance above a given flux density, expressed by a pdf, p(z), i.e.

$$\begin{array}{c}\hfill \mathrm{d}n=\overline{n}\stackrel{\sim }{p}(r)\mathrm{d}r=\overline{n}p(z)\mathrm{d}z,\end{array}$$(1)

where $\overline{n}$ denotes the average surface density of a flux density limited survey, r denotes the comoving radial distance and z the cosmological redshift. The functions $\stackrel{\sim }{p}(r)$ and p(z) denote the sample pdfs in comoving radial distance and redshift space, respectively. With an estimate of p(z) in hand, we can then infer statistical properties of the three dimensional large scale structure from the projected two dimensional information contained in a wide area continuum radio survey.

Previous studies of the clustering properties and redshift distribution of radio sources were done primarily for radio continuum surveys at frequencies around 1 GHz. The angular clustering properties of the VLA Faint Images of the Radio Sky at Twenty centimeters survey (FIRST; Becker et al. 1995) and the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) have been studied extensively (see e.g. Cress et al. 1996; Blake & Wall 2002; Overzier et al. 2003; Nusser & Tiwari 2015; Chen & Schwarz 2016). However, those studies have been limited by quite restricted knowledge of the redshift distribution of radio selected samples, as e.g. the Combined EIS-NVSS Survey of Radio Sources (CENSORS; Brookes et al. 2008) provided spectroscopic follow up of just 143 NVSS radio sources. Thus, extensive use has been made of optical galaxy redshift surveys, such as the 6 degree Field Galaxy Survey (6dFGS; Jones et al. 2004) in Mauch & Sadler (2007), covering the redshift range 0.003 < z < 0.3. Cross-matching of NVSS, FIRST and other radio surveys with catalogues from the Sloan Digital Sky Survey (SDSS; Eisenstein et al. 2001) produced matches for about a third of all radio objects, still for redshifts below one, see e.g. Kimball & Ivezić (2008), Donoso et al. (2009). A better understanding of the redshift distribution of also high redshift radio sources required deep fields with good multi-wavelength coverage, like the COSMOS field (Scoville et al. 2007; Smolčić et al. 2017). At frequencies well below 1 GHz and similar to the LoTSS frequency range, the TIFR GMRT Sky Survey alternative data release 1 (TGSS; Intema et al. 2017) provided the first opportunity for a wide survey area to estimate the redshift distribution of radio sources via a cross matching with SDSS quasars (Pâris et al. 2018) and to study the angular clustering properties at low radio frequencies (Dolfi et al. 2019). The p(z) obtained in these studies are different from our measurements in the sense that we make use of photometric measurements for about 95% of all observed radio sources.

In this work, we study the distribution of low-frequency radio sources as a function of redshift and flux density by combining observations from three deep survey regions. We infer the redshift distribution of the radio source sample by combining their individual posterior pdfs, which were derived from multi-wavelength observations of the three deep fields in Duncan et al. (2021)². We then make use of this flux density limited redshift distribution and the corresponding angular two-point correlation measurements (Siewert et al. 2020) from the wide field LoTSS data release 1 radio sources in the HETDEX field (LoTSS-DR1; Shimwell et al. 2019) to infer clustering properties such as the correlation length of these radio sources. We make use of the value added source catalogue of LoTSS-DR1, in which artefacts and multiple components of radio sources have been identified (Williams et al. 2019). Note that we use deep and wide field data for which artifacts of bright sources have largely been removed by a major effort of visual inspection by experts, in contrast to the more recent second data release of LoTSS (Shimwell et al. 2022). The LoTSS-DR1 value added source catalogue contains photometric redshifts for about 48% of all sources. Mainly the limited depth of the multi-wavelength data used in the LoTSS analysis in the HETDEX field (Duncan et al. 2019) gives rise to selection effects that prevent us from the direct application of the resulting redshift distribution on all radio sources.

The clustering measurements from LoTSS-DR1 are presented in Siewert et al. (2020). After carefully accounting for survey masks, systematic issues, and artefacts arising from multiple components of radio sources, they achieve clustering measures for radio sources that reasonably align with standard ΛCDM cosmology. A complementary clustering study that is based on the wider LoTSS-DR2 radio source catalogue (Shimwell et al. 2022) is presented in Hale et al. (2024). In Nakoneczny et al. (2024) the cross-correlation of LoTSS-DR2 radio sources with the cosmic microwave background is studied. The analysis of Hale et al. (2024) and Nakoneczny et al. (2024) fixes the flux density threshold to 1.5 mJy and is based on the method to infer the overall redshift distribution of LoTSS deep field sources presented and discussed in detail in this work, but also differs in order to make use of both spectroscopic (for 26% of all deep field sources) and photometric redshifts. In this work, we do not to make use of any spectroscopic information. The selection function for the available spectroscopic data of cross-matched radio sources in the three deep fields is unknown. Ideally they would be drawn from a random sample of radio sources. In order to avoid such biases, we stick to the photometric data, which are sampled homogeneously in each of the three fields. We assume that the remaining systematic uncertainties of the photometry are captured by the differences in multi-wavelength coverage and data quality between the three fields, which we account for by bootstrap sampling as described below.

Throughout, we assume the spatially flat Lambda cold dark matter (LCDM) model to convert the redshift of an object to a spatial distance (and vice versa) and use Ω_M = 0.317, in agreement with the Planck best-fit parameters (Planck Collaboration VI 2020). This work is structured as follows. In the next section we describe the LoTSS-DF-DR1 data and obtain the redshift distribution for flux density limited samples from the measured posteriors of photometric redshifts, as presented in Duncan et al. (2021). We describe our technique of weighted stacking of redshift pdfs of sources from the three aforementioned fields. In Sect. 3 we summarise and extend some of the results on the angular two-point correlation function from Siewert et al. (2020) with the wide field LoTSS-DR1 (Shimwell et al. 2019). We estimate the clustering length in Sect. 4. Finally, in Sect. 5 we present our conclusions.

2. Redshift distribution from LoTSS Deep Fields

We consider the task of obtaining a redshift distribution function from the LoTSS-DF-DR1 (Tasse et al. 2021; Sabater et al. 2021). Located in some of the best-studied northern extra-galactic survey fields – Boötes, European Large Area Infrared Survey field North 1 (ELAIS-N1 or EN1), and the Lockman Hole (LH) – the LoTSS Deep Fields data reach a rms sensitivity of ∼32, 20, 22 μJy/beam at a central frequency of 144 MHz for Boötes and LH, and at 146 MHz for EN1, respectively (Tasse et al. 2021; Sabater et al. 2021).

For the three deep fields multi-wavelength observations are available for different fractions of field of view. They cover the infrared, optical and X-ray and together allow us to identify and match 96% of the radio sources within about 26 square degrees of sky (Kondapally et al. 2021). In all three fields, the multi-wavelength matched aperture photometry used for source identification and photometric redshift analysis spans from the UV to mid-infrared, however the exact set of filters and their associated sensitivity varies from field to field (Kondapally et al. 2021, see also Fig. 1 of Duncan et al. 2021).

The photometric redshifts were obtained using a hybrid method that combines both template fitting and machine-learning estimates to produce a consensus redshift estimates and associated calibrated uncertainties. The full methodology is presented in Duncan et al. (2021), here we briefly summarise the implementation. Three different template based estimations are calculated using the EAZY software (Brammer et al. 2008) with three different template sets chosen to represent a range of different spectral energy distributions expected in the radio population, including both stellar only emission and combined stellar and active galactic nuclei (AGN) emission (Duncan et al. 2018a). The individual template fitting results are separately optimised using zero-point offsets calculated from the spectroscopic redshift sample in each field and the posterior redshift predictions calibrated such that they accurately represent the uncertainties in the estimates. Next, additional machine-learning estimates are produced using the Gaussian process redshift code, GPz (Almosallam et al. 2016a,b), with training and prediction performed separately for each field using the respective photometry and spectroscopic training samples.

Finally, the individual template and machine-learning estimates are then combined following the hierarchical Bayesian combination method presented in Dahlen et al. (2013), incorporating the additional improvements outlined in Duncan et al. (2018a,b). The consensus photometric redshift posteriors for an individual galaxy, p_i(z), are evaluated onto a grid based on the initial redshift steps used for template fitting, spanning from 0 ≤ z ≤ 7. A sample of photometric redshift pdfs for nine randomly selected sources, three from each of the three deep fields, is shown in Fig. 1. As the figure demonstrates, the posterior pdf of many sources has a well defined peak, e.g. the sources indicated by the green full line, the red line with dots, or the purple dashed line, while other posterior pdfs are multi modal, e.g. the sources shown by the blue and orange full lines. For other sources, like the ones indicated by the pink dashed line and the grey line with dots, it is clear that they are at z > 1, with a broad redshift distribution.

Fig. 1. Example photo-z redshift posterior distributions of some sources from each of the three LoTSS Deep Fields (coded by colour and line-style).

Table 1 shows the number of radio sources for various flux density thresholds per deep field and the fraction of sources with photometric and spectroscopic redshifts (not used in this work, but included in Hale et al. 2024; Nakoneczny et al. 2024). Note that these numbers do not include any quality assessments, besides the mere existence of the posterior photo-z estimate. The degree of completeness of the photometric redshifts decreases with increasing flux density from 96% below 0.5 mJy to 91% at 8 mJy. The brighter sources are almost exclusively AGN (Best et al. 2023), a population that extends to high redshifts where multi-wavelength data become incomplete. Introducing a flux density threshold of at least 0.5 mJy, the fraction of spectroscopic redshifts lies between 10% and 41%, depending on the field and its flux density threshold. Note that while in the EN1 field the overall fraction of spectroscopic redshifts is as low as 5% without a flux density threshold (besides the source detection criterion of a signal to noise ratio of 5), for a flux density limit of 8 mJy spectroscopic redshifts are available for 41% of all radio sources. For the Lockman Hole we observe a similar trend, but reach only a completeness of 20% at the highest flux density threshold. In contrast, also less complete at the highest flux densities, in the case of the Boötes field, the fraction of radio sources with spectroscopic redshifts varies just between 21% and 30%. Obviously, spectra have not been sampled in a homogeneous manner over the three deep fields. In contrast to the spectroscopic redshifts, the completeness level of photometric redshifts is not only significantly higher, but also shows less fluctuation between the three fields (the fluctuations are at most 4% at any given flux density threshold and at most 7% between all different flux density cuts in the same deep field).

The photometric redshift estimates for radio sources come with varying uncertainties in measurements (Benítez 2000; Brodwin et al. 2006; Duncan et al. 2021). As outlined above, the estimates are from a probability distribution function. For the purposes of reducing the full redshift posterior into a single photometric redshift for catalogues, Duncan et al. (2021) define the single photometric redshift value, z1_median, as the median of the primary peak in the p_i(z) above the 80% highest probability density credible interval (HPD CI; Wittman et al. 2016, see also Duncan et al. 2019 for a more detailed discussion). In the LoTSS-DF-DR1 release, the ‘best’ redshift, z_best, is then defined as the spectroscopic redshift if available, or z1_median otherwise. One could then compute the probabilistic distribution of all the sources as a function of redshift, p(z), using these ‘best’ estimate values for redshift, see Fig. 2. From this figure, the variation from field to field is evident, as well as likely nonphysical features or biases arising from the sensitivity limitations within the multi-wavelength dataset (e.g. the peak in redshift distribution at z ∼ 1.8).

Fig. 2. Weighted sum, p(z), using z_best value from the LoTSS Deep Fields for a flux threshold of 1 mJy. The numbers in the figures legend indicate the number of sources in each field used in the measurement.

However, Fig. 2 does not take into account the inherent uncertainties in the redshift measurement. In this work we therefore consider the full posterior pdf from photo-z estimates. Since the photometric redshift estimates are posterior pdfs that do incorporate all the available information about the redshift uncertainties, here we use a stacking approach to combine them. The advantage of such an approach is its simplicity, but it certainly makes strong implicit assumptions. The most important one that there is no correlation between the individual posterior pdfs, which is certainly not true as they all depend on the same systematic issues of a given set of multi-wavelength observations (see Malz & Hogg 2022 for a detailed discussion). Using estimates from three different fields with different multi-wavelength observations alleviates this problem, but does not solve it entirely.

In the following, we consider the redshift distribution as a function of the flux density threshold. Here we provide more details on the procedure already presented and applied in Hale et al. (2024) and Nakoneczny et al. (2024). Our procedure ensures that the stacked and weighted pdfs are properly normalised. The posterior pdfs of each source are stored at N_z = 701 redshifts z_i that are equally spaced between z₀ = 0 and z₇₀₁ = 7. They are normalised to unity using a trapezoidal integration rule. Our estimate of the redshift distribution in each field f is then

$$\begin{array}{c}\hfill p_f(z)=\frac{1}{N_f}\sum _{s=1}^{N_f}{p}_s(z).\end{array}$$(2)

Here, p_s(z) is the posterior pdf for source s in field f, N_f denotes the total number of sources with posterior pdf in field f. We obtain posterior pdfs for flux density limited samples for each field and a weighted sum is taken by combining each of these p_f(z), i.e.

$$\begin{array}{c}\hfill p(z)=\sum _{f=1}^3{w}_f{p}_f(z),\sum _{f=1}^3{w}_f=1.\end{array}$$(3)

Here, the weight w_f represents the fraction of sources in a field f.

The errors on the weighted sum p(z) are computed using the standard bootstrap resampling method. We make N_b = 50 random samples for different flux density thresholds for each field; this is done by applying the method of selection with repetition, i.e., a sample of N_f sources is formed from the original catalogue by selecting sources randomly from the original catalogues for each field with repetitions allowed. The weighted sum is computed from each of these, resulting in a set of bootstrap samples, {p_b(z)}, which are then used to compute uncertainties based on the empirical variance,

$$\begin{array}{c}\hfill \mathrm{\Delta }p(z)=\sqrt{\frac{\sum _{b=1}^{N_{\mathrm{b}}}{[p_b(z)-\overline{p}(z)]}^2}{N_{\mathrm{b}}-1}}·\end{array}$$(4)

Here, $\overline{p}(z)$ is the weighted sum obtained from the means of the bootstrap samples of the three fields.

Figure 3 presents the resulting p_f(z) and p(z) with an estimate of uncertainties for flux density thresholds of 1, 2, 4 and 8 mJy, respectively. We find reasonable agreement between all three fields and at all flux density thresholds. Some notable variations are seen at z < 0.5, which are likely due to real differences in the large scale structure at low redshifts (the effect of cosmic variance), but might also be influenced by differences in multi-wavelengths coverage or effects of aperture photometry for nearby sources. Another notable variation is observed in the redshift range between 1 and 2. These differences at truly cosmological distance are most likely due to the different multi-wavelength coverage of the three deep fields. By employing errors based on the bootstrap resampling method, we capture both the systematic differences of photo-z measurements and cosmic variance. A comparison with Fig. 2 shows that the pronounced peak in the redshift range 1–2 turns into an increased uncertainty on the pdfs in that range when the pdfs of individual sources are stacked.

Fig. 3. Redshift distribution for individual LoTSS Deep Fields, pf(z), and their weighted sum, p(z), using photo-z pdfs from LoTSS Deep Fields DR1 for flux density thresholds of 1, 2, 4 and 8 mJy, respectively (top left to bottom right panels). Uncertainties are estimated by bootstrap resampling.

In Fig. 4, we also compare our results to estimates of the redshift distribution for radio galaxies from COSMOS field (Smolčić et al. 2017) – after scaling the flux density S ∝ ν^α with an assumed spectral index of α = −0.8 and applying equivalent flux density thresholds, and the T-RECS simulation (Bonaldi et al. 2023). We find only reasonable agreement of our weighted sum p(z) with the p(z) estimated from these two references. The T-RECS simulation shows an excess of sources at 1 < z < 2 and a deficit at small redshifts compared to our results. The p(z) estimated from the COSMOS field is in agreement with our weighted sum p(z) except at around redshift values of 1. Note that the LoTSS Deep Field sample contains about an order of magnitude more sources than the COSMOS field at corresponding flux density thresholds and should therefore be less affected by cosmic variance.

Fig. 4. Comparison of p(z) from T-RECS medium deep field and COSMOS 3 GHz survey (scaled to 144 MHz) with LoTSS Deep Fields DR1 weighted sum p(z) for a flux density threshold of 1 mJy.

3. Angular two-point correlation from LoTSS-DR1

The angular two-point correlation function, w(θ), quantifies the angular clustering of extra-galactic sources (Peebles 1980). It measures the excess probability of finding a source in the vicinity of another source, separated by an angle θ. In case of Poisson distributed point sources this function would be zero. In this work we parameterise w(θ) by a simple power law,

$$\begin{array}{c}\hfill {\displaystyle w(\theta )=A{\left(\frac{\theta }{1\mathrm{deg}}\right)}^{1-\gamma }.}\end{array}$$(5)

This ansatz is motivated by a corresponding power-law ansatz for the three dimensional spatial correlation function, see the next section for more details.

Siewert et al. (2020) measured the angular two-point correlation function and fitted the amplitude of angular clustering A and the index γ for radio sources above different flux density thresholds LoTSS-DR1 radio sources from the HETDEX spring field (Shimwell et al. 2019). The basis for the angular clustering measurements is the LoTSS-DR1 value added source catalogue (Williams et al. 2019) containing 318 520 radio sources observed over 424 square degrees. The measurements of w(θ) make use of the optimal estimator originally defined by Landy & Szalay (1993). Measurements of w(θ) from Siewert et al. (2020) are reproduced in Fig. 5, where we also add a new measurement at S_min = 8 mJy. The estimated angular correlation, $\widehat{w}(\theta )$, is biased due to the integral constraint which arises due to the finite geometry of the survey. Therefore, we also report the estimated bias w_Ω, that we obtain from an iterative fit, i.e. $\widehat{w}(\theta )=A{(\theta /1\mathrm{deg})}^{1-\gamma }-w_{\mathrm{\Omega }}$ (see Appendix of Siewert et al. 2020 for details). The fit range was chosen in Siewert et al. (2020) as 0.2 deg <  θ < 2 deg, avoiding the effects of non-linear structures at the small scales and systematic flux density uncertainties between different pointing on the larger scales (the typical distance between two pointing is 1.7 deg). We consider the ‘mask 1’ measurement for flux density thresholds of 1, 2 and 4 mJy and the additional flux density threshold of 8 mJy, for which we follow the same analysis pipeline, but use only half the number of bins for the angular separation θ to retain a similar number of source pairs per bin. Fig. 5 presents the results for the 2, 4, 8 mJy flux density thresholds. We summarise the measurements from Siewert et al. (2020) and our new results in Table 2. We also quote the goodness of fit and number of radio sources after masking the survey area and flux density cut. We obtain the best goodness-of-fit for the 8 mJy sample, A = 7.69 ± 0.33 and γ − 1 = 0.89 ± 0.49. However, the rather small sample size leads to rather large uncertainties. The smallest uncertainties with a still acceptable goodness-of-fit are found for the 2 mJy sample, A = 5.11 ± 0.60 and γ − 1 = 0.74 ± 0.16, as discussed in Siewert et al. (2020).

Fig. 5. Absolute value of two-point angular correlation function for different fluxes and their corresponding fit lines (solid lines show the best-fit power-law and dashed lines show the fits after including the effect of the integral constraint) over the fit-range 0.2 deg <  θ < 2 deg with same colour codes. Open symbols for data points indicate negative values of w(θ). Data and fits for 8 mJy are new, for 2 and 4 mJy we take them from Siewert et al. (2020).

These measurements of A (reported at 1 deg) are higher to measurements in the literature: For NVSS, A = (1.45 ± 0.15)×10⁻³, γ − 1 = 1.05 ± 0.10 has been obtained by Blake et al. (2004) for S_1.4 GHz > 10 mJy and A = (1.0 ± 0.2)×10⁻³ for γ − 1 = 0.8 by Overzier et al. (2003) for similar flux density thresholds. In a more recent study Rana & Bagla (2019a,b) have measured A = (8.4 ± 0.5)×10⁻³ and γ − 1 = 0.77 ± 0.15 for TGSS at S_154 MHz > 100 mJy, which is in better agreement with our measurements, especially at our highest flux density threshold of 8 mJy. Note that Siewert et al. (2020) concluded that only LoTSS-DR1 results at and above 2 mJy should be used for cosmological analysis, as sources at lower flux densities still suffer from systematic issues that have not been understood well enough in LoTSS-DR1. Discussions of these potential systematic issues are commented on in Hale et al. (2024). The newly added data point for flux densities above 8 mJy is in agreement with the results at 2 and 4 mJy, however the significantly reduced number of radio sources at flux densities above 8 mJy reduces the statistical significance of the measurement. In Hale et al. (2024) the LoTSS-DR2 angular two-point correlation has been measured for a flux density threshold of 1.5 mJy, where for a fixed value of γ − 1 = 0.8, an amplitude of $A=(2.{88}_{-0.06}^{+0.07})\times {10}^{-3}$ was found.³

4. Clustering scale

Ignoring relativistic effects (Yoo 2010; Challinor & Lewis 2011; Bonvin & Durrer 2011), the relation between the typical comoving length scale of clustering, r₀, and the angular two-point correlation w(θ) is given by Limber’s equation (Limber 1954). One has to assume the statistical isotropy and homogeneity of the large-scale structure and a functional form for the spatial two-point correlation function, ξ(r), where r denotes the comoving distance. Often, a power law is assumed, ξ(r) = (r/r₀)^−γ, with γ > 0, see e.g. Peebles (1980). Increasing the clustering length r₀, implies an increasing correlation of any pair of objects at fixed distance r. This ansatz ignores the evolution of large scale structure. To take the evolution of galaxy clustering into account, a dependence on redshift must be introduced. A simple ansatz is to model this evolution as a power of 1 + z,

$$\begin{array}{c}\hfill \xi (r,z)={\left(\frac{r}{r_0}\right)}^{-\gamma }{(1+z)}^{\gamma -3-ϵ},\end{array}$$(6)

where ϵ parameterises the type of the galaxy clustering model (Groth & Peebles 1977; Overzier et al. 2003). ϵ = 0 describes the ‘stable clustering model’, which assumes that cosmic structures are gravitationally bound at small scales and do not evolve over the observed range in redshift. In this model, galaxy clusters neither expand nor contract with the Universe and have a correlation function which decreases with redshift; ϵ = γ − 3 parameterises the ‘comoving clustering model’, in which the large scale structures expand with the Universe and hence their correlation function remains fixed in comoving coordinates. In that model, cosmic large-scale structures would not (yet) be gravitationally bound. Finally, ϵ = γ − 1 parameterises the ‘linear growth model’, in which the clustering is described as per the linear perturbation theory (before the cosmological constant starts to dominate). However, the evolution of clustering properties is degenerate with the evolution of the galaxy clustering bias (the effect that the radio source density does not necessarily trace mass density).

We first state the result of Limber’s approximation, which holds for small angular scales and a Universe dominated by non-relativistic matter (see Simon 2007 for a detailed discussion of the range of validity of the approximation)

$$\begin{array}{c}\hfill w(\theta )=r_0^{\gamma }\sqrt{\pi }\frac{\mathrm{\Gamma }[(\gamma -1)/2]}{\mathrm{\Gamma }[\gamma /2]}{\theta }^{1-\gamma }{\int }_0^{\infty }d\overline{r}{\stackrel{\sim }{p}}^2(\overline{r}){\overline{r}}^{1-\gamma }{[1+z(\overline{r})]}^{\gamma -3-ϵ},\end{array}$$(7)

which we can write as w(θ) = A(r₀) (θ/1 deg)^1 − γ. Above, Γ[x] denotes the Gamma function. Based on a measurement of A and γ, this relation can be used to compute r₀ when a model for the evolution of galaxy clustering is assumed. In Eq. (7) $\overline{r}$ is the mean comoving radial distance of two sources separated by comoving distance r. The mean comoving distance corresponds to a redshift of $z=z(\overline{r})$ and $\stackrel{\sim }{p}(\overline{r})$ is the window function or pdf in comoving distance space for, in our case, the radio sources in the LoTSS deep fields. From observations one generally measures the window functions as a function of redshift p(z). Equation (7) needs to be modified accordingly. For this we need a relation between the comoving radial distance $\overline{r}$ and redshift z, for which one has to assume a cosmological model.

We consider the case where distances are given by a spatially isotropic, homogeneous metric and assume a flat LCDM model. The radial line-of-sight comoving distance for a flat Universe is then given by,

$$\begin{array}{c}\hfill {\displaystyle r(z)=\frac{c}{H_0}{\int }_0^z\frac{dz\prime }{E(z\prime )},}\end{array}$$(8)

where

$$\begin{array}{c}\hfill E(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{M}}{(1+z)}^3+1-{\mathrm{\Omega }}_{\mathrm{M}}},\end{array}$$(9)

with Ω_M denoting the dimensionless matter density of the present Universe and today’s Hubble rate H₀ = 100 h km/s/Mpc.

From the normalisation condition for the redshift distribution, ∫₀^∞p(z) dz = 1, we find

$$\begin{array}{c}\hfill \stackrel{\sim }{p}(r(z))=\frac{H_{0}E(z)}{c}p(z).\end{array}$$(10)

One can then re-write the integral term in relation (7) in terms of redshift z as,

$$\begin{array}{cc}\hfill I(\gamma ,ϵ,{\mathrm{\Omega }}_{\mathrm{M}})& ={\left(\frac{H_0}{c}\right)}^{-\gamma }{\int }_0^{\infty }\mathrm{d}\overline{r}{\stackrel{\sim }{p}}^2(\overline{r}){\overline{r}}^{1-\gamma }{(1+z)}^{\gamma -3-ϵ}\hfill \\ \hfill & {\displaystyle ={\int }_0^{\infty }\mathrm{d}zE(z)p^2(z){(1+z)}^{\gamma -3-ϵ}{\left({\int }_0^z\frac{dz\prime }{E(z\prime )}\right)}^{1-\gamma }.}\hfill \end{array}$$(11)

From this expression, given the values for the density parameters, an evolution model of galaxy clustering, and an observed redshift distribution of galaxies, one can easily compute the clustering length r₀. Its unit follows from the unit of the Hubble distance c/H₀, which is 3000 h⁻¹ Mpc. We obtain

$$\begin{array}{c}\hfill A=\sqrt{\pi }\frac{\mathrm{\Gamma }[(\gamma -1)/2]}{\mathrm{\Gamma }[\gamma /2]}I(\gamma ,ϵ,{\mathrm{\Omega }}_{\mathrm{M}}){\left(\frac{r_0{H}_0}{c}\right)}^{\gamma }{\left(\frac{\pi }{180}\right)}^{1-\gamma },\end{array}$$(12)

when we measure the angular separation in units of degrees. Thus the measured strength of the angular correlation depends on the cosmological model (H₀, Ω_M), the correlations length r₀, the exponent γ, as well as the evolution of clustering, described by ϵ of the given probe. Obviously, the cosmological model as well as the clustering properties of dark and baryonic matter should not depend on the flux density of the chosen sample. However, the specifics of the chosen probe might depend on the flux density cut and differ in clustering evolution.

For different flux thresholds and value of ϵ we compute the clustering length. The expectation and uncertainty of r₀ are calculated from the corresponding redshift distribution and the values of A and γ from the wide field LoTSS-DR1 as mentioned in Table 2. We construct a multivariate normal distribution using these values and the associated covariance matrix, to randomly select (A,γ) pairs from such a distribution. We then randomly select 100 distributions p_b(z), as inferred from LoTSS-DF-DR1 (see Sect. 2), and draw 1000 (A,γ) pairs for each of those p_b(z). Thus we sample 100 000 r₀ values for each case. Then, similar to the bootstrap error computation as described in Eq. (4) we compute expectation (mean) and uncertainty of r₀. The results are shown in Fig. 6 and in Table 2. We also show the median redshift of the four flux density limited redshift distributions from LoTSS-DF-DR1, which are between 0.9 and 1.

Fig. 6. Comoving correlation length for different flux density thresholds at the central LoTSS frequency of 144 MHz and clustering parameter ϵ. For comparison a selection of measurements from the literature are shown. The result of Peacock & Nicholson (1991) uses bright radio galaxies at z < 0.1 and does therefore not depend on the clustering evolution (ignoring any factors of 1 + z, they implicitly use the comoving clustering model), thus applies to any value of ϵ. The analysis of Mazumder et al. (2022) uses observations at 325 MHz and distinguishes AGNs and SFGs. We show their results for AGNs. The radio sources from Lindsay et al. (2014) (FIRST at 1.4 GHz) and Hale et al. (2024) (LoTSS-DR2 at 144 MHz) are mostly AGNs.

As shown in Table 2, the most precise measurement of r₀ is obtained for 2 mJy with (10.1 ± 2.6) h⁻¹ Mpc, (13.3 ± 2.1) h⁻¹ Mpc, and (15.1 ± 2.5) h⁻¹ Mpc for ϵ = γ − 3, 0, and γ − 1, respectively. The clustering length increases with the value of ϵ. We also find good agreement between the different flux thresholds for all three models of clustering evolution. For the comoving clustering model (ϵ = γ − 3) we find (9.5 ± 2.8) h⁻¹ Mpc and (13.3 ± 5.2) h⁻¹ Mpc, for 4 mJy and 8 mJy, respectively. In principle an inconsistent clustering model could be identified by wildly inconsistent clustering lengths for different flux density cuts, but that is apparently not the case here. However, we should not expect them to be identical, as the different flux density samples contain a different mix of AGNs and star forming galaxies (SFGs) (see Best et al. 2023), with almost all radio sources above 8 mJy being AGNs and an increasing (but small) fraction of SFGs as we lower the flux threshold. This is in line with the slight decrease of medium redshift with decreasing flux density threshold, as the first SFGs that are detected tend to be at smaller redshifts.

Several studies have measured the comoving clustering length for galaxies, including radio galaxies (see, e.g. Magliocchetti 2022). The vast literature on that subject suggests that r₀ varies depending on galaxy properties and environment. Especially luminous and old galaxies give rise to larger clustering length, in contrast to less luminous or blue galaxies. For nearby radio galaxies (z < 0.1) Peacock & Nicholson (1991) reported a clustering length of r₀ = (11.0 ± 1.2) h⁻¹ Mpc (as this measurement summarises the amount of clustering today, it can serve as a reference point for all three evolution models). Later Overzier et al. (2003) reported a value of r₀ = (14 ± 3) h⁻¹ Mpc for an analysis of Faranoff-Riley type II (FRII) radio galaxies, which contribute dominantly to the correlation length at high flux densities (S_1.4 GHz > 200 mJy) in NVSS and FIRST. At lower flux densities and for the mix of all radio sources, they find smaller clustering length of r₀ = (4 − 10) h⁻¹ Mpc, depending on the flux density threshold. Based on a simulated redshift distribution of FIRST sources at S_1.4 GHz > 1 mJy, Lindsay et al. (2014) find $r_0=8.{20}_{-0.42}^{+0.41}{h}^{-1}$ Mpc (scaling with a spectral index of α = −0.8 or −0.7 to the central LoTSS frequency we should compare to S > 6.2 mJy or 4.9 mJy, respectively). These numbers correspond to the analysis for the comoving clustering model (ϵ = γ − 3). Those findings have been confirmed by more recent studies (see, e.g. Magliocchetti 2022).

Inspecting our results for the comoving clustering model (see Fig. 6 and Table 2), we find that they are in good agreement among each other for all considered flux density thresholds and with the values that have been measured for radio galaxies at 1.4 GHz (see above). Our results are consistent with the trend observed at 1.4 GHz that increasing flux density thresholds seem to go along with stronger clustering. Combining our findings with results based on TGSS (Dolfi et al. 2019; Rana & Bagla 2019a,b), we see the same trend. In Hale et al. (2024) the analysis of LoTSS-DR2 reveals for γ − 1 = 0.8 a value of $r_0=(7.{32}_{-0.51}^{+0.59})h^{-1}$ Mpc at 1.5 mJy (note the different fit range: 0.03deg < θ < 1deg), also consistent with this trend. At 325 MHz, Mazumder et al. (2022) measured the clustering length from a deep observation of the LH for AGNs (with z_median = 1.02, very similar to our analysis) and SFGs separately at a flux density above 0.3 mJy, corresponding to about 0.6 mJy at the LoTSS frequencies. They find $r_0^{\mathrm{AGN}}=8.{30}_{-0.91}^{+0.96}{h}^{-1}$ Mpc, when assuming γ − 1 = 0.8. Our result is close to their finding.

Let us also investigate the stable clustering (ϵ = 0) and linear growth model (ϵ = γ − 1), with results for both of them presented in Fig. 6 and Table 2. We measure values for r₀ that are close to those measured for galaxy clusters (e.g. (24 ± 9) h⁻¹ Mpc from Bahcall & Soneira 1983; Postman et al. 1992). They show the same trends of more clustering for brighter objects, but provide clustering lengths that exceed the local reference measurement from Peacock & Nicholson (1991), which makes them less plausible, as our values probe the strength of clustering at a median redshift close to unity, and thus both models would suggest that the clustering of radio galaxies actually decreased since redshift of unity. While in principle this is expected for a LCDM model in the future (all not gravitationally bound structures will be diluted in the de Sitter future of the universe), the onset of acceleration is not far enough in the past to make those models plausible. It is interesting to note that the for models with larger ϵ the trend of more clustering for higher flux density threshold is more pronounced than for the comoving model.

Looking at our findings and measurements from the literature, the comoving clustering model can easily match the data for the local and radio loud AGNs, which implies that in fact there is no redshift dependence in ξ(r). This would mean that the galaxy clustering bias must be a function of redshift inversely proportional to the growth of the large scale structure. Indeed Hale et al. (2024) and Nakoneczny et al. (2024) find that such a bias model provides a better description of the LoTSS-DR2 data, compared to a redshift independent galaxy clustering bias.

Finally, we investigate the effect that using the improved posterior redshift distribution in Eq. (2), and shown in Fig. 3, has over obtaining the redshift distribution based on catalogued z_best values, for which an example is shown in Fig. 2. For the flux density threshold with the best statistics but still above the systematic limitations of LoTSS-DR1, namely S > 2 mJy, we find r₀ = (9.9 ± 2.4) h⁻¹ Mpc, (12.9 ± 1.9) h⁻¹ Mpc, and (14.3 ± 2.3) h⁻¹ Mpc for ϵ = γ − 3, 0, and γ − 1, respectively. Those are less than 1σ smaller than the correlation lengths measured by means of the full posterior distribution p(z) and shown in Fig. 6 and Table 2. For these estimates only the uncertainties in (A, γ) are taken into account. The correlation lengths based on z_best values tend to underestimates the high-z tail of the distribution, as AGNs are sometimes mistaken for SFGs at lower redshift, but the opposite happens less likely, as was shown by Duncan et al. (2021). Our method accounts for those systematic and thus results in a slightly larger value of the clustering length as moving objects at fixed angular distance to larger redshift increases their physical distance. Our findings seem to indicate that the measurement of the correlation length is not strongly depending on the assumptions made here and for a rather limited survey area of LoTSS-DR1 with its large uncertainties on A and γ. However, already with the significantly reduced statistical uncertainties of LoTSS-DR2, see Hale et al. (2024) who report uncertainties of Δr₀ ≈ 0.6 h⁻¹ Mpc, this is no longer the case.

5. Conclusions

In this work we obtained estimates of the redshift distribution p(z) of LoTSS Deep Field radio sources for various flux density limits and quantified their uncertainties Δp(z) by means of bootstrap resampling. We based our estimates on stacking the posterior pdfs for individual radio sources as described and determined in Duncan et al. (2021), which had made use of the good multi-wavelength coverage (Kondapally et al. 2021) of radio sources from three LoTSS Deep Fields (Tasse et al. 2021; Sabater et al. 2021). These had allowed Duncan et al. (2021) to obtain posterior pdfs for the photometric redshift of 96% of all radio sources in the survey region with good multi-wavelength information. After applying a flux density threshold, the photometric redshift completeness drops to 91%. We have implicitly assumed that the remaining 4–9% of radio sources follow the same distribution, which might result in an underestimation of the number of radio sources at redshifts above unity. By averaging over three different deep fields and over a total area of about 26 square degrees, we reduce cosmic variance and estimate the effects of systematic issues by bootstrap resampling of the data. We conclude that LoTSS radio sources above a flux density of 1 mJy have a median redshift of about unity.

We also used the inferred redshift distribution of the LoTSS deep fields and applied it on the wide field clustering data from LoTSS-DR1 for three different flux density thresholds, 2, 4 and 8 mJy (the results for 1 mJy are shown for completeness, but should not be trusted as they still suffer from systematic issues; see also the discussions in Siewert et al. 2020 and Hale et al. 2024). We find good consistency for the comoving clustering model in which the clustering structures have formed well before they are observed and probed by the survey, which indicates that halos hosting radio galaxies (which are rare compared to normal optical or infrared galaxies) are in most cases not gravitationally bound to each other and therefore expand with the Hubble flow. Our most precise result is that for 2 mJy, with a clustering length of r₀ = (10.1 ± 2.6) h⁻¹ Mpc. Our bootstrap analysis shows that the precision in this analysis is largely limited by the precision of the measured angular two-point correlation function rather than the redshift distribution of sources. This will change with larger LoTSS samples and significantly improved measurements for A and γ, as already clear from the LoTSS-DR2 analysis in Hale et al. (2024).

Our study of the redshift distribution of radio sources also allows us to estimate the size of the sampled comoving volume of the wide area LoTSS. As an estimate we weight the comoving volume of the Universe by the redshift distribution of LoTSS radio sources,

$$\begin{array}{c}\hfill V_{\mathrm{LoTSS}}=\int \mathrm{d}\mathrm{\Omega }\int \mathrm{d}z\frac{p(z)r^2(z)}{H(z)}·\end{array}$$(13)

Thus, LoTSS-DR2 (Shimwell et al. 2022), which covers about 1/8 of the full sky, probes ≈3.3 h⁻³ Gpc³. After LOFAR observing cycle 20 (which finishes in summer 2024) the coverage of LoTSS will allow us to increase this volume to about ≈10 h⁻³ Gpc³. This demonstrates the unique potential to measure the largest cosmic structure of combining wide area radio continuum surveys with multi-wavelength information. The WEAVE-LOFAR survey (Smith et al. 2016) aims at obtaining spectroscopic redshifts for all LoTSS sources above a flux density of 8 mJy. The clustering study of this work thus serves as a first reference point for much more detailed studies of the three dimensional large scale clustering based on radio selected spectroscopic redshifts, which, compared to optically and infrared selected samples, will provide an independent and complementary probe of the Universe at the largest scales.

Acknowledgments

We thank Maciej Bilicki for discussions and comments. NB and DJS acknowledge financial support by Deutsche Forschungsgemeinschaft (DFG) under grant RTG-1620 ‘Models of Gravity’. CLH acknowledges support from the Leverhulme Trust through an Early Career Research Fellowship and from the Hintze Family Charitable Foundation through the Oxford Hintze Centre for Astrophysical Surveys. 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 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster). 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. 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. SC acknowledges support from the Italian Ministry of University and Research (MUR) through PRIN 2022 ‘EXSKALIBUR – Euclid-Cross-SKA: Likelihood Inference Building for Universe’s Research’ and from the European Union – Next Generation EU. LOFAR data products were provided by the LOFAR Surveys Key Science project (LSKSP; https://lofar-surveys.org/) and were derived from observations with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) 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 parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The efforts of the LSKSP have benefited from funding from the European Research Council, NOVA, NWO, CNRS-INSU, the SURF Co-operative, the UK Science and Technology Funding Council and the Jülich Supercomputing Centre. This research made use of Astropy, a community developed core Python package for astronomy (Astropy Collaboration 2013) hosted at http://www.astropy.org/, matplotlib (Hunter 2007), NumPy (Walt et al. 2011), lmfit (Newville et al. 2016), TopCat (Taylor 2005), SciPy (Virtanen et al. 2020), h5py (Collete & Contributors 2014), TreeCorr (Jarvis et al. 2004) and Python language (van Rossum 1995).

References

All Tables

All Figures

	Fig. 1. Example photo-z redshift posterior distributions of some sources from each of the three LoTSS Deep Fields (coded by colour and line-style).
In the text

	Fig. 2. Weighted sum, p(z), using z_best value from the LoTSS Deep Fields for a flux threshold of 1 mJy. The numbers in the figures legend indicate the number of sources in each field used in the measurement.
In the text

	Fig. 3. Redshift distribution for individual LoTSS Deep Fields, pf(z), and their weighted sum, p(z), using photo-z pdfs from LoTSS Deep Fields DR1 for flux density thresholds of 1, 2, 4 and 8 mJy, respectively (top left to bottom right panels). Uncertainties are estimated by bootstrap resampling.
In the text

	Fig. 4. Comparison of p(z) from T-RECS medium deep field and COSMOS 3 GHz survey (scaled to 144 MHz) with LoTSS Deep Fields DR1 weighted sum p(z) for a flux density threshold of 1 mJy.
In the text

Fig. 5. Absolute value of two-point angular correlation function for different fluxes and their corresponding fit lines (solid lines show the best-fit power-law and dashed lines show the fits after including the effect of the integral constraint) over the fit-range 0.2 deg <  θ < 2 deg with same colour codes. Open symbols for data points indicate negative values of w(θ). Data and fits for 8 mJy are new, for 2 and 4 mJy we take them from Siewert et al. (2020).

In the text

Fig. 6. Comoving correlation length for different flux density thresholds at the central LoTSS frequency of 144 MHz and clustering parameter ϵ. For comparison a selection of measurements from the literature are shown. The result of Peacock & Nicholson (1991) uses bright radio galaxies at z < 0.1 and does therefore not depend on the clustering evolution (ignoring any factors of 1 + z, they implicitly use the comoving clustering model), thus applies to any value of ϵ. The analysis of Mazumder et al. (2022) uses observations at 325 MHz and distinguishes AGNs and SFGs. We show their results for AGNs. The radio sources from Lindsay et al. (2014) (FIRST at 1.4 GHz) and Hale et al. (2024) (LoTSS-DR2 at 144 MHz) are mostly AGNs.

In the text