© The Authors 2023

1. Introduction

When gas, dust, and stars accrete onto a supermassive black hole (SMBH) in the centre of a galaxy, collimated jets arise along the Kerr rotation axis that blast some of the infalling material into the intergalactic medium (IGM) (e.g., Blandford & Rees 1974). In this process, the ejecta dissolve into a relativistic plasma that drags along a magnetic field and glows in synchrotron light. The resulting luminous structure is called a radio galaxy (RG); the central black hole that has generated it an active galactic nucleus (AGN).

It is increasingly clear that RGs and their AGN play an important role in galaxy evolution and cosmology. By heating the interstellar medium (ISM) or even expelling it from their host galaxies through galactic superwinds, AGN quench star formation (e.g., Di Matteo et al. 2005). Beckmann et al. (2017) have shown that AGN-induced star formation quenching is most pronounced in massive galaxies. There is also compelling evidence that the accompanying RGs provide the energy necessary to stop (e.g., McNamara & Nulsen 2012; Yang et al. 2019) bremsstrahlung-mediated cooling flows (Fabian et al. 1984) in clusters of galaxies. In the absence of cooling flows, the intra-cluster medium (ICM) remains dilute and hot, and galaxies in the centres of clusters are denied infalling gas that could otherwise reignite star formation. Cosmological simulations that include this RG feedback indeed resolve (e.g., Croton et al. 2006) the overprediction of baryonic masses and luminosities of central cluster galaxies that early simulations found. Finally, RGs may be responsible for magnetising the IGM that pervades the filaments of the Cosmic Web (e.g., Vazza et al. 2017).

Despite the emerging picture that RGs trace quenched star formation, inhibit cooling flows, and magnetise filaments, our knowledge of them is far from complete. Concerning geometry, a major unknown is the exact connection between the morphology of RGs and the pressure field of the ambient IGM, especially in filaments and cluster outskirts. Another question is whether small and large RGs come from the same initial population, or whether their growth is driven by distinct physical processes. Finally, we do not know how large can RGs become, and, more generally, how many RGs there are of each length.

To test RG growth models that answer these and other questions, it is imperative to study the subpopulation of most spatially extreme RGs: the giant radio galaxies (GRGs). The defining feature of giants is that their proper lengths – when projected onto the celestial sphere – exceed some threshold l_p, GRG, which is canonically chosen as 0.7 Mpc or 1 Mpc. If l_p, GRG = 0.7 Mpc, then the preceding literature describes a total of 1281 giants.

In recent years, several studies have successfully searched for giants in systematic, wide-area surveys such as the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) and the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017). A combination of manual (i.e. visual) and automated searches (Solovyov & Verkhodanov 2011, 2014; Amirkhanyan 2016; Proctor 2016; Dabhade et al. 2017, 2020a) in the NVSS yielded 313 new giants (24% of the aforementioned literature population). Meanwhile, Dabhade et al. (2020b) discovered 225 new giants (17% of this same population) in the LoTSS DR1 (Shimwell et al. 2019), whose survey footprint is 80 times smaller than NVSS’s. Such searches have the advantage of introducing almost homogeneous selection effects throughout the survey footprint, which can potentially be modelled and thus corrected for during any subsequent statistical inference.

In this work this idea comes to fruition, by conducting a precision analysis of the intrinsic giant radio galaxy length distribution. To do so, we require two ingredients. First, in Sect. 2, we develop a statistical framework that allows one to answer probabilistic questions regarding both large samples of giants and individual specimina. Then, in Sect. 3, we describe our LoTSS DR2 (Shimwell et al. 2022) GRG search campaign and the trove of previously unknown giants it has yielded; moreover, we describe the assemblage of the most complete GRG catalogue to date. In Sect. 4, combining theory and data, we infer the tail index parameter that describes the intrinsic GRG length distribution, which constrains future models and simulations aimed at understanding RG growth. In Sect. 5, we discuss caveats of the current work and give recommendations for future extensions, before we present conclusions in Sect. 6.

We assume a concordance inflationary Λ cold dark matter cosmology with parameters 𝔐 from Planck Collaboration VI (2020); that is to say 𝔐 = (h = 0.6766, Ω_BM, 0 = 0.0490, Ω_M, 0 = 0.3111, Ω_Λ, 0 = 0.6889), where H₀ := h × 100 km s⁻¹ Mpc⁻¹. We define the spectral index α such that it relates to flux density F_ν at frequency ν as F_ν ∝ ν^α, and define giants using threshold l_p, GRG := 0.7 Mpc. Regarding the terminology, we use ‘angular length’ where others use ‘largest angular size’ (LAS), and ‘projected proper length’ where others use ‘largest linear size’ (LLS)¹.

2. Theory

To measure the intrinsic GRG length distribution, we must first establish a suitable statistical framework. In this section, we provide a summary of the theory developed in Appendix A. Following Occam’s razor, we construct a model with minimal assumptions that provides new insight into the GRG phenomenon and the detection biases inherent to systematic search campaigns.

2.1. Intrinsic proper length

Firstly, we assume that giants and non-giant RGs share a common length distribution². In particular, because power laws are abundant in Nature, we assume that the intrinsic proper length random variable (RV) L has a Pareto distribution with tail index ξ < −1 and support from l_min > 0 onwards. If an RV is Pareto distributed, then the relative occurrence of two possible outcomes equals their ratio raised to a power: the tail index ξ. In astrophysics, Pareto distributions describe the kinetic energies of freshly accelerated electrons in large-scale structure and supernova shocks (e.g., Kirk & Schneider 1987), the initial masses of main-sequence stars (e.g., Kroupa 2001), and the luminosities of gamma-ray bursts (e.g., Bloom et al. 2001), to name a few examples. Previous works (e.g., Andernach et al. 2021) have already hinted at the approximate validity of a Pareto distribution description for GRG lengths. By comparing our final model to observations, as discussed in Sect. 4 and visualised in Fig. 14, we demonstrate that this assumption is indeed a powerful approximation in the current case.

The probability density function (PDF) f_L : ℝ → ℝ_≥0 thus becomes

$$\begin{array}{c}\hfill f_L\left(l\right)=\{\begin{array}{cc}0\hfill & l\le l_{\min },\hfill \\ \frac{-(\xi +1)}{l_{\min }}{\left(\frac{l}{l_{\min }}\right)}^{\xi }\hfill & l>l_{\min }.\hfill \end{array}\end{array}$$(1)

We refer the reader to Appendix A.1 for a derivation of this expression and a demonstration of its connection to the literature’s most common parametrisation.

2.2. Projected proper length

From the distribution of intrinsic lengths and the assumption of random radio galaxy orientations, we now derive the distribution of projected lengths. This distribution is more easily compared to observations, which usually lack inclination angle information.

2.2.1. Distribution for RGs

To model length and orientation, we consider a vector $\overrightarrow{L}\in {\mathbb{R}}^3$ (of length L := ||L||₂) for each RG. In accordance with the IAU Solar System convention for positive poles, the unit vector $\widehat{L}\in {\mathbb{S}}^2$ marks the direction from which the central Kerr black hole is seen rotating in anticlockwise direction³. We define the inclination angle Θ as the angle between L and a vector parallel to the line of sight pointing towards the observer⁴. Observations that allow one to measure the orientation of the RG axis in 3D are time-intensive, and so usually only the RG length projected onto the plane of the sky is known.

Geometrically, we model RGs as line segments – as if they were ‘thin sticks’, with vanishing volumes – so that the projected proper length RV L_p relates to L and Θ through

$$\begin{array}{c}\hfill L_{\mathrm{p}}=Lsin\mathrm{\Theta }.\end{array}$$(2)

We assume that $\widehat{L}$ is drawn from a uniform distribution on 𝕊², so that the PDF f_Θ : [0, π]→ℝ_≥0 becomes $f_{\mathrm{\Theta }}\left(\theta \right)=\frac{1}{2}sin\theta$. Since L and Θ are independent, we find in Appendix A.2.1 that the PDF of L_p, f_Lp : ℝ → ℝ_≥0, is

$$\begin{array}{c}\hfill f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \mathrm{if}{l}_{\mathrm{p}}\le 0,\hfill \\ \frac{-(\xi +1)}{l_{\min }}\frac{l_{\mathrm{p}}}{l_{\min }}I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})\hfill & \mathrm{if}0<l_{\mathrm{p}}<l_{\min },\hfill \\ \frac{{(\xi +1)}^2}{l_{\min }}\frac{\sqrt{\pi }}{4}{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi }\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}\hfill & \mathrm{if}{l}_{\mathrm{p}}\ge l_{\min },\hfill \end{array}\end{array}$$(3)

where

$$\begin{array}{c}\hfill I(a,b):={\displaystyle {\int }_1^{\infty }\frac{{\eta }^a\mathrm{d}\eta }{\sqrt{{\eta }^2-b^2}}\mathrm{for}a<0,|b|<1.}\end{array}$$(4)

We note that for l_p > l_min, the projected proper length has a Pareto distribution with the same tail index as the intrinsic proper length distribution. We compare f_L and f_Lp in Fig. 1.

2.2.2. Distribution for giants

For giants specifically (i.e. RGs such that L_p > l_p, GRG, where l_p, GRG is some constant threshold; in this work, l_p, GRG := 0.7 Mpc), the projected proper length distribution becomes a Pareto distribution with tail index ξ again:

$$\begin{array}{c}\hfill f_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le l_{\mathrm{p},\mathrm{GRG}}\hfill \\ \frac{-(\xi +1)}{l_{\mathrm{p},\mathrm{GRG}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi }\hfill & \text{if}{l}_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}.\hfill \end{array}\end{array}$$(5)

In other words, for giants, projection retains the Paretianity of lengths. A measurement of the tail index of the projected length distribution is immediately also a measurement of the tail index of the intrinsic length distribution.

The survival function, which gives the probability that a GRG has a projected proper length exceeding l_p, takes on a particularly simple form:

$$\begin{array}{c}\hfill \mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})={\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi +1}.\end{array}$$(6)

The mean projected proper length of giants is the expectation value of L_p | L_p > l_p, GRG:

$$\begin{array}{c}\hfill \mathbb{E}[L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]=l_{\mathrm{p},\mathrm{GRG}}\frac{\xi +1}{\xi +2},\end{array}$$(7)

which is only defined when ξ < −2. For example, when ξ = −4, $\mathbb{E}[L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]=\frac{3}{2}{l}_{\mathrm{p},\mathrm{GRG}}$, which becomes 1.05 Mpc for l_p, GRG := 0.7 Mpc and 1.5 Mpc for l_p, GRG := 1 Mpc.

Appendix A.2.2 provides a derivation for all three expressions.

2.3. Deprojection factor

The deprojection factor, $D:=\frac{L}{L_{\mathrm{p}}}={sin}^{-1}\mathrm{\Theta }$, quantifies how much larger intrinsic lengths are compared with projected lengths. The PDF of D, f_D : ℝ → ℝ_≥0, is

$$\begin{array}{c}\hfill f_D\left(d\right)=\{\begin{array}{cc}0\hfill & \mathrm{if}d\le 1;\hfill \\ \frac{1}{d^2\sqrt{d^2-1}}\hfill & \mathrm{if}d>1.\hfill \end{array}\end{array}$$(8)

The mean deprojection factor $\mathbb{E}[D]=\frac{\pi }{2}$. Deprojection factors can become arbitrarily large under the current model, because projected lengths can become arbitrarily small. As discussed in Sect. 5.2, this is not a very realistic set-up. In reality, an RG’s projected length is bounded from below by its lobes, which have a non-vanishing volume and thus extend along all three spatial dimensions. Upon projection, the projected length therefore cannot shrink beyond some lower limit that depends on the lobe geometry. In Appendix A.3, we show that by enriching the conventional stick-like geometry with spherical lobes, deprojection factors indeed become bounded.

2.4. Intrinsic proper length posterior and its moments

Because an RG’s intrinsic length is more physically informative than its projected length, we ideally obtain the former. In this subsection, we quantify what a measurement L_p = l_p already reveals about L.

We first note that the projected length bounds the intrinsic length from below. The intrinsic length can be much larger, however, but this is improbable for two reasons: large lengths are rarer than small lengths, although how drastic this effect is depends on ξ; in addition, viewing directions with large inclination angles are uncommon. The best we can do is to construct a posterior distribution for L given L_p = l_p. This posterior has a concise analytic form. If l_p > l_min, which is the relevant case for giants, the distribution of L | L_p = l_p is

$$\begin{array}{c}\hfill f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)=\{\begin{array}{cc}0\hfill & \text{if}l\le l_{\mathrm{p}}\hfill \\ \frac{-\xi }{2^{1+\xi }\pi }\frac{{\mathrm{\Gamma }}^2(-\frac{\xi }{2})}{\mathrm{\Gamma }(-\xi )}\frac{1}{l_{\mathrm{p}}}\frac{1}{\sqrt{{\left(\frac{l}{l_{\mathrm{p}}}\right)}^2-1}}{\left(\frac{l}{l_{\mathrm{p}}}\right)}^{\xi -1}\hfill & \text{if}l>l_{\mathrm{p}}.\hfill \end{array}\end{array}$$(9)

For l ≫ l_p, $f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)\propto {\left(\frac{l}{l_{\mathrm{p}}}\right)}^{\xi -2}$: the posterior PDF tends to a power law in l with index ξ − 2. In Fig. 2, we visualise the posterior PDF for several values of ξ.

Clearly, to evaluate Eq. (9), one must choose l_p – however, the shape of the distribution is the same for all choices. We illustrate this by comparing the PDF for a comparatively small GRG (l_p = 1.0 Mpc; top panel) to the PDF for Alcyoneus⁵ (l_p = 5.0 Mpc; bottom panel).

The posterior mean is

$$\begin{array}{c}\hfill \mathbb{E}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\mathrm{p}}·\frac{-\xi }{2^{2\xi +3}\pi }\frac{{\mathrm{\Gamma }}^4(-\frac{\xi }{2})}{{\mathrm{\Gamma }}^2(-\xi )};\end{array}$$(10)

the posterior variance is

$$\begin{array}{c}\hfill \mathbb{V}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\mathrm{p}}^2(\frac{\xi }{\xi +1}-\frac{{\xi }^2}{2^{4\xi +6}{\pi }^2}\frac{{\mathrm{\Gamma }}^8(-\frac{\xi }{2})}{{\mathrm{\Gamma }}^4(-\xi )}).\end{array}$$(11)

Both mean and standard deviation scale linearly in l_p: the projection effect is a multiplicative noise source. In Table 1, we provide explicit values for various ξ.

Higher moments exist up to order ⌈ − ξ⌉; because the PDF f_{L | Lp = lp}(l) is strongly skewed, such moments do further specify the distribution.

It is important to note that it is formally incorrect to statistically deproject RGs by drawing samples from deprojection factor D and multiplying them with some measurement L_p = l_p, or even more crudely, by multiplying the latter with 𝔼[D]. The reason that renders such approaches invalid is that L_p = L sin Θ and D = sin⁻¹Θ are not independent RVs. We refer to Appendix A.4 for an explicit proof of this fact, and for derivations of this subsection’s expressions.

2.5. GRG inclination angle

Radio galaxies with jets that make a small angle with the plane of the sky are more likely to have a projected length exceeding l_p, GRG than those with jets that make a large angle with the plane of the sky. For this reason, the inclination angle distribution of giants is different from that of RGs: it is more peaked around θ = 90°. More precisely, the PDF f_{Θ | Lp > lp, GRG} : [0, π]→ℝ_≥0 of the GRG inclination angle Θ | L_p > l_p, GRG has the general form

$$\begin{array}{c}\hfill f_{\mathrm{\Theta }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\theta )=\frac{(1-F_L\left(\frac{l_{\mathrm{p},\mathrm{GRG}}}{sin\theta }\right))f_{\mathrm{\Theta }}\left(\theta \right)}{1-F_{L_{\mathrm{p}}}(l_{\mathrm{p},\mathrm{GRG}})}.\end{array}$$(12)

Under our Pareto distribution assumption for L, this concretises to

$$\begin{array}{c}\hfill f_{\mathrm{\Theta }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\theta )=\frac{2}{-(\xi +1)\sqrt{\pi }}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{sin}^{-\xi }\theta .\end{array}$$(13)

We note that f_{Θ | Lp > lp, GRG}(θ)∝sin^−ξθ; the factor in front serves only as a normalisation constant. The distribution is independent of the choice of l_p, GRG and depends on a single parameter: ξ. We visualise f_{Θ | Lp > lp, GRG}(θ) in Fig. 3.

Appendix A.5 contains a brief derivation.

2.6. GRG angular length

The model predicts the distribution of GRG angular lengths in the Local Universe up to comoving distance r_max. The GRG angular length RV Φ | L_p > l_p, GRG relates to the GRG projected proper length RV L_p | L_p > l_p, GRG and the comoving distance RV R as

$$\begin{array}{c}\hfill \mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}=L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}·\frac{1+z\left(R\right)}{R}.\end{array}$$(14)

(We note that this relation is valid only in a flat Friedmann–Lemaître–Robertson–Walker (FLRW) universe.) We also assume that the GRG number density is constant in the Local Universe. The PDF of Φ | L_p > l_p, GRG has useful analytic forms under two different idealisations.

In a Euclidean universe, z(R) = 0, and the minimal GRG angular length ${\varphi }_{\mathrm{GRG}}=\frac{l_{\mathrm{p},\mathrm{GRG}}}{r_{\max }}$. Then

$$\begin{array}{c}\hfill f_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )=\{\begin{array}{cc}0\hfill & \mathrm{if}\varphi \le {\varphi }_{\mathrm{GRG}}\hfill \\ -3\frac{\xi +1}{\xi +4}·\frac{1}{{\varphi }_{\mathrm{GRG}}}({\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{\xi }-{\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{-4})\hfill & \mathrm{if}\varphi >{\varphi }_{\mathrm{GRG}},\hfill \end{array}\end{array}$$(15)

which is valid as long as ξ ≠ −4.

In an expanding universe at low redshifts, the Hubble–Lemaître law $z(R)\approx \frac{R}{d_H}$ holds; the Hubble distance $d_H:=\frac{c}{H_0}$. In this case,

$$\begin{array}{cc}& f_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )\hfill \\ \hfill & =\{\begin{array}{cc}0\hfill & \text{if}\varphi \le {\varphi }_{\mathrm{GRG}}+\frac{l_{\mathrm{p},\mathrm{GRG}}}{d_H}\hfill \\ \frac{-3(\xi +1){\varphi }^{\xi }}{r_{\max }^3{l}_{\mathrm{p},\mathrm{GRG}}^{\xi +1}}{\int }_{\frac{1}{r_{\max }}+\frac{1}{d_H}}^{\frac{\varphi }{l_{\mathrm{p},\mathrm{GRG}}}}\frac{\mathrm{d}k}{k^{\xi +1}{(k-\frac{1}{d_H})}^4}\hfill & \text{if}\varphi >{\varphi }_{\mathrm{GRG}}+\frac{l_{\mathrm{p},\mathrm{GRG}}}{d_H}.\hfill \end{array}\hfill \end{array}$$(16)

Figure 4 shows GRG angular length PDFs under both idealisations.

The PDFs undergo a minor shift upon changing universe type but are otherwise similar. For most current-day applications, it will therefore be unnecessary to calculate an even more refined version of f_{Φ | Lp > lp, GRG}(ϕ). Appendix A.6 contains derivations and details.

2.7. Maximum likelihood estimation of the tail index

The GRG projected proper length distribution features just one parameter of physical interest: the tail index ξ. If observational selection effects are negligible, one can directly use maximum likelihood estimation (MLE) on GRG data to infer ξ. In particular, we consider a set of projected lengths {L_p, 1, ..., L_p, N}∼L_p | L_p > l_p, GRG from N giants. Appendix A.7 shows that the maximum likelihood estimate of ξ is the RV ξ_MLE, given by

$$\begin{array}{c}\hfill {\xi }_{\mathrm{MLE}}=-\frac{N}{{\sum }_{i=1}^{N}ln\frac{L_{\mathrm{p},i}}{l_{\mathrm{p},\mathrm{GRG}}}}-1.\end{array}$$(17)

2.8. Observed projected proper length

2.8.1. General considerations

In the preceding theory, we have ignored observational selection effects that favour some projected proper lengths over others. In practice, several such effects occur; the importance of each varies per survey and (G)RG search campaign within it. One of them is the bias against physically long RGs that the interferometer’s largest detectable angular scale can induce⁶. As a result, the projected proper length of an observed RG might not be adequately modelled through RV L_p. Instead, we must introduce a new RV L_p, obs.

We define the completeness C : ℝ_{> 0} × ℝ_{> 0} → [0,1] at (l_p, z_max) to be the fraction of all RGs with projected proper length l_p in the cosmological volume up to z_max that is detected in a particular RG search campaign. Then, assuming that the distribution of L_p does not evolve with redshift between z = z_max and z = 0 (i.e. ξ remains constant),

$$\begin{array}{c}\hfill C(l_{\mathrm{p}},z_{\max })=\frac{{\int }_0^{z_{\max }}{p}_{\mathrm{obs}}(l_{\mathrm{p}},z)r^2\left(z\right)E^{-1}\left(z\right)\mathrm{d}z}{{\int }_0^{z_{\max }}{r}^2\left(z\right)E^{-1}\left(z\right)\mathrm{d}z},\end{array}$$(18)

where p_obs(l_p, z) is the probability that an RG of projected proper length l_p at cosmological redshift z is detected through the campaign, and r(z) is the comoving radial distance at cosmological redshift z. In a flat FLRW universe, the dimensionless Hubble parameter E is

$$\begin{array}{c}\hfill E\left(z\right):=\frac{H\left(z\right)}{H_0}=\sqrt{{\mathrm{\Omega }}_{\mathrm{R},0}{(1+z)}^4+{\mathrm{\Omega }}_{\mathrm{M},0}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }},0}.\end{array}$$(19)

The PDF of the observed projected proper length RV L_p, obs becomes

$$\begin{array}{c}\hfill f_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p}}\right)=\frac{C\left(l_{\mathrm{p}}\right)f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)}{{\int }_0^{\infty }C\left(l_{\mathrm{p}}^{\prime }\right)f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}^{\prime }\right)\mathrm{d}{l}_{\mathrm{p}}^{\prime }},\end{array}$$(20)

where we suppress the z_max-dependence for succinctness. We note that multiplying p_obs(l_p, z) with an l_p- and z-independent factor affects the completeness C(l_p, z_max), but cancels in Eq. (20); f_Lp, obs will be independent of it. Finally, the PDF of the GRG observed projected proper length RV L_p, obs | L_p, obs > l_p, GRG is

$$\begin{array}{c}\hfill f_{L_{\mathrm{p},\mathrm{obs}}|L_{\mathrm{p},\mathrm{obs}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le l_{\mathrm{p},\mathrm{GRG}}\hfill \\ \frac{C\left(l_{\mathrm{p}}\right)f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)}{{\int }_{l_{\mathrm{p},\mathrm{GRG}}}^{\infty }C\left(l_{\mathrm{p}}^{\prime }\right)f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}^{\prime }\right)\mathrm{d}{l}_{\mathrm{p}}^{\prime }}\hfill & \text{if}{l}_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}.\hfill \end{array}\end{array}$$(21)

We derive these expressions in Appendix A.8.1.

2.8.2. Fuzzy angular length threshold

We provide a concrete example of an important observational selection effect in visual searches for GRG candidates in survey images. To cope with the sheer number of detectable RGs in modern surveys like the LoTSS, a natural criterion is to only add sources to a candidate list if they appear – by eye – to have an angular length larger than some threshold. However, it is hard to precisely assess the angular length of a candidate before actually measuring it; sometimes, a candidate with a smaller angular length than the threshold will feature in the list, while some candidates with a larger angular length than the threshold will not. This leads to the notion of a ‘fuzzy angular length threshold’, where the probability that an RG with angular length ϕ is observed through the visual search increases (e.g., linearly) from 0 to 1 between ϕ_min and ϕ_max:

$$\begin{array}{cc}& p_{\mathrm{obs},\mathrm{AL}}(l_{\mathrm{p}},z)=min\{max\{\frac{\varphi (l_{\mathrm{p}},z)-{\varphi }_{\min }}{{\varphi }_{\max }-{\varphi }_{\min }},0\},1\},\hfill \end{array}$$(22)

$$\begin{array}{cc}& \varphi (l_{\mathrm{p}},z)=\frac{l_{\mathrm{p}}(1+z)}{r\left(z\right)}.\hfill \end{array}$$(23)

See the top row of Fig. 5 for several examples of associated completeness curves C(l_p) and observed projected proper length PDFs. See Appendix A.8.2 for additional information.

2.8.3. Surface brightness limitations

Another important observational selection effect is due to a survey’s finite noise level. The noise determines the surface brightness threshold b_ν, th (typically comparable to the noise level itself) below which radio galaxy features remain visually undetected.

Fanaroff–Riley class II. We model the lobe surface brightness RV B_ν at the central observing frequency ν_obs as

$$\begin{array}{c}\hfill B_{\nu }=\frac{b_{\nu ,\mathrm{ref}}·S}{{(1+Z)}^{3-\alpha }}{\left(\frac{L}{l_{\mathrm{ref}}}\right)}^{\zeta },\end{array}$$(24)

where b_ν, ref is the median surface brightness of RGs of intrinsic proper length l_ref at cosmological redshift z = 0 and frequency ν_obs, and S is a lognormally distributed RV with median 1 and dispersion parameter σ_ref that captures the variability in surface brightness among this population of RGs. The denominator models the fact that surface brightness is not conserved with distance in an FLRW universe; Z is the cosmological redshift RV up to z = z_max and α is the typical lobe spectral index. The exponent ζ < 0 characterises the scaling between intrinsic proper length and surface brightness. (If RGs are self-similar, so that morphology does not predict length, one finds ζ = −2.) For this selection effect, the observing probability is

$$\begin{array}{cc}& p_{\mathrm{obs},\mathrm{SB}}(l_{\mathrm{p}},z)={\int }_{s_{\min }}^{\infty }\sqrt{1-{\left(\frac{s_{\min }}{s}\right)}^{-\frac{2}{\zeta }}}{f}_S\left(s\right)\mathrm{d}s,\hfill \end{array}$$(25)

$$\begin{array}{cc}& s_{\min }=\frac{b_{\nu ,\mathrm{th}}}{b_{\nu ,\mathrm{ref}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta }{(1+z)}^{3-\alpha },\hfill \end{array}$$(26)

$$\begin{array}{cc}& f_S\left(s\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{ref}}s}exp(-\frac{{ln}^{2}s}{2{\sigma }_{\mathrm{ref}}^2}).\hfill \end{array}$$(27)

We note that p_obs, SB does not depend on b_ν, th or b_ν, ref separately, but on their ratio only. See the bottom row of Fig. 5 for several examples of associated completeness curves C(l_p) and observed projected proper length PDFs.

Fanaroff–Riley class I. For FRI RGs, the assumption of a constant surface brightness beyond the core is inaccurate. The simplest correction in which FRI RGs retain a well-defined notion of length assumes a linearly decreasing surface brightness profile, which peaks at the core and goes to zero at the RG’s two endpoints. (A power-law profile does not work: in such case, the surface brightness only goes asymptotically to zero – but never actually reaches it.) In this case, we define RV B_ν to be the mean surface brightness along an RG’s jets, which can be regarded as a typical value for that RG. As B_ν again obeys Eq. (24), we find that the formulaic structure of p_obs, SB(l_p, z) is identical for FRI and FRII giants, except that FRI giants require a change

$$\begin{array}{c}\hfill b_{\nu ,\mathrm{th}}\to \frac{b_{\nu ,\mathrm{th}}}{2(1-\frac{l_{\mathrm{p},\mathrm{GRG}}}{l_{\mathrm{p}}})},\end{array}$$(28)

which affects p_obs, SB(l_p, z) through s_min. There is no change for l_p = 2l_p, GRG. Although the formulaic structure might be the same, the best-fit parameters can differ. For example, it is possible that σ_ref, FRI ≠ σ_ref, FRII or ζ_FRI ≠ ζ_FRII. See Appendix A.8.3 for derivations and numerical implementation considerations.

To include both aforementioned selection effects at the same time, a natural approximation is to assume that the observing probability simply factorises:

$$\begin{array}{c}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)\approx p_{\mathrm{obs},\mathrm{AL}}(l_{\mathrm{p}},z)·p_{\mathrm{obs},\mathrm{SB}}(l_{\mathrm{p}},z).\end{array}$$(29)

2.9. GRG number density

A central question in the field of radio galaxies is: how intrinsically rare are giants? By counting giants in a search campaign with well-understood selection effects, we can give an answer. More precisely, one can estimate the comoving number density of giants in the Local Universe, n_GRG, through the number of observed giants up to cosmological redshift z_max in a uniformly searched region of sky of solid angle Ω, N_GRG, obs(z_max,Ω). Then, under the standard assumption l_p, GRG > l_min, Appendix A.9 shows that

$$\begin{array}{cc}& n_{\mathrm{GRG}}=\frac{l_{\mathrm{p},\mathrm{GRG}}^{\xi +1}}{-(\xi +1)}·\frac{H_0}{c}·\hfill \\ \hfill & \frac{\frac{4\pi }{\mathrm{\Omega }}{N}_{\mathrm{GRG},\mathrm{obs}}(z_{\max },\mathrm{\Omega })}{{\int }_{l_{\mathrm{p},\mathrm{GRG}}}^{\infty }{l}_{\mathrm{p}}^{\xi }{\int }_0^{z_{\max }}{p}_{\mathrm{obs}}(l_{\mathrm{p}},z)4\pi r^2\left(z\right)E^{-1}\left(z\right)\mathrm{d}z\mathrm{d}{l}_{\mathrm{p}}}.\hfill \end{array}$$(30)

Although the appropriate p_obs(l_p, z) varies per search campaign, it is always possible to bound p_obs(l_p, z) from above – for example by 1. In such case, Eq. (30) bounds n_GRG from below.

2.10. GRG lobe volume-filling fraction

Because giants attain cosmological lengths, they might contribute to the energisation and magnetisation of Cosmic Web filaments in regions that smaller radio galaxies cannot reach. A key statistic that measures the enrichment of the Cosmic Web by giants is the volume-filling fraction (VFF) of their lobes. Assuming that lobes do not grow along with the expansion of the Universe, the proper VFF changes over cosmic time: VFF_GRG(z) = VFF_GRG(z = 0)⋅(1 + z)³, where

$$\begin{array}{cc}& {\mathrm{VFF}}_{\mathrm{GRG}}(z=0):=n_{\mathrm{GRG}}·\mathbb{E}[V|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]\hfill \\ \hfill & =n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}·L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}],\hfill \end{array}$$(31)

where V is the combined volume of the lobes and $\mathrm{{\rm Y}}:=\frac{V}{L^3}$ is a dimensionless RV that captures the diversity in radio galaxy lobe shapes. We find under self-similar growth

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)=n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}]·\mathbb{E}[L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}].\end{array}$$(32)

𝔼[Υ] can be estimated from observations, but one must be wary of selection effects. Unfortunately, 𝔼[L³ | L_p > l_p, GRG] does not exist for ξ ≥ −4, which is the regime supported by observations. A useful lower bound then is

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)>n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}]·{\mathbb{E}}^3[L|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}],\end{array}$$(33)

where

$$\begin{array}{c}\hfill \mathbb{E}[L|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]=l_{\mathrm{p},\mathrm{GRG}}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-1)\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})\mathrm{\Gamma }(-\frac{\xi }{2}+\frac{1}{2})}.\end{array}$$(34)

This is the mean intrinsic proper length of giants, which is only defined when ξ < −2.⁷

An alternative is to deviate slightly from our Pareto ansatz and truncate the GRG projected proper length distribution at some l_p, max; then

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)=n_{\mathrm{GRG}}·\mathbb{E}[{\mathrm{{\rm Y}}}_{\mathrm{p}}]·\frac{\xi +1}{\xi +4}·\frac{l_{\mathrm{p},\max }^{\xi +4}-l_{\mathrm{p},\mathrm{GRG}}^{\xi +4}}{l_{\mathrm{p},\max }^{\xi +1}-l_{\mathrm{p},\mathrm{GRG}}^{\xi +1}},\end{array}$$(35)

where ${\mathrm{{\rm Y}}}_{\mathrm{p}}:=\frac{V}{L_{\mathrm{p}}^3}$. See Appendix A.10 for a derivation and further details.

2.11. Unification model constraints from quasar and non-quasar giants

The unification model and its extensions (e.g., Hardcastle & Croston 2020) form an elegant family of hypotheses that aim to explain the observational diversity of active galaxies. It posits that active galaxies with quasars differ from those without quasars primarily because of differences in orientation of the dusty tori surrounding SMBHs. In particular, the central idea is that a quasar appears brighter to the observer than a non-quasar AGN because the axis of its dusty torus happens to be virtually parallel to the line of sight. As such, only quasars would offer an unobscured view of the luminous accretion disk surrounding the SMBH, whilst also beaming relativistic jet emission towards the observer. Using our statistical framework, we predict the general ramifications of the basic unification model on a GRG sample.

The basic unification model suggests to divide the radio galaxy population in two, distinguishing between RGs generated by AGN with quasar appearance (quasar RGs) and RGs generated by AGN without quasar appearance (non-quasar RGs). We assume that quasar RGs have inclination angles θ ≤ θ_max or θ ≥ 180° −θ_max whilst non-quasar RGs have θ_max < θ < 180° −θ_max⁸. If quasar RGs are more closely aligned with the line of sight than non-quasar RGs but are otherwise similar, fewer of them will satisfy l_p ≥ l_p, GRG and thus be classified as giants. Therefore, the quasar GRG fraction f_Q – the fraction of quasar giants in an actual GRG sample – constrains θ_max. We model f_Q as an RV:

$$\begin{array}{c}\hfill f_{\mathrm{Q}}:=\frac{N_{\mathrm{Q}}}{N};N_{\mathrm{Q}}\sim \mathrm{Binom}(N,p_{\mathrm{Q}}),\end{array}$$(36)

where the RV N_Q is the number of quasar giants in the sample, the constant N is the total number of giants in the sample and the parameter p_Q is the quasar GRG probability. Our framework predicts

$$\begin{array}{cc}& p_{\mathrm{Q}}:=\frac{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}},sin\mathrm{\Theta }\le sin{\theta }_{\max })}{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}})}\hfill \\ \hfill & =\frac{4\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{-(\xi +1)\sqrt{\pi }\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\displaystyle {\int }_0^{sin{\theta }_{\max }}\frac{x^{-\xi }\mathrm{d}x}{\sqrt{1-x^2}}.}\hfill \end{array}$$(37)

See Appendix A.11 for a derivation. Interestingly, as long as quasar giants and non-quasar giants are subject to the same selection effects, these selection effects do not affect p_Q.

Figure 6 shows that, for all relevant ξ, p_Q is a steeply and monotonically increasing function of θ_max. Thus, knowing ξ, one can use an empirical f_Q to determine p_Q and in turn θ_max.

Does one expect quasar giants to have a different distribution for L_p, obs | L_p, obs ≥ l_p, GRG than non-quasar giants? Interestingly, our framework allows us to prove that the inclination angle differences between the two classes affect their relative rarity, but not their observed projected proper length distribution. Under the basic unification model, quasar giants and non-quasar giants obey the same L_p, obs | L_p, obs ≥ l_p, GRG if they are subject to the same selection effects.

2.12. Extreme giants in a sample

The model can predict the occurrence of giants with extreme projected proper lengths in a sample of cardinality N. The probability that an observed GRG has a projected proper length exceeding l_p is

$$\begin{array}{c}\hfill p_{>l_{\mathrm{p}}}:=\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}>l_{\mathrm{p}}|L_{\mathrm{p},\mathrm{obs}}>l_{\mathrm{p},\mathrm{GRG}})=\frac{1-F_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p}}\right)}{1-F_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p},\mathrm{GRG}}\right)},\end{array}$$(38)

so that the number of observed giants with a projected proper length exceeding l_p is N_{> lp} ∼ Binom(N, p_{> lp}). Its expectation is 𝔼[N_{> lp}]=N ⋅ p_{> lp}. Furthermore, the probability that the sample contains at least one giant with projected proper length l_p or higher is

$$\begin{array}{c}\hfill \mathbb{P}(N_{>l_{\mathrm{p}}}\ge 1)=1-{(1-p_{>l_{\mathrm{p}}})}^N.\end{array}$$(39)

See Appendix A.12 for details. Figure 7 shows 𝔼[N_{> lp}] and ℙ(N_{> lp} ≥ 1) for various ξ. As an example, the case ξ = −3.0 predicts that a sample of N = 1000 giants with redshifts below z_max = 0.5 should contain almost ten giants of l_p > 5 Mpc, and still several of l_p > 6 Mpc. Such predictions are useful as they can be directly compared to elementary sample statistics.

3. Sample compilation and properties

To measure the intrinsic GRG length distribution, we also require a large sample of giants collected from a single survey through a systematic approach. This ensures approximately homogeneous selection effects, which we can correct for in subsequent analysis using the statistical framework of Sect. 2.

3.1. LoTSS DR2

The Low-Frequency Array (LOFAR; van Haarlem et al. 2013) is a powerful, Pan-European radio interferometer that features both (sub)arcsecond-scale resolution and sensitivity to degree-scale structures. Dabhade et al. (2020b) have already demonstrated that this combination is ideal for detecting giants: these authors found a record 225 new specimina in the LoTSS DR1, the first data release of the LOFAR’s Northern Sky survey at central frequency ν_obs = 144 MHz.⁹ Excitingly, the recent LoTSS DR2 (Shimwell et al. 2022) improves the data calibration and increases the survey footprint from 424 deg² to 5635 deg² – that is by more than a factor 13. By default, the LoTSS features imagery at 6″ and 20″ resolutions. To further facilitate the discovery of giants (among other goals), we reprocessed the LoTSS by subtracting compact sources and imaging at 60″ and 90″ resolution; more details are given in Oei et al. (2022a, and in prep.). This 60″ and 90″ imagery has turned out to be effective in highlighting jets and lobes of RGs of large angular and physical extent, whose surface brightnesses are low and which therefore have remained undetected in shallower surveys, and even in the LoTSS DR2 at higher resolutions. We demonstrate this fact in Fig. 8 by comparing the LoTSS DR2 at 6″ and 60″ for three giants whose discovery has relied on the lower-resolution images. After the serendipitous discovery of several such hitherto unknown giants in the 60″ and 90″ images, we decided to initiate a systematic, multi-resolution, visual GRG search through the area covered by LoTSS DR2 pipeline products as of September 2022¹⁰. This search comprised of a hundreds-of-hours-long inspection of the LoTSS maps at 6″ and 60″, alongside Pan-STARRS DR1 (Chambers et al. 2016) and SDSS DR9 (Ahn et al. 2012) maps, in Aladin Desktop 11.0 (Bonnarel et al. 2000).

Reliable automated search strategies do not yet exist for several reasons. Giants showcase a rich morphological variety (see Sect. 3.6) and are easily confused with other types of astrophysical sources (see Sects. 3.3 and 3.4). Moreover, the known population is too small to effectively apply supervised learning techniques. However, it appears possible to find giants in morphological outlier lists of unsupervised learning techniques such as self-organising maps (SOMs; Mostert et al. 2021). The efficacy of such techniques in GRG searches has not yet been quantified.

3.2. Angular length threshold

To limit the amount of manual work, we decided to search only for GRG candidates whose angular length exceeds some threshold. In Fig. 9, we show how the angular length ϕ of giants with six different projected proper lengths l_p varies as a function of cosmological redshift z. Because of the expansion of the Universe, giants have an arcminute-scale minimum angular length: if l_p, GRG = 0.7 Mpc, all giants obey ϕ > 1.3′¹¹. For the purpose of finding a GRG, it is therefore never useful to inspect a source with an angular length less than 1.3′. Our highest priority has been to find giants with z < 0.2, which lie in a volume for which the total matter density field is known or will be known in the coming years through the combined power of deep spectroscopic surveys and Bayesian inference frameworks, such as the Bayesian Origin Reconstruction from Galaxies (BORG; Jasche & Wandelt 2013; Jasche et al. 2015; Jasche & Lavaux 2019). In an upcoming publication, we combine a sample of low-redshift giants with the BORG to analyse the large-scale environments of giants (Oei et al., in prep.). Figure 9 shows that all giants with l_p > 1 Mpc at z < 0.2 have an angular length ϕ > 5′. For this reason, we have chosen 5′ as the angular length threshold of our search campaign. This choice has kept the visual inspection duration to order ∼10² h, while still enabling us to target all Mpc-exceeding giants in the Local Universe (z < 0.2). In practice, this threshold is ‘fuzzy’: it is hard to accurately estimate angular lengths by eye before performing an actual measurement, so that our list of GRG candidates does contain some with a smaller angular length than the specified threshold. Inversely, it will presumably lack some GRG candidates with an angular length exceeding the threshold.

3.3. GRG candidates in the radio

We first identified GRG candidates in the LoTSS at 6″ and 60″. These maps serve complementary roles. The 6″ images reveal the precise morphology of radio galaxy cores and jets, which are necessary to pinpoint the host galaxy. Figure 10 provides a representative sense of the host galaxy identification accuracy these data allow for when combined with modern optical surveys such as the DESI Legacy Imaging Surveys (Dey et al. 2019) DR9. In contrast, the 60″ images have such compact sources removed or highly suppressed, but better highlight diffuse structures, such as RG lobes. Being similar in morphology, we made sure not to interpret diffuse emission from low-redshift spiral galaxies and their circumgalactic media, or radio halos and relics in galaxy clusters, as RG lobe emission. We required that all new RGs feature a detection of at least two¹² lobes, or of at least one lobe and one jet oriented towards the lobe(s), at least one of the resolutions used in this work (6″, 20″, 60″, and 90″).

3.4. GRG candidates in the optical

To confirm that a radio structure really is a radio galaxy, we compared the radio images with optical images of the same sky region. If a patch of radio emission is indeed due to RG jets or lobes, the patch itself must have no codirectional galactic counterpart in the optical. If the radio emission is due to a low-redshift spiral galaxy or a galaxy cluster instead, a corresponding easily recognisable counterpart will exist. We also took care not to erroneously associate the lobes of two distinct RGs. For this reason, we were more cautious to associate a pair of lobes to a suspected host galaxy when, in the optical, one could discern other galaxies in the angular vicinity of the lobes that could have generated them instead.

The Pan-STARRS and SDSS images used for these purposes complement each other, as they differ in quality throughout the sky – and in particular around sources of high optical flux density. Neither consistently outperforms the other. Only Pan-STARRS covers the full Northern Sky and could thus always be relied upon.

3.5. Host galaxy identification

We also used the Pan-STARRS and SDSS maps for the identification of host galaxies. We collected host redshifts from the SDSS DR12 (Alam et al. 2015) and Gaia (Gaia Collaboration 2016) DR3 (Gaia Collaboration 2021) through automated VizieR queries, from the Galaxy List for the Advanced Detector Era (GLADE) 2.4 (Dálya et al. 2018), and from the DESI DR9 photometric redshift catalogue (Zou et al. 2022). If redshifts from multiple sources were available, we favoured SDSS over GLADE data, GLADE over Gaia data, and Gaia over DESI data. Similarly, we only adopted photometric redshifts if spectroscopic ones were not available.

For a small subset of RGs, a definite host galaxy could not be established beyond reasonable doubt, but a set of candidates containing the host galaxy could. In these cases, the lowest-redshift candidate provides a lower bound to the RG’s projected proper length¹³. If this lower bound exceeds l_p, GRG = 0.7 Mpc, then the actual projected proper length certainly does so too, and the RG can be classified as a GRG – despite persisting uncertainty concerning the identity of the host galaxy.

3.6. LoTSS DR2 GRG sample

Our search campaign has led to the identification of 2060 hitherto unknown giants. To establish novelty, we assembled a literature catalogue with all known giants as of September 2022, combining the catalogue of Dabhade et al. (2020a)¹⁴ with the giants discovered in Galvin et al. (2020), Ishwara-Chandra et al. (2020), Tang et al. (2020), Bassani et al. (2021), Brüggen et al. (2021), Delhaize et al. (2021), Masini et al. (2021), Kuźmicz & Jamrozy (2021), Andernach et al. (2021), Mahato et al. (2022), Gürkan et al. (2022), and Simonte et al. (2022). Fusing our sample with this literature catalogue, we obtain a final catalogue with N = 3341 giants.

Figure 11 shows the locations of all known giants in the sky. We list basic properties of the 50 projectively largest new discoveries in Table 2, and refer to Appendix F for access to these data for our entire sample. In Fig. 12 and Figs. B.1–B.3, we present images for discoveries with projected proper lengths l_p in the ranges 5.1–4 Mpc, 4–3 Mpc, 3–2 Mpc, and 2–0.7 Mpc, respectively. For each giant, we use the LoTSS resolution θ_FWHM ∈ {6″, 20″, 60″, 90″} that most clearly conveys the morphology through a single image. The selection reflects our sample’s diversity in shapes and sizes and provides a sense of the data quality.

3.6.1. Angular lengths

The angular length distribution of the newly found giants is as follows: the smallest ϕ = 1.5′, the median ϕ = 5.2′, the largest ϕ = 2.2°, and 80% of angular lengths fall within [3.4′,9.8′]. Thirteen of our discoveries – listed in Table 3 – are larger than the Moon in the sky (whose angular diameter varies over time, but here taken to be ϕ = 30′). Our search more than doubles the known number of such spectacular giants – from 10 to 23.

The GRG associated with NGC 2300 (see the middle-left panel of Fig. 13) is the giant with the largest angular length ever found, and the radio galaxy with the largest angular length in the Northern Sky¹⁵. It remains possible that the GRG has been generated by spiral galaxy NGC 2276 instead, with which elliptical galaxy NGC 2300 is interacting. However, this scenario seems unlikely, as only a fraction ∼10⁻³ of known giants are hosted by spirals. Its discovery emphasises that low-frequency interferometers like the LOFAR and the MWA, which are sensitive to degree-scale angular scales, are important to complete a low-redshift census of giant radio galaxies. Sky-wide, a simple extrapolation of our findings suggests that several (∼10¹) degree-scale angular length giants similar to the GRG of NGC 2300 still await discovery at the LoTSS DR2 depth.

3.6.2. Redshifts

The redshift distribution of the newly found giants is as follows: the lowest z = 0.00635 ± 6 × 10⁻⁵, the median z = 0.29, the highest z = 2.6394 ± 6 × 10⁻⁴, and 80% of redshifts fall within [0.12, 0.68]. Because of our focus on giants of large angular length (see Sect. 3.2), we have found only 36 giants beyond z > 1. One of these, the GRG at rank 13 of Table 2, is the largest secure giant found beyond redshift 1. It lies at z = 1.1 ± 0.1 and spans l_p = 3.9 ± 0.1 Mpc. Its host galaxy does not appear to contain a quasar.

3.6.3. Projected proper lengths

With l_p = 5.1 ± 0.2 Mpc and l_p = 4.99 ± 0.04 Mpc, our LoTSS DR2 sample contains the first two 5 Mpc–scale giants. We have presented a dedicated analysis of the latter GRG in Oei et al. (2022a). 11 discoveries have l_p ≥ 4 Mpc, 53 have 3 ≤ l_p < 4 Mpc, 291 have 2 ≤ l_p < 3 Mpc, 1215 have 1 ≤ l_p < 2 Mpc, and 490 have 0.7 ≤ l_p < 1 Mpc. The median l_p = 1.35 Mpc, and 80% of projected proper lengths fall within [0.82 Mpc, 2.29 Mpc].

3.6.4. Stellar and supermassive black hole masses

Following Oei et al. (2022a), we collected host stellar masses M_⋆ from Chang et al. (2015) and Salim et al. (2018), and estimated host SMBH masses M_• via SDSS DR12 stellar velocity dispersions (Alam et al. 2015) and the M-sigma relation of Kormendy & Ho (2013)’s Eq. (7). From all 3341 giants in our final catalogue, only 732 (22%) could be assigned a stellar mass in this way, and only 1115 (33%) an SMBH mass; for both quantities, our LoTSS DR2 sample accounts for four-fifths of the resulting subpopulation. Figure C.1 shows both M_⋆ and M_• in relation to projected proper length l_p.

The median M_⋆ = 3.4 × 10¹¹ M_⊙, and 80% of stellar masses fall within [1.8 × 10¹¹ M_⊙, 5.3 × 10¹¹ M_⊙]. We discover two giants whose hosts, J150329.07+374850.3 and J073505.24+415827.5, are the least massive known: both have a stellar mass M_⋆ = 4.8 × 10¹⁰ M_⊙. These are small giants, with l_p = 0.8 Mpc and l_p = 0.7 Mpc, respectively. The top panel of Fig. C.1 hints at a weak positive correlation between M_⋆ and l_p, which future work should confirm or reject.

The median M_• = 1.0 × 10⁹ M_⊙, and 80% of SMBH masses fall within [0.4 × 10⁹ M_⊙, 2.2 × 10⁹ M_⊙]. The SMBH masses of J123703.24+275819.5, J163659.07+541725.4, and J103129.54+502959.1 are the highest estimated yet, with M_• = 2 × 10¹¹ M_⊙, M_• = 5 × 10¹⁰ M_⊙, and M_• = 5 × 10¹⁰ M_⊙, respectively. The latter masses equal the theoretical maximum mass of accreting black holes of typical spin (King 2016). Curiously, although J103129.54+502959.1’s M_• is among the highest estimated SMBH mass of any GRG host, the GRG itself is relatively small: l_p = 0.81 ± 0.06 Mpc. Conversely, the bottom panel of Fig. C.1 shows that multi-Mpc giants can have hosts with SMBH masses that are two orders of magnitude lower than J103129.54+502959.1’s.

3.6.5. Spiral or lenticular host galaxies

Remarkably, although NGC 6185 is a spiral galaxy of Hubble–de Vaucouleurs class SAa (Jansen et al. 2000), it appears to have generated the giant shown in the top-left panel of Fig. 13. Such systems are exceedingly rare: not only are few RGs giants, but also virtually all giants have an elliptical galaxy as their host. With l_p = 2.54 ± 0.01 Mpc, this is the largest known spiral galaxy–hosted radio galaxy. Hitherto, the largest spiral galaxy–hosted giant in the literature has been J2345-0449 (Bagchi et al. 2014), with l_p = 1.6 Mpc. Given the favourably low redshift of z = 0.03430 ± 7 × 10⁻⁵, NGC 6185 and its enigmatic giant solicit a dedicated analysis (Oei et al. 2023). Besides NGC 6185, spiral or lenticular galaxies J080403.40+404809.3 and J091459.66+294348.8 (known alternatively as NGC 2789; see the bottom-right panel of Fig. 13) have also generated giants; these have projected proper lengths l_p = 1.1 Mpc and l_p = 0.9 Mpc, respectively. Morphological host classification through SDSS, Pan-STARRS, and DESI imagery is reliable only up to z ∼ 0.1–0.2, depending on viewing angle and various other factors. Our LoTSS DR2 sample contains 342 giants with definite hosts at z = 0.15 or below, among which are all 3 giants with spiral or lenticular hosts discussed here. It thus appears that the fraction of GRG hosts that is of such non-elliptical nature is ∼1%. A more detailed morphological characterisation of the sample appears possible, for example using data from Hart et al. (2016), but this is beyond the scope of the current work.

4. Results

After first building a statistical framework and then collecting a large sample of giants from a single survey through a systematic approach, we were ready to infer the intrinsic GRG length distribution. In particular, we aimed to establish whether the RG intrinsic proper length RV L is well described by a Pareto distribution, and if so, what its tail index ξ is. Subsequently, we inferred derived quantities, such as the comoving GRG number density in the Local Universe.

4.1. Giant radio galaxy length distribution

4.1.1. Empirical survival function

From our LoTSS DR2 GRG sample, we computed the empirical cumulative distribution function (ECDF) of the GRG observed projected proper length RV L_p, obs | L_p, obs > l_p, GRG (see Eq. (21)). We only included giants for which l_p itself – rather than a lower bound – is known, and set either z_max = 0.5 or z_max = 0.25; this retained 1473 or 811 giants for analysis, respectively¹⁶. The empirical survival function (ESF), which equals one minus the ECDF, is shown in dark grey in both panels of Fig. 14. Just as the ECDF approximates the CDF, the ESF approximates the survival function (SF). In this case, for any l_p, the ESF provides the probability that a randomly drawn LoTSS DR2 GRG will have a projected proper length exceeding l_p. If L_p, obs | L_p, obs > l_p, GRG were Paretian, its ESF (hereafter: ‘the’ ESF) would resemble a straight line in Fig. 14’s right panel – both axes have logarithmic scaling. However, the ESF clearly displays curvature.

4.1.2. Selection effects

Under the ansatz that L is Paretian, as proposed in Sect. 2.1, the aforementioned ESF’s curvature implies a significant role for observational selection effects. The reason is the following. The ansatz implies that the GRG projected proper length RV L_p | L_p > l_p, GRG is also Paretian (see Eq. (5)), as is L_p, obs | L_p, obs > l_p, GRG when selection effects are negligible (see Sect. 2.8.1 and set C(l_p) = 1). But in our case L_p, obs | L_p, obs > l_p, GRG is not Paretian: its ESF is curved. To avoid contradiction, we must relax at least one assumption. If the Pareto ansatz is held, then selection effects must be at play.

When selection effects are non-negligible, we must devise a procedure to disentangle them from the data if our ξ estimate is to be uncontaminated. To this end, we performed joint Bayesian inference with a model that includes both ξ and parameters that describe the selection effects.

In particular, we considered the roles of the fuzzy angular length threshold and surface brightness selection effects, introduced in Sects. 2.8.2 and 2.8.3, respectively. In Sect. 3.2, we explain that we have attempted to maintain a 5′ angular length threshold during our LoTSS DR2 GRG search. If we want to use Sect. 2.8.2 to correct for the bias against faraway and physically small giants that this threshold has imprinted onto our sample, we must estimate parameters ϕ_min and ϕ_max. A natural choice is to assume that they lie symmetrically around the intended angular length threshold of 5′ – but at what distance from it? We propose to consider this distance, $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$, as a yet unknown model parameter that we fitted to the data. Similarly, if we want to use Sect. 2.8.3 to correct for the bias against physically large giants that the limited depth of the LoTSS DR2 imprints onto our sample, we must estimate parameters b_ν, ref and σ_ref. Our approach was to, again, fit these parameters – possibly making use of any available prior knowledge.

The meaning of b_ν, ref depends on the choice of l_ref and ν_obs; we defined l_ref := 0.7 Mpc and used ν_obs = 144 MHz. We furthermore assumed b_ν, th = 1 × σ_Iν, with σ_Iν = 25 Jy deg⁻² being the typical LoTSS DR2 6″ RMS noise (Shimwell et al. 2022).

4.1.3. Prior

We exploited two sources of prior knowledge. Firstly, we attempted to directly estimate b_ν, ref and σ_ref by selecting from all LoTSS DR2 giants with l_p ≤ 1 Mpc a random subset of size 50 (10%). For these giants, we estimated the mean surface brightnesses of both lobes from the LoTSS DR2 6″ imagery, differentiating between the brighter and the fainter lobe. Because our goal was to estimate b_ν, ref, we attempted to undo cosmological and growth-induced surface brightness dimming assuming a universal lobe spectral index α = −1 and self-similar growth: ζ = −2. The resulting surface brightnesses correspond to z = 0, and to the epoch in each giant’s life when l = l_ref. The median of the corrected bright lobe mean surface brightnesses is b_ν, ref = 1.3 × 10³ Jy deg⁻², whilst the median of the corrected faint lobe mean surface brightnesses is b_ν, ref = 0.7 × 10³ Jy deg⁻². All lobes taken together, the median becomes b_ν, ref = 1.0 × 10³ Jy deg⁻². We performed maximum likelihood estimation assuming the surface brightness distribution is lognormal and found σ_ref = 1.3, again using all lobes. Note, however, that we have only used observed giants here, whilst b_ν, ref and σ_ref should correspond to the entire population of giants. As fainter giants will have preferentially fallen out, we might have overestimated b_ν, ref and underestimated σ_ref.

To probe whether we had overestimated b_ν, ref and underestimated σ_ref, we used data from Dabhade et al. (2020b) and a Monte Carlo approach. First, for their sample of 239 LoTSS DR1 giants, we computed total luminosity densities L_ν at rest-frame frequency ν = 144 MHz. (Given that LoTSS DR1 and DR2 noise levels are similar, this population is also representative of the LoTSS DR2.) In the top panel of Fig. 15, we show L_ν versus l_p for all 139 for which z < 0.6. To increase the range of lengths covered, we additionally show data on Alcyoneus (Oei et al. 2022a). Next, under the same assumptions of a constant spectral index and self-similar growth, we derived a simple luminosity density–surface brightness relationship that allows for back-and-forth conversion between the two – at least, given projected lengths and redshifts. The mean lobe surface brightness ⟨b_ν⟩ is proportional to the total luminosity density L_ν, and assuming a pair of spherical lobes

$$\begin{array}{c}\hfill ⟨b_{\nu }⟩=\frac{2f_{L_{\nu }}·L_{\nu }}{{\pi }^2·\mathbb{E}[D](\eta (f_l))·f_l^2·l_{\mathrm{p}}^2·{(1+z)}^{3-\alpha }}.\end{array}$$(40)

Here f_Lν is the fraction of the total luminosity density that belongs to the lobes, f_l is the fraction of the RG’s axis length that lies inside the lobes, and 𝔼[D](η(f_l)) is the mean deprojection factor as given by Eq. (A.29). The peak surface brightness b_ν, max relates to ⟨b_ν⟩ as $b_{\nu ,\max }=\frac{3}{2}⟨b_{\nu }⟩$. Appendix D contains derivations for both these results. For Alcyoneus, at z_Alc. = 0.25, f_Lν = 0.3 and f_l = 0.3 (Oei et al. 2022a). Assuming these parameter values, again in the top panel of Fig. 15, we show luminosity density–projected length pairs of RGs at z = z_Alc. whose peak (solid green curve) or mean (dash-dotted green curve) surface brightness equals the LoTSS DR2 noise level. Thus, each curve represents a family of borderline-detectable giants at Alcyoneus’s redshift. The other dash-dotted curves indicate similar barely detectable families, but for other redshifts. Without optimising any free parameters, the curves correctly predict that Alcyoneus’s lobes have surface brightnesses comparable to the 6″ LoTSS DR2 noise level and explain the absence of observations in the top-left corner of the figure. We conclude that Eq. (40) appears reasonable, but note that RGs may significantly differ in their values of f_Lν and f_l.¹⁷

Bolstered, we made use of Eq. (40) to Monte Carlo simulate – for particular values of b_ν, ref and σ_ref – luminosity density–projected length relationships as they appear to observers. The simulated giants have projected lengths adopted from the observed giants, randomly sampled redshifts up to z = 0.6 assuming a spatially constant GRG number density, and randomly sampled reference surface brightnesses (i.e. those for RGs at z = 0 that have intrinsic length l = l_ref) whose distribution is determined by b_ν, ref and σ_ref. We then used Eq. (24) to compute surface brightnesses as they would be observed, and retained only those giants whose surface brightness exceeds the LoTSS DR2 noise level. For these simulated detectable giants, we finally generated luminosity densities using Eq. (40), assuming wide uniform distributions f_Lν ∼ f_l ∼ Uniform(0.1, 0.9). The middle and bottom panels of Fig. 15 show results for b_ν, ref = 1750 Jy deg⁻² and b_ν, ref = 250 Jy deg⁻², respectively; we adopted σ_ref = 1.3 from our LoTSS DR2 GRG surface brightness measurements. The median luminosity density of the z < 0.6 LoTSS DR1 giants is L_ν = 1.1 × 10²⁶ W Hz⁻¹, that of the b_ν, ref = 1750 Jy deg⁻² simulated giants is L_ν = 0.6 × 10²⁶ W Hz⁻¹, and that of the b_ν, ref = 250 Jy deg⁻² simulated giants is L_ν = 0.2 × 10²⁶ W Hz⁻¹.¹⁸ Interestingly, the higher reference surface brightness median provides a better fit to the data. Even in case our Monte Carlo approach would predict luminosity densities that are biased low by a factor two, the higher median remains favoured.

In conclusion, it seems reasonable to suppose that our measurement b_ν, ref = 1.0 × 10³ Jy deg⁻² is not, or only mildly, biased high by selection effects. Still, we take a conservative approach and in setting priors we assume a 75% error on our measurement of b_ν, ref and a 50% error on our measurement of σ_ref. Thus, the priors for b_ν, ref and σ_ref – which we choose to be Gaussian – have 68% credible intervals [250 Jy deg⁻², 1750 Jy deg⁻²] and [0.65, 1.95]. We retain flat priors for ξ and $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$.

4.1.4. Inference

To compute the posterior distribution for ξ, $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$, b_ν, ref and σ_ref, we first brute-force evaluated the likelihood function over a regular grid that covers a total of 3.3 million parameter combinations¹⁹. For each proposed parameter quartet, we computed the PDF of L_p, obs | L_p, obs > l_p, GRG, and obtained the likelihood assuming that the LoTSS DR2 GRG projected proper lengths are IID draws from it. To obtain the PDF, we successively evaluated Eqs. (22), (25), (29), (18), (3), and (21), alongside their direct dependencies. This required the numerical evaluation of five integrals. Compared with using Riemann sums, we achieved substantial accuracy improvements at virtually no added numerical cost by approximating these integrals with the trapezoid rule and the composite Simpson’s rule.

We summarise the likelihood function in Table E.1 and Fig. E.1. To obtain the posterior, we simply multiplied the likelihood function by the prior and normalised the result.

In Table 4, for each parameter, we list the maximum a posteriori probability (MAP) estimate, alongside estimates for the posterior mean and standard deviation. In Fig. 16, we visualise all one- and two-dimensional posterior marginals, in which we mark the MAP (white circle) and the posterior mean (white cross). The joint marginal for ξ and b_ν, ref shows that these parameters have a strong negative correlation, indicating that with current data, the steep slope of the ESF at high l_p can equally be described with a steep intrinsic slope and mild surface brightness selection (i.e. ξ low and b_ν, ref high), or by a shallow intrinsic slope and strong surface brightness selection (i.e. ξ high and b_ν, ref low). We leave it up to future studies to break this degeneracy, either by using larger samples, by measuring b_ν, ref directly, or by improving survey sensitivities so that surface brightness selection effect modelling becomes superfluous altogether.

4.1.5. Goodness of fit

In both panels of Fig. 14, we compare the ESF and SF of L_p, obs | L_p, obs > l_p, GRG for z_max = 0.5, using the MAP parameters for the latter. The model appears able to produce a tight fit to the data. The mean and standard deviation of the ESF–SF residuals are 0.01% and 0.3%, whilst the mean and standard deviation of the absolute ESF–SF residuals are both 0.2%. Using a Kolmogorov–Smirnov test, we formally verified that our best parameters are indeed good parameters – in the sense that they represent a plausible model underlying the data. The Kolmogorov–Smirnov statistic is the maximum deviation between the ESF and SF, and equals 1% in our case. The p-value – the probability that an ESF–SF discrepancy of at least this magnitude would occur if the SF represents the true underlying distribution – is p = 99%. For any reasonable significance level, we do not reject the null hypothesis. The model, given our best parameters, indeed represents a possible description of the data. We conclude that the distribution of GRG intrinsic proper lengths, after correcting for selection effects, is consistent with a single Pareto distribution with tail index ξ = −3.5. We show the SF of this distribution in both panels of Fig. 14 (fading green dots). For low l_p, the observed slope is shallower (due to angular length selection), whilst for high l_p, the observed slope is steeper (due to surface brightness selection).

A quasi-Pareto distribution can arise naturally as the tail of a lognormal distribution (e.g., Malevergne et al. 2011), and there are reasons to believe that the entire radio galaxy length distribution is indeed approximately lognormal (Oei et al., in prep.). This provides an explanation of the approximately Paretian nature of the giant radio galaxy length distribution found in this section. The specific value of the tail index ξ is set by both the physics of radio galaxy growth and the distribution of radio galaxies over large-scale environments, the latter of which we measure in Oei et al. (in prep.). Our result ξ = −3.5 ± 0.5 is a new constraint for dynamical models such as those of Turner & Shabala (2015) and Hardcastle (2018).

4.2. Giant radio galaxy number density

If in addition to our discoveries, we know how many giants our search campaign has missed, then we can infer the true comoving GRG number density in the Local Universe. The posterior distribution over selection effect parameters $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$, b_ν, ref and σ_ref induces a probability distribution over the search completeness function C(l_p). C(l_p) denotes the probability that a giant of projected proper length l_p in comoving space up to z = z_max is detected through the search. We first generated parameter samples from our posterior using rejection sampling, and then used each to calculate a C(l_p) sample. We show the distribution over C(l_p) for z_max = 0.5 in Fig. 17. For small l_p, C is low as many giants drop out due to angular length selection; for large l_p, C is low as many giants drop out due to surface brightness selection. The completeness peaks around l_p ∼ 2 Mpc; however, even there the majority of giants remains undetected.

We inferred a probability distribution over the true comoving GRG number density n_GRG by combining Eqs. (18) and (30) with the LoTSS DR2 GRG catalogue and samples from Sect. 4.1.4’s posterior. The resulting skewed distribution, with mean and SD n_GRG = 4.6 ± 2.4(100 Mpc)⁻³ and 80% credible interval 3.1–6.7(100 Mpc)⁻³, is shown in Fig. 18. We note that, although the uncertainty in b_ν, ref induces a large uncertainty in C from l_p ∼ 1.5 Mpc onwards, the completeness uncertainty at large projected lengths does not substantially contribute to the uncertainty in n_GRG. This is because the GRG population is dominated by smaller giants, for which the completeness appears better constrained.

What picture arises regarding the abundance of giant radio galaxies in the Local Universe’s large-scale structure? If we model the Cosmic Web through comoving cubic unit cells (Oei et al. 2022b) with 50 Mpc sides, and each cubic unit cell contributes one cluster and three filaments, then a cube with 100 Mpc sides features a total of eight clusters and 24 filaments. For comparison, in a (100 Mpc)³ volume up to z_max = 0.5, the SDSS-III cluster catalogue of Wen et al. (2012) contains on average 11.2 clusters of any mass, and 4.5 clusters of mass M₂₀₀ > 10¹⁴ M_⊙. Since clusters contain ∼20% of giants (Oei et al., in prep.), we find the average number of giants per cluster to be ∼10⁻¹. If one assumes that filaments contain the remaining ∼80% of giants, and uses the fact that the average number of filaments per cluster is of order unity, it follows that the average number of giants per filament is also ∼10⁻¹. In all likelihood, most clusters and filaments do not currently contain a giant.

4.3. Giant radio galaxy lobe volume-filling fraction

Because giant radio galaxies enrich the IGM with hot plasma and magnetic fields far beyond the circumgalactic media of their hosts, they may provide a meaningful contribution to the heating and magnetisation of – in particular – the most rarefied parts of the filament IGM. By combining the GRG number density and the GRG jet power distribution (e.g., Dabhade et al. 2020a), one could estimate the instantaneous heating and magnetisation contributions directly. We recommend such analysis for future research.

We evaluated Eq. (35) to obtain an estimate of the fraction of the Local Universe’s proper volume that GRG lobes occupy. We used Alcyoneus as a reference giant, for which V = 2.5 ± 0.3 Mpc³ and l_p = 4.99 ± 0.04 Mpc (Oei et al. 2022a); this suggests 𝔼[Υ_p]≈2%. Future work should determine whether Alcyoneus’s case is typical, as observations, such as those shown in Figs. 12–13, suggest that giants exhibit a large variety of shapes – and thus total lobe volume–cubed length ratios. Interestingly, simulations by Krause et al. (2012) have found that these shapes also depend on environmental parameters such as ambient pressure and density. Truncating the GRG projected length distribution at l_p, max = 7 Mpc, so that its support is exactly an order of magnitude, Eq. (35) predicts ${\mathrm{VFF}}_{\mathrm{GRG}}(z=0)=3_{-1}^{+4}\times {10}^{-7}$.²⁰

Whether this result is sensitive to changes in l_p, max depends on ξ, with ξ = −3 being a special value under self-similar growth. In that case, small and large giants contribute equally to VFF_GRG: although large giants are rarer ($f_{L_{\text{p}}}\propto l_{\text{p}}^{-3}$), their larger lobe volumes ($V\propto l_{\text{p}}^3$) exactly compensate. For ξ < −3, small giants provide the dominant contribution to VFF_GRG and the choice of l_p, max can be irrelevant; for ξ > −3, large giants dominate and the choice of l_p, max always matters.

If we assume that giants occur in clusters and filaments only and use the fact that clusters and filaments comprise about 5% of the Local Universe’s volume (Forero-Romero et al. 2009), then the GRG lobe VFF within clusters and filaments specifically is ${\mathrm{VFF}}_{\mathrm{GRG}}(z=0)=5_{-2}^{+8}\times {10}^{-6}$. We conclude that, at each given moment, GRG lobes occupy just a small fraction of the WHIM and ICM. If the enrichment of the IGM by giants is to affect the WHIM and ICM on a large scale, mixing processes in the IGM are necessary and many galaxies must be able to form giants at some point in their evolution.

4.4. Unification model constraints from quasar and non-quasar giants

In Sect. 2.11, we have predicted general ramifications of the basic unification model on a GRG sample. We constrained and tested this model with our LoTSS DR2 GRG sample.

First, for all 3198 giants with definitively identified hosts, we queried the SDSS DR12 spectral class spCl: a Boolean label indicating whether or not the host contains a quasar. As many hosts have no SDSS DR12 spectrum, or even fall outside of SDSS DR12 coverage, we retrieved host classifications for just 1442 of these giants (45%). Of these classified giants, 318 are quasar giants (22%) (of which 45 (14%) are discoveries presented in this work) and 1124 are non-quasar giants (78%) (of which 876 (78%) are discoveries presented in this work). Therefore, the apparent LoTSS DR2 quasar GRG fraction $f_{\mathrm{Q}}=\frac{45}{45+876}=5\%$. However, spectral class labels are preferentially available for GRG hosts with higher optical flux densities, such as those at low redshifts or those containing quasars, because they are more probable spectroscopic targets. Through Fig. 19, we demonstrate that the fraction of observed GRG hosts with spectral class labels indeed decreases with redshift, whilst the fraction of quasar identifications increases. For each redshift interval, in white, we denote the fraction of quasar giants within the classified population. If spectral class labels would have been available for the non-classified observed populations as well, the quasar GRG fractions would probably have been lower. Assuming that all observed giants whose hosts have an unknown spectral class are non-quasar giants, we find quasar GRG fractions f_Q = 5%, 4%, 8%, 16%, and 28%, for redshift intervals 0–0.2, 0.2–0.4, 0.4–0.6, 0.6–0.8, and 0.8–1, respectively. The true quasar GRG fraction for a given redshift interval might still differ from the aforementioned observed quasar GRG fraction: namely, if selection effects make a given quasar GRG easier (or harder) to detect than a given non-quasar GRG. At higher redshifts, quasar giants certainly appear easier to detect than non-quasar giants, as hosts without quasars often become too faint to optically identify. For the lowest redshift intervals, this problem does not exist, and we therefore consider the observed quasar GRG fraction f_Q = 5% to be closest to the true one. A thorough analysis of the impact of selection effects on f_Q is a topic for future research.

Using Eq. (37), and assuming a quasar GRG probability p_Q = 5%, we found a maximum inclination angle ${\theta }_{\max }={39}_{-3}^{+2}°$. For p_Q ∼ Uniform  (4%,6%), the result remained the same. In conclusion, if the basic unification model considered in Sect. 2.11 is correct, then observations of giants predict that quasars are AGN seen along lines-of-sight that make an angle of at most ${\theta }_{\max }={39}_{-3}^{+2}°$ with the black hole rotation axis.

Finally, we tested whether the RV L_p, obs | L_p, obs ≥ l_p, GRG has the same distribution for quasar and non-quasar giants, as predicted by the unification model. In the top panel of Fig. 20, we show PDFs approximated through kernel density estimation (KDE) for newly discovered (LoTSS DR2) SDSS-classified quasar giants and non-quasar giants up to z_max = ∞. Despite the small number of quasar giants, N_Q = 45, the two-sample Kolmogorov–Smirnov (KS) test yielded a low p-value p < 1‰: we rejected the null hypothesis that these distributions stem from a single underlying one. However, this does not mean that the unification model hypothesis should immediately be rejected as well, because the selective availability of SDSS DR12 spectral class labels induces a severe selection effect²¹. This effect can be tempered by choosing a lower z_max, but for choices such as z_max = 0.25 and 0.5, just N_Q = 10 and 16 quasar giants remain – too few to extract meaningful information. The bottom panel of Fig. 20 again shows observed projected proper length KDE PDFs, but now for all SDSS-classified quasar giants and non-quasar giants up to z_max = ∞; in this case, N_Q = 318. Again, the two-sample KS test yields p < 1‰, but this time quasar giants appear smaller than non-quasar giants. Because we aggregated samples here that have different selection effects imprinted, it becomes hard to draw clear conclusions. Limiting z_max reduces the severity of most selection effects, but of course comes at the cost of reducing the sample size. Interestingly, the corresponding quasar GRG and non-quasar GRG projected length distributions do become more alike; for instance for z_max = 0.25 and 0.5, the two-sample KS test yields p = 4% and 1%, respectively. In conclusion, because quasar giants are intrinsically rare and the availability of spectral class labels is biased towards low redshifts and hosts containing quasars, it is challenging to robustly test the unification model with current GRG observations. We refrain from drawing final conclusions, and recommend a careful future analysis.

5. Discussion

5.1. Radio galaxy length definitions

How large are radio galaxies? Despite the simplicity of this question and more than half a century of research on radio galaxies, their intrinsic length distribution has not yet been rigorously characterised. In this work, we have carried out the first precision analysis of the tail of the radio galaxy intrinsic length distribution. Precision analyses tend to raise questions; firstly, whether the studied observable is well defined, and secondly whether one could conceive of more informative observables – that is to say those that make it easier to reveal underlying physical mechanisms. This work’s main observable is the radio galaxy projected proper length; we argue that it is neither well defined nor maximally informative.

5.1.1. The current length definition: survey-dependence

Contemporary research uses a survey-dependent definition for radio galaxy angular lengths, which then makes projected proper lengths survey-dependent too.

The angular length is canonically defined as the largest possible angular separation between two directions for which the RG’s specific intensity function I_ν, RG exceeds b_ν, th: some specified factor of order unity times the image noise σ_Iν. A complication is that not only I_ν, RG, but also σ_Iν varies with observing frequency; the latter because of observational factors such as (u,v)-coverage, bandpass, radio-frequency interference (RFI), ionospheric weather, the sky density of bright calibrators and the performance of calibration algorithms. I_ν, RG additionally depends on resolution, at least for point sources; σ_Iν additionally depends on resolution and integration time. As a result, both I_ν, RG and the concrete value of b_ν, th used in the angular length definition change from image to image. Each study thus far has therefore implicitly used a different definition for angular length, instead of a shared, absolute one. As the projected proper length follows from combining the angular length with the host redshift, it suffers from the same problem.

Whether the survey-dependence of the current length definition is problematic, depends on the radio galaxy. For archetypal FRII RGs, the angular length is roughly equal to the angular distance between the hotspots. These are an FRII RG’s brightest morphological components, and are thus the first to be picked up by a survey. In contrast, archetypal FRI RGs, which gradually fade with distance from the host, can have significantly larger angular lengths in surveys of higher sensitivity. The giants in the middle-right and bottom-left panel of Fig. 13 are good examples: these radio galaxies were known before the LoTSS DR2, but were not known to be giants; similarly, more sensitive surveys are poised to assign them even larger extents. If we are to move towards precision science, it therefore makes sense – at least for FRI RGs – to more explicitly recognise that the angular and projected proper lengths are functions of the observing frequency ν_obs and a surface brightness threshold b_ν, th. If catalogues would explicitly state for what combination (ν_obs,b_ν, th) they provide ϕ = ϕ(ν_obs,b_ν, th) and l_p = l_p(ν_obs,b_ν, th), it is possible to homogenise a collection of data sets by using universal angular and projected proper length definitions for all RGs²². If two angular lengths have been measured for the same RG, for instance ϕ₁ at (ν_obs, 1,b_{ν, th, 1}) and ϕ₂ at (ν_obs, 2,b_{ν, th, 2}), then we can estimate ϕ for any desired (ν_obs,b_ν, th) through interpolation. For example, the interpolation formula for a symmetric radio galaxy with jets or lobes of constant spectral index α, and with a specific intensity function contribution which fades to zero linearly with angular distance from the host, is

$$\begin{array}{c}\hfill \varphi =max\{{\varphi }_1+\frac{b_{\nu ,\mathrm{th}}-b_{\nu ,\mathrm{th},1}{\left(\frac{{\nu }_{\mathrm{obs}}}{{\nu }_{\mathrm{obs},1}}\right)}^{\alpha }}{b_{\nu ,\mathrm{th},2}{\left(\frac{{\nu }_{\mathrm{obs}}}{{\nu }_{\mathrm{obs},2}}\right)}^{\alpha }-b_{\nu ,\mathrm{th},1}{\left(\frac{{\nu }_{\mathrm{obs}}}{{\nu }_{\mathrm{obs},1}}\right)}^{\alpha }}({\varphi }_2-{\varphi }_1),0\}.\end{array}$$(41)

5.1.2. The ideal length definition: physical relevance

Alcyoneus, shown in Fig. 12’s top-right panel, is a 5 Mpc giant whose ageing lobes are revealed for the first time by the LoTSS DR2 (Oei et al. 2022a). Future image sensitivity improvements shall reveal more hitherto unseen, fading lobes around known RGs that formed in the aftermath of earlier AGN activity episodes. Could a large fraction of sufficiently old RGs turn out to be giants, once such sensitivity improvements start providing evidence of earlier and earlier AGN activity episodes?

Up to now, it has been informative to include all visible plasma in the angular length measurement. In future images, some of the visible plasma might be of such low pressure that it has become physically insignificant, in the sense that it does not affect the thermodynamics of the surrounding IGM anymore; we propose to exclude such plasma from a radio galaxy length. Practically, one option is to introduce an absolute threshold: for example to include plasma of pressure P ≥ 10⁻¹⁷ Pa only – this is the pressure of the warm–hot intergalactic medium (WHIM) in the ρ_BM ∼ 10 Ω_BM, 0ρ_c, 0 and T ∼ 5 × 10⁵ K regime. Another option is to introduce an environment-dependent pressure threshold; this would mean that a cluster RG sees its length measured against a higher pressure threshold than a filament RG, because its plasma becomes thermodynamically irrelevant sooner. A problem with (equipartition or minimum energy) pressure–based length definitions is that pressure is a derived quantity: the specific intensity only determines the product of pressure and line-of-sight length through the lobe.

5.2. Moving beyond line segment projection

In this work, we have adopted a classical approach to treating projection, by modelling radio galaxies as line segments. When a radio galaxy’s inclination angle θ is close to 0° or 180°, this approach predicts that the angular length ϕ and thus the projected proper length l_p vanishes. In reality however, because radio galaxy lobes have non-zero volumes, the RV L_p will never tend to zero. Figure 21 illustrates this point through three sources interpreted to be radio galaxies aligned closely with the line of sight. In each case, ϕ remains of arcminute scale.

A more realistic approach, outlined in Appendix A.3, moves beyond the simplistic line segment geometry by adding two lobes of radius R to the endpoints of the line segment, whose length is given by RV L. The ratio $\frac{L}{L_{\mathrm{p}}}$, which was unbounded under the classical approach, is now at most $1+\frac{1}{\eta }$, where $\eta :=\frac{2R}{L}$. (The classical approach simply is the limit η = 0.) We suggest a recalculation of this work’s theoretical and applied results under this more realistic RG geometry as a direction for future research. To obtain the applied results, one must either fix η as a hyperparameter, or include it as an additional model parameter. We note that even using conservatively low values of η, such as $\eta =\frac{1}{10}$, will represent an improvement in realism over η = 0.

5.3. Unmodelled selection effects

In this work, we have modelled both an angular length and a surface brightness selection effect. Several other plausible selection effects have not been included in the forward model, as we have judged each to be of minor importance. However, in unison, they could have a non-negligible influence on the distribution of the observed projected proper GRG length RV L_p, obs | L_p, obs > l_p, GRG. Some of their influence might have been absorbed by the parameters of the included selection effects – that is by $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$, b_ν, ref, and σ_ref – or, worse still, by ξ.

One of these unmodelled selection effects is that RGs whose axes are oriented almost parallel to the line of sight are more likely to be rejected from a sample than RGs whose axes are closer to the plane of the sky, as the former do not always have a characteristic double-lobe appearance. For instance, perhaps not all readers would regard our identification of the sources in the left column of Fig. 21 as RGs convincing. Nevertheless, as remarked in Sect. 2.11, conditioning L_p, obs | L_p, obs > l_p, GRG on inclination angle does not affect its distribution. As a result, this selection effect does not necessitate forward model modifications.

Furthermore, as shown in Sect. 4.4, there is a selection effect at play that favours the selection of quasar giants over non-quasar giants at high redshift, as host galaxies with quasars are more luminous in the optical and therefore have a better chance to be picked up in optical imagery. This selection effect will get less severe once deeper photometric surveys become available.

We encountered three more selection effects during our LoTSS DR2 GRG search. The larger an RG – and especially an FRII RG – becomes, the harder it is for an observer to identify its host galaxy, as an increasing number of plausible host candidates can lie interspersed in the strip of sky between the lobes. The severity of this effect, which diminishes the prevalence of the largest giants in a sample, depends on the balance chosen between avoiding false discoveries and avoiding rejections of true discoveries. Another selection effect runs against RGs in galaxy clusters. Such environments can contain multiple adjacent RGs, making it at times unclear which lobe belongs to which RG. When no confident double-lobe associations can be made, the RGs involved fail to make it into the sample. If galaxy clusters contain primarily smaller giants, this selection effect induces a bias against smaller giants. A final unmodelled bias runs against RGs at the end of their life cycle. Once the AGN stops launching jets for a prolonged period, it becomes hard to identify the host galaxy, which will no longer present as a bright, compact radio source. This effect preferentially deselects larger giants, which are even more likely to approach the end of their life than smaller giants.

5.4. Does the choice of prior matter?

We have taken a conservative approach to constraining the posterior distribution through the prior: we have left tail index ξ and angular length selection half-width $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$ fully unconstrained, and have adopted wide Gaussian priors for reference surface brightness parameters b_ν, ref and σ_ref, despite measuring them explicitly under assumptions. We provide summary statistics of the posterior in Table 4 and visualise its one- and two-parameter marginals in Fig. 16. Does our choice of prior significantly affect the inferences? To explore this question, we chose a different reasonable prior and compared results. One such prior is the fully uniform prior, which equivalises the posterior and the likelihood function. We provide analogous summary statistics of the likelihood function in Table E.1 and visualise analogous marginals in Fig. E.1. Reassuringly, no statistically significant parameter changes occur. In particular, for z_max = 0.5, ξ = −3.5 ± 0.5 becomes ξ = −3.4 ± 0.5 upon changing to the uniform prior; for z_max = 0.25, ξ = −3.5 ± 0.4 even remains the same. However, given the strong likelihood degeneracy between ξ and b_ν, ref apparent in Fig. E.1, more stringent priors on b_ν, ref are able to meaningfully shift ξ’s posterior mean. Such priors shall be appropriate only after studying the surface brightness properties of large radio galaxies – and the associated selection effect – in more detail.

5.5. Cosmological evolution of the GRG length distribution

In this work we have assumed that the parameter ξ, which fully characterises the intrinsic GRG length distribution under the ansatz of Paretianity, remains constant throughout the chosen redshift range [0, z_max]. To test this assumption, we have analysed our LoTSS DR2 GRG sample in Sect. 4 up to both z_max = 0.5 and z_max = 0.25. For z_max = 0.5, we included giants that existed in the last 5 Gyr of the Universe’s history, and found ξ = −3.5 ± 0.5. Meanwhile, for z_max = 0.25, we included giants that existed in the last 3 Gyr of the Universe’s history only, and found the very similar ξ = −3.5 ± 0.4. Thus, our analysis did not produce evidence that ξ evolves over cosmic time. However, given the large error bars, a modest time evolution cannot be excluded. Furthermore, the data sets that underlie these inferences are not disjoint: the 811 giants that inform the lower-maximum-redshift analysis make up 55% of the 1473 giants that inform the higher-maximum-redshift analysis.

Whether a time evolution of ξ is expected is presumably tied to whether giant radio galaxy growth varies with environmental density at a given epoch, because the combined effects of the Universe’s expansion and ongoing large-scale structure formation can similarly change environmental density. Combining our LoTSS DR2 GRG sample with Cosmic Web reconstructions to explore giant growth as a function of environmental density is the topic of a forthcoming work (Oei et al., in prep.). Interestingly, the recent exploration by Lan & Prochaska (2021) that compared the environments of giants and non-giants did not find significant differences.

The most straightforward model extension is to again assume that RG lengths are Pareto distributed with tail index ξ, but now ξ = ξ(z). This function’s first-degree Maclaurin polynomial, which provides the linearisation at the present day, is

$$\begin{array}{c}\hfill \xi (z)\approx \xi (z=0)+\frac{\mathrm{d}\xi }{\mathrm{d}z}(z=0)·z.\end{array}$$(42)

One would adopt ξ(z = 0) and $\frac{\mathrm{d}\xi }{\mathrm{d}z}(z=0)$ as the parameters of interest, replacing what used to be a constant ξ; the number of model parameters would thus increase by one. However, an attempt to infer the cosmic evolution of ξ appears promising only once the major selection effects are better constrained.

6. Conclusions

In this work, we have performed Bayesian inference on a LoTSS-derived sample of 2060 giant radio galaxy projected proper lengths, using a one-parameter model that assumes a spatially homogeneous, non-evolving population of radio galaxies with stick-like geometry, Pareto-distributed lengths, and isotropic inclination angles. Before fitting to data, we extended the forward model with two selection effects typical of contemporary manual GRG search campaigns. The best-fit survival function tightly reproduces the empirical one, leaving permille-scale absolute residuals. Having quantified the most important selection effects, we estimated the true comoving giant radio galaxy number density in the Local Universe.

1. We developed an analytical model through which statistical questions about radio galaxy (RG) lengths can be rigorously answered. In the current work, we applied this model to giant radio galaxies. We adopted the ansatz that the RG intrinsic proper length L, as measured in three spatial dimensions, is a random variable (RV) with a Pareto Type I distribution (i.e. a simple power-law distribution) characterised by tail index ξ. Next, by assuming that RGs have no preferential orientation with respect to the observer, we derived the distribution of the RG projected proper length L_p. By conditioning, one obtains the version relevant for giants, L_p | L_p > l_p, GRG (where we chose l_p, GRG := 0.7 Mpc). This RV is again Paretian, with the same tail index ξ. In summary, for giant radio galaxies, projection retains Paretianity. Finally, observers face selection effects; we modelled the observed projected proper length L_p, obs by considering an angular length threshold selection effect and a surface brightness selection effect. The angular length threshold selection effect assumes a linearly increasing angular-length-dependent probability of sample inclusion around a particular pre-defined threshold, meant to emulate manual visual searches that only target RGs of some angular length and above. The surface brightness selection effect assumes that giants are self-similar, and that their lobes have lognormally distributed surface brightnesses which must be above-noise to secure sample inclusion. The GRG observed projected proper length L_p, obs | L_p, obs > l_p, GRG again follows through conditioning. We assumed our data to be realisations of this RV.

2. The model also yielded explicit expressions for the (posterior) distribution of L | L_p = l_p. This allows one to deproject RGs in a statistical sense, providing the intrinsic proper length given the projected proper length in the limit of negligible selection effects. We also present practical expressions for the mean and variance of L | L_p = l_p. To unravel the driving factors that allow some RGs to become giants, most authors search for correlations between host or environmental physical parameters and the GRG projected length. However, if there is a causal chain that connects host or environmental parameters to GRG length, the connection will be to the intrinsic length; the observer’s vantage point does not play a role in the physics. Therefore, the projection effect merely serves as a multiplicative noise source. We suggest that future analyses should recognise the projection effect as such, and correlate host or environmental parameters with the intrinsic, rather than projected, proper length – using statistical deprojection.

3. Through a manual visual search of LoTSS DR2 pipeline products, which are part of the LOFAR’s Northern Sky survey at 144 MHz, we discovered a population of 2060 previously unknown giants. This is the largest single contribution to the literature yet, and increases the community-wide census by a factor 2.6. We present 11 discoveries with l_p ≥ 4 Mpc, 53 with 3 ≤ l_p < 4 Mpc, 291 with 2 ≤ l_p < 3 Mpc, 1215 with 1 ≤ l_p < 2 Mpc, and 490 with 0.7 ≤ l_p < 1 Mpc. Our study extends the known breadth of the giant radio galaxy phenomenon. Among the findings are both the giant hosted by J081956.41+323537.6 and Alcyoneus (Oei et al. 2022a), at l_p = 5.1 Mpc and l_p = 5.0 Mpc the projectively largest giants ever found. We discover that multi-Mpc radio galaxies can be generated before redshift 1, despite the Universe’s mean density being an order of magnitude higher, and by spiral galaxies, whose stellar masses are typically an order of magnitude lower than those of ellipticals. We discover giants whose hosts have a record-low stellar mass M_⋆ = 4.8 × 10¹⁰ M_⊙. We also discover giants whose hosts have a record-high supermassive black hole mass M_• ≳ 5 × 10¹⁰ M_⊙; interestingly, with l_p = 0.8 Mpc, one of these giants is relatively small. We more than double the number of known giants with angular lengths exceeding that of the Moon; one discovery, at 2.2°, is the angularly largest radio galaxy in the Northern Sky and the angularly largest giant overall. Excitingly, our LoTSS DR2 search has been far from exhaustive: many thousands of readily identifiable giants still await discovery in this public data set.

4. Using our LoTSS DR2 GRG sample up to z_max = 0.5, we generated a posterior distribution over ξ and three selection effect parameters. Our model provides an excellent fit to the data, with absolute residuals being on average 2‰. We inferred that the intrinsic proper length distribution of the largest radio galaxies resembles a Pareto distribution with tail index ξ = −3.5 ± 0.5. Our analysis did not yield evidence for an evolving ξ in the last 5 Gyr of cosmic time.

5. The selection effect parameters estimated through the posterior are far from nuisance parameters, as they allowed us to statistically undo the selection effects imprinted on our LoTSS DR2 GRG data. As a result, we could for the first time estimate the true comoving giant radio galaxy number density n_GRG in the Local Universe up to z_max. We relied on the crucial assumption that the surface brightness distribution of RGs with intrinsic length l_ref = 0.7 Mpc at z = 0 and frequency ν_obs = 144 MHz is unimodal – and lognormal in particular. We furthermore assumed a lobe spectral index α = −1 and self-similar growth. We found n_GRG(l_p, GRG = 0.7 Mpc, z_max = 0.5) = 5 ± 2(100 Mpc)⁻³. The implication is that giant radio galaxies are truly rare – not only from a current observational perspective, but also from a cosmological one. Current GRG lobes occupy just a few millionths of the IGM volume. At any given moment in time, most clusters and filaments – the building blocks of modern large-scale structure – do not harbour giants.

Giants embody the most extreme known mechanism by which galaxies can affect the Cosmic Web around them. Whereas this work has explored the geometric properties of giants, a thorough exploration of their Cosmic Web energisation and magnetisation potential is a future frontier. Excitingly, the interactions between giants and the ethereal intergalactic medium may also allow for new constraints on the thermodynamics in filaments.

---

Acknowledgments

M.S.S.L. Oei warmly thanks Frits Sweijen for coding the very useful https://github.com/tikk3r/legacystamps and for tireless ICT advice, and Aleksandar Shulevski for comments that have improved the manuscript. In dear memory of Nicoline. By staying strong despite life’s challenges, you showed what it means to be a giant. 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 Netherlands Organisation for Scientific Research (NWO). A. Drabent acknowledges support by the BMBF Verbundforschung under the grant 05A20STA. The LOFAR is the Low-Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT Foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS–INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF–NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; the Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; the Istituto Nazionale di Astrofisica (INAF), Italy. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of the University of Turin (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard–Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eötvös Loránd University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing–Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 11433005). The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

References

Appendix A: Framework derivations and details

A.1. Intrinsic proper length

Let f_L : ℝ → ℝ_≥0 be the PDF of the distribution of intrinsic (i.e. 3D) proper (i.e. not comoving) radio galaxy lengths L. We assume that f_L follows a power law between l_min and l_max:

$$\begin{array}{c}\hfill f_L\left(l\right)=\{\begin{array}{cc}{f}_L\left(l_{\mathrm{ref}}\right){\left(\frac{l}{l_{\mathrm{ref}}}\right)}^{\xi }\hfill & l\in [l_{\min },l_{\max }],\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$(A.1)

Then, when ξ ≠ −1,

$$\begin{array}{cc}\hfill 1& ={\displaystyle {\int }_{\mathbb{R}}{f}_L\left(l\right)\mathrm{d}l={\int }_{l_{\min }}^{l_{\max }}{f}_L\left(l_{\mathrm{ref}}\right){\left(\frac{l}{l_{\mathrm{ref}}}\right)}^{\xi }\mathrm{d}l}\hfill \\ \hfill & =\frac{f_L\left(l_{\mathrm{ref}}\right)l_{\mathrm{ref}}}{\xi +1}({\left(\frac{l_{\max }}{l_{\mathrm{ref}}}\right)}^{\xi +1}-{\left(\frac{l_{\min }}{l_{\mathrm{ref}}}\right)}^{\xi +1}).\hfill \end{array}$$(A.2)

This equation provides the normalisation factor f_L(l_ref) given l_min, l_max, ξ and an arbitrary choice for l_ref ≠ 0.

When l_max = ∞, L has a Pareto Type I distribution and we must have ξ < −1 for the integrals of Eq. A.2 to converge. Throughout the remaining analysis, for the sake of simplicity, we assume l_max = ∞ and ξ < −1 and choose l_ref = l_min. Then $f_L\left(l_{\mathrm{ref}}\right)=f_L\left(l_{\min }\right)=-\frac{\xi +1}{l_{\min }}$. The corresponding CDF is

$$\begin{array}{c}\hfill F_L(l)=\{\begin{array}{cc}0\hfill & \text{if}l\le l_{\min }\hfill \\ 1-{\left(\frac{l}{l_{\min }}\right)}^{\xi +1}\hfill & \text{if}l>l_{\min }.\hfill \end{array}\end{array}$$(A.3)

In reality, radio galaxies cannot become arbitrarily long (Hardcastle 2018): at some distance from the host, all energy carried by the jets will have been radiated away, used to perform work on the IGM, transferred to CMB photons through inverse Compton scattering, or converted into heat. However, at the moment of writing, the implied maximum length l_max remains ill-constrained. As we see throughout Appendix A, a major advantage of simple model assumptions is that explicit and thus insightful analytic expressions can be derived. Such easy-to-evaluate expressions complement the results of expensive numerical simulations, which aim to maximise realism; our aim to maximise insight is best served by setting l_max = ∞. From Eq. A.3, we see that for l_min = 0.7 Mpc and ξ = −3.5, the tail of the intrinsic proper length distribution (l > 5 Mpc) contains less than 1% of probability. Thus, for most framework applications, our choice l_max = ∞ is unproblematic — except when higher powers of L are involved. For example, the volume-filling fraction calculations of Appendix A.10 necessitate considering L³; for realistic ξ, results exist only for finite l_max.

Upon relabelling ξ → −α − 1, l_min → k, L → X and l → x, one obtains the literature’s most common form of the Pareto Type I PDF:

$$\begin{array}{c}\hfill f_X\left(x\right)=\{\begin{array}{cc}\frac{\alpha k^{\alpha }}{x^{\alpha +1}}\hfill & x\ge k,\hfill \\ 0\hfill & x<k.\hfill \end{array}\end{array}$$(A.4)

A.2. Projected proper length

A.2.1. Distribution for RGs

Let f_Lp : ℝ → ℝ_≥0 be the PDF of the distribution of projected proper radio galaxy lengths. The PDF f_Lp follows from the associated CDF F_Lp : ℝ → [0, 1] through differentiation; F_Lp and f_L relate through

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)& :=\mathbb{P}(L_{\mathrm{p}}\le l_{\mathrm{p}})={\displaystyle {\int }_{\mathbb{R}}\mathbb{P}(Lsin\mathrm{\Theta }\le l_{\mathrm{p}}|L=l)f_L\left(l\right)\mathrm{d}l}\hfill \\ \hfill & ={\displaystyle {\int }_{\mathbb{R}}\mathbb{P}(sin\mathrm{\Theta }\le \frac{l_{\mathrm{p}}}{l})f_L\left(l\right)\mathrm{d}l.}\hfill \end{array}$$(A.5)

We note that F_Lp(l_p) vanishes for l_p ≤ 0: f_L(l) vanishes when l is negative, whilst $\mathbb{P}(sin\mathrm{\Theta }\le \frac{l_{\mathrm{p}}}{l})$ vanishes when l is positive. Clearly, the interesting case is l_p > 0. We can differentiate between two cases: l_p ≤ l_min, and l_p > l_min. In the first case, considering that f_L has support from l_min onwards only, we have $\frac{l_{\mathrm{p}}}{l}\le 1$, and

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)& ={\displaystyle {\int }_{l_{\min }}^{\infty }(\mathbb{P}(\mathrm{\Theta }\le arcsin\frac{l_{\mathrm{p}}}{l})}\hfill \\ \hfill & +\mathbb{P}(\mathrm{\Theta }\ge \pi -arcsin\frac{l_{\mathrm{p}}}{l}))\frac{-(\xi +1)}{l_{\min }}{\left(\frac{l}{l_{\min }}\right)}^{\xi }\mathrm{d}l\hfill \\ \hfill & =-\frac{\xi +1}{l_{\min }}{\displaystyle {\int }_{l_{\min }}^{\infty }2\mathbb{P}(\mathrm{\Theta }\le arcsin\frac{l_{\mathrm{p}}}{l}){\left(\frac{l}{l_{\min }}\right)}^{\xi }\mathrm{d}l}\hfill \\ \hfill & =-\frac{\xi +1}{l_{\min }}{\displaystyle {\int }_{l_{\min }}^{\infty }(1-\sqrt{1-{\left(\frac{l_{\mathrm{p}}}{l}\right)}^2}){\left(\frac{l}{l_{\min }}\right)}^{\xi }\mathrm{d}l}\hfill \\ \hfill & =1+(\xi +1){\displaystyle {\int }_{l_{\min }}^{\infty }\sqrt{{\left(\frac{l}{l_{\min }}\right)}^2-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^2}{\left(\frac{l}{l_{\min }}\right)}^{\xi -1}\frac{\mathrm{d}l}{l_{\min }}}\hfill \\ \hfill & =1+(\xi +1){\displaystyle {\int }_1^{\infty }\sqrt{{\eta }^2-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^2}{\eta }^{\xi -1}\mathrm{d}\eta .}\hfill \end{array}$$(A.6)

In the second case, we split up the integral in two:

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)& ={\displaystyle {\int }_{l_{\min }}^{l_{\mathrm{p}}}}\mathbb{P}(sin\mathrm{\Theta }\le \frac{l_{\mathrm{p}}}{l})f_L\left(l\right)\mathrm{d}l\hfill \\ \hfill & +{\displaystyle {\int }_{l_{\mathrm{p}}}^{\infty }}\mathbb{P}(sin\mathrm{\Theta }\le \frac{l_{\mathrm{p}}}{l})f_L\left(l\right)\mathrm{d}l\hfill \\ \hfill & ={\displaystyle {\int }_{l_{\min }}^{l_{\mathrm{p}}}}{f}_L\left(l\right)\mathrm{d}l+{\displaystyle {\int }_{l_{\mathrm{p}}}^{\infty }}(1-\sqrt{1-{\left(\frac{l_{\mathrm{p}}}{l}\right)}^2})f_L\left(l\right)\mathrm{d}l\hfill \\ \hfill & =1-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi +1}+{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi +1}\hfill \\ \hfill & +(\xi +1){\displaystyle {\int }_{l_{\mathrm{p}}}^{\infty }}\sqrt{{\left(\frac{l}{l_{\min }}\right)}^2-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^2}{\left(\frac{l}{l_{\min }}\right)}^{\xi -1}\frac{\mathrm{d}l}{l_{\min }}\hfill \\ \hfill & =1+(\xi +1){\displaystyle {\int }_{\frac{l_{\mathrm{p}}}{l_{\min }}}^{\infty }}\sqrt{{\eta }^2-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^2}{\eta }^{\xi -1}\mathrm{d}\eta \hfill \\ \hfill & =1+(\xi +1)\frac{\sqrt{\pi }}{4}{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi +1}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}.\hfill \end{array}$$(A.7)

In summary,

$$\begin{array}{cc}& F_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le 0\hfill \\ 1+(\xi +1){\int }_1^{\infty }\sqrt{{\eta }^2-{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^2}{\eta }^{\xi -1}\mathrm{d}\eta \hfill & \text{if}0<l_{\mathrm{p}}<l_{\min }\hfill \\ 1+(\xi +1)\frac{\sqrt{\pi }}{4}{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi +1}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}\hfill & \text{if}{l}_{\mathrm{p}}\ge l_{\min }.\hfill \end{array}\hfill \end{array}$$(A.8)

Through differentiation,

$$\begin{array}{c}\hfill f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le 0\hfill \\ -\frac{\xi +1}{l_{\min }}\frac{l_{\mathrm{p}}}{l_{\min }}I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})\hfill & \text{if}0<l_{\mathrm{p}}<l_{\min }\hfill \\ \frac{{(\xi +1)}^2}{l_{\min }}\frac{\sqrt{\pi }}{4}{\left(\frac{l_{\mathrm{p}}}{l_{\min }}\right)}^{\xi }\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}\hfill & \text{if}{l}_{\mathrm{p}}\ge l_{\min }.\hfill \end{array}\end{array}$$(A.9)

where

$$\begin{array}{c}\hfill I(a,b):={\displaystyle {\int }_1^{\infty }}\frac{{\eta }^a\mathrm{d}\eta }{\sqrt{{\eta }^2-b^2}}\text{for}a<0,|b|<1.\end{array}$$(A.10)

Significantly, for l_p > l_min, the projected length distribution follows a power law in l_p with the same exponent as the power law for the intrinsic length distribution: $f_{L_{\text{p}}}\propto l_{\text{p}}^{\xi }$ (just as f_L ∝ l^ξ).²³ We compare f_L and f_Lp in Fig. 1.

A.2.2. Distribution for giants

To derive the projected length distribution for giants, we consider the distribution of the conditioned RV L_p | L_p > l_p, GRG:

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right):=& \mathbb{P}(L_{\mathrm{p}}\le l_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})\hfill \\ \hfill =& 1-\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})\hfill \\ \hfill =& 1-\frac{\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p}},L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})}{\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})}.\hfill \end{array}$$(A.11)

For l_p > l_p, GRG, this reduces to

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)=& 1-\frac{1-F_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)}{1-F_{L_{\mathrm{p}}}\left(l_{\mathrm{p},\mathrm{GRG}}\right)}.\hfill \end{array}$$(A.12)

Furthermore assuming l_p, GRG > l_min, we twice use the bottom expression of Eq. A.8 to obtain the final CDF expression. As before, the corresponding PDF follows through differentation. We find

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)=& \{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le l_{\mathrm{p},\mathrm{GRG}}\hfill \\ 1-{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi +1}\hfill & \text{if}{l}_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}.\hfill \end{array}\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill f_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)=& \{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le l_{\mathrm{p},\mathrm{GRG}}\hfill \\ -\frac{\xi +1}{l_{\mathrm{p},\mathrm{GRG}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi }\hfill & \text{if}{l}_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}.\hfill \end{array}\hfill \end{array}$$(A.14)

Thus, the associated survival function is

$$\begin{array}{c}\hfill \mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})={\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi +1}.\end{array}$$(A.15)

The mean projected proper length of giants follows from the PDF by direct computation:

$$\begin{array}{cc}\hfill \mathbb{E}[L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]:=& {\displaystyle {\int }_{-\infty }^{\infty }{f}_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(l_{\mathrm{p}}\right)·l_{\mathrm{p}}\mathrm{d}{l}_{\mathrm{p}}}\hfill \\ \hfill =& l_{\mathrm{p},\mathrm{GRG}}\frac{\xi +1}{\xi +2}.\hfill \end{array}$$(A.16)

A.3. Deprojection factor

Consider a radio galaxy (RG) with a projected proper length l_p. Let the inclination angle θ denote the angle between the RG’s central axis and the line of sight. The RG’s inclination angle, projected proper length and intrinsic proper length l relate via l_p = l sin θ. Switching to random variable (RV) notation by using capital letters, the intrinsic proper length (which is the most physically relevant quantity) follows from the projected proper length (which can be measured) and the inclination angle (which is typically unknown), according to

$$\begin{array}{c}\hfill L=\frac{1}{sin\mathrm{\Theta }}{L}_{\mathrm{p}}.\end{array}$$(A.17)

A.3.1. Without lobes

Calling the deprojection factor D := (sin Θ)⁻¹, we now calculate the distribution of D for $f_{\mathrm{\Theta }}\left(\theta \right)=\frac{1}{2}sin\theta ;\theta \in [0,\pi ]$. The result is a continuous univariate distribution without parameters supported on the semi-infinite interval (1,∞). Let F_D : ℝ → [0, 1] be the cumulative density function (CDF) of D. Then, for d > 1,

$$\begin{array}{cc}\hfill F_D\left(d\right)& :=\mathbb{P}(D\le d)=\mathbb{P}(sin\mathrm{\Theta }\ge \frac{1}{d})\hfill \\ \hfill & =\mathbb{P}(arcsin\frac{1}{d}\le \mathrm{\Theta }\le \pi -arcsin\frac{1}{d})\hfill \\ \hfill & =F_{\mathrm{\Theta }}(\pi -arcsin\frac{1}{d})-F_{\mathrm{\Theta }}(arcsin\frac{1}{d})\hfill \\ \hfill & =cosarcsin\frac{1}{d}=\sqrt{1-\frac{1}{d^2}}.\hfill \end{array}$$(A.18)

Meanwhile, F_D(d) = 0 for d ≤ 1. The quantile function $F_D^{-1}\text{}:[0,1]\to [1,\infty )$ follows from solving F_D(d) = p for d:

$$\begin{array}{c}\hfill F_D^{-1}\left(p\right)=\frac{1}{\sqrt{1-p^2}}.\end{array}$$(A.19)

Thus, the minimum factor is $F_D^{-1}p=0=1$, the median $F_D^{-1}(p=\frac{1}{2})=\frac{2}{\sqrt{3}}$, $F_D^{-1}(p=\frac{1}{2}\sqrt{2})=\sqrt{2}$, $F_D^{-1}(p=\frac{1}{2}\sqrt{3})=2$, and factors can grow arbitrarily large: $F_D^{-1}p\to \infty$ as p → 1. We conclude that half of all RGs have an intrinsic proper length more than $\frac{2}{\sqrt{3}}$ their projected proper length.

By differentiating F_D to d we obtain f_D : ℝ → ℝ_≥0, the probability density function (PDF) of D:

$$\begin{array}{c}\hfill f_D\left(d\right)=\{\begin{array}{cc}\frac{1}{d^2\sqrt{d^2-1}}\hfill & \mathrm{if}d>1;\hfill \\ 0\hfill & \mathrm{if}d\le 1.\hfill \end{array}\end{array}$$(A.20)

The mean of D is $\mathbb{E}\left(D\right)=\frac{\pi }{2}$; the variance of D is undefined, as the corresponding integral diverges. Because f_D(d) has no maximum, the mode is undefined too. The PDF and CDF of D are shown in the upper left and right panels of Fig. A.1, respectively.

A.3.2. With lobes

The simplest model that includes a pair of lobes approximates them as spheres of radius R, whose centres are connected by a line segment of length L. Regardless of the viewing angle, the spheres retain their size, in contrast to the line segment connecting them. Calling $\eta :=\frac{2R}{L}\in {\mathbb{R}}_{\ge 0}$, we have

$$\begin{array}{c}\hfill D:=\frac{L+2R}{Lsin\mathrm{\Theta }+2R}=\frac{1+\eta }{sin\mathrm{\Theta }+\eta }.\end{array}$$(A.21)

Again assuming $f_{\mathrm{\Theta }}\left(\theta \right)=\frac{1}{2}sin\theta ;\theta \in [0,\pi ]$, we repeat the derivation and find

$$\begin{array}{cc}\hfill F_D\left(d\right)& =\{\begin{array}{cc}0\hfill & \mathrm{if}d\le 1;\hfill \\ \sqrt{1-{(\frac{1+\eta }{d}-\eta )}^2}\hfill & \mathrm{if}1<d<\frac{1}{\eta }+1;\hfill \\ 1\hfill & \mathrm{if}d\ge \frac{1}{\eta }+1.\hfill \end{array}\hfill \end{array}$$(A.22)

$$\begin{array}{cc}\hfill f_D\left(d\right)& =\{\begin{array}{cc}0\hfill & \mathrm{if}d\le 1;\hfill \\ \frac{(1+\eta )(1+\eta -\eta d)}{d^2\sqrt{d^2-{(1+\eta -\eta d)}^2}}\hfill & \mathrm{if}1<d<\frac{1}{\eta }+1;\hfill \\ 0\hfill & \mathrm{if}d\ge \frac{1}{\eta }+1.\hfill \end{array}\hfill \end{array}$$(A.23)

The quantile function becomes

$$\begin{array}{c}\hfill F_D^{-1}\left(p\right)=\frac{1+\eta }{\sqrt{1-p^2}+\eta }.\end{array}$$(A.24)

The minimum factor remains $F_D^{-1}p=0=1$, but a maximum factor now exists: $F_D^{-1}(p=1)=\frac{1}{\eta }+1$. The median of D is

$$\begin{array}{c}\hfill F_D^{-1}(p=\frac{1}{2})=\frac{1+\eta }{\frac{1}{2}\sqrt{3}+\eta },\end{array}$$(A.25)

which tends to $\frac{2}{\sqrt{3}}$ for η → 0 (as before), and to 1 for η → ∞: projection ceases to be an appreciable effect when the lobes are much larger than the line segment connecting their centres. If η ≠ 0, D has finite support, and thus the mean, variance and higher moments exist. The n-th non-central moment is

$$\begin{array}{c}\hfill \mathbb{E}\left[D^n\right]={(1+\eta )}^2{I}_{n-2}\left(\eta \right)-\eta (1+\eta )I_{n-1}\left(\eta \right),\end{array}$$(A.26)

where

$$\begin{array}{c}\hfill I_n\left(\eta \right):={\displaystyle {\int }_1^{\frac{1}{\eta }+1}\frac{d^n\mathrm{d}d}{\sqrt{d^2-{(1+\eta -\eta d)}^2}}.}\end{array}$$(A.27)

In particular,

$$\begin{array}{c}\hfill I_{-1}\left(\eta \right)=\frac{\pi }{2(1+\eta )},I_0\left(\eta \right)=\frac{ln(\sqrt{1-{\eta }^2}+1)-ln\eta }{\sqrt{1-{\eta }^2}},\\ \hfill I_1\left(\eta \right)=\frac{1}{(1-\eta )\eta }(1-\frac{2{\eta }^2\mathrm{arcsinh}\sqrt{\frac{1-\eta }{2\eta }}}{\sqrt{1-{\eta }^2}}).\end{array}$$(A.28)

The expectation value of D is

$$\begin{array}{cc}\hfill \mathbb{E}\left[D\right]& ={(1+\eta )}^2{I}_{-1}\left(\eta \right)-\eta (1+\eta )I_0\left(\eta \right)\hfill \\ \hfill & =(1+\eta )(\frac{\pi }{2}-\frac{\eta }{\sqrt{1-{\eta }^2}}ln(\frac{\sqrt{1-{\eta }^2}}{\eta }+\frac{1}{\eta })).\hfill \end{array}$$(A.29)

The variance of D follows from combining

$$\begin{array}{c}\hfill \mathbb{E}\left[D^2\right]={(1+\eta )}^2{I}_0\left(\eta \right)-\eta (1+\eta )I_1\left(\eta \right)\end{array}$$(A.30)

and the identity 𝕍[D] = 𝔼[D²] − 𝔼²[D]. The mode remains undefined. For typical values of η, we show the PDF and CDF of D in the bottom left and right panels of Fig. A.1, respectively. Figure A.2 shows the median, mean and standard deviation of D as a function of η.

A.4. Intrinsic proper length posterior and its moments

Our next objective is to find the posterior PDF f_{L|Lp = lp}(l) through Bayes’ theorem:

$$\begin{array}{c}\hfill f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)=\frac{f_{L_{\mathrm{p}}|L=l}\left(l_{\mathrm{p}}\right)f_L\left(l\right)}{f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)}.\end{array}$$(A.31)

Our line of attack will be to first calculate the likelihood CDF F_Lp|L = l(l_p), and then through differentiation the likelihood PDF f_Lp|L = l(l_p). Of course, we always consider l > 0, and

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}|L=l}\left(l_{\mathrm{p}}\right)& :=\mathbb{P}(L_{\mathrm{p}}\le l_{\mathrm{p}}|L=l)=\mathbb{P}(lsin\mathrm{\Theta }\le l_{\mathrm{p}})\hfill \\ \hfill & =\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le 0\hfill \\ 2\mathbb{P}(\mathrm{\Theta }\le arcsin\frac{l_{\mathrm{p}}}{l})\hfill & \text{if}0<l_{\mathrm{p}}\le l\hfill \\ 1\hfill & \text{if}{l}_{\mathrm{p}}>l,\hfill \end{array}\hfill \end{array}$$(A.32)

where we have invoked Eq. A.17 and the fact that f_Θ(θ) is symmetric in $\theta =\frac{\pi }{2}$ for the case 0 < l_p < l. Concretising this case further yields

$$\begin{array}{cc}\hfill F_{L_{\mathrm{p}}|L=l}\left(l_{\mathrm{p}}\right)& =2F_{\mathrm{\Theta }}(arcsin\frac{l_{\mathrm{p}}}{l})=1-cosarcsin\frac{l_{\mathrm{p}}}{l}\hfill \\ \hfill & =1-\sqrt{1-{\left(\frac{l_{\mathrm{p}}}{l}\right)}^2}.\hfill \end{array}$$(A.33)

Through differentiation,

$$\begin{array}{c}\hfill f_{L_{\mathrm{p}}|L=l}\left(l_{\mathrm{p}}\right)=\{\begin{array}{cc}0\hfill & \text{if}{l}_{\mathrm{p}}\le 0\hfill \\ \frac{l_{\mathrm{p}}}{l^2}\frac{1}{\sqrt{1-{\left(\frac{l_{\mathrm{p}}}{l}\right)}^2}}\hfill & \text{if}0<l_{\mathrm{p}}<l\hfill \\ 0\hfill & \text{if}{l}_{\mathrm{p}}\ge l.\hfill \end{array}\end{array}$$(A.34)

Having found the likelihood, the posterior PDF follows directly through Bayes’ theorem and the PDFs computed hitherto.

In concreto, when 0 < l_p ≤ l_min,

$$\begin{array}{c}\hfill f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)=\{\begin{array}{cc}0\hfill & \text{if}l\le l_{\min }\hfill \\ \frac{1}{I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})}\frac{1}{l_{\mathrm{p}}}\frac{1}{\sqrt{{\left(\frac{l}{l_{\mathrm{p}}}\right)}^2-1}}{\left(\frac{l}{l_{\min }}\right)}^{\xi -1}\hfill & \text{if}l>l_{\min }.\hfill \end{array}\end{array}$$(A.35)

whereas for l_p > l_min,

$$\begin{array}{c}\hfill f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)=\{\begin{array}{cc}0\hfill & \text{if}l\le l_{\mathrm{p}}\hfill \\ -\frac{\xi }{2^{1+\xi }\pi }\frac{{\mathrm{\Gamma }}^2(-\frac{\xi }{2})}{\mathrm{\Gamma }(-\xi )}\frac{1}{l_{\mathrm{p}}}\frac{1}{\sqrt{{\left(\frac{l}{l_{\mathrm{p}}}\right)}^2-1}}{\left(\frac{l}{l_{\mathrm{p}}}\right)}^{\xi -1}\hfill & \text{if}l>l_{\mathrm{p}}.\hfill \end{array}\end{array}$$(A.36)

We note that in both cases, for l ≫ l_p, $f_{L|L_{\mathrm{p}}=l_{\mathrm{p}}}\left(l\right)\propto {\left(\frac{l}{l_{\mathrm{p}}}\right)}^{\xi -2}$: the posterior probability density follows a power law in l with exponent ξ − 2.

For l_p ≤ l_min, the posterior mean is

$$\begin{array}{c}\hfill \mathbb{E}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\min }\frac{I(\xi ,\frac{l_{\mathrm{p}}}{l_{\min }})}{I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})}.\end{array}$$(A.37)

The posterior variance 𝕍[L | L_p=l_p] follows from considering the second non-central moment:

$$\begin{array}{c}\hfill \mathbb{E}[L^2|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\min }^2\frac{I(\xi +1,\frac{l_{\mathrm{p}}}{l_{\min }})}{I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})},\end{array}$$(A.38)

so that

$$\begin{array}{cc}\hfill \mathbb{V}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]& =\mathbb{E}[L^2|L_{\mathrm{p}}=l_{\mathrm{p}}]-{\mathbb{E}}^2[L|L_{\mathrm{p}}=l_{\mathrm{p}}]\hfill \\ \hfill & =l_{\min }^2(\frac{I(\xi +1,\frac{l_{\mathrm{p}}}{l_{\min }})}{I(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})}-\frac{I^2(\xi ,\frac{l_{\mathrm{p}}}{l_{\min }})}{I^2(\xi -1,\frac{l_{\mathrm{p}}}{l_{\min }})}).\hfill \end{array}$$(A.39)

For l_p > l_min,

$$\begin{array}{cc}\hfill \mathbb{E}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]& =l_{\mathrm{p}}·\frac{-\xi }{2^{\xi +1}\pi }\frac{{\mathrm{\Gamma }}^2(-\frac{\xi }{2})}{\mathrm{\Gamma }(-\xi )}{\displaystyle {\int }_1^{\infty }\frac{{\eta }^{\xi }\mathrm{d}\eta }{\sqrt{{\eta }^2-1}}}\hfill \\ \hfill & =l_{\mathrm{p}}·\frac{-\xi }{2^{2\xi +3}\pi }\frac{{\mathrm{\Gamma }}^4(-\frac{\xi }{2})}{{\mathrm{\Gamma }}^2(-\xi )}.\hfill \end{array}$$(A.40)

Proceeding analogously,

$$\begin{array}{c}\hfill \mathbb{E}[L^2|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\mathrm{p}}^2\frac{\xi }{\xi +1},\end{array}$$(A.41)

so that

$$\begin{array}{c}\hfill \mathbb{V}[L|L_{\mathrm{p}}=l_{\mathrm{p}}]=l_{\mathrm{p}}^2(\frac{\xi }{\xi +1}-\frac{{\xi }^2}{2^{4\xi +6}{\pi }^2}\frac{{\mathrm{\Gamma }}^8(-\frac{\xi }{2})}{{\mathrm{\Gamma }}^4(-\xi )}).\end{array}$$(A.42)

In this case, both the mean and standard deviation of L | L_p = l_p are proportional to l_p. In the table below, we list the mean and standard deviation in multiples of l_p for several values of ξ. Since we assume ξ < −1, the mean and variance are guaranteed to exist. The existence of higher-order moments is ξ-dependent; in concreto, the highest defined order is ⌈ − ξ⌉.

We prove that L_p and D are not independent by contradiction. If we assume that L_p and D are independent, then 𝔼[L_pD]=𝔼[L_p]𝔼[D], or 𝔼[L]=𝔼[L]𝔼[sin Θ]𝔼[D] by the independence of L and sin Θ. In other words, if L_p and D are independent, then 1 = 𝔼[sin Θ]𝔼[D]. However, $\mathbb{E}[sin\mathrm{\Theta }]=\frac{\pi }{4}$ and $\mathbb{E}[D]=\frac{\pi }{2}$; because $1\ne \frac{{\pi }^2}{8}$, the initial assertion must be wrong.

A.5. GRG inclination angle

We derive the inclination angle distribution for giants. The probability that a GRG has inclination angle θ is

$$\begin{array}{cc}\hfill f_{\mathrm{\Theta }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\theta )\mathrm{d}\theta & =\frac{\mathbb{P}(\mathrm{\Theta }=\theta ,Lsin\mathrm{\Theta }>l_{\mathrm{p},\mathrm{GRG}})}{\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})}\hfill \\ \hfill & =\frac{f_{\mathrm{\Theta }}(\theta )\mathrm{d}\theta ·\mathbb{P}(L>\frac{l_{\mathrm{p},\mathrm{GRG}}}{sin\theta })}{\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})},\hfill \end{array}$$(A.43)

where we make use of the fact that the numerator’s joint probability factorises because L and Θ are independent RVs. We thus find

$$\begin{array}{c}\hfill f_{\mathrm{\Theta }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\theta )=\frac{(1-F_L\left(\frac{l_{\mathrm{p},\mathrm{GRG}}}{sin\theta }\right))f_{\mathrm{\Theta }}\left(\theta \right)}{1-F_{L_{\mathrm{p}}}(l_{\mathrm{p},\mathrm{GRG}})}.\end{array}$$(A.44)

Under the Paretian assumption for L, we have

$$\begin{array}{c}\hfill f_{\mathrm{\Theta }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\theta )=\frac{2}{-(\xi +1)\sqrt{\pi }}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{sin}^{-\xi }\theta .\end{array}$$(A.45)

A.6. GRG angular length

The GRG angular length RV Φ | L_p > l_p, GRG relates to the GRG projected proper length RV L_p | L_p > l_p, GRG and the comoving distance RV R as

$$\begin{array}{c}\hfill \mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}=L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}·\frac{1+z\left(R\right)}{R}.\end{array}$$(A.46)

The model predicts the distribution of GRG angular lengths in the Local Universe up to comoving distance r_max. If the GRG number density is constant in the Local Universe,

$$\begin{array}{cc}\hfill F_R\left(r\right)& =\{\begin{array}{cc}0\hfill & \text{if}r\le 0\hfill \\ {\left(\frac{r}{r_{\max }}\right)}^3\hfill & \text{if}0<r<r_{\max }\hfill \\ 1\hfill & \text{if}r\ge r_{\max }\hfill \end{array}\hfill \end{array}$$(A.47)

$$\begin{array}{cc}\hfill f_R\left(r\right)& =\{\begin{array}{cc}0\hfill & \text{if}r\le 0\hfill \\ \frac{3}{r_{\max }}{\left(\frac{r}{r_{\max }}\right)}^2\hfill & \text{if}0<r<r_{\max }\hfill \\ 0\hfill & \text{if}r\ge r_{\max }.\hfill \end{array}\hfill \end{array}$$(A.48)

Because GRG life cycles are shorter than the age of the Universe, L_p | L_p > l_p, GRG and R are independent RVs. In a non-expanding universe, z(R) = 0, and the distribution of Φ | L_p > l_p, GRG can be calculated analytically. From a well-known ratio distribution identity,

$$\begin{array}{cc}\hfill f_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )& ={\displaystyle {\int }_{-\infty }^{\infty }|r|·f_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi r)·f_R(r)\mathrm{d}r}\hfill \\ \hfill & ={\displaystyle {\int }_0^{r_{\max }}{f}_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi r)·3{\left(\frac{r}{r_{\max }}\right)}^3\mathrm{d}r.}\hfill \end{array}$$(A.49)

The integrand is non-zero only when ϕr ≥ l_p, GRG, suggesting a lower integration limit of $\frac{l_{\mathrm{p},\mathrm{GRG}}}{\varphi }$. The integral vanishes altogether when $\frac{l_{\mathrm{p},\mathrm{GRG}}}{\varphi }\ge r_{\max }$. Calling ${\varphi }_{\mathrm{GRG}}:=\frac{l_{\mathrm{p},\mathrm{GRG}}}{r_{\max }}$, we have f_{Φ | Lp > lp, GRG}(ϕ) = 0 for ϕ ≤ ϕ_GRG. For ϕ > ϕ_GRG,

$$\begin{array}{cc}\hfill f_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )& ={\displaystyle {\int }_{\frac{l_{\mathrm{p},\mathrm{GRG}}}{\varphi }}^{r_{\max }}-\frac{\xi +1}{l_{\mathrm{p},\mathrm{GRG}}}{\left(\frac{\varphi r}{l_{\mathrm{p},\mathrm{GRG}}}\right)}^{\xi }·3{\left(\frac{r}{r_{\max }}\right)}^3\mathrm{d}r}\hfill \\ \hfill & =-3\frac{\xi +1}{\xi +4}·\frac{1}{{\varphi }_{\mathrm{GRG}}}({\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{\xi }-{\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{-4}),\hfill \end{array}$$(A.50)

where we assume ξ ≠ −4. This PDF depends on ξ, l_p, GRG, and r_max only. The CDF follows from direct integration:

$$\begin{array}{cc}& F_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )=\hfill \\ \hfill & \{\begin{array}{cc}0\hfill & \text{if}\varphi \le {\varphi }_{\mathrm{GRG}}\hfill \\ 1+\frac{1}{\xi +4}(-3{\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{\xi +1}-(\xi +1){\left(\frac{\varphi }{{\varphi }_{\mathrm{GRG}}}\right)}^{-3})\hfill & \text{if}\varphi >{\varphi }_{\mathrm{GRG}},\hfill \end{array}\hfill \end{array}$$(A.51)

just as the mean:

$$\begin{array}{c}\hfill \mathbb{E}[\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]=3\frac{\xi +1}{\xi +4}(\frac{1}{\xi +2}+\frac{1}{2})·{\varphi }_{\mathrm{GRG}},\end{array}$$(A.52)

which exists only for ξ < −2. By solving

$$\begin{array}{c}\hfill \frac{\mathrm{d}{f}_{\mathrm{\Phi }}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}{\mathrm{d}\varphi }\left({\varphi }_{\mathrm{mode}}\right)=0\end{array}$$(A.53)

for the GRG angular length mode ϕ_mode, we find

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{mode}}={\left(\frac{-4}{\xi }\right)}^{\frac{1}{\xi +4}}·{\varphi }_{\mathrm{GRG}}.\end{array}$$(A.54)

There is no explicit expression for the associated median.

For a general non-Euclidean universe, no simple analytic form for the GRG angular length PDF appears to exist. We find an approximation valid at low redshifts by considering the Maclaurin polynomial of degree 1 for z(r), the relation between comoving distance and cosmological redshift. One finds $z(r)\approx \frac{H_0}{c}r$. As a result, for the Local Universe, Eq. 14 becomes

$$\begin{array}{cc}\hfill \mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}& \approx L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}·\frac{1+\frac{H_0}{c}R}{R}\hfill \\ \hfill & =L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}·(\frac{1}{R}+\frac{1}{d_H}),\hfill \end{array}$$(A.55)

where we use that the Hubble distance $d_H:=\frac{c}{H_0}$. From a well-known inverse distribution identity,

$$\begin{array}{c}\hfill f_{\frac{1}{R}}(k)=\{\begin{array}{cc}0\hfill & \text{if}k\le \frac{1}{r_{\max }}\hfill \\ \frac{3}{r_{\max }^3{k}^4}\hfill & \text{if}k>\frac{1}{r_{\max }},\hfill \end{array}\end{array}$$(A.56)

so that

$$\begin{array}{c}\hfill f_{\frac{1}{R}+\frac{1}{d_H}}(k)=\{\begin{array}{cc}0\hfill & \text{if}k\le \frac{1}{r_{\max }}+\frac{1}{d_H}\hfill \\ \frac{3}{r_{\max }^3{(k-\frac{1}{d_H})}^4}\hfill & \text{if}k>\frac{1}{r_{\max }}+\frac{1}{d_H}.\hfill \end{array}\end{array}$$(A.57)

Combining a well-known product distribution identity with the fact that L_p | L_p > l_p, GRG and $\frac{1}{R}+\frac{1}{d_H}$ are independent RVs, we find

$$\begin{array}{c}\hfill f_{\mathrm{\Phi }|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}(\varphi )=\frac{-3(\xi +1){\varphi }^{\xi }}{r_{\max }^3{l}_{\mathrm{p},\mathrm{GRG}}^{\xi +1}}{\displaystyle {\int }_{\frac{1}{r_{\max }}+\frac{1}{d_H}}^{\frac{\varphi }{l_{\mathrm{p},\mathrm{GRG}}}}\frac{\mathrm{d}k}{k^{\xi +1}{(k-\frac{1}{d_H})}^4}}\end{array}$$(A.58)

when $\varphi >{\varphi }_{\mathrm{GRG}}+\frac{l_{\mathrm{p},\mathrm{GRG}}}{d_H}$. When $\varphi \le {\varphi }_{\mathrm{GRG}}+\frac{l_{\mathrm{p},\mathrm{GRG}}}{d_H}$, f_{Φ | Lp > lp, GRG}(ϕ) = 0. It is easy to verify that in the Euclidean limit, $\frac{1}{d_H}\to 0$, this expression reduces to that of Eq. A.50.

A.7. Maximum likelihood estimation of the tail index

How can we estimate ξ from a sample of N giants? Let ξ_MLE be the maximum likelihood estimate (MLE) of ξ. This RV is a function of N IID RVs {L_p, 1, ..., L_p, N}∼L_p | L_p > l_p, GRG. Define the following likelihood and log-likelihood functions:

$$\begin{array}{cc}\hfill \mathcal{L}\left(\xi \right)& :={\displaystyle \prod _{i=1}^N{f}_{L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}\left(L_{\mathrm{p},\mathrm{i}}\right)}\hfill \end{array}$$(A.59)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mathcal{L}}\left(\xi \right)& :=ln\left(\mathcal{L}\left(\xi \right)l_{\mathrm{p},\mathrm{GRG}}^N\right)\hfill \\ \hfill & =Nln(-(\xi +1))+\xi {\displaystyle \sum _{i=1}^N}ln\frac{L_{\mathrm{p},\mathrm{i}}}{l_{\mathrm{p},\mathrm{GRG}}}.\hfill \end{array}$$(A.60)

We note the necessity to include a factor $l_{\text{p},\text{GRG}}^N$ in the definition of the log-likelihood to avoid a dimensionality error. The second derivative of $\stackrel{\sim }{\mathcal{L}}$ to ξ is

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\stackrel{\sim }{\mathcal{L}}}{\mathrm{d}{\xi }^2}=-\frac{N}{{(\xi +1)}^2}<0.\end{array}$$(A.61)

Thus, if there exists a solution to the equation $\frac{\mathrm{d}\stackrel{\sim }{\mathcal{L}}}{\mathrm{d}\xi }=0$, it must correspond to a global maximum of the likelihood and log-likelihood functions:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\stackrel{\sim }{\mathcal{L}}}{\mathrm{d}\xi }(\xi ={\xi }_{\mathrm{MLE}})=0,\text{or}{\xi }_{\mathrm{MLE}}& =-\frac{N}{{\sum }_{i=1}^{N}ln\frac{L_{\mathrm{p},\mathrm{i}}}{l_{\mathrm{p},\mathrm{GRG}}}}-1.\hfill \end{array}$$(A.62)

A.8. Observed projected proper length

A.8.1. General considerations

The model can be extended to incorporate observational selection effects. The relevant effects to consider vary from RG search campaign to RG search campaign, although some formulae apply in all cases. We derive these here.

To keep our extensions simple, we assume that the survey sensitivity is sufficient to detect RGs up to some redshift z_max only. In addition we assume that the projected proper length distribution does not evolve between z = z_max and z = 0. We let p_obs(l_p, z) denote the probability that an RG of projected proper length l_p at cosmological redshift z is detected during the campaign. Also, r is the radial comoving distance and n is the total RG number density, counting the intrinsic number of RGs (irrespective of length) per unit of comoving volume. The observed number of RGs with projected proper length between l_p and l_p + dl_p throughout a survey covering a solid angle Ω is dN_Lp, obs(l_p,Ω), where

$$\begin{array}{cc}& \mathrm{d}{N}_{L_{\mathrm{p}},\mathrm{obs}}(l_{\mathrm{p}},\mathrm{\Omega })\hfill \\ \hfill & =\frac{\mathrm{\Omega }}{4\pi }{\displaystyle {\int }_0^{z_{\max }}}n·f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)\mathrm{d}{l}_{\mathrm{p}}·p_{\mathrm{obs}}(l_{\mathrm{p}},z)4\pi r^2\left(z\right)\frac{\mathrm{d}r}{\mathrm{d}z}\mathrm{d}z.\hfill \end{array}$$(A.63)

The total number of RGs with projected proper length between l_p and l_p + dl_p throughout a survey with solid angle Ω is dN_Lp(l_p,Ω), where

$$\begin{array}{c}\hfill \mathrm{d}{N}_{L_{\mathrm{p}}}(l_{\mathrm{p}},\mathrm{\Omega })=\frac{\mathrm{\Omega }}{4\pi }{\displaystyle {\int }_0^{z_{\max }}n·f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)\mathrm{d}{l}_{\mathrm{p}}·4\pi r^2\left(z\right)\frac{\mathrm{d}r}{\mathrm{d}z}\mathrm{d}z.}\end{array}$$(A.64)

We define the completeness C(l_p, z_max) to be

$$\begin{array}{cc}\hfill C(l_{\mathrm{p}},z_{\max })& :=\frac{\mathrm{d}{N}_{L_{\mathrm{p}},\mathrm{obs}}(l_{\mathrm{p}},\mathrm{\Omega })}{\mathrm{d}{N}_{L_{\mathrm{p}}}(l_{\mathrm{p}},\mathrm{\Omega })}\hfill \\ \hfill & =\frac{{\int }_0^{z_{\max }}{p}_{\mathrm{obs}}(l_{\mathrm{p}},z)r^2\left(z\right)E^{-1}\left(z\right)\mathrm{d}z}{{\int }_0^{z_{\max }}{r}^2\left(z\right)E^{-1}\left(z\right)\mathrm{d}z}.\hfill \end{array}$$(A.65)

The completeness only depends on the function p_obs(l_p, z) and our choice of z_max.

Let the RV 𝒪 denote whether an RG picked at random within z < z_max is detected during the search campaign. We have 𝒪 ∼ Bernoulli(C(L_p,z_max)); the parameter that determines the distribution of 𝒪 is itself an RV. It immediately follows that ℙ(𝒪 = 1 | L_p=l_p) = C(l_p, z_max).

If L_p would be discrete,

$$\begin{array}{cc}\hfill \mathbb{P}(L_{\mathrm{p}}=l_{\mathrm{p}}|\mathcal{O}=1)& =\frac{\mathbb{P}(\mathcal{O}=1,L_{\mathrm{p}}=l_{\mathrm{p}})}{\mathbb{P}(\mathcal{O}=1)}\hfill \\ \hfill & =\frac{\mathbb{P}(\mathcal{O}=1|L_{\mathrm{p}}=l_{\mathrm{p}})\mathbb{P}(L_{\mathrm{p}}=l_{\mathrm{p}})}{{\sum }_{l_{\mathrm{p}}}\mathbb{P}(\mathcal{O}=1|L_{\mathrm{p}}=l_{\mathrm{p}})\mathbb{P}(L_{\mathrm{p}}=l_{\mathrm{p}})}.\hfill \end{array}$$(A.66)

Let L_p, obs be the observed projected proper length RV; that is L_p, obs := L_p | 𝒪 = 1. Its PDF is given by the continuous analogon of the preceding equation:

$$\begin{array}{c}\hfill f_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p}}\right)=\frac{C(l_{\mathrm{p}},z_{\max })f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)}{{\int }_0^{\infty }C(l_{\mathrm{p}},z_{\max })f_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)\mathrm{d}{l}_{\mathrm{p}}}.\end{array}$$(A.67)

Obviously, $F_{L_{\text{p},\text{obs}}}{l}_{\text{p}}={\int }_0^{l_{\text{p}}}{f}_{L_{\text{p},\text{obs}}}{l}_{\text{p}}^{\prime }\text{d}{l}_{\text{p}}^{\prime }$. The CDF F_{Lp, obs|Lp, obs > lp, GRG}(l_p) follows from an analogon of Eq. A.12.

We note that multiplying p_obs(l_p, z) with an l_p- and z-independent factor changes C(l_p, z_max) by the same factor, but leaves f_Lp, obs unaltered.

A.8.2. Fuzzy angular length threshold

Here we illustrate a simple extension. When performing visual searches for GRG candidates, a natural criterion is to only inspect sources with an angular length that is larger than some threshold. Researchers determine the threshold based on the amount of time available to them: lower thresholds will lead to more complete samples, but will take more time to collect. For humans, it is hard to estimate a source’s angular length precisely by eye; as a result, some of the GRG candidates included in the project’s GRG candidate catalogue will be sources with an angular length below the threshold, whilst others will be sources with an angular length above the threshold. We idealise this situation by asserting that sources with angular length ϕ_min or below are included in the catalogue with probability 0 (i.e. never), and that sources with an angular length ϕ_max or above are included in the catalogue with probability 1 (i.e. always). We assume a linear increase in probability as a function of ϕ for intermediate angular lengths: a source with angular length ϕ_min < ϕ < ϕ_max is included in the catalogue with probability

$$\begin{array}{c}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)=min\{max\{\frac{\varphi (l_{\mathrm{p}},z)-{\varphi }_{\min }}{{\varphi }_{\max }-{\varphi }_{\min }},0\},1\}.\end{array}$$(A.68)

In a flat Friedmann–Lemaître–Robertson–Walker (FLRW) universe, the angular length of an RG with projected proper length l_p at cosmological redshift z is

$$\begin{array}{c}\hfill \varphi (l_{\mathrm{p}},z)=\frac{l_{\mathrm{p}}(1+z)}{r\left(z\right)}.\end{array}$$(A.69)

A.8.3. Surface brightness limitations

Fanaroff–Riley class II If the surface brightness B_ν of the lobes is proportional to some negative power ζ of the GRG’s proper length L, that is

$$\begin{array}{c}\hfill B_{\nu }=b_{\nu ,\mathrm{ref}}{\left(\frac{L}{l_{\mathrm{ref}}}\right)}^{\zeta },\end{array}$$(A.70)

then B_ν is Pareto distributed — just like L:

$$\begin{array}{cc}\hfill F_{B_{\nu }}\left(b\right)& =\{\begin{array}{cc}0\hfill & \text{if}b\le 0\hfill \\ {\left(\frac{b}{b_{\nu ,\mathrm{ref}}}\right)}^{\frac{\xi +1}{\zeta }}\hfill & \text{if}0<b<b_{\nu ,\mathrm{ref}}\hfill \\ 1\hfill & \text{if}b\ge b_{\nu ,\mathrm{ref}},\hfill \end{array}\hfill \end{array}$$(A.71)

$$\begin{array}{cc}\hfill f_{B_{\nu }}\left(b\right)& =\{\begin{array}{cc}0\hfill & \text{if}b\le 0\hfill \\ \frac{\xi +1}{\zeta }\frac{1}{b_{\nu ,\mathrm{ref}}}{\left(\frac{b}{b_{\nu ,\mathrm{ref}}}\right)}^{\frac{\xi +1}{\zeta }-1}\hfill & \text{if}0<b<b_{\nu ,\mathrm{ref}}\hfill \\ 0\hfill & \text{if}b\ge b_{\nu ,\mathrm{ref}}.\hfill \end{array}\hfill \end{array}$$(A.72)

If RGs are self-similar in shape, then their lobe volumes are proportional to L³. RGs appear to retain constant lobe luminosity density over most of their lifetime, so that their lobe monochromatic emission coefficients (Rybicki & Lightman 1986) are proportional to L⁻³. As line-of-sight lengths through the lobes are proportional to L, surface brightness is proportional to L⁻². These arguments thus suggest ζ = −2.

It is a poor approximation to assume that all giants of proper length l_ref have the same surface brightness b_ν, ref: observations suggest a variability of several orders of magnitude. A better description is

$$\begin{array}{c}\hfill B_{\nu }=b_{\nu ,\mathrm{ref}}{\left(\frac{L}{l_{\mathrm{ref}}}\right)}^{\zeta }S,\end{array}$$(A.73)

where S is a lognormal RV whose median is 1. The PDF of S is then determined by parameter σ_ref:

$$\begin{array}{c}\hfill f_S\left(s\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{ref}}s}exp(-\frac{{ln}^{2}s}{2{\sigma }_{\mathrm{ref}}^2}).\end{array}$$(A.74)

The surface brightness CDF and PDF now become

$$\begin{array}{cc}\hfill F_{B_{\nu }}\left(b\right)& =\{\begin{array}{cc}0\hfill & \text{if}b\le 0\hfill \\ F_S\left(\frac{b}{b_{\nu ,\mathrm{ref}}}\right)+{\left(\frac{b}{b_{\nu ,\mathrm{ref}}}\right)}^{\frac{\xi +1}{\zeta }}{\int }_{\frac{b}{b_{\nu ,\mathrm{ref}}}}^{\infty }{s}^{-\frac{\xi +1}{\zeta }}{f}_S\left(s\right)\mathrm{d}s\hfill & \text{if}b>0,\hfill \end{array}\hfill \\ \hfill \\ \hfill f_{B_{\nu }}\left(b\right)& =\{\begin{array}{cc}0\hfill & \text{if}b\le 0\hfill \\ \frac{\xi +1}{\zeta }\frac{1}{b_{\nu ,\mathrm{ref}}}{\left(\frac{b}{b_{\nu ,\mathrm{ref}}}\right)}^{\frac{\xi +1}{\zeta }-1}{\int }_{\frac{b}{b_{\nu ,\mathrm{ref}}}}^{\infty }{s}^{-\frac{\xi +1}{\zeta }}{f}_S\left(s\right)\mathrm{d}s\hfill & \text{if}b>0.\hfill \end{array}\hfill \end{array}$$(A.75)

We note that B_ν is not exactly Pareto distributed anymore.

In a relativistic rather than Euclidean universe, surface brightness is not constant with distance. To describe RGs beyond z = 0, we introduce a final model refinement:

$$\begin{array}{c}\hfill B_{\nu }=\frac{b_{\nu ,\mathrm{ref}}·S}{{(1+Z)}^{3-\alpha }}{\left(\frac{L}{l_{\mathrm{ref}}}\right)}^{\zeta },\end{array}$$(A.76)

where the RV Z denotes cosmological redshift and α is the spectral index of the lobes. We interpret b_ν, ref ⋅ S as the (lognormally distributed) lobe surface brightness for RGs of intrinsic proper length l_ref at z = 0. In this case, the CDF and PDF of B_ν are most easily determined through sampling. To sample Z, we can first sample the comoving distance RV R instead, and subsequently use Z = z_𝔐(R). We stress that this conversion depends on cosmological parameters 𝔐.

To forward model a survey’s surface brightness selection effect, we must compute

$$\begin{array}{cc}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)& =\mathbb{P}(B_{\nu }>b_{\nu ,\mathrm{th}}|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z)\hfill \\ \hfill & =1-F_{B_{\nu }|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z}\left(b_{\nu ,\mathrm{th}}\right),\hfill \end{array}$$(A.77)

where b_ν, th > 0 is the surface brightness threshold. Typically, b_ν, th is comparable to the survey’s RMS noise. What is F_{Bν | Lp = lp,  Z = z}(b)? In the simplest case, devoid of S- and Z-dependence,

$$\begin{array}{cc}\hfill F_{B_{\nu }|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z}\left(b\right)& =\mathbb{P}(b_{\nu ,\mathrm{ref}}{\left(\frac{L_{\mathrm{p}}}{sin\mathrm{\Theta }·l_{\mathrm{ref}}}\right)}^{\zeta }\le b|L_{\mathrm{p}}=l_{\mathrm{p}})\hfill \\ \hfill & =\mathbb{P}({sin}^{-\zeta }\mathrm{\Theta }\le \frac{b}{b_{\nu ,\mathrm{ref}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta })\hfill \\ \hfill & =\mathbb{P}(sin\mathrm{\Theta }\le {\stackrel{\sim }{s}}^{-\frac{1}{\zeta }})\hfill \\ \hfill & =\{\begin{array}{cc}0\hfill & \text{if}b\le 0\hfill \\ 2F_{\mathrm{\Theta }}(arcsin\left({\stackrel{\sim }{s}}^{-\frac{1}{\zeta }}\right))\hfill & \text{if}0<b<b_{\nu ,\mathrm{ref}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{\zeta }\hfill \\ 1\hfill & \text{if}b\ge b_{\nu ,\mathrm{ref}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{\zeta }\hfill \end{array}\hfill \\ \hfill & =1-\sqrt{1-min\{{\stackrel{\sim }{s}}^{-\frac{2}{\zeta }}(max\{b,0\}),1\}},\hfill \end{array}$$(A.78)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{s}=\stackrel{\sim }{s}\left(b\right):=\frac{b}{b_{\nu ,\mathrm{ref}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta }.\end{array}$$(A.79)

Therefore,

$$\begin{array}{c}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)=\sqrt{1-min\{{\stackrel{\sim }{s}}^{-\frac{2}{\zeta }}\left(b_{\nu ,\mathrm{th}}\right),1\}}.\end{array}$$(A.80)

In the most refined case, for b > 0,

$$\begin{array}{cc}& F_{B_{\nu }|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z}\left(b\right)\hfill \\ \hfill & =\mathbb{P}(\frac{b_{\nu ,\mathrm{ref}}·S}{{(1+Z)}^{3-\alpha }}{\left(\frac{L_{\mathrm{p}}}{sin\mathrm{\Theta }·l_{\mathrm{ref}}}\right)}^{\zeta }\le b|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z)\hfill \\ \hfill & =\mathbb{P}({sin}^{-\zeta }\mathrm{\Theta }·S\le \frac{b}{b_{\nu ,\mathrm{ref}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta }{(1+z)}^{3-\alpha })\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }}\mathbb{P}({sin}^{-\zeta }\mathrm{\Theta }\le \frac{b}{b_{\nu ,\mathrm{ref}}}\frac{1}{s}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta }{(1+z)}^{3-\alpha })f_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }}\mathbb{P}(sin\mathrm{\Theta }\le {\left(\frac{\stackrel{\sim }{s}}{s}\right)}^{-\frac{1}{\zeta }})f_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & ={\displaystyle {\int }_0^{\stackrel{\sim }{s}}}{f}_S\left(s\right)\mathrm{d}s+{\displaystyle {\int }_{\stackrel{\sim }{s}}^{\infty }}2F_{\mathrm{\Theta }}(arcsin\left({\left(\frac{\stackrel{\sim }{s}}{s}\right)}^{-\frac{1}{\zeta }}\right))f_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & ={\displaystyle {\int }_0^{\stackrel{\sim }{s}}}{f}_S\left(s\right)\mathrm{d}s+{\displaystyle {\int }_{\stackrel{\sim }{s}}^{\infty }}{f}_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & -{\displaystyle {\int }_{\stackrel{\sim }{s}}^{\infty }}cosarcsin\left({\left(\frac{\stackrel{\sim }{s}}{s}\right)}^{-\frac{1}{\zeta }}\right)f_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & =1-{\displaystyle {\int }_{\stackrel{\sim }{s}}^{\infty }}\sqrt{1-{\left(\frac{\stackrel{\sim }{s}}{s}\right)}^{-\frac{2}{\zeta }}}{f}_S\left(s\right)\mathrm{d}s,\hfill \end{array}$$(A.81)

where

$$\begin{array}{c}\hfill \stackrel{\sim }{s}=\stackrel{\sim }{s}\left(b\right):=\frac{b}{b_{\nu ,\mathrm{ref}}}{\left(\frac{l_{\mathrm{p}}}{l_{\mathrm{ref}}}\right)}^{-\zeta }{(1+z)}^{3-\alpha }.\end{array}$$(A.82)

The minimum value of S for which an RG of projected length l_p at redshift z is detectable, is $s_{\min }(l_{\mathrm{p}},z):=\stackrel{\sim }{s}(b_{\nu ,\mathrm{th}})$. For brevity, we simply write s_min instead. We find

$$\begin{array}{cc}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)=& {\displaystyle {\int }_{s_{\min }}^{\infty }}\sqrt{1-{\left(\frac{s_{\min }}{s}\right)}^{-\frac{2}{\zeta }}}{f}_S\left(s\right)\mathrm{d}s.\hfill \end{array}$$(A.83)

The following approximation might facilitate the numerical evaluation of this integral. We note that for s ≫ s_min, the square root factor in the integral approaches 1. Now split up the original integral in two parts, where η governs the approximation’s accuracy:

$$\begin{array}{cc}\hfill p_{\mathrm{obs}}& \approx {\widehat{p}}_{\mathrm{obs}}\hfill \\ \hfill & ={\displaystyle {\int }_{s_{\min }}^{\eta s_{\min }}}\sqrt{1-{\left(\frac{s_{\min }}{s}\right)}^{-\frac{2}{\zeta }}}{f}_S\left(s\right)\mathrm{d}s+{\displaystyle {\int }_{\eta s_{\min }}^{\infty }}{f}_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & ={\int }_{s_{\min }}^{\eta s_{\min }}\sqrt{1-{\left(\frac{s_{\min }}{s}\right)}^{-\frac{2}{\zeta }}}{f}_S\left(s\right)\mathrm{d}s+1-F_S\left(\eta s_{\min }\right),\hfill \end{array}$$(A.84)

where we substitute numerically integrating to infinity for an evaluation of the CDF of the lognormally distributed RV S. The approximation error is bounded from above:

$$\begin{array}{cc}\hfill {\widehat{p}}_{\mathrm{obs}}-p_{\mathrm{obs}}& ={\displaystyle {\int }_{\eta s_{\min }}^{\infty }}(1-\sqrt{1-{\left(\frac{s_{\min }}{s}\right)}^{-\frac{2}{\zeta }}})f_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & <(1-\sqrt{1-{\left(\frac{s_{\min }}{\eta s_{\min }}\right)}^{-\frac{2}{\zeta }}}){\displaystyle {\int }_{\eta s_{\min }}^{\infty }}{f}_S\left(s\right)\mathrm{d}s\hfill \\ \hfill & =(1-\sqrt{1-{\eta }^{\frac{2}{\zeta }}})(1-F_S\left(\eta s_{\min }\right))\hfill \\ \hfill & <1-\sqrt{1-{\eta }^{\frac{2}{\zeta }}}.\hfill \end{array}$$(A.85)

Let us assume ζ = −2. For η = 100, the approximation error is at most 0.005, and for η = 1000, the approximation error is at most 0.0005.

Fanaroff–Riley class I The simplest correction in which FRI RGs retain a well-defined notion of length assumes a linearly decreasing surface brightness, from some value b_ν(0) at the core to zero at the RG’s two endpoints. For a symmetric FRI RG, the surface brightness profile along one of the jets is b_ν : ℝ_≥0 → ℝ_≥0, and depends on the projected proper distance from the core r_p as

$$\begin{array}{c}\hfill b_{\nu }(r_{\mathrm{p}})=\{\begin{array}{cc}{b}_{\nu }(0)(1-\frac{2r_{\mathrm{p}}}{l_{\mathrm{p}}})\hfill & \text{if}0\le r_{\mathrm{p}}<\frac{l_{\mathrm{p}}}{2}\hfill \\ 0\hfill & \text{if}{r}_{\mathrm{p}}\ge \frac{l_{\mathrm{p}}}{2}.\hfill \end{array}\end{array}$$(A.86)

Now we consider an FRI GRG, whose projected proper length l_p > l_p, GRG would only be observed in full in the absence of noise. In actual observations, this GRG is detected as a GRG if and only if

$$\begin{array}{c}\hfill b_{\nu }\left(\frac{l_{\mathrm{p},\mathrm{GRG}}}{2}\right)>b_{\nu ,\mathrm{th}},\text{or}{b}_{\nu }(0)(1-\frac{l_{\mathrm{p},\mathrm{GRG}}}{l_{\mathrm{p}}})>b_{\nu ,\mathrm{th}}.\end{array}$$(A.87)

Under our assumption of a linear surface brightness profile, the mean surface brightness along the jet axis $⟨b_{\nu }⟩=b_{\nu }(\frac{1}{2}·\frac{l_{\mathrm{p}}}{2})$. As this must be half of the surface brightness at the core,

$$\begin{array}{c}\hfill ⟨b_{\nu }⟩=b_{\nu }\left(\frac{l_{\mathrm{p}}}{4}\right)=\frac{b_{\nu }(0)}{2}.\end{array}$$(A.88)

Combining Eqs. A.87 and A.88, we find that the GRG is detected as such if

$$\begin{array}{c}\hfill ⟨b_{\nu }⟩>\frac{b_{\nu ,\mathrm{th}}}{2(1-\frac{l_{\mathrm{p},\mathrm{GRG}}}{l_{\mathrm{p}}})}.\end{array}$$(A.89)

Now regarding ⟨b_ν⟩ as an RV and recognising that it might behave exactly as in Eq. A.76, we find that the surface brightness selection effect for FRI giants may be modelled as

$$\begin{array}{c}\hfill p_{\mathrm{obs}}(l_{\mathrm{p}},z)=\mathbb{P}(B_{\nu }>\frac{b_{\nu ,\mathrm{th}}}{2(1-\frac{l_{\mathrm{p},\mathrm{GRG}}}{l_{\mathrm{p}}})}|L_{\mathrm{p}}=l_{\mathrm{p}},Z=z).\end{array}$$(A.90)

We see that the full formulaic structure of p_obs(l_p, z) is the same for FRI and FRII giants, except that for FRI giants a change

$$\begin{array}{c}\hfill b_{\nu ,\mathrm{th}}\to \frac{b_{\nu ,\mathrm{th}}}{2(1-\frac{l_{\mathrm{p},\mathrm{GRG}}}{l_{\mathrm{p}}})}\end{array}$$(A.91)

is necessary. There is no change for l_p = 2l_p, GRG.

A.9. GRG number density

A statistic of major interest is the number density of giants in the contemporary Universe. Let n_GRG be the comoving GRG number density, so that

$$\begin{array}{c}\hfill n_{\mathrm{GRG}}=n\mathbb{P}(L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}})=n(1-F_{L_{\mathrm{p}}}\left(l_{\mathrm{p},\mathrm{GRG}}\right)).\end{array}$$(A.92)

$$\begin{array}{c}\hfill N_{\mathrm{GRG},\mathrm{obs}}(\mathrm{\Omega },z_{\max }):={\displaystyle {\int }_{l_{\mathrm{p},\mathrm{GRG}}}^{\infty }}\mathrm{d}{N}_{L_{\mathrm{p},\mathrm{obs}}}(l_{\mathrm{p}},\mathrm{\Omega }).\end{array}$$(A.93)

After invoking Eq. A.63 and isolating n, we obtain

$$\begin{array}{c}\hfill n=\frac{\frac{4\pi }{\mathrm{\Omega }}{N}_{\mathrm{GRG},\mathrm{obs}}(\mathrm{\Omega },z_{\max })}{{\int }_{l_{\mathrm{p},\mathrm{GRG}}}^{\infty }{\int }_0^{z_{\max }}{f}_{L_{\mathrm{p}}}\left(l_{\mathrm{p}}\right)p_{\mathrm{obs}}(l_{\mathrm{p}},z)4\pi r^2\left(z\right)\frac{\mathrm{d}r}{\mathrm{d}z}\mathrm{d}z\mathrm{d}{l}_{\mathrm{p}}}.\end{array}$$(A.94)

Combining Eqs. A.8, A.9, A.92 and A.94 for l_p, GRG > l_min, we arrive at

$$\begin{array}{cc}& n_{\mathrm{GRG}}\hfill \\ \hfill & =-\frac{l_{\mathrm{p},\mathrm{GRG}}^{\xi +1}}{\xi +1}·\frac{\frac{4\pi }{\mathrm{\Omega }}{N}_{\mathrm{GRG},\mathrm{obs}}(\mathrm{\Omega },z_{\max })}{{\int }_{l_{\mathrm{p},\mathrm{GRG}}}^{\infty }{l}_{\mathrm{p}}^{\xi }{\int }_0^{z_{\max }}{p}_{\mathrm{obs}}(l_{\mathrm{p}},z)4\pi r^2\left(z\right)\frac{\mathrm{d}r}{\mathrm{d}z}\mathrm{d}z\mathrm{d}{l}_{\mathrm{p}}}.\hfill \end{array}$$(A.95)

From observations, we know Ω and can — for a given z_max — simply count N_GRG, obs(Ω,z_max). Moreover, we can fit ξ and the parameters that occur in p_obs(l_p, z) (e.g. $\frac{1}{2}({\varphi }_{\max }-{\varphi }_{\min })$, b_ν, ref, and σ_ref) to the ECDF of L_p, obs | L_p, obs > l_p, GRG. We note that n_GRG does not depend on l_min, which drops out through the division. However, n_GRG does depend on cosmological parameters through the relation between cosmological redshift z and radial comoving distance r.

A.10. GRG lobe volume-filling fraction

Assuming self-similar growth, the combined proper volume V of an RG’s lobes and its intrinsic proper length l obey V ∝ l³. The constant of proportionality varies per RG and depends on the shape of the lobes; we treat it as an RV $\mathrm{{\rm Y}}:=\frac{V}{L^3}$. Then the proper VFF of GRG lobes VFF_GRG(z) = VFF_GRG(z = 0)⋅(1 + z)³, where

$$\begin{array}{cc}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0):=& n_{\mathrm{GRG}}·\mathbb{E}[V|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]\hfill \\ \hfill =& n_{\mathrm{GRG}}·\mathbb{E}[\frac{V}{L^3}·L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]\hfill \\ \hfill =& n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}·L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}].\hfill \end{array}$$(A.96)

Assuming that RGs grow self-similarly, so that shape does not reveal length, Υ and L³ are conditionally independent given L_p > l_p, GRG: Υ ⫫ L³ | L_p > l_p, GRG. As a result,

$$\begin{array}{cc}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)& =n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]·\mathbb{E}[L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]\hfill \\ \hfill & =n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}]·\mathbb{E}[L^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}].\hfill \end{array}$$(A.97)

To obtain this last line, we once more exploit self-similarity: Υ | (L_p > l_p, GRG) = Υ. We can approximate 𝔼[Υ] by taking the mean of some $\frac{V}{l^3}$ deduced from observations. A technical complication arises from the fact that 𝔼[L³ | L_p > l_p, GRG] does not exist for ξ ≥ −4 under our model. This is an artefact of the Pareto distribution assumption for L, which unrealistically features support over an infinitely long part of the real line: {l ∈ ℝ | l > l_min}. This causes the expectation value integral to diverge for ξ ≥ −4. An approximation to VFF_GRG(z = 0) that works for all ξ is the lower bound

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)>n_{\mathrm{GRG}}·\mathbb{E}[\mathrm{{\rm Y}}]·{\mathbb{E}}^3[L|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}],\end{array}$$(A.98)

which follows from Jensen’s inequality. Here

$$\begin{array}{c}\hfill \mathbb{E}[L|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]=l_{\mathrm{p},\mathrm{GRG}}\frac{\mathrm{\Gamma }(-\frac{\xi }{2}-1)\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})\mathrm{\Gamma }(-\frac{\xi }{2}+\frac{1}{2})}.\end{array}$$(A.99)

Alternative approximation formulae, which use ${\mathrm{{\rm Y}}}_{\mathrm{p}}:=\frac{V}{L_{\mathrm{p}}^3}$, are

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)>n_{\mathrm{GRG}}·\mathbb{E}[{\mathrm{{\rm Y}}}_{\mathrm{p}}]·{\mathbb{E}}^3[L_{\mathrm{p}}|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}]\end{array}$$(A.100)

and

$$\begin{array}{c}\hfill {\mathrm{VFF}}_{\mathrm{GRG}}(z=0)>n_{\mathrm{GRG}}·\mathbb{E}[{\mathrm{{\rm Y}}}_{\mathrm{p}}]·m_{L_{\mathrm{p}}^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}},\end{array}$$(A.101)

where

$$\begin{array}{c}\hfill m_{L_{\mathrm{p}}^3|L_{\mathrm{p}}>l_{\mathrm{p},\mathrm{GRG}}}=l_{\mathrm{p},\mathrm{GRG}}^3·2^{\frac{3}{-(\xi +1)}}\end{array}$$(A.102)

is the median of the cubed projected proper length for giants. An advantage of these latter expressions is that there are more data available to estimate 𝔼[Υ_p] than there are to estimate 𝔼[Υ].

A.10. Unification model constraints from quasar and non-quasar giants

The quasar GRG probability p_Q is

$$\begin{array}{cc}\hfill p_{\mathrm{Q}}:=& \frac{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}},sin\mathrm{\Theta }\le sin{\theta }_{\max })}{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}})}\hfill \\ \hfill =& \mathbb{P}(sin\mathrm{\Theta }\le sin{\theta }_{\max }|L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}})\hfill \\ \hfill =& {\displaystyle {\int }_0^{sin{\theta }_{\max }}}\frac{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}}|sin\mathrm{\Theta }=x)f_{sin\mathrm{\Theta }}(x)\mathrm{d}x}{\mathbb{P}(L_{\mathrm{p},\mathrm{obs}}\ge l_{\mathrm{p},\mathrm{GRG}})}\hfill \\ \hfill =& {\displaystyle {\int }_0^{sin{\theta }_{\max }}}\frac{\mathbb{P}(L_{\mathrm{p}}\ge l_{\mathrm{p},\mathrm{GRG}}|sin\mathrm{\Theta }=x)f_{sin\mathrm{\Theta }}(x)\mathrm{d}x}{\mathbb{P}(L_{\mathrm{p}}\ge l_{\mathrm{p},\mathrm{GRG}})}\hfill \\ \hfill =& {\displaystyle {\int }_0^{sin{\theta }_{\max }}}\frac{(1-F_L\left(\frac{l_{\mathrm{p},\mathrm{GRG}}}{x}\right))f_{sin\mathrm{\Theta }}(x)\mathrm{d}x}{1-F_{L_{\mathrm{p}}}(l_{\mathrm{p},\mathrm{GRG}})}.\hfill \end{array}$$(A.103)

Now using Eqs. A.3 and A.8 and $f_{sin\mathrm{\Theta }}(x)=\frac{x}{\sqrt{1-x^2}}$ over the domain of integration,

$$\begin{array}{c}\hfill p_{\mathrm{Q}}=\frac{4\mathrm{\Gamma }(-\frac{\xi }{2}+1)}{-(\xi +1)\sqrt{\pi }\mathrm{\Gamma }(-\frac{\xi }{2}-\frac{1}{2})}{\displaystyle {\int }_0^{sin{\theta }_{\max }}}\frac{x^{-\xi }\mathrm{d}x}{\sqrt{1-x^2}}.\end{array}$$(A.104)

A.11. Extreme giants in a sample

An interesting feature of the model is its ability to predict the occurrence of giants with extreme projected proper lengths in a GRG sample of, say, size N. Now consider some l_p > l_p, GRG — what is the probability p_{> lp} that an observed GRG will have a projected proper length exceeding l_p? Proceeding as in the derivation of Eq. A.12, we find

$$\begin{array}{cc}\hfill p_{>l_{\mathrm{p}}}\left(l_{\mathrm{p}}\right):=& \mathbb{P}(L_{\mathrm{p},\mathrm{obs}}>l_{\mathrm{p}}|L_{\mathrm{p},\mathrm{obs}}>l_{\mathrm{p},\mathrm{GRG}})\hfill \\ \hfill =& \frac{1-F_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p}}\right)}{1-F_{L_{\mathrm{p},\mathrm{obs}}}\left(l_{\mathrm{p},\mathrm{GRG}}\right)}.\hfill \end{array}$$(A.105)

In the absence of selection effects, p_{> lp}(l_p) is given by Eq. A.15. The number of giants with extreme projected proper lengths N_{> lp} ∼ Binom(N, p_{> lp}(l_p)).

Interesting questions can be answered readily. For example, the probability that the sample contains at least one GRG with projected proper length l_p or larger, is

$$\begin{array}{c}\hfill \mathbb{P}(N_{>l_{\mathrm{p}}}\ge 1)=1-\mathbb{P}(N_{>l_{\mathrm{p}}}=0)=1-{(1-p_{>l_{\mathrm{p}}}\left(l_{\mathrm{p}}\right))}^N.\end{array}$$(A.106)

Appendix B: Additional images

In this appendix, as in Fig. 12, we show newly discovered giants. These cover the projected length range l_p ∈ [0.7 Mpc, 4 Mpc).

Appendix C: Stellar and supermassive black hole masses

In Fig. C.1, we present the SDSS-derived relations between host stellar mass and projected proper length, and between host supermassive black hole mass and projected proper length, for all giants in our final catalogue. We obtained the data as in Section 3.7 of Oei et al. (2022a).

Appendix D: Surface brightness prior

Consider a radio galaxy at cosmological redshift z of intrinsic proper length l bounded by spherical lobes of radius R and spectral index α. If a fraction f_l of the radio galaxy’s central axis lies within the lobes, then

$$\begin{array}{c}\hfill R=\frac{1}{4}{f}_l·l.\end{array}$$(D.1)

If a fraction f_Lν of the radio galaxy’s total luminosity density L_ν comes from the lobes, then the monochromatic emission coefficient

$$\begin{array}{c}\hfill j_{\nu }=\frac{L_{\nu }}{4\pi }\frac{f_{L_{\nu }}}{2}\frac{1}{V},\end{array}$$(D.2)

where the lobe volume $V=\frac{4}{3}\pi R^3$. A formula of practical value features projected proper length l_p instead of l. Given the approximate nature of our approach, we therefore simply assume l ≈ 𝔼[D](η(f_l)) ⋅ l_p, with the deprojection factor expectation 𝔼[D] given in Eq. A.29; $\eta (f_l)=\frac{f_l}{2-f_l}$. The maximum surface brightness of the lobes as seen by an observer is

$$\begin{array}{c}\hfill b_{\nu ,\max }=\frac{3f_{L_{\nu }}·L_{\nu }}{{\pi }^2·\mathbb{E}[D](\eta (f_l))·f_l^2·l_{\mathrm{p}}^2·{(1+z)}^{3-\alpha }},\end{array}$$(D.3)

valid for the line of sight that pierces through a lobe along a diameter. The average line of sight length d within a lobe is smaller than d_max = 2R, though, and given by

$$\begin{array}{c}\hfill ⟨d⟩:=\frac{{\int }_0^{R}2\sqrt{R^2-x^2}2\pi x\mathrm{d}x}{\pi R^2}=\frac{4}{3}R.\end{array}$$(D.4)

Therefore, the mean surface brightness of the lobes as seen by an observer is

$$\begin{array}{c}\hfill ⟨b_{\nu }⟩=\frac{⟨d⟩}{d_{\max }}{b}_{\nu ,\max }=\frac{2}{3}{b}_{\nu ,\max }.\end{array}$$(D.5)

Appendix E: Likelihood function

In Table E.1, we present maximum likelihood and likelihood mean and standard deviation estimates for the inference described in Section 4.1.4; one may compare the results to those in Table 4. In Fig. E.1, we visually summarise the likelihood function; one may compare to Fig. 16.

The strong degeneracy between ξ and b_ν, ref, directly apparent from the central two-parameter marginal in the leftmost column of Fig. E.1, translates to a ridge of essentially constant likelihood that extends from ξ ≈ −4 to ξ ≈ −2. Compared to a non-degenerate case, this makes the maximum likelihood parameters both intrinsically less meaningful and more prone to numerical approximation errors.

Appendix F: Properties of newly discovered giants

Table 2 provides properties of the 50 projectively longest giants discovered during this work’s LoTSS DR2 search campaign. We share these data, alongside those for the other 2010 (98%) giants in our sample, in Flexible Image Transport System (FITS) format through the Centre de Données astronomiques de Strasbourg (CDS).

For our final catalogue, which also includes literature giants, please contact the authors.

All Tables

All Figures