© The Authors 2024

1. Introduction

Supermassive Kerr black holes are key building blocks of the Universe, on galactic and cosmological scales alike. During episodes of baryon accretion, they turn into active galactic nuclei (AGN), launching winds and jets that warm and rarefy the interstellar medium (e.g., Fabian 2012; King & Pounds 2015). These energy flows generally suppress the formation of new stars, although local star formation enhancement can occur within expanding kiloparsec-radius rings (e.g., Dugan et al. 2017). Meanwhile, jet-mediated AGN outflows – also known as radio galaxies (RGs) – can have a vast, megaparsec-scale reach, protruding from both the galaxy and its circumgalactic medium. The synchrotron radiation from RGs illuminates their jets, lobes, cocoons, and threads (e.g., Ramatsoku et al. 2020), and dominates the known extragalactic radio sky. Because the behaviour of RGs is linked to the physics of black hole accretion and galactic winds, the pressure field and magnetisation history (e.g., Vazza et al. 2017) of the warm–hot intergalactic medium (warm–hot IGM, or WHIM), as well as shocks (e.g., Nolting et al. 2019a), vorticity (e.g., Nolting et al. 2019b), and cooling flows (e.g., Fabian et al. 1984) in the intracluster medium, a precise understanding of the RG phenomenon is indispensable to modern astrophysics.

A key goal of the study of RGs is to identify the main factors that determine their dynamics and to formulate models (e.g., Scheuer 1974; Turner & Shabala 2015; Hardcastle 2018) that describe them quantitatively. One way to investigate the growth of RGs is to search for particularly large examples and to analyse what internal or external traits set them apart from the rest. Following this logic, Oei et al. (2022) presented and studied Alcyoneus, a giant radio galaxy (GRG, giant RG, or simply giant) whose proper length component in the plane of the sky l_p = 4.99 ± 0.04 Mpc. In general, giants are members of the RG population for which l_p ≥ l_p, GRG – where the latter is some fixed megaparsec-scale threshold – and that therefore rank among the largest RGs in existence. However, despite being one of the largest RGs known, Alcyoneus is not particularly luminous, and its host galaxy does not feature a particularly massive central black hole or stellar population – at least, when compared to other giants and their hosts (Oei et al. 2022). The question naturally arises whether external properties, rather than those internal to the host galaxy, are the most important drivers of RG growth.

It is well established that the IGM resists the growth of RGs by forcing jets to convert a part of their kinetic energy into work spent to form lobe cavities (e.g., Hardcastle & Croston 2020). The work needed to free up a cavity is the product of its volume and the local pressure¹. Following this line of reasoning, an RG with fixed intrinsic properties should reach a larger end-of-life extent in more tenuous (and colder) Cosmic Web (CW) environments. It has therefore been proposed that the astonishing growth of giants might be explained by their presumptive tendency to reside in tenuous parts of the Cosmic Web.

Several previous works have investigated the role of the enveloping Cosmic Web density field on GRG growth. In all cases, the authors traced Cosmic Web environments through three-dimensional galaxy positions estimated with photometric or spectroscopic redshifts. The pioneering work of Subrahmanyan et al. (2008) presented a case study of GRG MSH 05-22 and its Cosmic Web environment as traced by 6dF data (Jones et al. 2004). Using a sample of 12 giants and environments traced by 2dF/AAOmega data (Sharp et al. 2006), Malarecki et al. (2015) concluded that the lobes of giants grow in directions that avoid denser regions of the Cosmic Web. In the most comprehensive study yet, Lan & Prochaska (2021) used a sample of 110 giants and environments traced by DESI Legacy Imaging Surveys Data Release (DR) 9 data (Dey et al. 2019). They did not find evidence that giants occur in more dilute environments than non-giant RGs. A limitation of this latest study is the use of photometric redshifts to probe the Cosmic Web density field.

In order to determine decisively how GRG growth and the Cosmic Web relate, analyses with both more giants and more accurate Cosmic Web reconstructions seem necessary. In this work, we present major developments in both regards. First, the Low-Frequency Array (LOFAR; van Haarlem et al. 2013) Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) has made possible the discovery of thousands of previously unknown giants in its Northern Sky imagery at the observing frequency ν_obs = 144 MHz and at resolutions θ_FWHM ∈ {6″, 20″, 60″, 90″}. In particular, the joint search efforts of Dabhade et al. (2020) in the LoTSS DR1 (Shimwell et al. 2019), those of Simonte et al. (2022) in the LoTSS Boötes Deep Field, and those of Oei et al. (2023a) in the LoTSS DR2 (Shimwell et al. 2022) have tripled the total number of known giants, which now stands at ∼3 × 10³. The number of known giants in the part of the Local Universe covered by the Sloan Digital Sky Survey (SDSS; York et al. 2000) DR7 (Abazajian et al. 2009) has even quintupled – a fact whose relevance becomes clear in Sect. 2.3. Second, the last two decades have seen the development of principled, physics-based Bayesian inference techniques through which the three-dimensional total (i.e., baryonic plus dark) matter density field of the Cosmic Web can be reconstructed (e.g., Wandelt et al. 2004; Kitaura & Enßlin 2008; Jasche et al. 2010; Jasche & Kitaura 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2019). These techniques make use of the fact that the statistical behaviour of the Early Universe’s total matter density field is known, as is the dominant process by which this field evolved over cosmic time: gravity. By simulating gravity acting on Early Universe density fields and comparing the evolved fields to the observed spatial distribution of galaxies, the late-time density field can be inferred. These late-time density fields subsequently enable megaparsec-scale density measurements with uncertainties around individual (giant) RGs. We shall use these measurements as an improved probe of the Cosmic Web density field.

Section 2 presents the data: radio galaxy observables and late-time density field reconstructions via which we probed the influence of the Cosmic Web on the growth of giants. In Sect. 3, we explain how we combined these data to determine megaparsec-scale densities and large-scale structure type probability distributions for giants and general RGs. Section 4 presents these results, alongside a quantification of the relation between RG radio luminosity and Cosmic Web density. We finally evaluate evidence for the claim that the Cosmic Web affects RG growth. In Sect. 5, we explain how we tested the reliability of RG Cosmic Web density measurements and dynamical classifications and discuss caveats of and promising future extensions to the current work, before we present conclusions in Sect. 6.

For consistency with Oei et al. (2022, 2023a), and Mostert et al. (2024), we assumed a flat, inflationary ΛCDM model with parameters from Planck Collaboration VI (2020): h = 0.6766, Ω_BM, 0 = 0.0490, Ω_M, 0 = 0.3111, and Ω_Λ, 0 = 0.6889, where Ω_DM, 0 ≔ Ω_M, 0 − Ω_BM, 0 = 0.2621 and H₀ ≔ h × 100 km s⁻¹ Mpc⁻¹. With ‘Local Universe’, we refer to the spherical region of space observed to have redshift z < z_max ≔ 0.16. All reported redshifts are heliocentric. Terminology-wise, we strictly distinguish an RG (a radio-bright structure of plasma and magnetic fields, consisting of a core, jets, hotspots, lobes, a cocoon, and collimated synchrotron threads) from the host galaxy that has generated it. As in our previous works, we define giants to be RGs with projected proper lengths l_p ≥ l_p, GRG ≔ 0.7 Mpc². We define the spectral index α so that it relates to flux density F_ν at frequency ν as F_ν ∝ ν^α. Under this convention, and when synchrotron self-absorption is negligible, radio spectral indices are typically negative.

2. Data

To compare the Cosmic Web environments of Local Universe giants to those of Local Universe radio galaxies, we combined a GRG catalogue, an RG catalogue, and a Cosmic Web density field reconstruction. All three data sets are publicly available.

2.1. Giant radio galaxies

As our source of giants we used the catalogue aggregated by Oei et al. (2023a), which contains 3341 giants with l_p ≥ l_p, GRG ≔ 0.7 Mpc³. This catalogue is intended to be complete up to (and including) September 2022. Thanks to the steradian-scale Northern Sky coverage of the LoTSS, combined with its arcsecond-scale resolution and sensitivity up to degree scales, LoTSS-discovered giants (Dabhade et al. 2020; Simonte et al. 2022; Oei et al. 2023a) dominate the catalogue.

In their LoTSS DR2 manual search for giants, Oei et al. (2023a) used an angular length threshold of $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })=5^{\prime }$ to limit the duration of their search to a manageable few hundred hours. At the same time, this 5′-threshold ensured that the search would yield most giants with sufficient surface brightness in the Local Universe⁴. This has been by design: Oei et al. (2023a) aimed to build a surface brightness–limited, but otherwise complete census of Local Universe giants in the LoTSS DR2 footprint, with the intent of localising them within the Cosmic Web density field reconstructions of Sect. 2.3.

As elaborated upon in Sect. 3.1, we retained 281 giants in the part of the Local Universe where Cosmic Web analysis is possible. Of these, 260 have spectroscopic redshifts; only these giants could be reliably localised. We show their sky locations in Fig. 1. LoTSS DR2 discoveries make up 208 of the final 260 giants (80%). Through six example LoTSS DR2 giants, Fig. 2 provides the reader a sense of the quality of the radio imagery underpinning this work, from which angular lengths ϕ have been inferred, a sense of the reliability of our host galaxy identification, from which spectroscopic redshifts z_s have inferred, and a sense of the morphological and surface brightness diversity of the objects under consideration in this study.

Fig. 1. Mollweide view of the sky showing the locations of all giants and LoTSS DR1 RGs in the Local Universe for which we inferred Cosmic Web densities and dynamical states. The background shows the Milky Way at 150 MHz (Zheng et al. 2017), on which we overlaid the LoTSS DR1 footprint (hatched dark rectangle). All RGs are drawn as lemniscates of Bernoulli; we did not attempt to portray realistic shapes or position angles. The colours represent redshifts z ∈ (0, zmax ≔ 0.16), whilst the diameters are proportional to projected proper lengths lp ∈ ( ∼ 1 kpc, 4.6 Mpc). Giants are translucent. Upon zooming in, the reader can appreciate the wide variety of observed radio galaxy sizes.

2.2. General radio galaxies

In order to determine whether the Cosmic Web environments of giants are exceptional, we must create a sample of reference Cosmic Web environments. For this reason, we also localised within the Cosmic Web a sample of general RGs, selected without regard for their length.

In particular, our starting point for this sample is the radio-bright active galactic nucleus (RLAGN) sub-sample described in Hardcastle et al. (2019), which contains 23 344 of the 318 520 sources (7%) in the LoTSS DR1 value-added catalogue of Williams et al. (2019). To construct their sub-sample, Hardcastle et al. (2019) combined multiple criteria that separate RLAGN from star-forming galaxies (SFGs). At the low redshifts considered in this work (z < z_max ≔ 0.16), most RLAGN are identified through SDSS spectroscopy, and we therefore expect high RLAGN completeness and minimal contamination from SFGs.

As elaborated upon in Sect. 3.1, we retained 1870 LoTSS DR1 RGs in the part of the Local Universe where Cosmic Web analysis is possible. Of these, 1443 have spectroscopic redshifts; only these RGs could be reliably localised. We show also their sky locations in Fig. 1. The LoTSS DR1 image quality is very similar to the LoTSS DR2 image quality, which Fig. 2 illustrates.

Fig. 2. LoTSS DR2 cutouts at central observing frequency νobs = 144 MHz and resolution θFWHM = 6″, centred around giant-generating BCGs of Local Universe clusters. All giants shown are discoveries of Oei et al. (2023a) for which no previous images have been published. Each cutout covers a solid angle of 15′×15′. Contours signify 3, 5, and 10 sigma-clipped standard deviations above the sigma-clipped median. For scale, we show the stellar Milky Way disk (with a diameter of 50 kpc) generated using the Ringermacher & Mead (2009) formula, alongside a 3 times inflated version. Each DESI Legacy Imaging Surveys DR9 (g, r, z) inset shows the central 3′×3′ region.

2.3. Cosmic Web late-time density field

To localise RGs within the Cosmic Web, we used data products from the Bayesian Origin Reconstruction from Galaxies (BORG; Jasche & Wandelt 2013) SDSS run (Jasche et al. 2015). The BORG SDSS uses second-order Lagrangian perturbation theory (2LPT; Bouchet et al. 1995) to forward model structure formation and evaluates the plausibility of a proposed structure formation history by comparing its late-time density field to the three-dimensional positions of galaxies in the SDSS DR7 Main Galaxy Sample (MGS; Abazajian et al. 2009). As a result, the BORG SDSS provides a late-time total⁵ matter density field posterior for the part of the Local Universe covered by the SDSS DR7 footprint. The posterior is represented by a Hamiltonian Monte Carlo (HMC; Duane et al. 1987) Markov chain of approximately ten thousand samples after the burn-in phase. Each sample covers the same volume of (750 Mpc h⁻¹)³ extent with a cubical grid of 256³ voxels. Thus, the side length of a BORG SDSS voxel $L=\frac{1}{256}·750\mathrm{Mpc}{h}^{-1}\approx 2.9\mathrm{Mpc}{h}^{-1}$. The late-time posterior can be compactly summarised by taking the mean and standard deviation (SD) of the samples on a per-voxel basis. In this work, we use both the individual samples and the summary cubes – the latter of which we call the BORG SDSS mean and SD. However, by using these summary cubes, information contained in higher moments of single-voxel posteriors and inter-voxel correlations remains unused.

Leclercq et al. (2015) have extended the BORG SDSS by calculating, for each sample, the Cosmic Web classification as stipulated by the T-web definition (Hahn et al. 2007). The result is a probabilistic classification with a marginal distribution for each voxel. As one aspect of a broader information theoretic analysis, Leclercq et al. (2016) demonstrated how these classifications can be used to predict galaxy properties, such as g − r colours. We used these classifications to characterise the dynamical environments of luminous giants and general RGs.

Both the BORG SDSS posterior and a posterior from another BORG run, the BORG 2M++ (Jasche & Lavaux 2019), have been used before to relate properties of active galaxies to the density of the enveloping Cosmic Web (Frank et al. 2016; Porqueres et al. 2018). Our work is the first to relate properties of radio galaxies to BORG (or BORG-like) Cosmic Web reconstructions.

3. Methods

As pointed out in Sect. 1, the number of known giants has increased substantially in recent years. For example, within the Local Universe covered by the SDSS DR7, the manual search of Oei et al. (2023a) alone has quintupled the number of known giants with spectroscopic redshifts – from 52 to 260. Meanwhile, the BORG SDSS now offers the first physically principled, probabilistic reconstruction of the total matter density field over this volume. In this section, we explain how we combined both advances by localising giants and general RGs within the BORG SDSS.

3.1. BORG SDSS localisation procedure

For each RG in our two samples (be it a giant or a general RG), we first transformed its right ascension, declination, and redshift into a vector with comoving coordinates r ≔ [x,y,z]^⊤ following Jasche et al. (2015)’s coordinate system convention. This transformation is cosmology-dependent. The BORG SDSS adopts the cosmological parameters of Jasche et al. (2015) to convert SDSS DR7 MGS redshifts into radial comoving distances and subsequently infer the structure formation history of the Local Universe. In order to obtain valid RG localisations, it was imperative that we used the same conversion between redshift and radial comoving distance. We therefore adopted the Jasche et al. (2015) cosmology for this particular step. Afterwards, we simply associated each RG to the BORG SDSS voxel nearest to r.

After associating RGs to voxels, we checked whether the late-time BORG SDSS density field at these locations is sufficiently constrained. More precisely, we evaluated whether – at each RG’s voxel – the survey response operator of Jasche et al. (2015)’s lowest r-band absolute magnitude bin (−21.00 < M_r < −20.33), R⁰, equals or exceeds some threshold $R_{\text{min}}^0$. If indeed $R^0(\mathit{r})\ge R_{\text{min}}^0$, we retained the RG for Cosmic Web analysis; if not, we discarded it. In this work, we chose $R_{\text{min}}^0=0.1$. Higher choices for $R_{\text{min}}^0$ rid the sample of comparatively uninformative (that is, more prior-dominated) density and dynamical state measurements, but come at the obvious cost of a reduced sample size.

In Fig. 3, we show six example localisations of giants in the Local Universe. We pinpointed the giants in the upper two panels to galaxy clusters, the giants in the middle two panels to filaments, and the giants in the bottom two panels to sheets. Next, in order to obtain density distributions for each retained RG, we explored two methods.

Fig. 3. Example localisations of giant radio galaxies within the large-scale structure of the Local Universe. The top row shows two cluster giants, the middle row shows two filament giants, and the bottom row shows two sheet giants. For each giant, we show a slice of constant Cartesian comoving z through the late-time BORG SDSS posterior mean total matter density field and a LoTSS DR2 6″ image at νobs = 144 MHz (inset). Outside of the SDSS DR7–constrained volume, the posterior mean tends to the Universe’s late-time mean total matter density ${\overline{\rho }}_0$. The locations of the giants are marked by white circles.

3.1.1. Fixed voxel method

In the ‘fixed voxel method’, we considered for each RG the marginal posterior density distribution at its voxel – that is the posterior density distribution for that voxel in isolation – despite the fact that the BORG SDSS provides reconstructions with complex inter-voxel density correlations. This method is the simplest of the two methods we have used.

3.1.2. Flexible voxel method

The fixed voxel method has at least two disadvantages. One issue is the fact that the BORG SDSS has inferred the Cosmic Web with a limited set of bright SDSS DR7 galaxies, which (amongst other factors) causes reconstruction uncertainty. In practice, this means that a given cluster or filament may morph and wiggle around in different BORG SDSS samples. If one does not wiggle around the voxel to sample from accordingly (but sticks with the same voxel all the time), one regularly samples outside of the cluster or filament in which the RG of interest resides. This, of course, biases the inferred densities low.

On top of Cosmic Web reconstruction uncertainty, host galaxy peculiar motion uncertainty leads to additional difficulty in the determination of an RG’s Cosmic Web density. In Appendix A, we show that spectroscopic redshift uncertainties generally lead to sub-voxel localisation uncertainty, even after taking low-mass galaxy cluster–like peculiar motion into account. For high-mass galaxy clusters however, peculiar motion can cause multi-voxel localisation errors. We also show that photometric redshift uncertainties cause localisation uncertainties of 10¹–10² Mpc, or up to tens of voxels, that are unworkably large. Cosmic Web localisation therefore only seems possible for RGs with spectroscopically detected hosts.

To counteract the fixed voxel method’s tendency to sample densities outside of the clusters and filaments in which the RGs truly reside, we propose a flexible voxel method. In this method, the voxel considered hitherto serves as a reference voxel, around which we search for the most likely correct voxel to sample from. More precisely, we iterate over 1000 BORG SDSS samples equally spaced within the MCMC, and consider for each sample all voxels in a sphere of radius 5 Mpc h⁻¹ around the reference voxel. We then simply adopt the voxel with the highest density as the most likely correct voxel for that sample. This procedure encapsulates our prior belief that the massive ellipticals which give rise to RGs are more likely to occur in a high-density region than in a low-density region of the same volume.

The search radius of 5 Mpc h⁻¹ chosen here is arbitrary to some degree. Clearly, for the flexible voxel method to be any different from the fixed voxel method, this radius must exceed the side length of a single voxel, 2.9 Mpc h⁻¹. For larger radii, we are able to correct for larger peculiar motion errors, and thus provide more accurate densities for RGs located in massive galaxy clusters. At the same time, for larger radii, we are at risk of straying too far from the reference voxel; for example, this could lead to sampling cluster-like densities for RGs that actually reside in an adjacent filament. Because the majority of RGs appears to reside in filaments, and because the gravity solver of the BORG SDSS significantly limits the usefulness of cluster densities (see Sect. 5), we chose a relatively ‘small’ search radius of 5 Mpc h⁻¹ – less than two voxels in each direction.

In Fig. 4, we compare the mean Cosmic Web densities inferred through the two methods for all RGs considered in this work. Flexible voxel–based relative densities are typically a factor two higher. The methods agree most often at the lowest redshifts, where the BORG SDSS Cosmic Web reconstructions are least uncertain.

Fig. 4. Comparison between Cosmic Web densities of giants and LoTSS DR1 RGs in the Local Universe, inferred via two variations of our BORG SDSS–based approach. For each RG, we show the mean density measured through Sect. 3.1.1’s fixed voxel method (horizontal axis) and through Sect. 3.1.2’s flexible voxel method (vertical axis).

3.2. BORG SDSS localisation in practice

As described in Sect. 3.1, for each giant and general RG in our samples, we measured the marginal posterior density RV. Of 281 giants that lie within the part of the BORG SDSS volume where $R^0(\mathit{r})\ge R_{\text{min}}^0$, there are 260 with a spectroscopic redshift (93%), which are therefore suitable for Cosmic Web analysis. Of these, 208 (80%) are LoTSS DR2 discoveries (Oei et al. 2023a). In exactly the same way, of the 1870 LoTSS DR1 RGs that lie within the constrained BORG SDSS volume, we retained 1443 specimina with spectroscopic redshifts (77%), which we selected for Cosmic Web analysis.

At the BORG SDSS resolution of 2.9 Mpc h⁻¹ per voxel side, the baryonic matter density field approximately equals the dark matter density field scaled down by a factor $\frac{{\mathrm{\Omega }}_{\mathrm{BM},0}}{{\mathrm{\Omega }}_{\mathrm{DM},0}}$. The BORG SDSS does not distinguish between these two fields. Instead, it provides the sum of the baryonic and dark matter density field: the total matter density field. In this article, we exclusively mention relative total matter densities – total matter densities divided by today’s cosmic mean total matter density Ω_M, 0ρ_c, 0, where ρ_c, 0 is today’s critical density. To be less verbose, we refer to just ‘relative densities’, and write 1 + δ, where δ denotes overdensity in the usual sense.

4. Results

In this section, we present the first empirical distributions of the Cosmic Web density and dynamical state around luminous giants and general RGs. We also determine the RG radio luminosity–Cosmic Web density relation, and test whether it can cause the density distribution discrepancy between luminous giants and general RGs.

4.1. Cosmic Web density distributions

The relative density posterior distributions of individual RGs resemble lognormal distributions, irrespective of whether the fixed or flexible voxel method was used to determine them. Figure 5 demonstrates this through example marginals for a typical luminous giant and a typical LoTSS DR1 RG. Indeed modelling 1 + Δ_GRG, obs | 1 + Δ_GRG = 1 + δ ∼ Lognormal(μ, σ²) and 1 + Δ_RG, obs | 1 + Δ_RG = 1 + δ ∼ Lognormal(μ, σ²), we can succinctly summarise each RG’s measured density distribution with two parameters. We provide these, for 50 out of 260 giants, in Table 1. (For access to such data for all giants, and for similar data on LoTSS DR1 giants, see the table’s footnote).

Fig. 5. Distributions for the measured relative total matter density RVs 1 + ΔGRG, obs | 1 + ΔGRG = 1 + δ and 1 + ΔRG, obs | 1 + ΔRG = 1 + δ for an individual giant (top row; solid), and for an individual RG (bottom row; hatched). Fixed voxel method densities (left column; green) are lower than flexible voxel method densities (right column; orange). The MLE-fitted lognormal PDFs (solid curves) demonstrate that the distributions are almost lognormal.

By aggregating just the means of these marginal relative density distributions, we could analyse the distributions for our observed populations as a whole. In green, Fig. 6 shows kernel density estimated (KDE) relative density distributions for both luminous giants (top panels) and the broader population of RGs (bottom panels) in the Local Universe. The panels in the left column represent the fixed voxel method, whilst the panels in the right column represent the flexible voxel method. More precisely, these KDE distributions approximate the distributions of the observed GRG relative density RV 1 + Δ_GRG, obs and the observed RG relative density RV 1 + Δ_RG, obs.

Fig. 6. Probability density functions (PDFs) of the relative total matter density RVs of 260 Local Universe giants (top row), and of 1443 Local Universe RGs (bottom row), determined through the fixed voxel method (left column) and flexible voxel method (right column). These RVs correspond to the density field smoothed to a scale of 2.9 Mpc h−1. We also show PDFs of gamma-distributed RVs with parameters obtained via MLE, both before (solid lines) and after (dash-dotted lines) a heteroskedasticity correction. We warn that these distributions are affected by surface brightness selection, and thus represent known populations only.

We sought to summarise these distributions parametrically. After testing various two-parameter distributions for continuous, non-negative RVs (such as the gamma distribution and the lognormal distribution), the KDE distributions of 1 + Δ_RG, obs in the bottom panels – which are based on all 1443 selected LoTSS DR1 RGs – appeared best approximated by a gamma distribution⁶. Maximum likelihood estimation (MLE) suggested 1 + Δ_RG, obs ∼ Γ(k = 2.1, θ = 1.8) for the fixed voxel method and 1 + Δ_RG, obs ∼ Γ(k = 3.6, θ = 1.8) for the flexible voxel method. Similarly, fitting a gamma distribution to 1 + Δ_GRG, obs through MLE gave 1 + Δ_GRG, obs ∼ Γ(k = 2.5, θ = 1.9) for the fixed voxel method and 1 + Δ_GRG, obs ∼ Γ(k = 4.7, θ = 1.6) for the flexible voxel method, although the latter does not provide a tight fit. These gamma distributions constitute practical, two-parameter representations of the underlying data and are drawn as solid lines in Fig. 6.

The mean of a gamma-distributed RV 1 + Δ ∼ Γ(k, θ) is 𝔼[1 + Δ]=kθ. Thus, for the fixed voxel method, 𝔼[1 + Δ_GRG, obs] = 4.8 and 𝔼[1 + Δ_RG, obs] = 3.8. Similarly, for the flexible voxel method, 𝔼[1 + Δ_GRG, obs] = 7.5 and 𝔼[1 + Δ_RG, obs] = 6.5. Naively, it appears that we can conclude that, in a statistical sense, giants occupy denser regions of the Cosmic Web than radio galaxies in general. Using a two-sample Kolmogorov–Smirnov (KS) test, we formally tested the null hypothesis that the observed giant and general RG relative density distributions of Fig. 6 share a common underlying distribution. For both the fixed and flexible voxel methods, the p-value p ≲ 10⁻⁶. The null hypothesis is thus rejected (for typical significance levels). These distributions come with at least two major caveats, however.

4.1.1. Heteroskedasticity

First, BORG SDSS density measurements are heteroskedastic: higher densities have higher measurement errors than lower densities (see e.g., Fig. 6 of Jasche et al. 2015). Heteroskedasticity causes the observed distributions to differ from the intrinsic distributions; not only by widening them (as also occurs in the more familiar homoskedastic setting), but also by systematically shifting the distributions towards lower densities. Appendix C presents a simple method to infer distributions corrected for this effect. We thus distinguish between 1 + Δ_GRG, obs and 1 + Δ_GRG, and between 1 + Δ_RG, obs and 1 + Δ_RG; in both cases, the latter RV is heteroskedasticity-corrected. As the effect induces minor shifts only, we assumed that the heteroskedasticity-free distributions are also approximately gamma. We then applied the method of Appendix C to both our GRG and RG data. Having no closed-form MLE expressions at our disposal, we obtained the MLE parameters of the heteroskedasticity-free gamma distributions by evaluating the likelihood function over an exhaustive grid of (k, θ) values. For the fixed voxel method, this yielded 1 + Δ_GRG ∼ Γ(k = 3.7, θ = 1.3) and 1 + Δ_RG ∼ Γ(k = 3.0, θ = 1.3); 𝔼[1 + Δ_GRG] = 4.8 and 𝔼[1 + Δ_RG] = 3.9. For the flexible voxel method, this yielded 1 + Δ_GRG ∼ Γ(k = 6.7, θ = 1.1) and 1 + Δ_RG ∼ Γ(k = 4.8, θ = 1.4); 𝔼[1 + Δ_GRG] = 7.4 and 𝔼[1 + Δ_RG] = 6.7. We draw the corresponding PDFs as grey, dash-dotted lines in Fig. 6.

4.1.2. Surface brightness selection

The second – and plausibly most important – caveat is the imprint of selection effects, and in particular the surface brightness selection effect induced by the LoTSS noise level σ_Iν = 25 Jy deg⁻² at θ_FWHM = 6″. This effect causes RGs whose lobe surface brightnesses are below the noise level to evade sample inclusion. Lobe surface brightness depends on the RG’s projected proper length l_p, morphological parameters ϝ_l and ϝ_lν, radio luminosity l_ν, redshift z, and lobe spectral index α. Equation (40) of Oei et al. (2023a) presents an approximate formula for the mean lobe surface brightness B_ν:

$$\begin{array}{c}\hfill B_{\nu }({\nu }_{\mathrm{ref}})=\frac{2·{ϝ}_{l_{\nu }}·l_{\nu }({\nu }_{\mathrm{ref}})}{{\pi }^2·\mathbb{E}[D](\eta ({ϝ}_l))·{ϝ}_l^2·l_{\mathrm{p}}^2·{(1+z)}^{3-\alpha }},\end{array}$$(1)

where ϝ_lν is the fraction of the total radio luminosity that belongs to the lobes, ϝ_l is the fraction of the RG’s axis that lies inside the lobes, $\eta ({ϝ}_l)=\frac{{ϝ}_l}{2-{ϝ}_l}$, and 𝔼[D] is the mean deprojection factor as given by Eq. (A.29) of Oei et al. (2023a):

$$\begin{array}{c}\hfill \mathbb{E}[D](\eta )=(1+\eta )(\frac{\pi }{2}-\frac{\eta }{\sqrt{1-{\eta }^2}}ln(\frac{\sqrt{1-{\eta }^2}}{\eta }+\frac{1}{\eta })).\end{array}$$(2)

To explicitly demonstrate that surface brightness selection biases our Local Universe GRG and RG samples towards short lengths and high radio luminosities, we approximated B_ν(ν_ref ≔ 150 MHz) for all LoTSS DR1 RGs with z < 0.2, assuming ϝ_l = ϝ_lν = 0.3 (Oei et al. 2022) and α = −0.7. Figure 7 shows the results. As mature RGs grow, they trace out a curve approximately parallel to the golden lines (which represent constant radio luminosities), directed towards the bottom-right of the plot. If their lives last long enough, they therefore reach a length – which depends on their radio luminosity – beyond which they disappear into the noise of the LoTSS. In particular, our calculations suggest that the only giants that can be detected through the LoTSS are those with l_ν(ν_ref)≳10²⁴ W Hz⁻¹.⁷ Indeed, the lowest–radio luminosity giant in the LoTSS DR1 GRG sample of Dabhade et al. (2020) – which consists of 239 giants – has l_ν = 1 × 10²⁴ W Hz⁻¹. Clearly, if RGs with l_ν ∼ 10²³ W Hz⁻¹ can also form giants, then this subpopulation will be missing in LoTSS GRG samples. Equivalently, the giants that are included in our LoTSS GRG sample necessarily have radio luminosities l_ν ≳ 10²⁴ W Hz⁻¹. It is overwhelmingly likely that the GRG sample’s radio luminosities are biased high. Concerningly, Croston et al. (2019) have shown that, for RGs at z < 0.4, radio luminosity and environmental richness correlate positively. This gives rise to the possibility that the trend shown in Fig. 6, where giants appear to reside in denser Cosmic Web environments than RGs in general, is a result of surface brightness selection: a GRG sample biased high in radio luminosity will be biased high in Cosmic Web density.

Fig. 7. Mean lobe surface brightness estimates Bν at νref = 150 MHz for LoTSS DR1 RGs at z < 0.2, shown as a function of projected proper length lp. During the majority (∼90%) of their lifetime, RGs grow while maintaining a roughly constant radio luminosity. If, additionally, growth is shape-preserving, surface brightness decreases quadratically with lp. The seven golden lines therefore denote approximate evolutionary tracks of RGs at z = 0 with end-of-life radio luminosities lν ∈ {1021 W Hz−1, …, 1027 W Hz−1}. The region of parameter space where RGs can become detectable giants, lν ≥ lν, min ≔ 1024 W Hz−1, is shaded gold. The region of parameter space occupied by giants, lp ≥ lp, GRG ≔ 0.7 Mpc, is shaded green. The region of parameter space inaccessible to the LoTSS, Bν < σIν = 25 Jy deg−2, is shaded red.

4.2. Radio luminosity–Cosmic Web density relation

Interestingly, our data allow us to revisit the radio luminosity–environmental richness correlation that Croston et al. (2019) have found. In Fig. 8, we plot all 1443 selected LoTSS DR1 RGs in Cosmic Web density–radio luminosity parameter space. By eye, it appears that for higher radio luminosities, the Cosmic Web density distribution shifts upwards. To quantify the apparent tendency that more luminous radio galaxies occupy denser regions of the Cosmic Web, we extended the gamma distribution fitting discussed in Sect. 4.1. In particular, we assumed that the RG relative density RV at some fixed radio luminosity, 1 + Δ_RG | L_ν = l_ν, is gamma distributed:

$$\begin{array}{c}\hfill 1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }\sim \mathrm{\Gamma }(k(\surd ),\theta (\surd )).\end{array}$$(3)

Fig. 8. PDFs of the observed RG relative total matter density RV given a radio luminosity at νref = 150 MHz, using MLE parameter values for the model described in Sect. 4.2. The black contours denote CDF values. We overplot all 1443 selected LoTSS DR1 RGs (dots), with those above the empirical median density coloured red, and those below coloured blue. We used flexible voxel method densities. For fixed voxel method densities, see Fig. D.1.

Instead of assuming constant parameters k and θ, we assumed that k and θ depend on radio luminosity:

$$\begin{array}{c}\hfill \surd =\surd (l_{\nu }){\log }_{10}\left(\frac{l_{\nu }}{l_{\nu ,\mathrm{ref}}}\right).\end{array}$$(4)

The natural minimal extension to assuming that k and θ remain constant as √ changes, is to consider Taylor polynomial approximations of degree one (around √ = 0 – i.e., Maclaurin polynomials):

$$\begin{array}{c}\hfill k(\surd )\approx k(0)+k^{\prime }\left(0\right)·\surd ,\theta (\surd )\approx \theta (0)+{\theta }^{\prime }\left(0\right)·\surd .\end{array}$$(5)

However, this model does not allow for a direct comparison to our data, which are samples from 1 + Δ_RG, obs | L_ν = l_ν rather than from 1 + Δ_RG | L_ν = l_ν. Therefore, we had to extend the model by applying heteroskedastic measurement errors. In line with Appendix C, we assumed

$$\begin{array}{c}\hfill 1+{\mathrm{\Delta }}_{\mathrm{RG},\mathrm{obs}}|1+{\mathrm{\Delta }}_{\mathrm{RG}}=1+\delta \sim \mathrm{Lognormal}(\mu ,{\sigma }^2),\end{array}$$(6)

where the parameters μ = μ(a, b) and σ² = σ²(a, b) are given by

$$\begin{array}{cc}\hfill \mu & =ln(1+\delta )-\frac{1}{2}ln(1+a{(1+\delta )}^{b-2});\hfill \\ \hfill {\sigma }^2& =ln(1+a{(1+\delta )}^{b-2}).\hfill \end{array}$$(7)

The parameters a and b characterise the heteroskedasticity. By considering the means and SDs of the relative densities of all selected LoTSS DR1 RGs, we found a = 0.4 and b = 1.1 for the fixed voxel method and a = 0.4 and b = 0.9 for the flexible voxel method.

In order to obtain a fit that is valid over many orders of magnitude in radio luminosity, rather than only over the order of magnitude in which most data are concentrated, we calculated the likelihood in such a way that all decades in radio luminosity that contain at least 20 RGs receive equal weight. We simply determined the likelihood function, and thus the posterior distribution for a flat prior, through a grid search over the four parameters k(0), θ(0), k′(0), and θ′(0). We defined l_ν, ref ≔ 10²³ W Hz⁻¹. Figure 9 visualises this flat-prior posterior. In addition, Fig. 8 shows percentile scores of 10%, 50%, 90%, 95%, and 99% corresponding to the MAP model. Figures D.1 and D.2 demonstrate that, although the densities are lower, the results remain qualitatively the same upon switching to the fixed voxel method. We summarise the posterior for both methods in Table 2.

Fig. 9. Posterior distribution over k(0), θ(0), k′(0), and θ′(0), based on selected Local Universe LoTSS DR1 RGs. We show all two-parameter marginals of the likelihood function, with contours enclosing 50% and 90% of total probability. We mark the MLE (grey dot) and the maximum a posteriori (MAP; grey cross). The one-parameter marginals again show the MLE (dash-dotted line), a mean-centred interval of standard deviation–sized half-width (hashed region), and a median-centred 90% credible interval (shaded region). We used flexible voxel method densities. For fixed voxel method densities, see Fig. D.2.

As values of k′(0) and especially θ′(0) tend to be positive, we find evidence for a positive radio luminosity–environmental richness relation, in line with the results of Croston et al. (2019). More quantitatively, we call the relation positive (at l_ν = l_ν, ref) when

$$\begin{array}{c}\hfill {\frac{\mathrm{d}\mathbb{E}[1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }]}{\mathrm{d}\surd }|}_{\surd =0}=k(0){\theta }^{\prime }(0)+k^{\prime }(0)\theta (0)>0.\end{array}$$(8)

We combined the posterior with Eq. (9) to obtain the posterior probability that the radio luminosity–Cosmic Web density relation is positive – that is to say that more luminous RGs tend to live in denser regions of the Cosmic Web. For both the fixed and flexible voxel methods, this probability exceeds 99%.

4.2.1. Surface brightness selection

However, we must again be mindful of a possible surface brightness selection effect. RGs of a given radio luminosity may be able to grow larger in more tenuous Cosmic Web environments, as radio luminosity traces jet power (in close-to-linear fashion; e.g., Willott et al. 1999; Shabala & Godfrey 2013; Hardcastle 2018). If this is indeed the case, then the lobe surface brightnesses of RGs in tenuous environments are likely lower than those of RGs with the same radio luminosity in dense environments. As a result, of all RGs with the same radio luminosity, those in tenuous environments are more likely to be missing from observational samples. This selection effect, that favours higher relative densities, holds for all radio luminosities, and could only explain the trend of Fig. 8 if the effect becomes more severe for higher radio luminosities. On the one hand, RGs with higher radio luminosities at a given relative density are expected to grow longer. On the other hand, as is clear from Fig. 7, the projected proper length range over which RGs retain detectable lobe surface brightnesses increases with radio luminosity. In other words, at higher radio luminosities, larger RGs remain detectable.

Whether the net result of these counteracting effects is less or more bias towards higher relative densities as radio luminosity increases, depends on the scaling between radio luminosity and growth. Self-similar growth predicts (Kaiser & Alexander 1997) that an RG’s projected proper length scales with jet power Q_jet, Cosmic Web density at the host galaxy 1 + δ_g, and time since birth t as

$$\begin{array}{c}\hfill l_{\mathrm{p}}\propto Q_{\mathrm{jet}}^{\frac{1}{5-\beta }}·{(1+{\delta }_{\mathrm{g}})}^{\frac{-1}{5-\beta }}·t^{\frac{3}{5-\beta }},\end{array}$$(9)

where −β is the exponent of the local Cosmic Web density profile:

$$\begin{array}{c}\hfill (1+\delta )(r)\propto (1+{\delta }_{\mathrm{g}})·r^{-\beta }.\end{array}$$(10)

Under the additional assumption of equipartition between the relativistic electron–positron kinetic energy density and the magnetic field energy density, and upon neglecting electron–positron energy losses (such as adiabatic, synchrotron, and inverse Compton losses), the RG’s radio luminosity obeys (Shabala & Godfrey 2013)

$$\begin{array}{c}\hfill l_{\nu }\propto Q_{\mathrm{jet}}^{\frac{5+p}{6}}{l}_{\mathrm{p}}^{3-\frac{(4+\beta )(5+p)}{12}}.\end{array}$$(11)

Here, −p is the exponent of the initial post-acceleration electron–positron kinetic energy distribution n(γ)∝γ^−p, valid for Lorentz factors between some minimum and maximum values γ_min and γ_max. Solving for Q_jet, substituting the result in Eq. (10), and solving for l_p, yields – at fixed 1 + δ_g and t –

$$\begin{array}{c}\hfill l_{\mathrm{p}}\propto l_{\nu }^{\frac{12}{36+(5+p)(6-3\beta )}}.\end{array}$$(12)

For a constant-density environment (i.e., β = 0) and a strong-shock spectral index p = 2, one finds $l_{\text{p}}\propto l_{\nu }^{2/13}$. For an environment where density falls off quadratically with distance, β = 2, and therefore $l_{\text{p}}\propto l_{\nu }^{1/3}$. Meanwhile, Eq. (1) predicts that the maximum detectable projected proper length for RGs with a given morphology, at a given redshift, and with given lobe spectral index, scales as $l_{\mathrm{p},\max }\propto \sqrt{l_{\nu }}$. These arguments suggest that l_p, max increases more rapidly with l_ν than l_p does. Thus, amongst a family of RGs with varying radio luminosity but seen at the same age (i.e., time after birth) and with the same morphology and Cosmic Web density, the most luminous is most easily detectable. This provides some reassurance that the trend seen in Fig. 8 is real: surface brightness selection at a fixed relative density seems to become less severe as radio luminosity increases.

However, there are at least two caveats here. First, more luminous RGs may live longer than less luminous RGs, as jointly suggested by the radio luminosity–stellar mass relation and the stellar mass–radio AGN occurrence relation of Sabater et al. (2019). As a result, in a snapshot of the RG population at a given instant, more luminous RGs may typically be older than less luminous RGs. Figure 10 illustrates this point⁸. If more luminous RGs indeed tend to be older than less luminous RGs, then the relation between projected length and radio luminosity at a fixed Cosmic Web density for a snapshot population will be steeper than the one of Eq. (13), which assumes equal ages. Secondly, we warn for the possibility – which remains subject of debate (Mingo et al. 2019) – that the RG shape distribution at a fixed Cosmic Web density may change with radio luminosity.

Fig. 10. Monte Carlo simulation of a population of nine galaxies (thick stripes) with successive active and quiescent phases. We mark the active phases (thin stripes). Each galaxy spends a different fraction ϝa of the time in the active phase. We model the time spent in both phases as draws from exponential distributions: ta ∼ Exp(λa) and tq ∼ Exp(λq). We assume that all galaxies share the same quiescent phase rate parameter λq = 10 Gyr−1, although this may not be entirely correct (Turner & Shabala 2015). This assumption fixes the active phase rate parameter: ${\lambda }_{\mathrm{a}}=({ϝ}_{\mathrm{a}}^{-1}-1){\lambda }_{\mathrm{q}}$. On human time scales, we can access only a snapshot view of this population. During the present snapshot (vertical line), seven galaxies generate RGs, and could be included in an RG sample. RG ages at the snapshot (hatched parts of the thin green stripes) generally increase with ϝa.

4.2.2. Predicting the density distribution of observed giants by treating them as luminous general RGs

We then tested whether Fig. 6’s higher densities for observed giants than for observed general RGs could be explained by these giants’ higher radio luminosities and Fig. 8’s radio luminosity–Cosmic Web density relation. We first calculated the RG relative total matter density RV for radio luminosities of l_ν, min and higher. By applying standard PDF identities, we found

$$\begin{array}{cc}& f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }\ge l_{\nu ,\min }}(1+\delta )=\hfill \\ \hfill & \frac{{\int }_{l_{\nu ,\min }}^{\infty }{f}_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }}(1+\delta )·f_{L_{\nu }}(l_{\nu })\mathrm{d}{l}_{\nu }}{{\int }_{l_{\nu ,\min }}^{\infty }{f}_{L_{\nu }}(l_{\nu })\mathrm{d}{l}_{\nu }}=\hfill \\ \hfill & \frac{{\int }_{l_{\nu ,\min }}^{\infty }{f}_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }}(1+\delta )·f_{L_{\nu }|L_{\nu }\ge l_{\nu ,\min }}(l_{\nu })\mathrm{d}{l}_{\nu }}{{\displaystyle {\int }_{l_{\nu ,\min }}^{\infty }}{f}_{L_{\nu }|L_{\nu }\ge l_{\nu ,\min }}(l_{\nu })\mathrm{d}{l}_{\nu }}.\hfill \end{array}$$(16)

To proceed, one must specify f_{Lν | Lν ≥ lν, min}(l_ν) – up to a constant, at least. We modelled L_ν | L_ν ≥ l_ν, min ∼ Pareto(l_ν, min, λ), where λ is the tail index or rate parameter of the associated exponential distribution⁹. For both l_ν, min = 10²³ W Hz⁻¹ and l_ν, min = 10²⁴ W Hz⁻¹, we found λ_MLE = 0.6 (based on 1 × 10³ and 3 × 10² RGs, respectively). For higher l_ν, min, λ_MLE increases somewhat, but becomes less reliable because of decreasing sample sizes. We therefore stuck with λ_MLE = 0.6. We then evaluated Eq. (17) for l_ν, min = 10²⁵ W Hz⁻¹ and added heteroskedastistic noise in order to compare our prediction to observations. We propagated the correlated uncertainties on the parameters k(0), θ(0), k′(0), and θ′(0), which together determine f_{1 + ΔRG | Lν = lν}, by rejection sampling from the joint posterior distribution. We show the resulting mean prediction with standard deviation half-width (pink) in Fig. 11.

Fig. 11. Density distribution predicted for observed giants combining Sect. 4.2’s RG radio luminosity–Cosmic Web density relation and a radio luminosity cut-off lν, min = 1025 W Hz−1 (pink); in addition, we show the KDE density distributions for giants with gamma distribution MLE fit (solid green) and for general RGs (hatched green) – all as in the left column of Fig. 6.

Although the prediction is based solely on observations of general RGs – and so is built without knowledge of the observed GRG distribution (solid green KDE) – it appears to fit it better than a direct MLE gamma fit to the GRG data (solid green curve) does. However, a better fit is not necessarily also a good fit. We therefore tested whether the prediction is consistent with our sample of Cosmic Web densities for observed giants. A KS test suggested an answer in the affirmative, as it yielded a p-value of 0.7.

Section 4.1 presented a paradox: known giant radio galaxies reside in denser Cosmic Web environments than their smaller kin, in apparent violation of the long-standing hypothesis that giants should emerge predominantly in the dilute Cosmic Web. The modelling in this section suggests that the higher densities found for known giants than for general RGs can be explained by a combination of these giants’ higher radio luminosities and the general relation between RG radio luminosity and Cosmic Web density. Fainter giants – which might inhabit the more dilute Cosmic Web – might exist, but are currently not observationally accessible.

4.3. Radio galaxy number densities

After accreting baryons, supermassive black holes (SMBHs) can launch jets that give rise to radio galaxies. Although it is clear that SMBH jets and RGs play an important role in galaxy evolution and cosmology, the physics of jet launching is not yet fully understood. Moreover, because SMBHs are astronomical unit–sized, galaxies are kiloparsec-sized (ratio 10⁸–10⁹), and the Cosmic Web is megaparsec-sized (ratio 10¹¹–10¹²), it is not yet possible to build cosmological simulations in which a realistic interplay between SMBHs, their host galaxies, and the enveloping Cosmic Web naturally arises. Thus, for now, major advances in our understanding of RGs on the cosmological scale must come from observations instead of from simulations. It is therefore of considerable interest to observationally constrain the occurrence of RGs as a function of Cosmic Web density.

Enticingly, our relative total matter density measurements – in combination with the BORG SDSS HMC Markov chain samples in their entirety – allow for a determination of the RG number density as a function of Cosmic Web density: n_RG(1 + δ). To see why, we first emphasise that the PDF f_1 + ΔRG, which we have approximated through observations and visualised in the bottom row of Fig. 6, provides a key ingredient for n_RG(1 + δ). It shows that if a volume of cosmological extent is surveyed, one will obtain more RGs with moderate (filament-like) densities than with low (void-like) or high (cluster-like) densities. We remark that in such a survey sample the number of RGs whose densities fall within a given interval depends not only on the RG number density at the interval’s densities, but also on the prevalence of environments with such densities. As a result, if a sample contains more filament-inhabiting RGs than, say, cluster-inhabiting RGs, one need not necessarily conclude that the RG number density in filaments is higher than in clusters; in fact, the RG number density in the latter environment could well be higher, provided that clusters are sufficiently rare.

Conveniently, the density fields of the BORG SDSS HMC Markov chain samples enable one to quantify the rarity of each Cosmic Web environment through a PDF f_1 + ΔCW, which corresponds to the Cosmic Web relative density RV 1 + Δ_CW. This RV represents the relative total matter density at a randomly chosen point in the Local Universe (and at the BORG SDSS resolution of 2.9 Mpc h⁻¹). We determine the distribution of 1 + Δ_CW simply by binning the density fields of a few hundred Markov chain samples. As a result, f_1 + ΔCW is best determined at low densities, which correspond to the most common environments, and least so at high densities, which are rare¹⁰. As we expect f_1 + ΔCW to be smooth, we applied a Savitzky–Golay filter (of polynomial order 1) to the binned data; however, this step is not essential.

Appendix E shows that n_RG(1 + δ), the RG number density as would arise in an environment of constant density 1 + δ, obeys

$$\begin{array}{c}\hfill \frac{n_{\mathrm{RG}}(1+\delta )}{{\overline{n}}_{\mathrm{RG}}}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(19)

Here, ${\overline{n}}_{\mathrm{RG}}$ is the cosmic mean RG number density. Thus, a simple point-wise division of the PDFs of 1 + Δ_RG and 1 + Δ_CW yields the RG number density at a given Cosmic Web density relative to the cosmic mean value.

Unfortunately, as discussed in Sect. 4.1, our measurements of f_1 + ΔRG suffer from surface brightness selection. To avoid selection effects, we formulated a more specific version of Eq. (20). In particular, we considered the number density of RGs per unit of radio luminosity as a function of radio luminosity and relative density: _RG(l_ν, 1 + δ) – again, relative to a cosmic mean value: ${\overline{}}_{\mathrm{RG}}(l_{\nu })$. Appendix E shows that this ratio equals

$$\begin{array}{c}\hfill \frac{{}_{\mathrm{RG}}(l_{\nu },1+\delta )}{{\overline{}}_{\mathrm{RG}}(l_{\nu })}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(20)

This radio luminosity–dependent relative number density is, as before, equal to a point-wise division of PDFs. With respect to the right-hand side (RHS) of Eq. (20), only the PDF in the numerator has changed. Section 4.2 explains how we modelled this PDF. It depends on the four parameters k(0), k′(0), θ(0), and θ′(0), for which the data shown in Fig. 8 induce the posterior distribution shown in Fig. 9. We evaluated the RHS of Eq. (21) by rejection sampling from the posterior, calculating f_{1 + ΔRG|Lν = lν} for each sampled quartet (k(0),k′(0),θ(0),θ′(0)), and dividing the resulting PDFs point-wise by f_1 + ΔCW. In Fig. 12, for four observationally accessible radio luminosities, we visualise the posterior mean alongside a range that is two standard deviations wide.

Fig. 12. Number densities of RGs with a given radio luminosity lν, as a function of Cosmic Web density at the 2.9 Mpc h−1 scale. The number density functions reflect conditions in the Local Universe, and are given relative to the cosmic mean number density at lν. The solid curves reflect the posterior mean, whilst the shaded areas around the mean denote −1 to +1 posterior standard deviation ranges. We used flexible voxel method densities.

As expected, RG number densities generally increase with Cosmic Web density. The more luminous RGs become, the more pronounced the contrast between their number density in clusters and their cosmic mean number density is. Our results indicate that high-luminosity RGs, with l_ν ∈ 10²⁵–10²⁷ W Hz⁻¹, attain number densities in clusters that are ∼10² times higher than average. Interestingly, for radio luminosities l_ν ∈ 10²¹–10²³ W Hz⁻¹, we obtain statistically significant, preliminary evidence that the RG number density peaks at a cluster-like density that depends positively on l_ν¹¹. Beyond this density, these low-luminosity RGs become less prevalent per unit of volume, possibly because the underlying galaxy population starts generating more luminous RGs instead.

We stress that these results are tentative, as both f_{1 + ΔRG|Lν = lν} and f_1 + ΔCW suffer from systematic errors at higher densities. These systematic errors are not reflected in the uncertainties of Fig. 12. In the case of f_{1 + ΔRG|Lν = lν}, the gamma distribution approximation loses validity in the tail: beyond some point, it will consistently under- or overestimate the probability density. Our LoTSS DR1 sample contains only a modest number (∼10²) of RGs at relative densities 1 + δ ≳ 10 – especially when selecting those around a given value of l_ν only (in which case just ∼10¹ RGs remain). As a result, in its tail the gamma PDF takes on the character of an extrapolation function rather than of a fitting function, and should consequently be treated with caution. In the future, larger RG samples can resolve this situation. In the case of f_1 + ΔCW, the 2LPT model of structure formation used by the BORG SDSS loses validity at cluster densities. As a result, this work’s estimate of f_1 + ΔCW underpredicts the prevalence of cluster-like densities. Newer BORG runs, such as the BORG 2M++, employ an enhanced gravity solver (Jasche & Lavaux 2019), relieving this limitation.

Although preliminary in nature, the inferences shown in Fig. 12 could in principle be used to calibrate cosmological simulations that feature RG feedback. The fact that our number densities are given with respect to a cosmic mean value is not problematic, because such a cosmic mean number density is straightforward to calculate from simulation snapshots. However, apart from making sure that the RG number densities in simulated environments of various densities are correctly proportioned, it is also important to calibrate the absolute level of RG feedback in simulations. Interestingly, there are variations on Eq. (20) for which an absolute number density can be explicitly calculated with current-day data. One can consider, for example, non-giant radio galaxies (i.e., l_p < l_p, GRG ≔ 0.7 Mpc) in the Local Universe with radio luminosities l_ν ≥ 10²⁴ W Hz⁻¹. Figure 7 suggests that this RG subpopulation is already fully accessible by the LoTSS. Calling these RGs luminous non-giant radio galaxies (LNGRGs), we have

$$\begin{array}{c}\hfill n_{\mathrm{LNGRG}}(1+\delta )={\overline{n}}_{\mathrm{LNGRG}}·\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{LNGRG}}}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(21)

We calculated ${\overline{n}}_{\mathrm{LNGRG}}$ by counting all LNGRGs in the LoTSS DR1 footprint up to z_max, finding N_LNGRG = 172. Because of calibration problems and higher noise levels towards the edges of the footprint, a fraction of the footprint is – artificially – devoid of LNGRGs. We thus applied a correction factor to N_LNGRG of 1.2, yielding N_LNGRG = 200. Assuming that N_LNGRG ∼ Poisson(λ_LNGRG), maximum likelihood estimation simply suggests ${\widehat{\lambda }}_{\mathrm{LNGRG},\mathrm{MLE}}=200$. We combined this MLE Poisson distribution for N_LNGRG with the comoving volume V = 16 × 10⁶ Mpc³ to obtain a probability distribution for ${\overline{n}}_{\mathrm{LNGRG}}$. We find ${\overline{n}}_{\mathrm{LNGRG}}=12\pm 1{(100\mathrm{Mpc})}^{-3}$. Next, we approximated f_{1 + ΔLNGRG} ≈ f_{1 + ΔRG | Lν ≥ 10²⁴ W Hz⁻¹}, and calculated the latter through Eq. (17) and Sect. 4.2’s subsequent procedure. Taking f_1 + ΔCW as before, we obtained a probability distribution over LNGRG number density functions, visualised in Fig. 13.

Fig. 13. Number density of luminous non-giant radio galaxies in the Local Universe, as a function of Cosmic Web density. We define these as RGs with 150 MHz radio luminosities lν ≥ 1024 W Hz−1, projected proper lengths lp < lp, GRG ≔ 0.7 Mpc, and redshifts z < zmax ≔ 0.16. The Cosmic Web densities encompass baryonic and dark matter, are defined on a 2.9 Mpc h−1 scale, and are given relative to the cosmic mean. We used flexible voxel method densities here. We show the posterior mean with a shaded −1 to +1 posterior standard deviation range. We also show five individual realisations (grey).

As it should, the LNGRG number density function follows the same trend as the functions of Fig. 12; the novelty here is the prediction of absolute number densities. We reiterate that systematic errors at cluster-like densities render only the left side of the plot reliable. It predicts that, in filaments, n_LNGRG ∼ 10⁻¹ (10 Mpc)⁻³, or n_LNGRG ∼ 10⁻⁴ Mpc⁻³. As major filaments have volumes V ∈ 10²–10³ Mpc³, only one in ten major filaments harbours an LNGRG.

4.4. Cosmic Web dynamical state distributions

For each RG in the two samples, we obtained a probability distribution over T-web environment classes from the BORG SDSS extensions of Leclercq et al. (2015). The T-web scheme classifies the Cosmic Web at each point on the basis of its local gravitational dynamics. On an intuitive level, the scheme counts the number of dimensions along which the orbit of a test particle released at rest with respect to the environment is stable. More formally, at any location the Hessian of the gravitational potential, or tidal tensor T, has either three, two, one, or zero positive eigenvalues; the environment is interpreted correspondingly as a cluster, filament, sheet, or void. Leclercq et al. (2015) have determined the T-web environment for each voxel of each BORG SDSS sample. By iterating over samples and counting the frequencies of the four environment classes on a per-voxel basis, we obtained a marginal environment class distribution p = (p₁, p₂, p₃, p₄) for each voxel. In Fig. 14, we show how each of p’s components is distributed for luminous giants. We note the following.

Fig. 14. Bayesian Cosmic Web reconstructions allow one to determine, for each giant, a probability distribution p over the four T-web environment classes (clusters, filaments, sheets, and voids). Here we summarise these distributions: each panel shows how a component of p ≔ (p1, p2, p3, p4) is itself distributed.

For 99% of our giants, the most likely environment class is more than 50% probable. For 20% of giants, it is more likely than not that they inhabit a cluster (i.e., p₁ > 50%); for 68% of giants, it is more likely than not that they inhabit a filament (i.e., p₂ > 50%); and for 11% of giants, it is more likely than not that they inhabit a sheet (i.e., p₃ > 50%). Finally, for 0% of giants, it is more likely than not that they inhabit a void (i.e., p₄ > 50%).

Some giants inhabit sheets (e.g., those in the bottom row of Fig. 3), but concluding this in individual cases rarely has a high degree of certainty: for just 0.5% of giants, it is very likely that they inhabit a sheet (i.e., p₃ > 80%).

A cluster environment is very unlikely (i.e., p₁ < 20%) for a clear majority of 63% of giants, whilst a filament environment is very unlikely (i.e., p₂ < 20%) for only 9% of giants; the filament environment hypothesis can thus seldom be excluded. A sheet environment is very unlikely (i.e., p₃ < 20%) for 71% of giants; finally, a void environment is nearly always very unlikely (i.e., p₄ < 20%) – this holds for 99% of giants.

The corresponding distributions for LoTSS DR1 RGs (not shown) are largely similar. This is also clear from Fig. 15, in which we compare the environment class distributions of luminous giants, LoTSS DR1 RGs, and the Local Universe in its entirety. Evidently, luminous giants and general RGs are not scattered uniformly throughout the Cosmic Web, but trace cluster and filament environments. Remarkably, nearly a quarter of all luminous giants occurs in clusters, which have a volume-filling fraction (VFF) of just 1–2%. It is easy to derive that the ratio between a GRG environment class probability (left pie) and the corresponding VFF (right pie) yields the GRG number density of that class relative to the cosmic mean GRG number density. Of course, the analogous statement for RGs also holds. In Table 3, we present relative number densities for luminous giants and general RGs in the Local Universe. Clusters are strongly (G)RG overdense; voids are strongly (G)RG underdense.

Fig. 15. Dynamical environment classification of giant radio galaxies (left), radio galaxies in general (centre), and the Local Universe in its entirety (right; from Table 3 of Leclercq et al. 2015). We define clusters, filaments, sheets, and voids in the T-web sense. These distributions depend on the scale to which the Cosmic Web density field is smoothed; in this case, the smoothing scale is 2.9 Mpc h−1. If giants and RGs were scattered uniformly throughout the Cosmic Web, their distributions would be similar to that of the volume-filling fractions (VFFs). Instead, observed RGs – and observed giants in particular – favour cluster and filament environments.

The left and central pie charts of Fig. 15 can be regarded as prior distributions over Cosmic Web dynamical states for luminous giants and RGs in or near the Local Universe. This fact can be useful whenever direct BORG SDSS environment characterisation is not possible. For a randomly picked luminous giant, Fig. 15’s left pie chart specifies the prior probability distribution over dynamical states to be p = (23%,59%,17%,1%). Interestingly, there is often additional information available to constrain the dynamical state. For example, multiwavelength imagery demonstrates that Alcyoneus (Oei et al. 2022), one of the longest known RGs, does not inhabit a galaxy cluster; in addition, the number of luminous galaxies in its vicinity excludes the possibility of a void environment¹². Combining this information with the prior, we obtain – after rounding – the posterior probability distribution p = (0%,80%,20%,0%). We thus estimate the probability that Alcyoneus inhabits a filament to be 80%. Of course, this line of reasoning generalises to any luminous giant outside of galaxy clusters that has massive galactic neighbours in its megaparsec-scale vicinity.

5. Discussion

5.1. Cosmic Web density accuracy

In this work, we have determined megaparsec-scale total matter densities around RGs in the Local Universe, and used these to compare the environments of luminous giants to those of general RGs. A natural question to ask is how accurate these densities are, and to which degree they can be used in other analyses, such as in the fitting of dynamical models to individual RGs (e.g., Hardcastle 2018). To quantify accuracy, one would ideally know the ground-truth densities that the inferred densities are meant to approach. Problematically, however, no ground-truth densities are known – especially for filaments, the environment type which most RGs seem to inhabit. For Local Universe clusters though, approximate masses are known. We exploited this fact to test the accuracy of our densities for cluster RGs. We immediately remark that our density accuracy is likely to vary with density, so that the accuracy of cluster RG densities might not be informative of the accuracy of filament RG densities. This is because the 2LPT gravity solver of the BORG SDSS has particularly limited validity in the strongly non-linear cluster regime.

To test the cluster RG density accuracy, we cross-matched our samples with Wen & Han (2015)’s SDSS DR12–based galaxy cluster catalogue. To specify each cluster’s location, this catalogue provides the right ascension, declination, and redshift of the brightest cluster galaxy (BCG). First, we selected all clusters in the Local Universe (z < z_max). Next, we cross-matched in an angular sense, and considered an RG matched with a BCG if their angular separation is less than 3″. We thus retained 60 BCG giants and 113 BCG LoTSS DR1 RGs. Figure 2 shows six members (10%) of the BCG giants sample. The perturbed jet and lobe shapes seen in the radio, and the massive ellipticals and crowded environments seen in the optical, together visually demonstrate the adequacy of our BCG giant identification. We could identify LoTSS DR1 BCG RGs with equal confidence.

As a quality check of the cluster masses, we tested whether they are accurate enough to contain fingerprints of known physical effects. In particular, we tested whether clusters with RG-generating BCGs are more massive than other Local Universe clusters. The AGN feedback paradigm suggests that this is the case: in low-mass clusters, AGN feedback converts cool cores to non-cool cores within a gigayear, shutting down further AGN feedback for many gigayears to come; by contrast, high-mass clusters maintain stable AGN feedback cycles with gigayear-scale periods (Nobels et al. 2022), and thus their BCGs have a higher probability of being observed while generating RGs. The median mass of clusters with RG-generating BCGs is M₅₀₀ = 1.1 × 10¹⁴ M_⊙, while for other clusters it is M₅₀₀ = 1.0 × 10¹⁴ M_⊙. A Kolmogorov–Smirnov test yielded a p-value p = 7%, suggesting to reject the null hypothesis of equality in distribution only for high significance levels (e.g., α = 10%). We also tested whether clusters with giant-generating BCGs are more massive than other clusters. The median mass of clusters with giant-generating BCGs is M₅₀₀ = 1.2 × 10¹⁴ M_⊙, while for other clusters it is M₅₀₀ = 1.0 × 10¹⁴ M_⊙. A KS test yielded p = 5%, suggesting to reject equality in distribution only for rather high significance levels (e.g., α = 5% or α = 10%).

As a first quality check of our densities, we tested whether our BCG giants have higher densities than our non-BCG giants. The median relative densities for BCG giants and non-BCG giants are 1 + δ = 8.5 and 1 + δ = 7.2, respectively. A KS test yielded p = 0.2%, indeed suggesting to reject equality in distribution. Similarly, we tested whether our LoTSS DR1 BCG RGs have higher densities than our LoTSS DR1 non-BCG RGs. The median relative densities for BCG RGs and non-BCG RGs are 1 + δ = 8.4 and 1 + δ = 5.9, respectively. A KS test yielded p ≪ 0.1%, again suggesting to reject equality in distribution.

Figure 16 shows the relation between cluster mass and relative density. Although it is qualitatively clear that BCG RGs in more massive clusters have higher inferred densities, the densities are all lower than one would expect from spreading out the cluster mass over a BORG SDSS voxel. (For reference, a cluster of mass M ∼ 10¹⁴ M_⊙ should have a 2.9 Mpc h⁻¹–scale relative density 1 + δ ∼ 40.) This systematic discrepancy cannot be attributed to errors in the cluster masses, as Wen & Han (2015) have made sure that their estimates are unbiased. The alternative is to conclude that our densities are biased low in the cluster regime. This is unsurprising, as the assumptions of the BORG SDSS’s 2LPT gravity solver are violated in clusters. Our densities are least accurate for high-mass clusters, where they can be off by an order of magnitude, and most accurate for low-mass clusters, where they can be less than a factor two too low. It is unclear whether, and if so how, these accuracy characterisations can be extrapolated to the filament density regime. The increased applicability of 2LPT to the filament density regime suggests that the accuracy will increase. We therefore expect our filament densities to be correct up to a small factor of order unity.

Fig. 16. Comparison between environmental estimates for luminous giants (large dots) and LoTSS DR1 RGs (small dots) generated by BCGs in the Local Universe. We compare cluster masses M500 from Wen & Han (2015) with relative densities 1 + δ from the BORG SDSS, and colour the dots based on their cluster probability p1 (relative to the average cluster probability ${\overline{p}}_1$ for luminous giants and LoTSS DR1 RGs). For three randomly chosen giants and three randomly chosen RGs, we visualise uncertainties. Generally, as expected, RG-generating BCGs claimed to be in more massive clusters by Wen & Han (2015) have higher BORG SDSS densities and T-web cluster probabilities. The contours indicate fractions of the density expected through M500, and show that our relative densities fall short significantly in the cluster regime – by a factor that ranges from less than two to more than an order of magnitude. Densities of low-mass clusters are less biased than densities of high-mass clusters. We used flexible voxel method densities here.

Until Bayesian Cosmic Web reconstructions are further developed, we urge caution in using our inferred densities as precision estimates: in filament environments, they should be trusted up to a factor of order unity only; in cluster environments, for the moment, other sources (such as cluster catalogues) appear to provide more accurate density estimates. If one would insist in using the cluster densities, then Fig. 16 suggests that one should correct them by multiplying by a factor of ∼5.

5.2. Cosmic Web dynamical state accuracy

Similar questions of accuracy can be raised about the Cosmic Web dynamical state distributions for luminous giants and RGs presented in the left and central pie charts of Fig. 15. To test the cluster GRG and cluster RG percentages, we again used the cluster catalogue of Wen & Han (2015).

Cluster giants can be subdivided into BCG giants and non-BCG cluster giants. In Sect. 5.1, we determined that 60 out of 260 giants (23%) are BCG giants¹³. We defined non-BCG cluster giants as giants generated by non-BCG galaxies for which the comoving distance to the BCG is less than some megaparsec-scale threshold. We obtained approximate comoving coordinates from right ascensions, declinations, and spectroscopic redshifts, assuming vanishing peculiar motion. In order to allow for a fair comparison to the BORG SDSS, we set the megaparsec-scale threshold to a BORG SDSS voxel side length: 2.9 Mpc h⁻¹. Although much larger than a galaxy cluster virial radius, this threshold does ensure that cluster galaxies with significant peculiar motion – a standard deviation σ_vp = 10²–10³ km s⁻¹ is typical – have a reasonable chance to be linked to their cluster. For this choice of threshold, we found 21 non-BCG cluster giants. We thus find a total of 81 cluster giants (31%), which is more than expected based on the BORG SDSS measurement (${23}_{-3}^{+3}\%$) shown in Fig. 15’s GRG pie chart. This suggests that the BORG SDSS underpredicts the occurrence of cluster giants. Indeed, the strongly underpredicted densities of galaxy clusters discussed in Sect. 5.1 also leave their mark on the accuracy of our dynamical states. For example, the giant shown in the top-left panel of Fig. 2 and tabulated on the tenth row of Table 1, clearly occupies a cluster-like environment with an estimated M₅₀₀ = 1.3 × 10¹⁴ M_⊙, but is erroneously assigned an 80% chance to inhabit a filament.

Similarly, cluster RGs can be subdivided into BCG RGs and non-BCG cluster RGs. In Sect. 5.1, we determined that 113 out of 1443 selected LoTSS DR1 RGs (8%) are BCG RGs. By defining non-BCG cluster RGs in a way analogous to non-BCG cluster giants, we found 109 non-BCG cluster RGs¹⁴. We thus find a total of 222 cluster RGs (15%), consistent with the BORG SDSS measurement (${16}_{-1}^{+1}\%$) shown in Fig. 15’s RG pie chart.

To assess the accuracy of the filament GRG and RG fractions shown in Fig. 15, we cross-matched our giants and selected LoTSS DR1 RGs with the galaxy group and cluster catalogue of Tempel et al. (2017), which provides estimates of M₂₀₀. Wen & Han (2015) define clusters through M₅₀₀ ≥ 0.6 × 10¹⁴ M_⊙, which translates to M₂₀₀ ≥ 0.9 × 10¹⁴ M_⊙: a typical conversion is M₅₀₀ = 0.7 M₂₀₀ (e.g., Pierpaoli et al. 2003). To separate groups from clusters, we therefore defined groups through M₂₀₀ < 0.9 × 10¹⁴ M_⊙. As with our cluster tests, we associated RGs to groups on the basis of estimated comoving coordinates, assuming vanishing peculiar motion. In the case of groups, this assumption is more reasonable. We again used a 2.9 Mpc h⁻¹ association threshold. We associated 155 out of 260 luminous giants (60%) and 903 out of 1443 general RGs (63%) with galaxy groups. These percentages provide lower bounds to the occurrence of filament giants and RGs, as there are also filament giants and RGs that do not reside in galaxy groups. It is therefore likely that the filament GRG and RG fractions of Fig. 15 (${59}_{-3}^{+3}\%$ and ${60}_{-1}^{+1}\%$, respectively) are underestimates.

Finally, as a probe of sheet and void environments, we explored the occurrence of giants and RGs for which the closest galaxy group or cluster is more than 10 Mpc away. We identified 18 out of 260 luminous giants (7%) and 115 out of 1443 general RGs (8%) as occupying sheets and voids. This suggests that the combined sheet and void GRG and RG fractions of Fig. 15 (${18}_{-2}^{+2}\%$ and ${24}_{-1}^{+1}\%$, respectively) are significant overestimates.

Rounded to an accuracy of 10%, our ancillary environment analysis suggests an environment class probability distribution for luminous giants of p = (30%,60%,10%), representing clusters, filaments, and a combination of sheets and voids, respectively. For general RGs, it suggests p = (20%,70%,10%). Both distributions are in fair agreement with the BORG SDSS measurement p = (20%,60%,20%). The ancillary analysis reinforces the idea that, in the Local Universe, luminous giants detectable at LoTSS-like image qualities occupy denser environments than general RGs.

5.3. Cosmic Web density resolution

All Cosmic Web density and dynamical state measurements presented in this work correspond to a density field with a resolution, or smoothing scale, of 2.9 Mpc h⁻¹. In fact, given the granular nature of matter, one must always adopt a smoothing scale in order to define an informative notion of density¹⁵. Ideally, this scale is chosen such that the resulting density field varies smoothly over the physical scales of interest – avoiding both unwanted stochastic variations over space and time (as would occur when the smoothing scale is too small), and excessive blurring of the physical phenomena under study (as would occur when the smoothing scale is too large). In the case of this work, the smoothing scale of the Cosmic Web density field could be considered too large. It is set by the limitations of the BORG SDSS 2LPT gravity solver, the ambition to cover the entire Local Universe in the direction of the SDSS MGS footprint, and – most practically – the relation between the number of voxels and the numerical cost of generating the HMC Markov chain. As galaxy clusters have typical radii of ∼1 Mpc h⁻¹, they are not resolved by the BORG SDSS; the same therefore certainly holds for galaxy groups. Thus, even if the gravitational dynamics of the reconstructions were exact, the utility of the BORG SDSS in improving dynamical model fits to individual RGs is limited: the smoothing scale of the density field is larger than the scale required to resolve (beta model–like) density profiles within filaments and clusters (0.1–1 Mpc h⁻¹). At best, our densities reflect the correct multi-megaparsec-scale averages, but these are still much lower than the, say, central 1 Mpc–scale densities. This is because the massive galaxies that usually host RGs tend to occupy the bottom of their local gravitational potential wells, where the large-scale density is typically higher than everywhere in the vicinity. Relying on simulations, Oei et al. (2023b) have estimated the relation between the 2.9 Mpc h⁻¹–scale relative total matter density 1 + δ (as used throughout this work), and the 1 Mpc–scale relative IGM density 1 + δ_IGM. For galaxies with a stellar mass M_⋆ = 10¹¹ M_⊙, their Fig. 11 suggests 1 + δ_IGM = 10 ⋅ (1 + δ)^0.75. For luminous giants, the mode of the 2.9 Mpc h⁻¹–scale relative total matter density RV is 1 + δ = 3.5 for the fixed voxel method and 1 + δ = 6.3 for the flexible voxel method (see Fig. 6); the corresponding host-centred 1 Mpc–scale relative IGM density is 1 + δ_IGM = 23 for the fixed voxel method and 1 + δ_IGM = 40 for the flexible voxel method. This relation in principle allows one to perform Bayesian inference of the distributions of the 1 Mpc–scale Cosmic Web density around RG hosts – possibly even conditioned on RG radio luminosity, as in Sect. 4.2 – by converting parametrised high-resolution distributions to the lower BORG SDSS resolution, applying heteroskedastic measurement errors, and finally comparing the resulting distribution to data. We leave this to future work.

Concerning resolution, we finally warn that care must be taken when calculating (AGN feedback–related) quantities that relate quadratically to density, such as the local bremsstrahlung emissivity: these tend to be underestimated if computed from low-resolution densities. To see why, we considered ρ, the density at some natural smoothing scale, alongside ⟨ρ⟩, the density averaged over some larger-than-natural smoothing scale. We assumed that only the latter is known. Next, we picked a quantity y = aρ², which we aimed to compute at the larger smoothing scale: ⟨y⟩=⟨aρ²⟩=a⟨ρ²⟩. Naively, one would be tempted to proceed in calculating ⟨y⟩ by assuming ⟨ρ²⟩≈⟨ρ⟩², as to use ⟨y⟩≈a⟨ρ⟩². However, by doing so, one risks underestimating ⟨y⟩, because ⟨ρ²⟩≥⟨ρ⟩² by Jensen’s inequality. In Appendix F, we explicitly calculated the ratio between ⟨ρ²⟩ and ⟨ρ⟩² for sheets and voids, filaments, and clusters. For sheets and voids, this ratio is at most $\frac{4}{3}$; for filaments, it is at most a number of order unity; for clusters, it can exceed a thousand.

5.4. Spectroscopic redshift selection

To localise RGs in the Cosmic Web and subsequently measure their relative total matter densities, we required them to have a spectroscopic redshift. By selecting on spectroscopic redshift availability, we implicitly selected on both Cosmic Web density and radio luminosity. The key reason is that spectroscopic redshifts are preferentially available for galaxies with high optical flux densities – so galaxies that are nearby¹⁶, and galaxies that have high optical luminosities. Galaxies with high optical luminosities tend to inhabit dense environments and, when active, tend to generate RGs with high radio luminosities.

To demonstrate that spectroscopic redshift selection indeed affects our RG radio luminosity–Cosmic Web density relation, Fig. 17 revisits the same Local Universe LoTSS DR1 RG data as shown in Fig. 8, but now with the Pan-STARRS r-band Kron flux densities F_ν of the host galaxies indicated. The top panel of Fig. 17 shows KDEs of 1 + Δ_RG, obs for hosts with F_ν < 1 nJy, and for hosts with F_ν ≥ 1 nJy. Hosts with low optical flux densities tend to inhabit more tenuous environments than hosts with high optical flux densities. By requiring spectroscopic redshifts, we have discarded RGs with only photometric redshifts, whose hosts almost exclusively have F_ν < 1 nJy. As a result, our estimates of 1 + Δ_RG, obs, 1 + Δ_RG, and 1 + Δ_RG | L_ν = l_ν are biased high. The severity of the bias depends on the percentage of RGs that have been selected out. As discussed in Sect. 3.2, 23% of all LoTSS DR1 RGs within the constrained BORG SDSS volume lack a spectroscopic redshift. Our sample of luminous giants is much less affected: only 7% of the population within the constrained BORG SDSS volume lack spectroscopy. Thus, correcting for spectroscopic redshift selection would shift the Cosmic Web density distribution for general RGs to somewhat lower densities, while it would not appreciably shift the corresponding distribution for luminous giants. The result would be an even bigger discrepancy between the Cosmic Web density distributions of luminous giants and general RGs than established in this work, strengthening our results.

Fig. 17. Observations underpinning the RG radio luminosity–Cosmic Web density relation of Fig. 8, but with Pan-STARRS r-band Kron flux densities Fν of the host galaxies indicated. We show uncertainty ellipses whose semi-major and semi-minor axes represent one standard deviation along each dimension. In the top panel, we show KDE Cosmic Web density distributions for Fν < 1 nJy, and for Fν ≥ 1 nJy. These suggest that spectroscopic redshift selection biases the observed Cosmic Web density distributions high.

Spectroscopic redshift selection also leaves an imprint on the sample’s radio luminosities. The top panel of Fig. 18 shows that the radio luminosities of the 1870 LoTSS DR1 RGs in the constrained BORG SDSS volume are distributed differently between the subpopulations with and without spectroscopic redshifts. In particular, the subpopulation with spectroscopic redshifts is more spread out over radio luminosity compared to the subpopulation without: both very low and very high values are more probable¹⁷. To understand why, we turn to the bottom panel of Fig. 18. The LoTSS 6″ noise level, the inverse square law, and cosmological surface brightness dimming together determine an RG radio luminosity detection threshold that increases with redshift. This threshold induces the hatched region, which is inaccessible to LoTSS-like radio observations. As a straightforward consequence, RGs with low radio luminosities can be detected only at low redshifts. It is not immediately clear whether the host galaxies of low–radio luminosity RGs at low redshifts have low or high optical flux densities. On the one hand, low radio luminosities suggest low optical luminosities: both forms of spectral luminosity correlate positively with stellar mass (as evinced by, for example, Fig. 2 of Mahajan et al. 2018; Sabater et al. 2019) and therefore¹⁸ positively with eachother. However, one must again be cautious of selection effects, which can inject strong but spurious luminosity–luminosity correlations in observations (Singal et al. 2019). On the other hand, observing low-redshift galaxies is easier, again because of the inverse square law and cosmological surface brightness dimming. The bottom panel of Fig. 18, and in particular its lower left corner, shows that the correlation between radio and optical luminosity is not strong enough to deny the observational advantages of low redshifts: low–radio luminosity RGs at low redshifts have high optical flux densities. As a result, spectroscopic redshifts are more common for such RGs. Similarly, the correlation between radio and optical luminosity explains that spectroscopic redshifts are more common for high–radio luminosity RGs.

Fig. 18. Overview of radio luminosity selection resulting from spectroscopic redshift selection, radio surface brightness limitations, and the correlation between radio and optical luminosity. The top panel compares radio luminosity distributions at νref = 150 MHz of LoTSS DR1 RGs in the constrained part of the BORG SDSS, distinguishing between those with spectroscopic redshifts (green) and those without (grey) (We perform KDE on samples of log10(Lν(νref)⋅(W Hz−1)−1) rather than of Lν(νref), and thus the probability densities on the vertical axis are dimensionless.). Radio galaxies with spectroscopic redshifts dominate (77%) the population. The bottom panel shows the two subpopulations in radio luminosity–redshift space, with RGs without spectroscopic redshifts marked (grey circles). We show uncertainty ellipses whose semi-major and semi-minor axes represent half a standard deviation along each dimension. One retrieves the top panel by collapsing the data on the horizontal axis.

This indirect radio luminosity selection effect is not expected to change the radio luminosity–Cosmic Web density relation inferred in this work and shown in Figs. 8 and 9 (as well as in Figs. D.1 and D.2). This is because the inferred relation concerns the Cosmic Web density distribution at given radio luminosity. Changes in the number of RGs available in each decade of radio luminosity – and in particular an increase in the number of RGs with 10²² W Hz⁻¹ < L_ν(ν = 150 MHz) < 10²³ W Hz⁻¹ – do not materially affect our analysis, as we already performed data rebalancing that downweights our sample’s most populated decades¹⁹.

5.5. Future outlook

Excitingly, there appear to be many opportuntities to improve our understanding of the large-scale environments of radio galaxies through Cosmic Web reconstructions, both with current and near-future data.

In this work, we have used LoTSS DR1 Local Universe RLAGN with spectroscopic redshifts from Hardcastle et al. (2019). With the release of the LoTSS DR2 radio data (Shimwell et al. 2022) and the associated optical crossmatching results (Hardcastle et al. 2023), it will become possible to redo our analysis on a footprint that is more than 13 times larger. Doing so would increase the size of the general RG sample by an order of magnitude, from ∼10³ to ∼10⁴. In turn, this would allow us to revisit Sect. 4.2’s RG radio luminosity–Cosmic Web density relation, but with all of Fig. 8’s six radio luminosity decades fully evenly weighted. If the positive correlation between radio luminosity and Cosmic Web density is genuine, it will emerge more clearly upon weighting evenly. Such a LoTSS DR2 sample might even be large enough to meaningfully investigate possible differences in the radio luminosity–Cosmic Web density relation for different types of RGs, such as those generated by quasars, and those generated by other AGN. In the near future, the WEAVE–LOFAR Survey (Smith et al. 2016) will generate ∼10⁶ spectra of LOFAR-selected galaxies, which will help alleviate the effects of spectroscopic redshift selection described in Sect. 5.4. Finally, the availability of optical spectra for all AGN in the current analysis and for those in similar follow-ups means that a comprehensive classification into HERGs and LERGs appears possible. In turn, this classification would allow for studies into the relationships between AGN accretion mode, RG growth, and the surrounding Cosmic Web. As these AGN reside in the Local Universe, many could be further characterised by crossmatching with X-ray catalogues.

Another promising extension of our work would be to use Cosmic Web reconstructions that are complementary to or improve upon the BORG SDSS. Notably, the BORG 2M++ (Jasche & Lavaux 2019) offers one such possibility: it provides density reconstructions for the entire sky and at a higher spatial resolution, and features more accurate selection functions, bias modelling, and gravitational dynamics than the BORG SDSS. However, its reconstructions are limited to z_max = 0.1, whereas those of the BORG SDSS extend to z_max = 0.16. Using the BORG SDSS and the BORG 2M++ in conjunction would allow for an extension of the number of giants for which a density can be determined: the catalogue of giants by Oei et al. (2023a), introduced in Sect. 2, contains 83 giants with spectroscopic redshifts z_s < 0.1 that lie outside of the constrained BORG SDSS volume. Adding these to the 260 giants for which we have already determined densities would increase the number of BORG giants by 32%. Perhaps more importantly, there are 89 giants for which both BORG SDSS and BORG 2M++ density estimates are possible. Comparing these densities allows for a quality assessment of the BORG SDSS density measurements at large. In particular, given its better gravity solver, we expect BORG 2M++ cluster densities to perform much better in Sect. 5.1’s Cosmic Web density–cluster mass test. In addition to the BORG 2M++, additional BORG runs are on their way. Besides, there exist comparable Cosmic Web reconstruction frameworks, such as COSMIC BIRTH (e.g., Kitaura et al. 2021). In any case, thanks to the combination of new spectroscopic galaxy surveys and advances in computing power, high-fidelity Bayesian Cosmic Web reconstructions will become available for an increasingly large fraction of the Local Universe and beyond. We expect that Cosmic Web density 1 + δ will therefore become a standard RG observable, alongside properties such as projected proper length l_p and radio luminosity l_ν.

Through the work of Leclercq et al. (2015), we have derived probabilistic Cosmic Web environment classifications, based on the T-web definition, for 260 luminous giants and 1443 LoTSS DR1 RGs. However, Leclercq et al. (2015) provides additional environmental classifications for the BORG SDSS volume that allow us to measure and understand the Cosmic Web environments of giants and other RGs in complementary ways. Future work could explore and compare RG environments in the T-web, DIVA, ORIGAMI, and LICH sense (Leclercq et al. 2016, 2017).

6. Conclusion

Using Bayesian Cosmic Web reconstructions of the Local Universe, we compared the large-scale environments of giant radio galaxies with those of the radio galaxy population in general. In particular, we measured multi-megaparsec-scale Cosmic Web densities around the hosts of ∼10² luminous giants and ∼10³ general RGs. This reveals that the currently observable population of giants inhabits denser regions of the Cosmic Web than general RGs – contradicting the popular hypothesis that giants primarily form in the dilute Cosmic Web. Our interpretation is that high jet powers, as implied by the high radio luminosities of the known population of giants, enable these RGs to overcome the IGM’s resistance to their megaparsec-scale growth.

Next, we quantified, for the first time, the relation between radio luminosity, a proxy for jet power, and Cosmic Web density. Our radio luminosity–conditioned Cosmic Web density distributions reinforce the idea that AGN that generate weak jets primarily reside in the dilute Cosmic Web, while AGN that generate powerful jets primarily reside in the dense Cosmic Web. We subsequently showed that the radio luminosity–Cosmic Web density relation, which we have inferred using our sample of general RGs only, can accurately explain the density distribution observed for luminous giants.

Additional evidence corroborates that luminous giants form more often in the dense Cosmic Web than RGs in general. Classifying Cosmic Web environments as clusters, filaments, sheets, or voids on the basis of the local gravitational field, we find that luminous giants occur more often in clusters than RGs in general. We also find evidence that clusters with giant-generating BCGs are more massive than other clusters; the evidence for the analogous claim for clusters with RG-generating BCGs is weaker. We now summarise our main findings in more detail.

1. Using the BORG SDSS, we have performed the most physically principled measurements yet of the Cosmic Web environments of giant radio galaxies and radio galaxies in general. Our characterisations require spectroscopic redshifts, and are currently confined to the Local Universe (z < z_max ≔ 0.16). In particular, we have determined 2.9 Mpc h⁻¹–scale total matter densities and gravitational environment classifications for 260 luminous giants, of which 208 (80%) have recently been discovered through the LoTSS DR2, and for 1443 LoTSS DR1 RGs.

2. While the marginal probability distributions of the Cosmic Web density around individual RGs are approximately lognormal, the probability distribution of the Cosmic Web density for the RG population as a whole is approximately gamma. We provide a physical argument that motivates why the gamma distribution could have arisen. The probability distribution of the Cosmic Web density for luminous giants also appears roughly gamma, but favours higher densities. Heteroskedastic measurement errors, which are a general feature of our densities, spuriously shift the population distributions towards lower densities. We demonstrate a forward modelling method to correct for heteroskedasticity.

3. We quantified the radio luminosity–Cosmic Web density relation for general RGs. We find that the 2.9 Mpc h⁻¹–scale Cosmic Web density distribution at a given 150 MHz radio luminosity is well described by 1 + Δ_RG | L_ν = l_ν ∼ Γ(k, θ), where k = 4.8 + 0.2 · √, θ = 1.4 + 0.02 · √, and √:= log₁₀(l_ν · (10²³ W Hz⁻¹)⁻¹). This shows that more luminous RGs tend to live in denser regions of the Cosmic Web. Treating giants as ordinary RGs with l_ν ≥ 10²⁵ W Hz⁻¹, we used this relation to predict the Cosmic Web density distribution for luminous giants. The prediction is consistent with the observed distribution for giants: a Kolmogorov–Smirnov test yields a p-value of 0.7. Whether less luminous giants – assuming they exist – also form more often in the dense Cosmic Web, remains to be seen when surveys with increased surface brightness sensitivity commence. If such giants obey the same radio luminosity–Cosmic Web density relation as other RGs, this will not be the case.

4. We show that our methodology enables the inference of radio galaxy number densities as a function of Cosmic Web density. Whether or not the RG number density is a strong function of Cosmic Web density, depends on the radio luminosity considered. In clusters, the number densities of high-luminosity RGs (e.g., 10²⁵–10²⁷ W Hz⁻¹) are ∼10² times higher than on average; in the same environment type, the number densities of low-luminosity RGs (e.g., 10²¹–10²³ W Hz⁻¹) are just ∼10¹ times higher than on average. Furthermore, we obtain tentative evidence that the number densities of low-luminosity RGs peak in low-mass clusters, before decreasing again. A possible explanation could be that, as cluster mass increases, the underlying population of active galaxies starts favouring the generation of high-luminosity RGs over low-luminosity RGs.

5. We used the BORG SDSS T-web classification to generate probability distributions over four Cosmic Web structure types – clusters, filaments, sheets, and voids – for both luminous giants and general RGs. Luminous giants inhabit clusters more often (${23}_{-3}^{+3}\%$ versus ${16}_{-1}^{+1}\%$). Independently, by crossmatching our RGs with a cluster catalogue, we find that the former percentage presumably is an underestimate, while confirming the latter percentage. The same crossmatching procedure reveals tentative evidence (p = 5%) that clusters with giant-generating BCGs are more massive than other clusters. At the same time, evidence that clusters with RG-generating BCGs are more massive than other clusters is weaker (p = 7%), even though the sample involved is almost twice as big.

Although RG environment characterisations with Bayesian Cosmic Web reconstructions are still far from perfect, we have demonstrated that the current generation of reconstructions already allows one to address open questions in radio galaxy research. Given the pitfalls of the BORG SDSS, every measured Cosmic Web density should be treated with skepsis, and we warn in particular that the densities around cluster RGs are strongly underestimated. Nevertheless, given the ever-increasing number of galaxies with spectroscopic redshifts, a growing list of technical refinements, and the widespread availability of computing power, exciting contemporary and near-future opportunities exist to expand and improve upon the results of this work. It is well possible that Cosmic Web density will soon become a staple quantity that helps to unravel the physics of radio galaxies.

Acknowledgments

M.S.S.L. Oei, R.J. van Weeren, and A. Botteon acknowledge support from the VIDI research programme with project number 639.042.729, which is financed by the Dutch Research Council (NWO). A. Botteon acknowledges financial support from the European Union – Next Generation EU. M.S.S.L. Oei warmly thanks Rafaël Mostert, Erik Osinga, and Wendy Williams for helpful discussions and the sharing of data. 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. Oei et al. (2023a)’s LoTSS DR2 GRG data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/672/A163. The LoTSS DR1 RG data with optical identifications from the value-added catalogue of Williams et al. (2019) are available at https://lofar-surveys.org/dr1_release.html; those interested in Hardcastle et al. (2019)’s RLAGN subsample should contact Martin J. Hardcastle. The BORG SDSS data release is publicly available at https://github.com/florent-leclercq/borg_sdss_data_release. Individual BORG SDSS HMC Markov chain samples are not available in the public domain; those interested should contact Florent Leclercq. The cluster catalogue of Wen & Han (2015) is available at the CDS.

References

Appendix A: Cosmic Web localisation accuracy with spectroscopic and photometric redshifts

In this appendix, we show that RGs with spectroscopic redshifts can be reliably localised within Cosmic Web reconstructions, while RGs with only photometric redshifts generally cannot.

To determine the accuracy of our radio galaxy localisations in the Cosmic Web, we Monte Carlo simulated a radial comoving distance distribution for each RG and compared its dispersion to the BORG SDSS voxel length. Each RG has a peculiar velocity v_p with respect to us. We treated v_p as a zero-mean Gaussian random variable (RV): v_p ∼ 𝒩(0, σ_vp), and chose a standard deviation σ_vp representative of conditions in low-mass galaxy clusters: σ_vp = 100 km s⁻¹. Similarly, we treated the measured redshift z_m as an RV, and again assumed Gaussianity: z_m ∼ 𝒩(μ_zm,σ_zm). Our catalogue provides the parameters μ_zm and σ_zm for each RG. Using Eqs. A.1 from left to right, we calculated the relative peculiar velocity RV β_p, the peculiar velocity redshift RV z_p, and the cosmological redshift RV z_c:

$$\begin{array}{c}\hfill {\beta }_{\mathrm{p}}\frac{v_{\mathrm{p}}}{c};z_{\mathrm{p}}=\sqrt{\frac{1+{\beta }_{\mathrm{p}}}{1-{\beta }_{\mathrm{p}}}}-1;z_{\mathrm{c}}=\frac{1+z_{\mathrm{m}}}{1+z_{\mathrm{p}}}-1.\end{array}$$(A.1)

Finally, we calculated the radial comoving distance RV r = r(z_c) = ||r||₂. The distributions of r are approximately Gaussian: for each of the 248 giants in the BORG SDSS volume with a spectroscopic redshift and an associated error, we performed the Shapiro–Wilk test on 1000 Monte Carlo samples and rejected the null hypothesis that r is Gaussian for one at significance level α = 0.01.²⁰ For these giants, 95% of the standard deviations σ_r are between 1.00 Mpc h⁻¹ and 1.12 Mpc h⁻¹; the sample median is 1.06 Mpc h⁻¹. By contrast, the 21 giants in the BORG SDSS volume with only a photometric redshift and an associated error have a σ_r that ranges from 20 Mpc h⁻¹ to 110 Mpc h⁻¹; for details, see Table A.1. To localise a giant in the BORG SDSS, precision up to the scale of a voxel, which are approximately 2.9 Mpc h⁻¹ long, is necessary – but sub-voxel localisation is redundant. We conclude that giants with a spectroscopic redshift can be tied to an individual voxel, even when taking peculiar velocity into account; on the other hand, giants with only a photometric redshift can be mislocalised more than a typical filament length, and are therefore not subjectable to a precise environment analysis. Of course, these conclusions extend to any radio galaxy with an associated host, or indeed to any galaxy.²¹

As only few (i.e. ∼10²) giants are currently known in the volume reconstructed by the BORG SDSS, it is worthwhile to perform host galaxy spectroscopic follow-up for those with photometric redshifts only. We provide a list with targets in Table A.1.

Appendix B: Cosmic Web density distribution: the gamma Ansatz

In this work, we have modelled several RG Cosmic Web density RVs, such as 1 + Δ_RG, obs, 1 + Δ_RG, and 1 + Δ_RG | L_ν = l_ν, as gamma variates. This choice is largely driven by data (and a preference for simplicity) rather than by theory, although the requirements of continuity and a strictly positive support have a physical basis. Gamma distributions adequately fit 1 + Δ_RG, obs (see the bottom row of Fig. 6), and their use in modelling 1 + Δ_RG | L_ν = l_ν is justified by the RG prediction of 1 + Δ_GRG, obs, which matches observations (see Fig. 11). We post hoc theorised why the gamma distribution arises in the current context.

The late-time density in a voxel is, of course, equal to the late-time mass in the voxel divided by its (fixed) volume of (2.9 Mpc h⁻¹)³. Therefore, up to a dimensionful constant, the late-time density and the late-time mass have the same probability distribution. The late-time mass M in the voxel equals the sum of the mass aggregated over cosmic time. If we are to consider a simplistic treatment of mass aggregation, or structure formation, we should pick the most relevant large-scale structure regime. The central pie chart of Fig. 15 suggests that most observed RGs have filament environments. In addition, Pasini et al. (2021) suggest that within filaments, most RGs inhabit galaxy groups, which grow by merging with other groups. Therefore, we considered a proto-filament – a massive structure that extends essentially along a single dimension – on which we considered a Poisson point process with constant spatial rate. The Poisson points represent peaks in the density field where galaxy groups arise. As with any homogeneous Poisson process, the distances between the points are exponentially distributed. It is now helpful to partition the filament in cells, each associated to a single Poisson point. Assuming that the filament has a constant cross-sectional area and an approximately constant density, the volumes and masses of the cells are simply proportional to the distances between the points, and are thus also exponentially distributed. Finally, if a fixed fraction of the mass in each cell collapses into the cell’s galaxy group, then the group masses within the filament – at early times – are, again, exponentially distributed. A group that exists at late times has built up its mass M by aggregating the masses of early-time groups. Assuming a constant temporal rate of early-time group aggregation, there exists a typical number N of early-time groups that contribute to the mass of a given late-time group: $M={\displaystyle \sum _{i=1}^N}{M}_i$. We then invoked the fact that the sum of a fixed number of independent and identically distributed RVs with an exponential distribution is an RV with a gamma distribution. Thus, as the M_i are exponentially distributed, M is gamma distributed. This, in turn, implies that the late-time densities of the voxels in which late-time groups fall, are gamma distributed also. As late-time RGs trace late-time groups, the late-time RG density distribution should be approximately gamma.

In the future, by crossmatching RG catalogues with galaxy group catalogues with accurate masses, the mass distribution of RG-hosting groups in the Local Universe can be sampled and tested against the gamma distribution hypothesis.

Appendix C: Modelling relative baryon density measurement heteroskedasticity

We must describe the measured RG relative baryon density given a true RG relative baryon density with a continuous probability distribution that has the positive half-line as its support. Preferably, one uses a simple analytical prescription.

Motivated by the quality of the fits shown in Fig. 5, we assumed 1 + Δ_RG, obs | 1 + Δ_RG = 1 + δ ∼ Lognormal(μ, σ²), and likewise for giants. The parameters μ and σ² are determined by assumptions on the RV’s mean and variance. In particular, we assumed that measured RG relative densities are unbiased estimators of the underlying true RG relative densities and that their variances are the same when the underlying true RG relative densities are the same:

$$\begin{array}{cc}\hfill \mathbb{E}[1+{\mathrm{\Delta }}_{\mathrm{RG},\mathrm{obs}}|1+{\mathrm{\Delta }}_{\mathrm{RG}}=1+\delta ]& =1+\delta ;\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill \mathbb{V}[1+{\mathrm{\Delta }}_{\mathrm{RG},\mathrm{obs}}|1+{\mathrm{\Delta }}_{\mathrm{RG}}=1+\delta ]& =g(1+\delta ),\hfill \end{array}$$(C.2)

where g(1+δ) is some non-negative function characterising the heteroskedasticity. For the current work, we adopted a two-parameter power-law variance model g(1+δ) ≔ a(1+δ)^b. From basic identities of the lognormal distribution, we found for this choice of g

$$\begin{array}{cc}\hfill \mu & =ln(1+\delta )-\frac{1}{2}ln(1+a{(1+\delta )}^{b-2});\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill {\sigma }^2& =ln(1+a{(1+\delta )}^{b-2}).\hfill \end{array}$$(C.4)

Using the fact that

$$\begin{array}{cc}& f_{1+{\mathrm{\Delta }}_{\mathrm{RG},\mathrm{obs}}}(1+\delta )=\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }}{f}_{1+{\mathrm{\Delta }}_{\mathrm{RG},\mathrm{obs}}|1+{\mathrm{\Delta }}_{\mathrm{RG}}=x}(1+\delta )·f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}}\left(x\right)\mathrm{d}x,\hfill \end{array}$$(C.5)

we could deduce the distribution of 1 + Δ_RG given some parametrisation of it and a procedure like MLE. We chose 1 + Δ_RG ∼ Γ(k,θ) on the basis that 1 + Δ_RG, obs appears to be much better described by a gamma distribution than by a lognormal distribution and that the effect of heteroskedasticity appears minor. We determined a and b from data.²²

Appendix D: Radio luminosity–Cosmic Web density relation: Fixed voxel method

In Sect. 4.2, we have presented a quantification of the relation between RG radio luminosity and Cosmic Web density. In particular, we inferred a probability distribution for the 2.9 Mpc h⁻¹–scale total matter density around the hosts of Local Universe radio galaxies with a given 150 MHz radio luminosity – which can range from 10²¹ to 10²⁷ W Hz⁻¹. The MLE model shown in Fig. 8 and the one- and two-dimensional posterior marginals shown in Fig. 9, are based on densities measured through Sect. 3.1.2’s flexible voxel method. In this appendix, through Figs. D.1 and D.2, we present the analogous results based on densities measured through Sect. 3.1.1’s fixed voxel method. Both methods reveal a positive scaling between RG radio luminosity and Cosmic Web density.

Fig. D.1. PDFs of the observed RG relative total matter density RV given a radio luminosity at νref = 150 MHz, using MLE parameter values for the model described in Sect. 4.2. The black contours denote CDF values. We overplot all 1443 selected LoTSS DR1 RGs (dots), with those above the empirical median density coloured red, and those below coloured blue. We used fixed voxel method densities. For flexible voxel method densities, see Fig. 8.

Fig. D.2. Posterior distribution over k(0), θ(0), k′(0), and θ′(0), based on selected Local Universe LoTSS DR1 RGs. We show all two-parameter marginals of the likelihood function, with contours enclosing 50% and 90% of total probability. We mark the MLE (grey dot) and the MAP (grey cross). The one-parameter marginals again show the MLE (dash-dotted line), a mean-centred interval of standard deviation–sized half-width (hashed region), and a median-centred 90% credible interval (shaded region). We used fixed voxel method densities. For flexible voxel method densities, see Fig. 9.

Appendix E: Relative number density derivation

We let the RV 1 + Δ_CW represent the 2.9 Mpc h⁻¹–scale relative density at a randomly chosen point in the contemporary Cosmic Web, with f_1 + ΔCW being its PDF. Similarly, we let the RV 1 + Δ_RG represent the 2.9 Mpc h⁻¹–scale relative density at a randomly chosen RG in the contemporary Cosmic Web, with f_1 + ΔRG being its PDF. We considered a cosmologically sized, comoving volume of extent V, in which a total of N_RG radio galaxies exist. The subvolume in which the relative density is between 1 + δ and 1 + δ + dδ has extent

$$\begin{array}{c}\hfill \mathrm{d}V=V{f}_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )\mathrm{d}\delta .\end{array}$$(E.1)

The number of RGs in this subvolume is

$$\begin{array}{c}\hfill \mathrm{d}{N}_{\mathrm{RG}}=n_{\mathrm{RG}}(1+\delta )\mathrm{d}V,\end{array}$$(E.2)

where n_RG(1 + δ) is the number density of RGs that would arise in an environment of constant density 1 + δ. In practice, one never encounters environments of constant density; the average RG number density in the part of the Cosmic Web where the relative density is between 1 + δ₁ and 1 + δ₂ equals

$$\begin{array}{c}\hfill {\overline{n}}_{\mathrm{RG}}(1+{\delta }_1,1+{\delta }_2)=\frac{{\int }_{{\delta }_1}^{{\delta }_2}{n}_{\mathrm{RG}}(1+\delta )f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )\mathrm{d}\delta }{{\int }_{{\delta }_1}^{{\delta }_2}{f}_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )\mathrm{d}\delta }.\end{array}$$(E.3)

As a result, the cosmic mean RG number density ${\overline{n}}_{\mathrm{RG}}$ is given by

$$\begin{array}{c}\hfill {\overline{n}}_{\mathrm{RG}}{\displaystyle {\int }_{-1}^{\infty }}{n}_{\mathrm{RG}}(1+\delta )f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )\mathrm{d}\delta .\end{array}$$(E.4)

The probability that an RG has a relative density between 1 + δ and 1 + δ + dδ, f_1 + ΔRG(1 + δ) dδ, is

$$\begin{array}{c}\hfill f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}}(1+\delta )\mathrm{d}\delta =\frac{\mathrm{d}{N}_{\mathrm{RG}}}{N_{\mathrm{RG}}}=\frac{n_{\mathrm{RG}}(1+\delta )V{f}_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )\mathrm{d}\delta }{N_{\mathrm{RG}}}.\end{array}$$(E.5)

Combining $N_{\mathrm{RG}}={\overline{n}}_{\mathrm{RG}}V$ with Eq. E.5 leads to

$$\begin{array}{c}\hfill \frac{n_{\mathrm{RG}}(1+\delta )}{{\overline{n}}_{\mathrm{RG}}}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(E.6)

In other words, through point-wise division of the PDFs of 1 + Δ_RG and 1 + Δ_CW, both of which we determined in this work, we obtained the RG number density that would arise at a given density, up to a constant. This constant is (the reciprocal of) the cosmic mean RG number density; the LHS of Eq. E.6 can thus be interpreted as the relative RG number density at 1 + δ – that is, the number density relative to the cosmic mean value. Fully analogously, the relative GRG number density at 1 + δ is

$$\begin{array}{c}\hfill \frac{n_{\mathrm{GRG}}(1+\delta )}{{\overline{n}}_{\mathrm{GRG}}}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{GRG}}}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(E.7)

By multiplying both sides of Eq. E.6 by f_{Lν|1 + ΔRG = 1 + δ}(l_ν), using the identity f_X|Y = y(x)f_Y(y) = f_Y|X = x(y)f_X(x), and dividing both sides by f_Lν(l_ν), we obtained

$$\begin{array}{c}\hfill \frac{n_{\mathrm{RG}}(1+\delta )·f_{L_{\nu }|1+{\mathrm{\Delta }}_{\mathrm{RG}}=1+\delta }(l_{\nu })}{{\overline{n}}_{\mathrm{RG}}·f_{L_{\nu }}(l_{\nu })}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(E.8)

The numerator and denominator at the left-hand side are both number density densities: they denote radio galaxy numbers per unit of comoving volume and unit of radio luminosity. As their physical dimensions are distinct from those of ordinary number densities, we denoted them differently. Calling the numerator _RG(l_ν, 1 + δ) and the denominator ${\overline{}}_{\mathrm{RG}}(l_{\nu })$, we wrote²³

$$\begin{array}{c}\hfill \frac{{}_{\mathrm{RG}}(l_{\nu },1+\delta )}{{\overline{}}_{\mathrm{RG}}(l_{\nu })}=\frac{f_{1+{\mathrm{\Delta }}_{\mathrm{RG}}|L_{\nu }=l_{\nu }}(1+\delta )}{f_{1+{\mathrm{\Delta }}_{\mathrm{CW}}}(1+\delta )}.\end{array}$$(E.9)

Appendix F: Disparity between mean squared density and squared mean density

Several physical quantities that relate to AGN and RGs, such as the bremsstrahlung emissivity around the host galaxy, scale with the square of IGM density ρ. In Sect. 5.3, we argue qualitatively that such quantities are likely to be underestimated if calculated using low-resolution densities ⟨ρ⟩, such as the ones offered by the BORG SDSS. In this appendix, we explicitly calculate how much bigger ⟨ρ²⟩ is than ⟨ρ⟩² – in the sense of a multiplicative Jensen’s gap. Our calculations demonstrate a strong dependence on Cosmic Web environment.

When the true density field varies over scales much larger than a voxel side length L, such as in sheets or voids, the density field variation within a voxel resembles a gradient from (1−ϵ)⟨ρ⟩ to (1+ϵ)⟨ρ⟩, where ϵ ∈ [0, 1]. If the planes of constant density are parallel to two of the voxel’s faces, one finds

$$\begin{array}{c}\hfill \frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}=\frac{{ϵ}^2}{3}+1,\text{so that}1\le \frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}\le \frac{4}{3}.\end{array}$$(F.1)

We show the multiplicative Jensen’s gap of Eq. F.1 in Fig. F.1 (blue curve).

Fig. F.1. Multiplicative discrepancy between the mean squared density ⟨ρ2⟩ and the squared mean density ⟨ρ⟩2 in voxelised Cosmic Web reconstructions, for density profiles of clusters (red), filaments (orange), and sheets or voids (blue). The length scale R ∼ L, where L is the voxel side length; here, L = 2.9 Mpc h−1. The asymptote (grey) corresponds to the ratio $\frac{4}{3}$.

For a filament whose baryon density is modelled with an isothermal single-β model (Cavaliere & Fusco-Femiano 1976, 1978) with $\beta =\frac{2}{3}$ and core radius r_c, one finds

$$\begin{array}{c}\hfill \frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}=\frac{{\left(\frac{R}{r_{\mathrm{c}}}\right)}^4}{(1+{\left(\frac{R}{r_{\mathrm{c}}}\right)}^2){ln}^2(1+{\left(\frac{R}{r_{\mathrm{c}}}\right)}^2)},\end{array}$$(F.2)

where we considered a cylinder with length L and radius R. Here, R is such that the area of a cylindrical section perpendicular to the axis equals L², the area of a voxel face. Thus, $R=\frac{1}{\sqrt{\pi }}L$. We show the multiplicative Jensen’s gap of Eq. F.2 in Fig. F.1 (orange curve).

When the true density field varies over scales much smaller than a voxel side length L, such as in clusters, the ratio between the mean of the squared density and the square of the mean density can be much larger than 1. In particular, for a Navarro–Frenk–White (NFW) profile (Navarro et al. 1996) with scale radius r_s,

$$\begin{array}{c}\hfill \frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}=\frac{1-{\left(\frac{r_{\mathrm{s}}}{r_{\mathrm{s}}+R}\right)}^3}{9{\left(\frac{r_{\mathrm{s}}}{R}\right)}^3{(ln\frac{r_{\mathrm{s}}}{r_{\mathrm{s}}+R}+1-\frac{r_{\mathrm{s}}}{r_{\mathrm{s}}+R})}^2}.\end{array}$$(F.3)

where we considered a sphere with radius R. Here, R is such that the volume of the sphere is equal to the volume of the voxel: $R={\left(\frac{3}{4\pi }\right)}^{\frac{1}{3}}L$. We show the multiplicative Jensen’s gap of Eq. F.3 in Fig. F.1 (red curve). In the limit $\frac{r_{\mathrm{s}}}{R}\to \infty$, $\frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}\to \frac{4}{3}$.

All Tables

All Figures

Fig. 1. Mollweide view of the sky showing the locations of all giants and LoTSS DR1 RGs in the Local Universe for which we inferred Cosmic Web densities and dynamical states. The background shows the Milky Way at 150 MHz (Zheng et al. 2017), on which we overlaid the LoTSS DR1 footprint (hatched dark rectangle). All RGs are drawn as lemniscates of Bernoulli; we did not attempt to portray realistic shapes or position angles. The colours represent redshifts z ∈ (0, zmax ≔ 0.16), whilst the diameters are proportional to projected proper lengths lp ∈ ( ∼ 1 kpc, 4.6 Mpc). Giants are translucent. Upon zooming in, the reader can appreciate the wide variety of observed radio galaxy sizes.

In the text

Fig. 2. LoTSS DR2 cutouts at central observing frequency νobs = 144 MHz and resolution θFWHM = 6″, centred around giant-generating BCGs of Local Universe clusters. All giants shown are discoveries of Oei et al. (2023a) for which no previous images have been published. Each cutout covers a solid angle of 15′×15′. Contours signify 3, 5, and 10 sigma-clipped standard deviations above the sigma-clipped median. For scale, we show the stellar Milky Way disk (with a diameter of 50 kpc) generated using the Ringermacher & Mead (2009) formula, alongside a 3 times inflated version. Each DESI Legacy Imaging Surveys DR9 (g, r, z) inset shows the central 3′×3′ region.

In the text

Fig. 3. Example localisations of giant radio galaxies within the large-scale structure of the Local Universe. The top row shows two cluster giants, the middle row shows two filament giants, and the bottom row shows two sheet giants. For each giant, we show a slice of constant Cartesian comoving z through the late-time BORG SDSS posterior mean total matter density field and a LoTSS DR2 6″ image at νobs = 144 MHz (inset). Outside of the SDSS DR7–constrained volume, the posterior mean tends to the Universe’s late-time mean total matter density ${\overline{\rho }}_0$. The locations of the giants are marked by white circles.

In the text

	Fig. 4. Comparison between Cosmic Web densities of giants and LoTSS DR1 RGs in the Local Universe, inferred via two variations of our BORG SDSS–based approach. For each RG, we show the mean density measured through Sect. 3.1.1’s fixed voxel method (horizontal axis) and through Sect. 3.1.2’s flexible voxel method (vertical axis).
In the text

Fig. 5. Distributions for the measured relative total matter density RVs 1 + ΔGRG, obs | 1 + ΔGRG = 1 + δ and 1 + ΔRG, obs | 1 + ΔRG = 1 + δ for an individual giant (top row; solid), and for an individual RG (bottom row; hatched). Fixed voxel method densities (left column; green) are lower than flexible voxel method densities (right column; orange). The MLE-fitted lognormal PDFs (solid curves) demonstrate that the distributions are almost lognormal.

In the text

Fig. 6. Probability density functions (PDFs) of the relative total matter density RVs of 260 Local Universe giants (top row), and of 1443 Local Universe RGs (bottom row), determined through the fixed voxel method (left column) and flexible voxel method (right column). These RVs correspond to the density field smoothed to a scale of 2.9 Mpc h−1. We also show PDFs of gamma-distributed RVs with parameters obtained via MLE, both before (solid lines) and after (dash-dotted lines) a heteroskedasticity correction. We warn that these distributions are affected by surface brightness selection, and thus represent known populations only.

In the text

Fig. 7. Mean lobe surface brightness estimates Bν at νref = 150 MHz for LoTSS DR1 RGs at z < 0.2, shown as a function of projected proper length lp. During the majority (∼90%) of their lifetime, RGs grow while maintaining a roughly constant radio luminosity. If, additionally, growth is shape-preserving, surface brightness decreases quadratically with lp. The seven golden lines therefore denote approximate evolutionary tracks of RGs at z = 0 with end-of-life radio luminosities lν ∈ {1021 W Hz−1, …, 1027 W Hz−1}. The region of parameter space where RGs can become detectable giants, lν ≥ lν, min ≔ 1024 W Hz−1, is shaded gold. The region of parameter space occupied by giants, lp ≥ lp, GRG ≔ 0.7 Mpc, is shaded green. The region of parameter space inaccessible to the LoTSS, Bν < σIν = 25 Jy deg−2, is shaded red.

In the text

Fig. 8. PDFs of the observed RG relative total matter density RV given a radio luminosity at νref = 150 MHz, using MLE parameter values for the model described in Sect. 4.2. The black contours denote CDF values. We overplot all 1443 selected LoTSS DR1 RGs (dots), with those above the empirical median density coloured red, and those below coloured blue. We used flexible voxel method densities. For fixed voxel method densities, see Fig. D.1.

In the text

Fig. 9. Posterior distribution over k(0), θ(0), k′(0), and θ′(0), based on selected Local Universe LoTSS DR1 RGs. We show all two-parameter marginals of the likelihood function, with contours enclosing 50% and 90% of total probability. We mark the MLE (grey dot) and the maximum a posteriori (MAP; grey cross). The one-parameter marginals again show the MLE (dash-dotted line), a mean-centred interval of standard deviation–sized half-width (hashed region), and a median-centred 90% credible interval (shaded region). We used flexible voxel method densities. For fixed voxel method densities, see Fig. D.2.

In the text

Fig. 10. Monte Carlo simulation of a population of nine galaxies (thick stripes) with successive active and quiescent phases. We mark the active phases (thin stripes). Each galaxy spends a different fraction ϝa of the time in the active phase. We model the time spent in both phases as draws from exponential distributions: ta ∼ Exp(λa) and tq ∼ Exp(λq). We assume that all galaxies share the same quiescent phase rate parameter λq = 10 Gyr−1, although this may not be entirely correct (Turner & Shabala 2015). This assumption fixes the active phase rate parameter: ${\lambda }_{\mathrm{a}}=({ϝ}_{\mathrm{a}}^{-1}-1){\lambda }_{\mathrm{q}}$. On human time scales, we can access only a snapshot view of this population. During the present snapshot (vertical line), seven galaxies generate RGs, and could be included in an RG sample. RG ages at the snapshot (hatched parts of the thin green stripes) generally increase with ϝa.

In the text

	Fig. 11. Density distribution predicted for observed giants combining Sect. 4.2’s RG radio luminosity–Cosmic Web density relation and a radio luminosity cut-off lν, min = 1025 W Hz−1 (pink); in addition, we show the KDE density distributions for giants with gamma distribution MLE fit (solid green) and for general RGs (hatched green) – all as in the left column of Fig. 6.
In the text

Fig. 12. Number densities of RGs with a given radio luminosity lν, as a function of Cosmic Web density at the 2.9 Mpc h−1 scale. The number density functions reflect conditions in the Local Universe, and are given relative to the cosmic mean number density at lν. The solid curves reflect the posterior mean, whilst the shaded areas around the mean denote −1 to +1 posterior standard deviation ranges. We used flexible voxel method densities.

In the text

Fig. 13. Number density of luminous non-giant radio galaxies in the Local Universe, as a function of Cosmic Web density. We define these as RGs with 150 MHz radio luminosities lν ≥ 1024 W Hz−1, projected proper lengths lp < lp, GRG ≔ 0.7 Mpc, and redshifts z < zmax ≔ 0.16. The Cosmic Web densities encompass baryonic and dark matter, are defined on a 2.9 Mpc h−1 scale, and are given relative to the cosmic mean. We used flexible voxel method densities here. We show the posterior mean with a shaded −1 to +1 posterior standard deviation range. We also show five individual realisations (grey).

In the text

	Fig. 14. Bayesian Cosmic Web reconstructions allow one to determine, for each giant, a probability distribution p over the four T-web environment classes (clusters, filaments, sheets, and voids). Here we summarise these distributions: each panel shows how a component of p ≔ (p1, p2, p3, p4) is itself distributed.
In the text

Fig. 15. Dynamical environment classification of giant radio galaxies (left), radio galaxies in general (centre), and the Local Universe in its entirety (right; from Table 3 of Leclercq et al. 2015). We define clusters, filaments, sheets, and voids in the T-web sense. These distributions depend on the scale to which the Cosmic Web density field is smoothed; in this case, the smoothing scale is 2.9 Mpc h−1. If giants and RGs were scattered uniformly throughout the Cosmic Web, their distributions would be similar to that of the volume-filling fractions (VFFs). Instead, observed RGs – and observed giants in particular – favour cluster and filament environments.

In the text

Fig. 16. Comparison between environmental estimates for luminous giants (large dots) and LoTSS DR1 RGs (small dots) generated by BCGs in the Local Universe. We compare cluster masses M500 from Wen & Han (2015) with relative densities 1 + δ from the BORG SDSS, and colour the dots based on their cluster probability p1 (relative to the average cluster probability ${\overline{p}}_1$ for luminous giants and LoTSS DR1 RGs). For three randomly chosen giants and three randomly chosen RGs, we visualise uncertainties. Generally, as expected, RG-generating BCGs claimed to be in more massive clusters by Wen & Han (2015) have higher BORG SDSS densities and T-web cluster probabilities. The contours indicate fractions of the density expected through M500, and show that our relative densities fall short significantly in the cluster regime – by a factor that ranges from less than two to more than an order of magnitude. Densities of low-mass clusters are less biased than densities of high-mass clusters. We used flexible voxel method densities here.

In the text

Fig. 17. Observations underpinning the RG radio luminosity–Cosmic Web density relation of Fig. 8, but with Pan-STARRS r-band Kron flux densities Fν of the host galaxies indicated. We show uncertainty ellipses whose semi-major and semi-minor axes represent one standard deviation along each dimension. In the top panel, we show KDE Cosmic Web density distributions for Fν < 1 nJy, and for Fν ≥ 1 nJy. These suggest that spectroscopic redshift selection biases the observed Cosmic Web density distributions high.

In the text

Fig. 18. Overview of radio luminosity selection resulting from spectroscopic redshift selection, radio surface brightness limitations, and the correlation between radio and optical luminosity. The top panel compares radio luminosity distributions at νref = 150 MHz of LoTSS DR1 RGs in the constrained part of the BORG SDSS, distinguishing between those with spectroscopic redshifts (green) and those without (grey) (We perform KDE on samples of log10(Lν(νref)⋅(W Hz−1)−1) rather than of Lν(νref), and thus the probability densities on the vertical axis are dimensionless.). Radio galaxies with spectroscopic redshifts dominate (77%) the population. The bottom panel shows the two subpopulations in radio luminosity–redshift space, with RGs without spectroscopic redshifts marked (grey circles). We show uncertainty ellipses whose semi-major and semi-minor axes represent half a standard deviation along each dimension. One retrieves the top panel by collapsing the data on the horizontal axis.

In the text

Fig. D.1. PDFs of the observed RG relative total matter density RV given a radio luminosity at νref = 150 MHz, using MLE parameter values for the model described in Sect. 4.2. The black contours denote CDF values. We overplot all 1443 selected LoTSS DR1 RGs (dots), with those above the empirical median density coloured red, and those below coloured blue. We used fixed voxel method densities. For flexible voxel method densities, see Fig. 8.

In the text

Fig. D.2. Posterior distribution over k(0), θ(0), k′(0), and θ′(0), based on selected Local Universe LoTSS DR1 RGs. We show all two-parameter marginals of the likelihood function, with contours enclosing 50% and 90% of total probability. We mark the MLE (grey dot) and the MAP (grey cross). The one-parameter marginals again show the MLE (dash-dotted line), a mean-centred interval of standard deviation–sized half-width (hashed region), and a median-centred 90% credible interval (shaded region). We used fixed voxel method densities. For flexible voxel method densities, see Fig. 9.

In the text

	Fig. F.1. Multiplicative discrepancy between the mean squared density ⟨ρ2⟩ and the squared mean density ⟨ρ⟩2 in voxelised Cosmic Web reconstructions, for density profiles of clusters (red), filaments (orange), and sheets or voids (blue). The length scale R ∼ L, where L is the voxel side length; here, L = 2.9 Mpc h−1. The asymptote (grey) corresponds to the ratio $\frac{4}{3}$.
In the text