Constraining the giant radio galaxy population with machine learning and Bayesian inference
Mostert & Oei et al.
Giant radio galaxy constraints from machine learning and Bayesian inference

Rafaël I.J. MostertThese authors contributed equally to this article.E-mail<EMAIL_ADDRESS>oei@strw.leidenuniv.nl1,2 Martijn S.S.L. Oei[1][2]1,3 B. Barkus4 Lara Alegre5 Martin J. Hardcastle6Kenneth J. Duncan5Huub J.A. Röttgering1Reinout J. van Weeren1 Maya Horton6

Leiden Observatory, Leiden University, Einsteinweg 55, 2333 CA Leiden, The Netherlands ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands Cahill Center for Astronomy and Astrophysics, California Institute of Technology, 1216 E California Blvd, CA 91125 Pasadena, USA School of Physical Sciences, The Open University, Walton Hall, Milton Keynes, MK7 6AA, UK Institute for Astronomy, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK

Large-scale sky surveys at low frequencies, such as the LOFAR Two-metre Sky Survey (LoTSS), allow for the detection and characterisation of unprecedented numbers of giant radio galaxies (GRGs, or `giants', of at least $l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$ long).
This, in turn, enables us to study giants in a cosmological context.
A tantalising prospect of such studies is a measurement of the contribution of giants to cosmic magnetogenesis.
However, this measurement requires en masse radio–optical association for well-resolved radio galaxies and a statistical framework to infer GRG population properties.

By automating the creation of radio–optical catalogues, we aim to significantly expand the census of known giants.
With the resulting sample and a forward model that takes into account selection effects, we aim to constrain their intrinsic length distribution, number density, and lobe volume-filling fraction (VFF) in the Cosmic Web.

We combined five existing codes into a single machine learning (ML)–driven pipeline that automates radio source component association and optical host identification for well-resolved radio sources.
We created a radio–optical catalogue for the entire LoTSS Data Release 2 (DR2) footprint and subsequently selected all sources that qualify as possible giants.
We combined the list of ML pipeline GRG candidates with an existing list of LoTSS DR2 crowd-sourced GRG candidates and visually confirmed or rejected all members of the merged sample.
To infer intrinsic GRG properties from GRG observations, we developed further a population-based forward model and constrained its parameters using Bayesian inference.

Roughly half of all GRG candidates that our ML pipeline identifies indeed turn out to be giants upon visual inspection, whereas the success rate is 1 in 11 for the previous best giant-finding ML technique in the literature.
We confirm $5,647$ previously unknown giants from the crowd-sourced LoTSS DR2 catalogue and $2,597$ previously unknown giants from the ML pipeline.
Our confirmations and discoveries bring the total number of known giants to at least $11,585$.
Our intrinsic GRG population forward model provides a good fit to the data.
The posterior indicates that the projected lengths of giants are consistent with a curved power law probability density function whose initial tail index $\xi(l_\mathrm{p,GRG}) = -2.8\pm0.2$ changes by $\Delta\xi = -2.4 \pm 0.3$ over the interval up to $l_\mathrm{p} = 5\ \mathrm{Mpc}$.
We predict a comoving GRG number density $n_\mathrm{GRG} = 13 \pm 10\ (100\ \mathrm{Mpc})^{-3}$, close to a recent estimate of the number density of luminous non-giant radio galaxies.
With the projected length distribution, number density, and additional assumptions, we derive a present-day GRG lobe VFF $\mathcal{V}_\mathrm{GRG-CW}(z=0) = 1.4 \pm 1.1 \cdot 10^{-5}$ in clusters and filaments of the Cosmic Web.

We present a state-of-the-art ML-accelerated pipeline for finding giants, whose complex morphologies, arcminute extents, and radio-emitting surroundings pose challenges.
Our data analysis suggests that giants are more common than previously thought.
More work is needed to make GRG lobe VFF estimates reliable, but tentative results imply that it is possible that magnetic fields once contained in giants pervade a significant (${\gtrsim}10\%$) fraction of today's Cosmic Web.

§ INTRODUCTION

Recent radio Stokes $I$ imaging and rotation measure observations show that filaments of the Cosmic Web are magnetised <cit.> with $B \sim 10^0$–$10^2\ \mathrm{nG}$ <cit.>.
However, the origin of these magnetic fields remains highly uncertain.
In a primordial magnetogenesis scenario <cit.>, the seeds of intergalactic magnetic fields can be traced to the Early Universe.
This scenario is not problem-free: primordial magnetic fields that arise before the end of inflation are typically too weak to match observations, while fields that arise after inflation (but before recombination) typically have coherence lengths that are too small.
Alternatively, in an astrophysical magnetogenesis scenario, the seeds of intergalactic magnetic fields are predominantly spread by energetic astrophysical phenomena in the more recent Universe, such as radio galaxies (RGs) and supernova-driven winds <cit.>.
In this latter scenario, giant radio galaxies (GRGs, or `giants') may play a significant role in the magnetisation of the intergalactic medium (IGM), as their associated jets can carry magnetic fields of strength $B \sim 10^2\ \mathrm{nG}$ from host galaxies to cosmological, megaparsec-scale distances <cit.>.

Efforts to measure the contribution of giants to astrophysical magnetogenesis in filaments of the Cosmic Web have only recently begun, with the advent of systematically processed, sensitive, low-frequency sky surveys such as the Low Frequency Array <cit.> Two-metre Sky Survey <cit.>.
By carrying out both a manual search for giants in LoTSS DR2 [Shimwell et al., 2022] pipeline products and a rigorous statistical analysis, Oei et al., 2023 inferred a key statistic: the volume-filling fraction (VFF) of GRG lobes within clusters and filaments of the Local Universe, $\mathcal{V}_\mathrm{GRG-CW}(z = 0)$.
However, much uncertainty remains as to its value, which is set by the GRG number density, GRG length distribution, and GRG length–lobe volume relation.

As the number of observed radio galaxies rapidly increases with decreasing angular length, the time required to manually associate radio source components and identify host galaxies in the optical logically increases.
Machine learning (ML)–based methods have the potential to massively accelerate the detection of various radio source classes by complementing and eventually replacing manual methods <cit.>.
The potential for detecting giants was demonstrated by Dabhade et al., 2020, who visually inspected the $1,600$ ML-identified GRG candidates of Proctor, 2016 and thereby discovered $151$ giants.
By combining into a single pipeline multiple ML-based and rule-based algorithms that automate both the radio component association and the optical host galaxy identification, we aim to improve upon the $9\%$ precision of Proctor, 2016's ML predictions.

As part of the present study, we constructed a LoTSS DR2 GRG sample of unparalleled size, by combining results from a visual search by astronomers [Oei et al., 2023], a visual search by citizen scientists [Hardcastle et al., 2023], and a ML-accelerated search (this article; Sect. <ref>).
With a definitive LoTSS DR2 GRG sample in hand, we refined the Bayesian forward model presented in Oei et al., 2023, and finally constrained several key geometric quantities pertaining to giants.

In Sect. <ref>, we briefly review, generalise, and develop the statistical GRG geometry theory of Oei et al., 2023.
In Sect. <ref>, we introduce the LoTSS DR2 data in which we searched for giants.
In Sect. <ref>, we describe the methods that we used to build our definite LoTSS DR2 GRG sample, and explain how we used the theory of Sect. <ref> in practice to infer GRG quantities of interest.
In Sect. <ref>, we present our findings regarding the projected proper length distribution for giants, their comoving number density, and their instantaneous lobe VFF in clusters and filaments of the Cosmic Web.
In Sect. <ref>, we discuss caveats of the present work, compare our results with previous results, and propose directions for future work, before we conclude in Sect. <ref>.

We assume a flat $\mathrm{\Lambda CDM}$ model with
parameters from Planck Collaboration et al., 2020: $h=0.6766$, $\Omega_\mathrm{BM,0}=0.0490$, $\Omega_\mathrm{M,0}=0.3111$, and $\Omega_\mathrm{\Lambda,0}=0.6889$,
where $H_0 := h \cdot 100$ $\mathrm{km}$ $\mathrm{s}^{-1}$ $\mathrm{Mpc}^{-1}$.
We define giants as radio galaxies with a projected proper[In Cosmic Web filament environments, where giants appear most common [Oei] et al., 2023], lobes may expand along the Hubble flow, rendering their proper and comoving extents different. To avoid ambiguity, we stress that our projected lengths are proper instead of comoving. A less precise synonym for `projected proper length' often found in the literature is `largest linear size'. length $l_\mathrm{p} \geq l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$.
We define the spectral index $\alpha$ such that it relates to flux density $F_\nu$ at frequency $\nu$ as $F_\nu \propto \nu^\alpha$; under this convention, $\alpha < 0$ at radio frequencies for which synchrotron self-absorption is negligible.

§ THEORY

To infer the intrinsic length distribution, number density, and lobe VFF of giants, we used a Bayesian forward modelling approach that incorporates selection effects.
We adopt the framework described in Oei et al., 2023, but generalise a few key formulae.
Furthermore, in a change that allows for the extraction of tighter parameter constraints from the data, we now predict joint projected proper length–redshift histograms rather than projected proper length distributions.

§.§ RG total and projected proper lengths

The central geometric quantity predicted by models of radio galaxy evolution <cit.> is, simply, the RG's intrinsic proper length $l$.
Once the probability distribution of the intrinsic proper length random variable (RV) $L$ is known, one can estimate other geometric quantities of interest, such as the VFF of RG lobes in the Cosmic Web.
However, for the vast majority of observed RGs only a projected proper length $l_\mathrm{p}$ is available, as accurate measurements of jet inclination angles $\theta$ are currently challenging.
In order to fit statistical models to data from surveys such as LoTSS DR2, models should therefore predict the distribution of the projected proper length RV $L_\mathrm{p}$.

§.§ GRG projected proper length: General

We now show, first without adopting a specific parametric form for the distribution of $L$, how the cumulative density function (CDF) and probability density function (PDF) of the GRG projected proper length RV $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ can be calculated.
In particular, we suppose that $L$ has support from some length $l_\mathrm{min} \geq 0$ onwards.
Then $L_\mathrm{p} = L \sin{\Theta}$, where $\Theta$ is the inclination angle RV.
Assuming that, at least on cosmological scales, all RG orientations in three dimensions are equally likely,[This assumption is admissible, because even if the relative orientations of RGs and filaments are not random <cit.>, the uniformity of filament orientations on large scales leads to uniform RG orientations on large scales.] the CDF of $L_\mathrm{p}$ relates to the PDF of $L$ via

\begin{align}
F_{L_\mathrm{p}}(l_\mathrm{p}) =
\begin{cases}
0 & \text{if } l_\mathrm{p} \leq 0;\\
1 - \int_{l_\mathrm{min}}^\infty \sqrt{1 - \left(\frac{l_\mathrm{p}}{l}\right)^2}\ f_L(l)\ \mathrm{d}l & \text{if } 0 < l_\mathrm{p} \leq l_\mathrm{min};\\
1 - \int_{l_\mathrm{p}}^\infty \sqrt{1 - \left(\frac{l_\mathrm{p}}{l}\right)^2}\ f_L(l)\ \mathrm{d}l & \text{if } l_\mathrm{p} > l_\mathrm{min}.\\
\end{cases}
\label{eq:projectedLengthCDF}
\end{align}

We note that, in the usual scenario of $l_\mathrm{min} = 0$, the second case disappears.
Equation <ref> generalises Eq. A.8 from Oei et al., 2023; its derivation closely follows the one presented there.

The CDF of the GRG projected proper length RV $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ is

\begin{align}
F_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p}) =
\begin{cases}
0 & \text{if } l_\mathrm{p} < l_\mathrm{p,GRG};\\
1 - \frac{\int_{l_\mathrm{p}}^\infty \sqrt{1 - \left(\frac{l_\mathrm{p}}{l}\right)^2}\ f_L(l)\ \mathrm{d}l}{\int_{l_\mathrm{p,GRG}}^\infty \sqrt{1 - \left(\frac{l_\mathrm{p,GRG}}{l}\right)^2}\ f_L(l)\ \mathrm{d}l} & \text{if } l_\mathrm{p} \geq l_\mathrm{p,GRG}.
\end{cases}
\end{align}

This result follows from combining Eq. <ref> and Eq. A.12 from Oei et al., 2023.[We have also assumed that $l_\mathrm{p,GRG} > l_\mathrm{min}$, which is the obvious case to consider.]
As PDFs follow from CDFs by differentiation, we find that the PDFs of $L_\mathrm{p}$ and $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ are related by

\begin{align}
f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p}) =
\begin{cases}
0 & \text{if } l_\mathrm{p} < l_\mathrm{p,GRG};\\
\frac{f_{L_\mathrm{p}}(l_\mathrm{p})}{\int_{l_\mathrm{p,GRG}}^\infty \sqrt{1 - \left(\frac{l_\mathrm{p,GRG}}{l}\right)^2}\ f_L(l)\ \mathrm{d}l} & \text{if } l_\mathrm{p} \geq l_\mathrm{p,GRG}.
\end{cases}
\label{eq:projectedLengthPDFGRG}
\end{align}

We note that, throughout the support of $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$, $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$ and $f_{L_\mathrm{p}}(l_\mathrm{p})$ are directly proportional – the quantity in the denominator of Eq. <ref> is merely a normalisation constant.[This is an example of a more general rule: for any RV $X$,

\begin{align}
f_{X\ \vert\ X \geq y}(x) = \begin{cases}
0 & \text{if }x < y;\\
\frac{f_X(x)}{1 - F_X(y)} & \text{if } x \geq y.
\end{cases}
\end{align}

To find $f_{L_\mathrm{p}}(l_\mathrm{p})$ if $l_\mathrm{p} > l_\mathrm{min}$, it is helpful to perform a change of variables.
By defining $\eta \coloneqq \frac{l}{l_\mathrm{p}}$, we rewrite

\begin{align}
F_{L_\mathrm{p}}(l_\mathrm{p}) = 1 - l_\mathrm{p}\int_1^\infty \sqrt{1 - \frac{1}{\eta^2}}\ f_L(l_\mathrm{p} \eta)\ \mathrm{d}\eta\ \ \ \ \text{if } l_\mathrm{p} > l_\mathrm{min}.
\end{align}

This form has the advantage that – within the integral – $l_\mathrm{p}$ occurs only in the integrand, whereas the form of Eq. <ref> features $l_\mathrm{p}$ in both the integrand and in the lower integration limit.
By differentation,

\begin{align}
f_{L_\mathrm{p}}(l_\mathrm{p}) = &-\int_1^\infty \sqrt{1 - \frac{1}{\eta^2}}\ f_L(l_\mathrm{p} \eta)\ \mathrm{d}\eta\nonumber\\
&-l_\mathrm{p}\int_1^\infty \sqrt{1 - \frac{1}{\eta^2}}\ \frac{\mathrm{d}f_L(l_\mathrm{p}\eta)}{\mathrm{d}l_\mathrm{p}}\ \mathrm{d}\eta\ \ \ \ \text{if } l_\mathrm{p} > l_\mathrm{min}.
\end{align}

To arrive at concrete expressions for the GRG projected proper length PDF of Eq. <ref>, we must choose a specific parametric form for the distribution of $L$ or $L_\mathrm{p}$.

§.§ GRG projected proper length: Curved power law

Oei et al., 2023 show that models that assume a Paretian tail for the RG intrinsic proper length distribution, and that include angular and surface brightness selection effects, can tightly reproduce the observed GRG projected proper length distribution.
The PDF of a Pareto-distributed RV is a simple power law, which is fully specified by a lower cut-off $l_\mathrm{min}$ and a tail index $\xi$.
However, there is a good reason to believe that the true GRG projected proper length PDF deviates from simple power law behaviour.
The true RG projected proper length PDF $f_{L_\mathrm{p}}$ will peak around a value set by the typical jet power, environment, lifetime, and inclination angle (amongst other properties).
Below this value, $f_{L_\mathrm{p}}$ will necessarily be an increasing function of $l_\mathrm{p}$; above this value, $f_{L_\mathrm{p}}$ will be a decreasing function.[This line of reasoning implicitly assumes that the distribution of $L_\mathrm{p}$ is unimodal.]
As giants embody the large-length tail of the distribution of $L_\mathrm{p}$, it is likely that the slope of $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$ first becomes more negative (and later becomes less negative) as $l_\mathrm{p}$ increases.

To remain close to the seemingly effective Pareto assumption of Oei et al., 2023, we assume in this work that, at least for $l_\mathrm{p} \geq l_\mathrm{p,GRG}$, the RG projected proper length PDF is a curved power law:

\begin{align}
f_{L_\mathrm{p}}(l_\mathrm{p}) \propto \left(\frac{l_\mathrm{p}}{l_\mathrm{p,GRG}}\right)^{\xi(l_\mathrm{p})}\ \ \ \ \text{if } l_\mathrm{p} \geq l_\mathrm{p,GRG},
\label{eq:RGProjectedProperLength}
\end{align}

where the exponent

\begin{align}
\xi(l_\mathrm{p}) \coloneqq \xi(l_\mathrm{p,1}) + \frac{l_\mathrm{p} - l_\mathrm{p,1}}{l_\mathrm{p,2} - l_\mathrm{p,1}}\left(\xi(l_\mathrm{p,2}) - \xi(l_\mathrm{p,1})\right)
\label{eq:exponentFunction}
\end{align}

is a linear function of $l_\mathrm{p}$.
As long as $l_\mathrm{p,1} \neq l_\mathrm{p,2}$, both projected proper length constants can be chosen arbitrarily; however, $l_\mathrm{p,1} \coloneqq l_\mathrm{p,GRG}$ seems to be a natural choice.
Adopting this choice, and defining $\Delta\xi \coloneqq \xi(l_\mathrm{p,2}) - \xi(l_\mathrm{p,1})$, leads to the final exponent formula

\begin{align}
\xi(l_\mathrm{p}) = \xi(l_\mathrm{p,GRG}) + \frac{l_\mathrm{p} - l_\mathrm{p,GRG}}{l_\mathrm{p,2} - l_\mathrm{p,GRG}}\Delta\xi.
\label{eq:exponentFunction2}
\end{align}

We adopted $\xi(l_\mathrm{p,GRG})$ and $\Delta\xi$ as two parameters of our model.
We furthermore chose $l_\mathrm{p,2} \coloneqq 5\ \mathrm{Mpc}$, which is close to the largest currently known radio galaxy projected proper length [Oei et al., 2022, Oei et al., 2023].
Being the first-order Taylor polynomial of an arbitrary function $\xi(l_\mathrm{p})$ at $l_\mathrm{p,GRG}$, Eq. <ref> represents a natural generalisation of the constant tail index assumption of Oei et al., 2023.
In particular, if model parameter $\Delta\xi = 0$, we recover the earlier Paretian model.

By the same reasoning as before, we find that if the RG projected proper length PDF is a curved power law for $l_\mathrm{p} \geq l_\mathrm{p,GRG}$, then the GRG projected proper length PDF is also a curved power law over this range:

\begin{align}
f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p}) \propto \left(\frac{l_\mathrm{p}}{l_\mathrm{p,GRG}}\right)^{\xi(l_\mathrm{p})}\ \ \ \ \text{if } l_\mathrm{p} \geq l_\mathrm{p,GRG}.
\label{eq:GRGProjectedProperLengthCurved}
\end{align}

The factors required to normalise $f_{L_\mathrm{p}}(l_\mathrm{p})$ and $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$ can be obtained numerically.

Whereas Oei et al., 2023 parametrised $f_L(l)$ and derived $f_{L_\mathrm{p}}(l_\mathrm{p})$ and $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$, we now parametrise $f_{L_\mathrm{p}}(l_\mathrm{p})$ and derive only $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$.
It is possible to start modelling at the level of $f_L(l)$, also in the context of curved power law PDFs, but the resulting expressions for $f_{L_\mathrm{p}}(l_\mathrm{p})$ and $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$ become tedious and rather uninsightful.
For simplicity, we therefore choose to parametrise $f_{L_\mathrm{p}}(l_\mathrm{p})$; we explore the alternative set-up in Appendix <ref>.

§.§ GRG observed projected proper length

Equation <ref> describes a distribution of GRG projected proper lengths in the absence of observational selection effects.
Unfortunately, this distribution cannot be directly tested against GRG samples obtained from surveys, which are always affected by selection.
For a thorough description and derivation of selection effect modelling in the context of our framework, we refer the reader to Sect. 2.8 and Appendix A.8 of Oei et al., 2023; here, we shall only briefly introduce the expressions that we require.

A key result, adopted from Eq. 21 of Oei et al., 2023, is that the GRG observed projected proper length RV $L_\mathrm{p,obs}\ \vert\ L_\mathrm{p,obs} \geq l_\mathrm{p,GRG}$ can be expressed as

\begin{align}
&f_{L_\mathrm{p,obs}\ \vert\ L_\mathrm{p,obs} \geq l_\mathrm{p,GRG}}(l_\mathrm{p}) =\begin{cases}
0 & \text{if } l_\mathrm{p} < l_\mathrm{p, GRG};\\
\frac{C\left(l_\mathrm{p},z_\mathrm{max}\right)f_{L_\mathrm{p}}\left(l_\mathrm{p}\right)}{\int_{l_\mathrm{p,GRG}}^\infty C\left(l_\mathrm{p}',z_\mathrm{max}\right)f_{L_\mathrm{p}}\left(l_\mathrm{p}'\right)\ \mathrm{d}l_\mathrm{p}'} & \text{if } l_\mathrm{p} \geq l_\mathrm{p, GRG},
\end{cases}
\label{eq:GRGObservedProjectedProperLength}
\end{align}

where $C$ is the completeness function.
More precisely, $C(l_\mathrm{p}, z_\mathrm{max})$ denotes the fraction of all RGs with projected proper length $l_\mathrm{p}$ in the volume up to cosmological redshift $z_\mathrm{max}$ that is detected and identified through the survey considered – in this work, this will be LoTSS DR2.
The repeated factors in numerator and denominator reveal that, in order to compute $f_{L_\mathrm{p,obs}\ \vert\ L_\mathrm{p,obs} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$, we need to know $f_{L_\mathrm{p}}(l_\mathrm{p})$ up to a constant only (and on $l_\mathrm{p} \geq l_\mathrm{p,GRG}$ only).
More concerningly, we also see that selection effects that reduce the completeness by the same factor for all $l_\mathrm{p} \geq l_\mathrm{p,GRG}$ leave no imprint on $f_{L_\mathrm{p,obs}\ \vert\ L_\mathrm{p,obs} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$.
Therefore, such selection effects cannot be constrained by a GRG observed projected proper length analysis alone.

Under the assumption that the RG projected proper length PDF $f_{L_\mathrm{p}}(l_\mathrm{p})$ does not evolve between redshifts $z = z_\mathrm{max}$ and $z = 0$, the completeness function becomes

\begin{align}
C\left(l_\mathrm{p}, z_\mathrm{max}\right) = \frac{\int_0^{z_\mathrm{max}}p_\mathrm{obs}\left(l_\mathrm{p},z\right) r^2\left(z\right) E^{-1}\left(z\right)\ \mathrm{d}z}{\int_0^{z_\mathrm{max}}r^2\left(z\right)E^{-1}\left(z\right)\ \mathrm{d}z},
\label{eq:completeness}
\end{align}

where the observing probability $p_\mathrm{obs}(l_\mathrm{p},z)$ is the probability that an RG of projected proper length $l_\mathrm{p}$ at redshift $z$ is detected by a survey and its subsequent analysis steps (such as the machine learning pipeline considered in this work), $r$ denotes comoving radial distance, and $E(z)$ is the dimensionless Hubble parameter[In a flat Friedmann–Lemaître–Robertson–Walker universe, the dimensionless Hubble parameter $E$ is

\begin{align}
E\left(z\right) \coloneqq& \frac{H\left(z\right)}{H_0} = \sqrt{\Omega_\mathrm{R,0} \left(1+z\right)^4 + \Omega_\mathrm{M,0} \left(1+z\right)^3 + \Omega_\mathrm{\Lambda,0}}.
\end{align}
The appropriate form of $p_\mathrm{obs}(l_\mathrm{p},z)$ is determined by the selection effects relevant to the survey of interest and its analysis.

In this work, we consider GRG lobe surface brightness (SB) selection, which at present renders many members of the GRG population undetectable, and selection by limitations of our analysis, which causes in principle detectable giants to evade identification.
We described the former effect parametrically, and determined the latter effect empirically.
The effects yield functions $p_\mathrm{obs,SB}(l_\mathrm{p},z)$ and $p_\mathrm{obs,ID}(l_\mathrm{p},z)$, respectively, which then combine to form a single observing probability function through

\begin{align}
p_\mathrm{obs}(l_\mathrm{p},z) = p_\mathrm{obs,SB}(l_\mathrm{p},z) \cdot p_\mathrm{obs,ID}(l_\mathrm{p},z).
\end{align}

§.§.§ Selection effects: Surface brightness limit

RG lobes whose SBs are lower than some threshold $b_{\nu,\mathrm{th}}$, which typically equals the survey noise level $\sigma$ times a factor of order unity, cannot be detected.
Following Sect. 2.8.3 of Oei et al., 2023, we modelled SB selection by assuming that the lobe SBs $B_\nu(\nu,l,z)$ at $\nu = \nu_\mathrm{obs}$ of RGs of intrinsic proper length $l = l_\mathrm{ref}$ residing at redshift $z = 0$ are lognormally distributed.
We parametrised $B_\nu(\nu_\mathrm{obs}, l_\mathrm{ref}, 0) = b_{\nu,\mathrm{ref}}\ S$, where $b_{\nu,\mathrm{ref}}$ is the median lobe SB and $S$ is a lognormally distributed RV with median $1$ and dispersion parameter $\sigma_\mathrm{ref}$.
The observing probability due to SB selection then is

\begin{align}
p_\mathrm{obs,SB}\left(l_\mathrm{p},z\right) &= \int_{s_\mathrm{min}}^\infty\sqrt{1 - \left(\frac{s_\mathrm{min}}{s}\right)^{-\frac{2}{\zeta}}}f_S\left(s\right) \mathrm{d}s;\label{eq:probabilityObservedSB}\\
s_\mathrm{min} &= \frac{b_{\nu,\mathrm{th}}}{b_{\nu,\mathrm{ref}}}\left(\frac{l_\mathrm{p}}{l_\mathrm{ref}}\right)^{-\zeta}\left(1+z\right)^{3-\alpha};\\
f_S\left(s\right) &= \frac{1}{\sqrt{2\pi}\sigma_\mathrm{ref} s}\exp{\left(-\frac{\ln^2s}{2\sigma_\mathrm{ref}^2}\right)}.
\end{align}

Here, $\alpha$ is the typical RG lobe spectral index, which we assumed fixed at $\alpha = -1$.
The exponent $\zeta$ determines how the SB distribution scales with projected proper length $l_\mathrm{p}$.

In a departure from Oei et al., 2023, we did not fix $\zeta = -2$, but rather left $\zeta$ a free parameter which we fitted to the data.
Deviations from $\zeta = -2$ occur in at least two cases: when giant growth is not shape-preserving, and if the radio luminosity distributions of giants of different $l_\mathrm{p}$ are distinct.
Dynamical models of RGs in general predict that both cocoons <cit.> and lobes <cit.> change shape over time, and in a jet power–dependent way.
There remains considerable uncertainty as to how shapes change throughout the giant phase: axial ratio–like measures generally show that RG lobes become more elongated during growth, but this trend could possibly reverse for giants, whose lobes might protrude from the clusters and filaments in which they are born.
Simulations suggest that, for such protrusions, the usual constant power-law profile assumptions for the ambient baryon density and temperature break down <cit.>.
If lobes of giants widen over time, then $\zeta$ would decrease.
The second case occurs if the end-of-life lengths of RGs increase with jet power, so that the subpopulation that survives up to some $l_\mathrm{p}$ has its jet power distribution – and therefore its radio luminosity distribution – shifted upwards with respect to subpopulations at smaller $l_\mathrm{p}$.
This effect, which appears plausible given models <cit.>, would increase $\zeta$.
At present, it seems hard to predict the net result on $\zeta$ of these counteracting effects.

§.§.§ Selection effects: Non-identification

Every present-day survey search method (such as visual inspection by scientists, visual inspection by citizen scientists, and ML-based approaches) will fail to identify some giants that are in principle identifiable (in the sense that they lie above the detection threshold set by the noise).
For automated approaches, such as the ML-based approach presented in this work, identification can become more challenging for larger angular lengths $\phi$: one reason being the increased number of unrelated, interloping radio sources that cover the solid angle occupied by the RG.
We call the probability that an identifiable RG is indeed identified – and therefore becomes part of the final sample – $p_\mathrm{obs,ID}$.

Say we have $M$ methods to search for giants in the same survey.
Let $\mathcal{G} = \{g_1, g_2, ..., g_N\}$ be the set of all identified giants (so that $|\mathcal{G}| = N$), and let $\mathcal{G}_i \subseteq \mathcal{G}$ be the subset identified by method $i$.
Figure <ref> illustrates the set-up.
The projected proper length and cosmological redshift of giant $g$ are $l_\mathrm{p}(g)$ and $z(g)$, respectively.
To determine the identification probability $p_{\mathrm{obs,ID},i}(l_\mathrm{p},z)$ for method $i$, we first assume it to be of logistic form

\begin{align}
p_{\mathrm{obs,ID},i}(l_\mathrm{p},z) = \frac{1}{1 + \exp{(-(\beta_{0,i} + \beta_{l_\mathrm{p},i} \cdot l_\mathrm{p} + \beta_{z,i} \cdot z))}}.
\end{align}

We obtain best-fit parameters $\hat{\beta}_{0,i}$, $\hat{\beta}_{l_\mathrm{p},i}$, and $\hat{\beta}_{z,i}$ by performing binary logistic regression with two explanatory variables on the data set $\mathcal{D}_i$, where

\begin{align}
\mathcal{D}_i \coloneqq \left\{\left(\left[l_\mathrm{p}(g), z(g)\right], \mathbb{I}(g \in \mathcal{G}_i)\right)\ \vert\ g \in \bigcup_{j =1, j \neq i}^M \mathcal{G}_j\right\}.
\end{align}

The first element of each pair in $\mathcal{D}_i$ is a point in projected length–redshift space, whilst the second element is 0 or 1: $\mathbb{I}$ denotes the indicator function.
Qualitatively, $\mathcal{D}_i$ stores for each giant in the union of all GRG subsets except $\mathcal{G}_i$ its projected length–redshift coordinates, together with the success or failure of its identification through method $i$.

The implicit assumption here is that all $ g \in \bigcup_{j =1, j \neq i}^M \mathcal{G}_j$ are typical examples of identifiable giants at the relevant projected proper length and redshift.
We caution that this might not be true: giants with a peculiar morphology, or those lying in parts of the sky where optical identification is hard (e.g. towards the Galactic Plane or crowded regions of large-scale structure), may be identifiable in a radio surface brightness sense, but will nonetheless evade sample inclusion more often than other giants.
As a result, giants that do end up in a sample – such as $\bigcup_{j =1, j \neq i}^M \mathcal{G}_j$ – will have more regular morphologies than giants in general and will lie in regions of the sky where optical identification is easier than for giants in general.
Typically, such giants are also more likely to be found by method $i$, and as a result our approach will probably render $p_{\mathrm{obs,ID},i}$ biased high.

Given a set of $M$ functions $\{p_{\mathrm{obs,ID},i}(l_\mathrm{p},z)\ \vert\ i \in \{1, 2, ..., M\}\}$, several possibilities exist to combine them into a single $p_\mathrm{obs,ID}(l_\mathrm{p},z)$.
At the minimum, $p_\mathrm{obs,ID}(l_\mathrm{p},z)$ is given by a point-wise maximum:

\begin{align}
p_\mathrm{obs,ID}(l_\mathrm{p},z) = \max_{i \in \{1, 2, ..., M\}} p_{\mathrm{obs,ID},i}(l_\mathrm{p},z),
\label{eq:identificationProbabilityCombinationMaximum}
\end{align}

which is appropriate if methods tend to find the same identifiable giants – as in our case.[
In case methods tend to find independent subsets of identifiable giants,

\begin{align}
p_\mathrm{obs,ID}(l_\mathrm{p},z) = 1 - \prod_{i=1}^M (1 - p_{\mathrm{obs,ID},i}(l_\mathrm{p},z)).
\end{align}

We note that it is possible to design methods that find subsets of identifiable giants that have even less overlap than independent subsets have.]

[fill opacity=0.5, draw opacity = 0.25]
[CadetBlue!7!white] (0,0) circle(3.6*);
[CadetBlue!14!white] (0,-.3*) circle(3.2*);

[CadetBlue!50!white](0,.5*) circle(1.85*);
[CadetBlue!50!white] (-0.9899*,-1.49899*) circle(1.45*);
[CadetBlue!50!white] (0.9899*,-1.49899*) circle(1.45*);

[CadetBlue!21!white](0,.5*) circle(1.8*);

[color=CadetBlue!7!white, rotate=-90, postaction=decorate, decoration=text along path, raise=-8pt, text align=align=center, text color=CadetBlue!50!white, text=||All, reverse path] (0,0) circle (3.6*);
[color=CadetBlue!14!white, rotate=-90, postaction=decorate, decoration=text along path, raise=-8pt, text align=align=center, text color=CadetBlue!75!white, text=||Identifiable, reverse path] (0.3*,0) circle (3.2*);
[color=CadetBlue!50!white, line width=0pt, rotate=-90, postaction=decorate, decoration=text along path, raise=-8pt, text align=align=center, text color=CadetBlue!100!white, text=||Identified, reverse path] (-0.5*,0) circle (1.85*);

[CadetBlue!30!white] (-0.9899*,-1.49899*) circle(1.4*);
[CadetBlue!30!white] (0.9899*,-1.49899*) circle(1.4*);
[dashed, CadetBlue!50!white, thin] (0,.5*) circle(1.8*);
[dashed, CadetBlue!50!white, thin] (-0.9899*,-1.49899*) circle(1.4*);
[dashed, CadetBlue!50!white, thin] (0.9899*,-1.49899*) circle(1.4*);

at (0,-3.2*) (G) $\mathcal{G} = \mathcal{G}_1 \cup \mathcal{G}_2 \cup \mathcal{G}_3$;
at (0,.5*) (D) $\mathcal{G}_1$;
at (-0.9899*,-1.49899*) (E) $\mathcal{G}_2$;
at (0.9899*,-1.49899*) (F) $\mathcal{G}_3$;

Schematic of a three-method search for giants.
Of all giants in the survey footprint up to $z = z_\mathrm{max}$, only those for which the lobe surface brightness at the observing frequency $\nu_\mathrm{obs}$ is above detection threshold $b_{\nu,\mathrm{th}}$ are identifiable.
$\mathcal{G}$ denotes the actually identified set of giants.
$\mathcal{G}_1$, $\mathcal{G}_2$, and $\mathcal{G}_3$ are the subsets identified by each method individually.
As an example, we shade $\mathcal{G}_2 \cup \mathcal{G}_3$, which has overlap with $\mathcal{G}_1$, and which can be used to measure $p_\mathrm{obs,ID,1}(l_\mathrm{p},z)$.

§.§ GRG number density

The preceding theory allows us to find the intrinsic, comoving number density of giants, $n_\mathrm{GRG}$, if we know the observed number of giants within a solid angle of extent $\Omega$ and in the volume up to $z_\mathrm{max}$, $N_\mathrm{GRG,obs}(\Omega, z_\mathrm{max})$.
We assume that, up to this redshift, $n_\mathrm{GRG}$ remains constant.
We note that we cannot calculate $n_\mathrm{GRG}$ using Eq. 30 from Oei et al., 2023: this equation assumes $\xi(l_\mathrm{p,1}) = \xi(l_\mathrm{p,2})$.
We derive a more general expression by first noting that the number of giants observed within a solid angle of extent $\Omega$ in the volume up to $z_\mathrm{max}$ and with projected proper lengths between $l_\mathrm{p}$ and $l_\mathrm{p} + \mathrm{d}l_\mathrm{p}$ is

\begin{align}
\mathrm{d}N_\mathrm{GRG,obs}(l_\mathrm{p},\Omega,z_\mathrm{max}) &= \frac{\Omega}{4\pi}\ n_\mathrm{GRG}\ f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})\ \mathrm{d}l_\mathrm{p}\ \cdot\nonumber\\
&\int_0^{z_\mathrm{max}} p_\mathrm{obs}(l_\mathrm{p},z)\ 4\pi r^2(z)\ \frac{\mathrm{d}r}{\mathrm{d}z}\ \mathrm{d}z.
\label{eq:numberOfGiantsDifferential}
\end{align}

\begin{align}
N_\mathrm{GRG,obs}(\Omega, z_\mathrm{max}) = \int_{l_\mathrm{p,GRG}}^\infty \mathrm{d}N_\mathrm{GRG,obs}(l_\mathrm{p},\Omega,z_\mathrm{max}),
\end{align}

we find, by isolating $n_\mathrm{GRG}$, that

\begin{align}
\label{eq:GRGNumberDensity}
&n_\mathrm{GRG}(l_\mathrm{p,GRG}, z_\mathrm{max}) = \frac{H_0}{c}\ \frac{4\pi}{\Omega}\ N_\mathrm{GRG,obs}(\Omega, z_\mathrm{max})\ \cdot\\
&\left(\int_{l_\mathrm{p,GRG}}^\infty f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p}) \int_0^{z_\mathrm{max}} p_\mathrm{obs}(l_\mathrm{p},z)\ 4\pi r^2(z)\ E^{-1}(z)\ \mathrm{d}z\ \mathrm{d}l_\mathrm{p}\right)^{-1}.\nonumber
\end{align}

This expression is valid also beyond the context of power law or curved power law PDFs $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$.
We remark that $n_\mathrm{GRG}$ can depend sensitively on $l_\mathrm{p,GRG}$, the projected proper length used to define giants.
In contrast to the approach of Oei et al., 2023, in this work we do not calculate $n_\mathrm{GRG}$ in a step following inference of the framework's parameters, but rather include it as a parameter to be constrained during inference.

§.§ GRG lobe volume-filling fraction

To constrain the contribution of giants to astrophysical magnetogenesis, we wish to know the volume-filling fraction of their lobes in clusters and filaments of the Cosmic Web, $\mathcal{V}_\mathrm{GRG-CW}$.
This quantity may have changed over cosmic time; in this work, we calculate it at the present day.
First, we make the approximation that all GRG lobes (in the Local Universe) lie in clusters and filaments.
In addition, we model the general RG relation between the two-lobe proper volume RV $V$ and the projected proper length RV $L_\mathrm{p}$ as a power law with scatter:

\begin{align}
V = V_\mathrm{GRG} \cdot \left(\frac{L_\mathrm{p}}{l_{\mathrm{p,GRG}}}\right)^\gamma \cdot X,
\label{eq:VFFModel}
\end{align}

where $V_\mathrm{GRG}$ is the mean two-lobe proper volume of an RG with a projected proper length $l_\mathrm{p,GRG}$, and $\gamma$ is a constant exponent.
Furthermore, $X$, which we take to be independent of $L_\mathrm{p}$, is a non-negative, dimensionless RV with a mean of unity and an otherwise arbitrary distribution.
In Sect. <ref>, we present observations that indicate that this model is reasonable.
Under this model,

\begin{align}
\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}] = \frac{V_\mathrm{GRG}}{l_\mathrm{p,GRG}^\gamma} \cdot \mathbb{E}[L_\mathrm{p}^{\gamma}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}].
\end{align}

Finally, if GRG lobes are sufficiently small or giants sufficiently rare (or both), the probability that there exist overlapping GRG lobes will be low.
If, indeed, GRG lobes do not overlap, $\mathcal{V}_\mathrm{GRG-CW} \propto \mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}]$ and $\mathcal{V}_\mathrm{GRG-CW} \propto n_\mathrm{GRG}$, so that

\begin{align}
\mathcal{V}_\mathrm{GRG-CW}(z) = &\frac{\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}](z)\ n_\mathrm{GRG}(z)\ (1+z)^3}{\mathcal{V}_\mathrm{CW}(z)},
\label{eq:VFF}
\end{align}

where $\mathcal{V}_\mathrm{CW}$ denotes the VFF of clusters and filaments combined.

§.§ GRG angular lengths

An object's angular length $\phi$, projected proper length $l_\mathrm{p}$, and cosmological redshift $z$ are related through

\begin{align}
\phi(l_\mathrm{p},z) = \frac{l_\mathrm{p}(1+z)}{r(z)}.
\end{align}

Due to the expansion of the Universe, there exists a minimum angular length for objects of a given projected proper length.
If one defines giants as RGs with projected proper lengths $l_\mathrm{p} \geq l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$, as in this work, then all giants have an angular length $\phi \geq 1.3'$ [Oei et al., 2023].
This fact has important consequences for GRG search campaigns.
At the LoTSS resolution of $\theta_\mathrm{FWHM} = 6''$, it implies that giants are always resolved and span at least 13 resolution elements.
Therefore, to model the detectability of giants at this resolution, one must consider their surface brightness (profiles), rather than their flux densities.

§.§ Inference

Finally, we describe how the framework's six free parameters $\vec{\theta} \coloneqq [\xi(l_\mathrm{p,GRG}), \Delta\xi, b_{\nu,\mathrm{ref}}, \sigma_\mathrm{ref},\zeta, n_\mathrm{GRG}]$ can be inferred from a data set containing a projected length and redshift for each observed giant.
In particular, we consider a rectangle in projected proper length–cosmological redshift parameter space, within which our model assumptions are expected to hold.
We partition this rectangle into $N_\mathrm{b}$ equiareal bins of width $\Delta l_\mathrm{p}$ and height $\Delta z$.
We denote the coordinates of bin $i$'s centre as $(l_{\mathrm{p},i}, z_i)$.

We binned the data to obtain a two-dimensional histogram.
The number of giants found in bin $i$, $N_i$, is an RV with a Poisson distribution: $N_i \sim \mathrm{Poisson}(\lambda_i)$.
Its expectation $\lambda_i$ depends on the model parameters $\vec{\theta}$.
Assuming that the $\{N_i\}$ are independent, the log-likelihood becomes

\begin{align}
\ln{\mathcal{L}(\{N_i\}\ \vert\ \vec{\theta})} = \sum_{i=1}^{N_\mathrm{b}} N_i \ln{\lambda_i(\vec{\theta})} - \lambda_i(\vec{\theta}) - \ln{(N_i!)}.
\label{eq:logLikelihood}
\end{align}

The last term on the right-hand side of Eq. <ref> is the same for all $\vec{\theta}$, and need not be calculated if one is interested in $\mathcal{L}$ up to a global constant only.[If one includes this term, it only needs to be calculated once. For numerical stability, it is helpful to note that

\begin{align}
-\sum_{i=1}^{N_\mathrm{b}} \ln{(N_i!)} = -\sum_{i=1}^{N_\mathrm{b}}\sum_{j=2}^{N_i}\ln{j}.
\end{align}
Following Eq. <ref>, but avoiding integration over $z$ and assuming narrow bins in both dimensions, we approximate

\begin{align}
\lambda_i \approx n_\mathrm{GRG} V_i \cdot f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_{\mathrm{p},i}) \Delta l_\mathrm{p} \cdot p_\mathrm{obs}(l_{\mathrm{p},i},z_i).
\end{align}

The volume in which the giants of bin $i$ fall, $V_i$, is

\begin{align}
V_i = \Omega r^2(z_i)\Delta r_i,\ \text{with}\ \Delta r_i = \frac{c}{H_0}\frac{\Delta z}{E(z_i)}.
\end{align}

Appendix <ref> details a particularly efficient trick to compute the likelihood for a range of $n_\mathrm{GRG}$, whilst leaving the other parameters fixed.
By multiplying the likelihood function with a prior distribution, for which we chose the uniform distribution, we obtained a posterior distribution over $\vec{\theta}$ – up to a constant.

§ DATA

We applied our automated radio–optical catalogue creation methods to all Stokes $I$ maps from LoTSS DR2 [Shimwell et al., 2022].[LoTSS DR2 is publicly available at <https://lofar-surveys.org/dr2_release.html>.]
The LoTSS DR2 observations cover the $120$–$168\ \mathrm{MHz}$ frequency range, have a $6''$ resolution, a median RMS sensitivity of $83\ \mathrm{\mu Jy\ beam}^{-1}$, and a flux density scale uncertainty of approximately $10\%$.
The observations are split into a region centred at 12h45m $+44\degree 30'$ and a region centred at 1h00m $+28\degree 00'$; both avoid the Galactic Plane.
These regions span $4,178$ and $1,457$ square degrees respectively, and together cover $27\%$ of the Northern Sky.
The observations consist of $841$ partly overlapping pointings with diameters of $4.0\degree$.
The vast majority of the pointings were observed for 8h, all within the May 2014–February 2020 time frame.

The ML pipeline presented in Sect. <ref> does not only rely on LoTSS DR2 Stokes $I$ maps, but also on an infrared–optical source catalogue.
This catalogue contains the positions, magnitudes, and colours of unWISE [Schlafly et al., 2019] infrared sources and of DESI Legacy Imaging Surveys DR9 [Dey et al., 2019] optical sources.

To discover as many giants as possible, we supplemented our ML pipeline's sample of GRG candidates with the GRG candidates from the value-added LoTSS DR2 catalogue [Hardcastle et al., 2023].
For angularly extended ($\phi > 15''$) radio components (and thus all giants), the radio source component association and most of the host galaxy identification for the value-added LoTSS DR2 catalogue were performed via a public project named `Radio Galaxy Zoo: LOFAR' on Zooniverse.
Zooniverse is an online citizen science platform for crowd-sourced visual inspection.[The Zooniverse website is <https://zooniverse.org>.]
We will refer to the value-added LoTSS DR2 catalogue as the `RGZ catalogue' and to the GRG candidates in that catalogue as the `RGZ GRG candidates'.
The detailed source component information provided by the RGZ catalogue allowed us to homogenise the angular length estimates of the ML pipeline GRG candidates and the RGZ GRG candidates (see Sect. <ref>).
After visual confirmation, we supplemented our GRG sample with literature GRG samples (see Sect. <ref>).

§ METHODS

Overview of our approach, which consists of two stages.
In the first stage we built a GRG sample, and in the second stage we inferred the properties of the intrinsic GRG population using a forward model.
The brackets indicate the different parts of our approach and mention the sections containing the corresponding details.

To derive the projected length distribution, number density, and lobe VFF for the intrinsic population of giants, we followed a two-stage approach.
In the first stage, we gathered all giants that we detected in the LoTSS DR2 Stokes $I$ images using our automatic ML pipeline and added all other giants that we found in the RGZ catalogue.
We re-evaluated and homogenised the source size estimates over the combined GRG sample, and manually inspected the plausibility of the associated radio source components and optical/infrared host galaxy.
Finally, we merged this GRG sample with the other GRG samples from the literature.
In the second stage, we search for the most likely parameters for the forward model presented in Sect. <ref> that describe the GRG observed projected proper length distribution and the selection effects of the merged GRG sample.
Figure <ref> shows an overview of our approach.

§.§ Detecting radio emission

We started out with the publicly available calibrated LoTSS DR2 Stokes $I$ images [Shimwell et al., 2022].
For each of the $841$ pointings, we ran the PyBDSF radio blob detection software [Mohan & Rafferty, 2015] using the same parameters as used in LoTSS DR2 [Shimwell et al., 2022] – notably, this means that we used a $5\sigma$ detection threshold.
Appendix <ref> provides the full list of PyBDSF parameters and their values.

The output we generated consists of a list of radio blobs with their location and properties.
PyBDSF can decompose each radio blob it detects, into one or more 2D Gaussians. For each radio blob, we also saved the corresponding list of Gaussians.
These Gaussians function as a source model for each radio blob and will be used in later steps in the ML pipeline.[However, as we discuss in Sect. <ref> these source models are not always adequate for extended, well-resolved radio sources.]

§.§ Calculating radio-to-optical/infrared likelihood ratios

For radio sources, the location of the host galaxy on the sky is close to the flux density–weighted centre of the radio source.[In this context, `close' refers to angular distances comparable to the survey's resolution.][The flux density–weighted centre of a radio component is the sum of the product of position and flux density for each pixel in the area where PyBDSF found significant emission, divided by the area's total flux density. This calculation is performed twice: once for determining the right ascension and once for determining the declination of the centroid. The flux density–weighted centre of a multi-component radio source is the flux density–weighted average of the components' individual flux density–weighted centres.]
The likelihood ratio method, which exploits this idea, quantifies the likelihood that a source in one observing band is the correct counterpart to a source in another observing band <cit.>.
Williams et al., 2019 used this method to cross-match the unresolved – and some resolved – radio sources of LoTSS DR1 to a combined catalogue of infrared and optical sources.
More specifically, the infrared sources came from AllWISE <cit.>, whilst the optical sources came from the Panoramic Survey Telescope and Rapid Response System 1 <cit.> DR1 $3\pi$ steradian survey.
The likelihood ratio function that Williams et al., 2019 used is a function of the angular distance between the flux density–weighted centre of the radio source and the flux density–weighted centre of the optical or infrared source, the magnitude of the optical or infrared source, and the colour of the optical or infrared source.
The likelihood ratio function also takes into account uncertainties in each of these three dependencies.

We adopted the same procedure as detailed by Williams et al., 2019 to cross-match our simple radio sources (where `simple' is to be understood as in Sect. <ref>) to a combined catalogue of infrared and optical sources.
The infrared sources came from unWISE [Schlafly et al., 2019], and the optical sources were now taken from the DESI Legacy Imaging Surveys DR9 [Dey et al., 2019], which boasts deeper imagery than Pan-STARRS1 DR1 used for LoTSS DR1.
The unWISE [Schlafly et al., 2019] and DESI Legacy Imaging Surveys DR9 source catalogues are used for LoTSS DR2 cross-matching more generally [Hardcastle et al., 2023].
Per pointing, we applied the likelihood ratio method to the full list of radio blobs and to the full list of Gaussians.
For both the blobs and the Gaussians, we stored the identifier of the optical or infrared source that produced the highest likelihood ratio, alongside this highest likelihood ratio itself.

§.§ Sorting radio emission with a gradient-boosting classifier

Most radio sources that consist of a single radio blob (mostly unresolved or barely resolved radio sources) can be cross-matched using the likelihood ratio method.
However, some resolved radio sources, and certainly most resolved giants (Sect. <ref>), consist of multiple radio blobs as parametrised by PyBDSF, and therefore require radio blob association and cannot be cross-matched using the likelihood ratio alone.
To separate the simple from the complex radio blobs in LoTSS DR1, a considerable amount of visual inspection was applied [Williams et al., 2019].
For LoTSS DR2, Alegre et al., 2022 trained a gradient-boosting classifier <cit.> to classify radio blobs as either `simple' or `complex' based on the properties of the radio blobs, the properties of the Gaussians fitted to these blobs, the likelihood ratios for each, and the distance to and properties of the nearest neighbours.

We adopted the procedure of Alegre et al., 2022 and use their trained GBC to separate the simple radio blobs from those that require radio component association beyond PyBDSF's capabilities and/or optical host identification beyond the scope of the likelihood ratio method of Sutherland & Saunders, 1992.
We expect most giants to fall in the latter case.

§.§ Associating radio emission into radio sources

We proceeded with automatic radio source component association for the complex radio blobs.
Following the procedure laid out by Mostert et al., 2022, for each of these radio blobs, we created a $300'' \times 300''$ LoTSS DR2 image cutout centred on the radio blob.
Next, a fast region-based convolutional neural network <cit.>, adapted and trained for this purpose by Mostert et al., 2022, was applied to these cutouts to predict which (if any) other radio blobs – whether they be complex or simple – form a single physical structure with the central radio blob.
For example, the two lobes of an RG, each represented by a radio blob, might be associated together to form a single physical radio source.
Due to the fixed $300'' \times 300''$ image size for which the Fast R-CNN was trained, we expect most radio sources that are associated in our pipeline to have an angular length $\phi < 424''$.[If predicted associations from neighbouring cutouts have an overlapping radio blob, the associations will be merged. For example: in cutout 1 lobe A and core B are associated and in cutout 2 core B and lobe C are associated, then the set (lobe A, core B, and lobe C) will enter the catalogue as a single radio source, thereby creating the possibility of detecting radio sources with angular length $\phi > 424''$.]

The result is a radio source catalogue in which some of the radio blobs have been merged, and a component catalogue that lists for each radio blob to which radio source it belongs.
The radio and the component catalogue were completed by appending to them the remaining list of simple radio blobs.

§.§ Identifying host galaxies in the optical and infrared

Barkus et al., 2022 created a method for identifying the optical or infrared host of an extended radio source.
The method described by Barkus et al., 2022 takes the radio morphology into account by drawing a ridgeline along the regions of high flux density.
The method continues with the application of the likelihood ratio method to quantify which pairs of host galaxy candidates and radio sources are a plausible match.
The likelihood ratio $LR$ used in this context follows Eq. 1 of Sutherland & Saunders, 1992, with the slight simplification of having the latter's dependence on two angular offsets replaced by a dependence on a radial angular offset only:

\begin{align}
\label{eq:LR}
LR = \frac{q(m,c) f(r)}{n(m,c)},
\end{align}

where $q(m,c)$ is a prior on the magnitude $m$ and colour $c$ of the optical host, $f(r)$ is a function of the angular offset between the optical centroid and the radio centroid, and $n(m,c)$ normalises for the number density of optical sources with a certain magnitude $m$ and colour $c$ in the catalogue used for the cross-matching.

To adapt the likelihood ratio for use in the case of extended radio sources, Barkus et al., 2022 implemented the different components of the ratio as follows.
For $n(m,c)$, Barkus et al., 2022 estimated the probability density over $m$ and $c$ for a distribution of $50,000$ randomly sampled sources from a combined Pan-STARRS–AllWISE catalogue in the region of the sky that overlapped with LoTSS DR1.
For $q(m,c)$, they estimated the probability over $m$ and $c$ for sources from the combined Pan-STARRS–AllWISE catalogue that were manually selected to be the most likely optical/near-infrared host for a sample of $950$ radio sources with angular length $\phi > 15''$.
For both $n(m,c)$ and $q(m,c)$, the AllWISE W1 magnitudes were used for $m$, the Pan-STARRS i-band magnitudes minus the AllWISE W1 magnitudes were used for colour $c$, and the PDF was formed using a 2D kernel density estimator <cit.> with a Gaussian kernel and a bandwidth of $0.2$.
For extended asymmetric or bent radio galaxies, the optical host is not likely to be found at the radio centroid.
Therefore, Barkus et al., 2022 proposed that $f(r)$ should be a function of both the distance between the radio centroid and the optical source $r_\mathrm{opt,centroid}$ and the smallest distance between the optical source and a ridgeline fitted to the radio source $r_\mathrm{opt,ridge}$.

\begin{equation}
f(r)=f_\mathrm{ridge}(r_\mathrm{opt,ridge})\cdot f_\mathrm{centroid}(r_\mathrm{opt,centroid}),
\end{equation}

\begin{equation} \label{eq:fridge}
f_\mathrm{ridge}(r_\mathrm{opt,ridge})=\frac{1}{2 \pi \sigma_r^2} e^{\frac{-r_\mathrm{opt,ridge}^2}{2\sigma_r^2}},
\end{equation}

\begin{equation} \label{eq:fcentroid}
f_\mathrm{centroid}(r_\mathrm{opt,centroid})=\frac{1}{2 \pi \sigma_\mathrm{c}^2} e^{\frac{-r_\mathrm{opt,centroid}^2}{2\sigma_\mathrm{c}^2}},
\end{equation}

where $\sigma_r^2=\sigma_\mathrm{opt}^2 + \sigma_\mathrm{radio}^2 + \sigma_\mathrm{astr}^2$.
We fixed the astrometric uncertainty $\sigma_\mathrm{astr} = 0.2''$.
The optical position uncertainties $\sigma_\mathrm{opt}$ are taken from the optical catalogue (generally ${\sim}0.1''$), the radio position uncertainty $\sigma_\mathrm{radio}$ is fixed to $3''$, and the uncertainty in the centroid position $\sigma_\mathrm{c}$ is empirically estimated at $0.2$ times the length of the considered radio source.
For the $30$ optical sources closest to the radio ridgeline, Barkus et al., 2022 calculated the $LR$ and considered the source with the highest $LR$ to be the most likely host galaxy.

We used the method of Barkus et al., 2022 but made three minor adaptations.
First, we introduce explicit regularisation for $q(m,c)$ and $n(m,c)$.
As the PDF estimates for $q(m,c)$ and $n(m,c)$ are 2D KDEs over sampled $(m,c)$-distributions, the parts of the $(m,c)$–parameter space that are sparsely sampled can lead to probabilities that are effectively zero when the realistic theoretical probability should be small but non-zero.
Through the $q/n$-fraction in Eq. <ref>, the resulting values of $LR$ in the sparsely sampled parts of the $(m,c)$–parameter space blow up to unrealistic large values or collapse to almost 0 (see Fig. <ref>). In practice, these unsampled parts of parameter space are almost never visited by new sources for which we calculate $LR$.
Even so, we add a constant factor to the KDE estimate of $q$ and $n$ to get more robust $LR$ values (see Fig. <ref>) and to express the model uncertainties in our functions of $q$ and $n$.
Using 10-fold cross-validation, we empirically select the bandwidths for the KDEs leading to $q$ and $n$ to be $0.4$.
Second, we propose an alternate form of $f(r)$.
For giants, $f(r)$ is rarely dominated by errors in the position of the optical source or that of the radio source.
As $r_\mathrm{opt,centroid}$ and $r_\mathrm{opt,ridge}$ are slightly correlated, multiplication of $f_\mathrm{ridge}(r_\mathrm{opt,ridge})$ and $f_\mathrm{centroid}(r_\mathrm{opt,centroid})$ under-estimates the chance of low values of $r_\mathrm{opt,centroid}$ or $r_\mathrm{opt,ridge}$.
Therefore, we combine $r_\mathrm{opt,centroid}$ and $r_\mathrm{opt,ridge}$ into a single parameter $r_\mathrm{mean}$ that is the mean of the two distance parameters.
Furthermore, we observe that the empirical distributions of $r_\mathrm{opt,centroid}$, $r_\mathrm{opt,ridge}$ and $r_\mathrm{mean}$ for a sample of radio sources with angular length $>1'$ $A_\mathrm{radio,opt}$ for which optical counterparts were determined via visual inspection do not follow a normal distribution as assumed by Barkus et al., 2022 but rather a lognormal distribution (see Fig. <ref>).
Instead of estimating the values of the different error components (astrometric error, error in optical position, error in radio position) we use the empirical values of the distribution of $f(r)$ for the sources in $A_\mathrm{radio,opt}$; see Appendix <ref> for details.
Third, we replaced the Pan-STARRS1 DR1 catalogue (from which colour $c$ was derived) with the DESI Legacy Imaging Surveys DR9 catalogue, as the latter goes up to an $i$-band magnitude of 24.

We applied the modified ridgeline method to all radio sources in our pipeline catalogue with angular lengths larger than $1'$ and brighter than $10\ \mathrm{mJy}$.
We limit the ridgeline procedure to these sources to save time, as the procedure takes multiple seconds per radio source.

After detecting the host galaxies of our radio sources, we checked for spectroscopic redshifts from SDSS <cit.>, or if not available, for photometric redshifts from DESI <cit.>.
The SDSS catalogue also provides velocity dispersions and a quasar flag.
The DESI catalogue includes a flag (fclean) that indicates whether the optical source used in photometric redshift estimation is free from blending and image artefacts.
The catalogue also includes a column (pstar) that estimates how likely it is that the optical source is a star based on its colours.
In both the ML pipeline and RGZ catalogues, we only retained sources for which fclean $=1$ and pstar $\leq 0.2$.

§.§ Reassessing angular source lengths

Summary of the angular length re-evaluation for radio source ILTJ130738.79+270355.1.
Panels A–D show the initial cutout, the removal of neighbouring sources, the masking of emission outside a convex hull based on the old angular length, and the emission that is left after masking all emission below thrice the local noise $\sigma$.
The red line segments delineate the convex hull of the left-over emission, and the red points indicate the furthest removed points in this convex hull.
The great-circle distance between these points is the $3\sigma$ angular length.

Next, we proceeded to reassess the angular source lengths, for both the radio catalogue created using the ML pipeline and the  column reported by the RGZ catalogue described by Hardcastle et al., 2023.
The angular source lengths in the ML catalogue and the  column in the RGZ catalogue are the full width at half maximum (FWHM) of the combined Gaussian components that make up a source, if the source is only composed of a single radio blob.
If the radio source is composed of multiple radio blobs, the reported size is the distance between the two furthest removed points on a convex hull that encloses the FWHMs of the blobs that make up the radio source.
However, in the literature, the length of a GRG is often reported to be the maximum distance between the signal of a radio source that exceeds three times the image noise $\sigma$.

To get the $3\sigma$ angular lengths, we applied five steps to all sources in both catalogues with a reported angular length $\phi > 1'$.
First, we created a square image cutout with a width and height equal to $1.5$ times the old angular source length.
Second, we mask all neighbouring radio emission.
Third, we mask all emission outside an ellipse with a major axis that is the old source length, and a minor axis that is 1.1 times the old source width or a quarter of the old source length if that value is bigger.
These numbers are a result of the observation that, with respect to the $3\sigma$ angular lengths, the old lengths were almost always significantly overestimated, while the source width tended to be underestimated.
Fourth, we mask all remaining emission that is below three times the local noise.
Fifth, we fitted a convex hull around the remaining emission and determined the distance between the points on this convex hull that were farthest apart.
See Fig. <ref> for an illustrative example.

The entire process from source detection (Sect. <ref>) to source list with optical identifications and updated angular lengths (this subsection) took roughly half an hour to one hour per LoTSS DR2 pointing, depending on the detected number of sources.
Each pointing can be processed independently, which allowed us to spread the processing of all $841$ LoTSS DR2 pointings over 5 nodes of a heterogeneous computer cluster with 80 physical CPU cores in total for three to four days.

Finally, for both the ML pipeline and RGZ catalogues, we calculate the projected proper lengths using the $3\sigma$ angular lengths and the redshift estimates corresponding to each source, and discard all sources that do not meet the $l_\mathrm{p} \geq l_\mathrm{p,GRG}$ criterion.
For the ML pipeline catalogue, we discarded all internally duplicate GRG candidates using a $1'$ cone search.
The RGZ catalogue did not contain any internal duplicates.
That left us with $7,001$ GRG candidates in the ML pipeline catalogue and $7,044$ GRG candidates in the RGZ catalogue.

§.§ Inspecting GRG candidates visually

The following step we took in the creation of our GRG sample was a manual visual inspection of all our GRG candidates.
For the RGZ GRG candidates, as described by Hardcastle et al., 2023, at least five different volunteers already inspected the radio and corresponding optical emission, and in most cases the candidates identified in this way were reinspected by a professional astronomer.
The purpose of our manual visual inspection was therefore to exclude only those sources where either the radio component association or the host identification was obviously incorrect.
For each GRG candidate, a single expert looked at a panel showing the candidate with its neighbouring sources masked and most neighbouring emission masked (akin to panel C in Fig. <ref>) and a panel showing the candidate in its wider context (akin to panel A in Fig. <ref>); additionally, the location of the optical host was indicated.
We sorted the candidates into three categories: candidates that looked reasonable, candidates that clearly missed (or included too many) significant radio components, and candidates that showed a very unlikely host galaxy location.
For the ML pipeline GRG candidates, we initially followed the same procedure as for the RGZ GRG candidates.
To speed up the visual inspection, we skipped the $4,272$ ML pipeline GRG candidates that were verified RGZ giants.
After inspecting the ML pipeline GRG candidates once, we subjected all that were not rejected to a second round of visual inspection.
The second round was aided by inspecting LoTSS DR2 radio contours over a Legacy Survey DR9 $(g,r,z)$ image cube, where sources from the combined optical–infrared catalogue within the field of view were highlighted.

For the RGZ catalogue, we judged $6,550$ ($93\%$) GRG candidates to be without issues, $389$ ($6\%$) to have radio component issues, and $105$ ($1\%$) to have been assigned an unlikely host galaxy.
For the $5,864$ (unique) ML pipeline GRG candidates, we judged $2,722$ ($47\%$) candidates to be without issues, $1,963$ ($33\%$) to have radio component issues, and $1,179$ ($20\%$) to have been assigned an unlikely host galaxy.
Radio component association issues for the ML-identified candidates occur because the association method leverages an object detection neural network with rectangular bounding boxes to capture the radio components [Mostert et al., 2022].
The large extent of these sources causes many unrelated (fore- and background) radio sources to appear in the rectangular bounding box that encompasses the candidate, increasing the likelihood of erroneous component associations.
Future ML radio association methods should consider using instance segmentation instead of object detection with a rectangular bounding box.[For the distinction between object detection and instance segmentation, see Lakshmanan et al., 2021.]
Of the $6,550$ RGZ giants, $5,647$ do not appear in previous literature and are thus new discoveries.
Of the $2,722$ ML pipeline giants that are not RGZ giants, $2,597$ are new discoveries.

Qualitatively, from the visual inspection, we noticed that the verified ML pipeline GRG sample contained more symmetric giants with colinear jets, while the verified RGZ GRG sample contained more giants with complex, bent structures indicative of interaction with the IGM.
The ML pipeline did also detect giants with complex structures, but was often unable to fully separate them from all neighbouring unrelated emission.
An in-depth comparison between the ML pipeline and RGZ GRG samples is beyond the scope of this work.
Figures <ref> and <ref> each show six examples of previously unknown giants found through our ML-based approach.
Through cutouts covering $3' \times 3'$, Fig. <ref> shows angularly compact giants; through cutouts covering $6' \times 6'$, Fig. <ref> shows more angularly extended specimina.

LoTSS DR2 cutouts at central observing frequency $\nu_\mathrm{obs} = 144\ \mathrm{MHz}$ and resolution $\theta_\mathrm{FWHM} = 6''$, centred around the hosts of newly discovered giants.
Each cutout covers a solid angle of $3' \times 3'$.
Contours signify 3, 5, and 10 sigma-clipped standard deviations above the sigma-clipped median.
For scale, we show the stellar Milky Way disk (with a diameter of $50\ \mathrm{kpc}$) generated using the Ringermacher & Mead, 2009 formula, alongside a 3 times inflated version.
Each DESI Legacy Imaging Surveys DR9 $(g,r,z)$ inset shows the central $1' \times 1'$ square region.
As all giants obey $\phi \geq 1.3'$, they must – if not oriented along one of the square's diagonals – necessarily protrude from this region.
Rowwise from left to right, from top to bottom, these giants are ILTJ000212.45+222116.2, ILTJ001115.77+220316.6, ILTJ001350.25+324530.8, ILTJ001831.84+322247.7, ILTJ003025.90+334729.2, and ILTJ003534.45+221937.8.

LoTSS DR2 cutouts at central observing frequency $\nu_\mathrm{obs} = 144\ \mathrm{MHz}$ and resolution $\theta_\mathrm{FWHM} = 6''$, centred around the hosts of newly discovered giants.
Each cutout covers a solid angle of $6' \times 6'$.
Contours signify 3, 5, and 10 sigma-clipped standard deviations above the sigma-clipped median.
For scale, we show the stellar Milky Way disk (with a diameter of $50\ \mathrm{kpc}$) generated using the Ringermacher & Mead, 2009 formula, alongside a 3 times inflated version.
Each DESI Legacy Imaging Surveys DR9 $(g,r,z)$ inset shows the central $1' \times 1'$ square region.
As all giants obey $\phi \geq 1.3'$, they must – if not oriented along one of the square's diagonals – necessarily protrude from this region.
Rowwise from left to right, from top to bottom, these giants are ILTJ002943.72+295700.3, ILTJ003010.58+170948.6, ILTJ003521.87+233625.9, ILTJ003712.91+284436.8, ILTJ004002.30+252550.9, and ILTJ235802.49+331838.5.

With $11,585$ unique giants, we present the largest catalogue of large-scale galactic feedback to the Cosmic Web.
The RGZ (orange) and ML pipeline (red) samples are strictly confined to the LoTSS DR2 area, while the sample by Oei et al., 2023 extends to yet-to-be-released LoTSS pointings processed with the DR2 pipeline.

§.§ Merging RGZ, ML pipeline, and literature samples

To complete our GRG sample, we iteratively added giants from the literature, going from the newest to the oldest publication.
This approach follows from the assumption that newer publications are generally based on more sensitive and higher-resolution observations, leading to more accurate angular length measurements.
In an effort to avoid having duplicate giants in the final sample, we only added giants when their host galaxies were more than $10''$ away from all host galaxies of already aggregated giants.

The joint RGZ–ML pipeline sample contains $9,272$ giants.
We added $1,471$ out of the $2,193$ giants presented by Oei et al., 2023,
$43$ out of the $69$ giants presented by Simonte et al., 2022,
$62$ out of the $62$ giants presented by Gürkan et al., 2022,
$165$ out of the $263$ giants presented by Mahato et al., 2022,
$178$ out of the $178$ giants presented by Andernach et al., 2021,
$0$ out of the $1$ giants presented by Masini et al., 2021,
$2$ out of the $2$ giants presented by Delhaize et al., 2021,
$1$ out of the $2$ giants presented by Bassani et al., 2021,
$1$ out of the $4$ giants presented by Tang et al., 2020,
$390$ out of the $694$ giants presented by Dabhade et al., 2020, and
$0$ out of the $6$ giants presented by Ishwara-Chandra et al., 2020.
These additions result in a final catalogue with $11,585$ unique giants.
This is the first catalogue of giants to contain more than $10^4$ specimina.

Figure <ref> shows a Mollweide view of the sky with the locations of both the newly confirmed giants and the giants from the literature.
Almost all known giants stay clear of the Galactic Plane, where radio emission from the Milky Way – of which we show the specific intensity function at $\nu_\mathrm{obs} = 150\ \mathrm{MHz}$ in greyscale [Zheng et al., 2017] – makes calibration and imaging harder.
In addition, optical host identification is much harder near the Galactic Plane.
The default field of view set-up of both our ML pipeline (Sect. <ref>) and of RGZ favours the discovery of giants with angular lengths of a few arcminutes at most.
By contrast, the GRG search campaign of Oei et al., 2023 featured a `fuzzy' ${\sim}5'$ lower threshold to allow for an exhaustive manual search with an interactive and dynamic field of view <cit.>.
Figure <ref> demonstrates that these design choices lead to GRG samples with markedly different angular length distributions.

Observed distributions of angular length $\phi$, showing that our three LoTSS DR2 search methods target different ranges of $\phi$.
The largest angular lengths detected by Oei et al., 2023, RGZ, and the ML pipeline are $132'$, $43'$, and $8'$ respectively, but we limit the horizontal axis to $12'$ for interpretability.
The vertical line marks the minimum angular length that giants can attain: $\phi_\mathrm{GRG}(l_\mathrm{p,GRG} = 0.7\ \mathrm{Mpc}) = 1.3'$.

Our sample of RGZ giants (orange squares) and ML pipeline giants (red squares) effectively complements the sample of giants with large angular lengths (blue dots) from the manual search of Oei et al., 2023.
The remaining giants (green pluses) are from earlier literature, as specified in Sect. <ref>.

[]Statistics of the GRG samples that we discovered, confirmed, or aggregated.
From left to right, the columns provide the number of giants in each sample, $N$, and the $10$th, the median,  and the $90$th percentile of the angular length $\phi$, redshift $z$, and projected proper length $l_{\mathrm{p}}$.

Sample     $N$      $\phi_\mathrm{10th}$ $(')$     $\phi_\mathrm{median}$ $(')$     $\phi_\mathrm{90th}$ $(')$      $z_\mathrm{10th}$     $z_\mathrm{median}$     $z_\mathrm{90th}$      $l_\mathrm{p,10th}$ $(\mathrm{Mpc})$     $l_\mathrm{p,median}$ $(\mathrm{Mpc})$     $l_\mathrm{p,90th}$ $(\mathrm{Mpc})$

ML pipeline     $2,722$      $1.50$     $1.97$     $3.09$     $0.44$     0.87     $1.28$      $0.72$     $0.87$     $1.29$

RGZ     $6,550$      $1.56$     $2.19$     $4.40$     $0.31$     0.75     $1.19$      $0.73$     $0.91$     $1.57$

Known giants     $11,585$      $1.57$     $2.33$     $5.70$     $0.23$     0.72     $1.19$      $0.73$     $0.94$     $1.68$

As a result, the samples complement each other: the sample of Oei et al., 2023 is more complete at lower redshifts and higher projected lengths, while the RGZ and ML pipeline samples are more complete at higher redshifts and lower projected lengths.
Figure <ref> demonstrates this point, while Table <ref> presents the corresponding statistics of the GRG samples.

For comparison of the $3\sigma$ lengths of the ML pipeline and RGZ giants to those in other surveys, we inform the reader that the central frequency and the average surface brightness threshold of the observations that we use are $\nu_\mathrm{obs} = 144\ \mathrm{MHz}$ and  $b_{\nu,\mathrm{th}} = 25\ \mathrm{Jy\ deg}^{-2}$ respectively.

§.§ Estimating Bayesian model parameters

After having refined our statistical GRG framework (Sect. <ref>), and after having assembled the largest sample of giants yet (Sects. <ref>–<ref>), we combined both advances to perform inference of the length distribution, number density, and lobe volume-filling fraction of giants.

Given that our goal has been to infer properties of the full population of giants, rather than just of those currently observed, we included two main selection effects in our forward modelling.
As detailed in Sect. <ref>, we parametrised surface brightness selection with three parameters, which are free parameters of the model.
As detailed in Sect. <ref>, a second cause of selection is the imperfect operation of our three LoTSS DR2 search methods, all of which fail to identify a significant fraction of giants with lobe surface brightnesses above the survey noise level.
We modelled this identification selection $p_{\mathrm{obs,ID}}$ with a set of logistic functions, regressed to GRG data.
We now provide details of this process.

§.§.§ Identification probability functions

To estimate $p_\mathrm{obs,ID}(l_\mathrm{p},z)$ from data, we first selected all giants detected by the joint efforts of our machine learning pipeline, RGZ, and the manual, visual search of Oei et al., 2023.
Next, we retained only those giants that are located in regions of the sky that have been scanned by all three searches.
This overlap region in principle corresponds to the full LoTSS DR2 coverage – were it not for the fact that the search of Oei et al., 2023 skipped over the LoTSS DR1, which had already been scanned by Dabhade et al., 2020.
Therefore, the actual overlap region amounts to the LoTSS DR2 coverage with a spherical quadrangle removed, whose minimum and maximum right ascensions are $\alpha_\mathrm{min} = 160\degree$ and $\alpha_\mathrm{max} = 230\degree$ and whose minimum and maximum declinations are $\delta_\mathrm{min} = 45\degree$ and $\delta_\mathrm{max} = 56\degree$.
Appendix <ref> provides an explicit decomposition of our assumed LoTSS DR2 coverage – and therefore implicitly of the overlap region – in terms of disjoint spherical quadrangles.

Some of the retained giants have been detected only in the combined RGZ–ML search, others have been detected only in the Oei et al., 2023 search, and yet others have been detected in both.
We note that, had it operated flawlessly, the combined RGZ–ML search would have detected all sources claimed by Oei et al., 2023 (or at least those that are genuine giants – which should be the vast majority).
Therefore, by mapping the (in)ability of the RGZ–ML search to detect the giants of Oei et al., 2023 as a function of $l_\mathrm{p}$ and $z$, one can estimate the RGZ–ML search's identification probability function, $p_\mathrm{obs,ID,1}(l_\mathrm{p},z)$.
More precisely, for each giant detected by Oei et al., 2023, we evaluated whether it was also detected in the RGZ–ML search, and stored a corresponding Boolean (that is to say, either $1$ or $0$).
We show these Booleans, at the $(l_\mathrm{p},z)$ coordinates of the giants they belong to, as yellow (representing $1$) and blue (representing $0$) dots in the top-left panel of Fig. <ref>.
Viewing the Boolean at $(l_\mathrm{p},z)$ as a realisation of a Bernoulli RV with success probability $p = p_\mathrm{obs,ID,1}(l_\mathrm{p},z)$, we recognise the inference of the identification probability function as a binary logistic regression problem with two explanatory variables.
The background of Fig. <ref>'s top-left panel shows the corresponding best fit.

By symmetry, this approach can be reversed to estimate the Oei et al., 2023 search's identification probability function, $p_\mathrm{obs,ID,2}(l_\mathrm{p},z)$.
Therefore, for each giant detected in the RGZ–ML search, we evaluated whether it was also detected by Oei et al., 2023, and stored a corresponding Boolean.
In the same way as before, we show these Booleans in the middle-left panel of Fig. <ref>.
The panel's background shows the best logistic fit.

We combine the two identification probability functions, $p_{\mathrm{obs,ID,1}}(l_\mathrm{p},z)$ and $p_{\mathrm{obs,ID,2}}(l_\mathrm{p},z)$, in point-wise fashion as to obtain a single function $p_{\mathrm{obs,ID}}(l_\mathrm{p},z)$.
To do so, we follow the minimal combination rule of Eq. <ref>.

We remark that, by giving each Boolean in these logistic regressions an equal weight, the resulting functions are tuned to fit crowded regions of projected length–redshift parameter space best – at the expense of accuracy in sparser regions.
To increase the accuracy of the functions for the parameter space at large, we performed a simple rebalancing step.
First, we calculated the mean number density in the parameter space given by $l_\mathrm{p} \in [0.7, 5\ \mathrm{Mpc}] \times [0, 0.5] \ni z$.
We then selectively subsampled the data in crowded regions, following the rule that the number density in each bin of width $0.5\ \mathrm{Mpc}$ and height $0.05$ should not exceed twice the mean number density of the entire parameter space.
We show the rebalanced data, alongside refitted logistic models, in the right column of Fig. <ref>.
We report the rebalanced model coefficients in Table <ref>, and treat them as constants during the Bayesian inference.

Overview of our determination of the probability to identify giants in the LoTSS DR2 with above-noise surface brightnesses, as a function of projected length and redshift – through Radio Galaxy Zoo: LOFAR and our machine learning pipeline (top row), through the search of Oei et al., 2023 (middle row), and through these methods in unison (bottom row).
Each of the upper four panels shows a binary logistic regression following the theory of Sect. <ref> and the practical considerations of Sect. <ref>.
The left column shows results from all available data, whilst the right column shows results from rebalanced data.
In our Bayesian inference, we used the latter results.

§.§.§ Inference in practice

In this work, we constrained the parameters of Sect. <ref>'s GRG population model via a projected length–redshift histogram.
From our most extensive sample of giants, we selected those with $0.7\ \mathrm{Mpc} \eqqcolon l_\mathrm{p,GRG} < l_\mathrm{p} < 5.1\ \mathrm{Mpc}$ and $0 < z < z_\mathrm{max} \coloneqq 0.5$ that lie in the LoTSS DR2 coverage as specified in Appendix <ref>.
We did not include the giants from Oei et al., 2023 for which only a lower bound to the redshift is known.
This selection retained $2,685$ out of $11,585$ giants.
We used these giants to fill a histogram with bins of width $\Delta l_\mathrm{p} = 0.1\ \mathrm{Mpc}$ and $\Delta z = 0.02$.
We did not systematically explore the effect of these bin size parameters on the resulting inference.
However, the smaller one chooses the bins, the higher the numerical cost will be.
On the other hand, if the bins are chosen much larger than the typical scales over which the underlying observed projected length–redshift distribution[With the `underlying' observed projected length–redshift distribution, we mean the observed projected length–redshift distribution one would obtain in the limit of an infinite number of observed giants.] varies, then some ability to extract parameter constraints will be lost.

To compute the posterior distribution over the six parameters $\vec{\theta} = [\xi(l_\mathrm{p,GRG}), \Delta\xi, b_\mathrm{\nu,ref}, \sigma_\mathrm{ref}, \zeta, n_\mathrm{GRG}]$, we assumed a uniform prior and brute-force evaluated the likelihood function over a regular grid that covers a total of $2.1 \cdot 10^9$ parameter combinations.[This approach is feasible by virtue of the low numerical cost of each likelihood function evaluation.
Its main advantage is its simplicity: there are no parameters to tune that govern the method's convergence behaviour.
Once the model is expanded to include more parameters, or when selection effects with higher numerical cost are incorporated, more efficient (though more complicated) methods such as Markov chain Monte Carlo or nested sampling will become necessary.]
In doing so, we applied the Poissonian likelihood trick described in Appendix <ref>, which sped up our computations by one to two orders of magnitude.
Table <ref> provides the parameter ranges for which we evaluated the likelihood (which coincide with their prior distribution ranges), alongside all of the model's constants and their assumed values.
Because each likelihood function evaluation can be computed independently of the others, the problem is fully parallelisable.
In practice, we distributed the ${\sim}10^4$ core-hours Python calculation over ${\sim}1500$ virtual cores, which were spread across ${\sim}20$ nodes of a computer cluster.
Next, we generated samples from the posterior distribution using rejection sampling <cit.>.
We subsequently used these samples to calculate probability distributions for derived quantities.[To calculate probability distributions over quantities that are a function of the parameters, such as the Local Universe GRG lobe VFF, $\mathcal{V}_\mathrm{GRG-CW}(z = 0)$, or the joint search completeness function $C$, we could in principle evaluate these quantities for each parameter combination of the aforementioned grid and weigh each grid point's result by the associated likelihood (or, equivalently, posterior probability).
However, some derived quantities are costly to compute, so that excessive evaluations should be avoided.]

Parameters and constants of GRG population forward model from Sect. <ref> alongside their prior ranges and values, as used in the Bayesian inference presented in Sect. <ref>.
The first six constants serve to define the quantitative meaning of the parameters and set the scope of the analysis.
The other eleven constants are not arbitrary: they affect the likelihood function and posterior distribution for a given set of parameter definitions and scope.

Parameter     Uniform prior range     Explanation
$\xi(l_\mathrm{p,1} = l_\mathrm{p,GRG})$          $[-3.5, -2]$     Sect. <ref>

$\Delta\xi$     $[-3.5, -1.5]$     Sect. <ref>

$b_{\nu,\mathrm{ref}}$     $[1, 100] \cdot \mathrm{Jy\ deg^{-2}}$     Sect. <ref>

$\sigma_\mathrm{ref}$     $[0.5, 2]$     Sect. <ref>

$\zeta$     $[-0.5, 0]$     Sect. <ref>

$n_\mathrm{GRG}$     $[0, 50] \cdot (100\ \mathrm{Mpc})^{-3}$     Sect. <ref>

Constant     Value     Explanation
$l_\mathrm{p,GRG}$          $0.7\ \mathrm{Mpc}$     Sect. <ref>

$l_\mathrm{p,1}$     $0.7\ \mathrm{Mpc}$     Sect. <ref>

$l_\mathrm{p,2}$     $5\ \mathrm{Mpc}$     Sect. <ref>

$l_\mathrm{ref}$     $0.7\ \mathrm{Mpc}$     Sect. <ref>

$\nu_\mathrm{obs}$     $144\ \mathrm{MHz}$     Sect. <ref>

$z_\mathrm{max}$     $0.5$     Sect. <ref>

$\alpha$     $-1$     Sect. <ref>

$b_{\nu,\mathrm{th}}$     $25\ \mathrm{Jy\ deg^{-2}}$     Sect. <ref>

$\beta_{0,1}$     $-1.0$     Sect. <ref>

$\beta_{0,2}$     $-1.0$     Sect. <ref>

$\beta_{l_\mathrm{p},1}$     $-0.1\ \mathrm{Mpc}^{-1}$     Sect. <ref>

$\beta_{l_\mathrm{p},2}$     $2.4\ \mathrm{Mpc}^{-1}$     Sect. <ref>

$\beta_{z,1}$     $2.8$     Sect. <ref>

$\beta_{z,2}$     $-6.4$     Sect. <ref>

$\Delta l_\mathrm{p}$     $0.1\ \mathrm{Mpc}$     Sect. <ref>

$\Delta z$     $0.02$     Sect. <ref>

$\Omega$     $1.62\ \mathrm{sr}$     Appendix <ref>

§ RESULTS

By combining an unparalleled sample of giant radio galaxies with a rigorous forward model, we have produced a posterior distribution over parameters that characterise the intrinsic population of giants.
Figure <ref> summarises the posterior over parameter hexads $\vec{\theta} = [\xi(l_\mathrm{p,GRG}), \Delta\xi, b_{\nu,\mathrm{ref}}, \sigma_\mathrm{ref}, \zeta, n_\mathrm{GRG}]$ by means of its one- and two-dimensional marginal distributions.

Likelihood function over $\vec{\theta} = [\xi(l_\mathrm{p,GRG}), \Delta\xi, b_{\nu,\mathrm{ref}}, \sigma_\mathrm{ref}, \zeta, n_\mathrm{GRG}]$, based on 2,685 projected lengths and redshifts of giants up to $z_\mathrm{max} = 0.5$.
We show all two-parameter marginals of the likelihood function, with contours enclosing 50% and 90% of total probability.
We mark the maximum likelihood estimate (MLE) values (grey dot) and the likelihood mean values (grey cross).
The one-parameter marginals again show the MLE (dash-dotted line), a mean-centred interval of standard deviation–sized half-width (hashed region), and a median-centred 90% credible interval (shaded region).

In this section, we analyse our newfound parameter constraints.

§.§ Length distribution of giant radio galaxies

Radio galaxies enrich the IGM with magnetic fields, but giants – given their megaparsec-scale reach – appear uniquely capable of seeding the more remote regions of the Cosmic Web.
Consequently, scientific interest in quantifying the length distribution of giants has arisen from the possibility that giants contribute significantly to cosmic magnetogenesis.
The question at hand is deceivingly simple: how common are giants of various lengths?

As pointed out by Oei et al., 2023, due to selection effects, the observed projected length distribution is not a reliable estimate of the true projected length distribution.
Worse still, the relevant selection effects might not be quantitatively known a priori, requiring joint inference of the length distribution, and the selection effect parameters.
Oei et al., 2023 performed such a joint inference, and found that their data were consistent with an underlying population of giants with Pareto-distributed lengths, characterised by tail index $\xi = -3.4 \pm 0.5$.
In the current work, we have relaxed the assumption of perfect Paretianity, and explore whether the data are consistent with a curved power law PDF for the GRG projected proper length RV $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$.
The marginals of Fig. <ref> suggest that they are – in fact, the data strongly favour curvature, with a tail index at $l_\mathrm{p,1} \coloneqq l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$ of $\xi(l_\mathrm{p,GRG}) = -2.8 \pm 0.2$ and a total increase in tail index up to $l_\mathrm{p,2} \coloneqq 5\ \mathrm{Mpc}$ of $\Delta\xi = -2.4 \pm 0.3$.
Given the small relative uncertainty on the latter value, our data appear inconsistent with perfect Paretianity ($\Delta \xi = 0$).
We note that our notion of `data' is different from that in Oei et al., 2023: not only do we use more than a thousand additional giants, we also make more effective use of their redshift information.
For further discussion, see Sect. <ref>.

It remains an open question whether giants can be understood as part of the ordinary radio galaxy population, or whether they evolve through qualitatively different physical processes.
As pointed out in Sect. 4.1.5 of Oei et al., 2023, a curved power law PDF for $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ is consistent with a scenario in which giants share a broader length continuum with smaller radio galaxies.
More specifically, if the broader radio galaxy length distribution is approximately lognormal, as appears justifiable on statistical grounds, then $\xi$ should decrease throughout the distribution's right tail – that is, throughout the GRG range.
Future research should determine whether such a unified non-giant RG–GRG scenario is also quantitatively consistent with the decrease in $\xi$ we have inferred here.
In addition, our inferences of $\xi(l_\mathrm{p,GRG})$ and $\Delta\xi$ are important in constraining Sect. <ref>'s GRG lobe volume-filling fraction.

§.§ Number density of giant radio galaxies

The extent to which giants have contributed to cosmic magnetogenesis depends on their intrinsic number density – which need not necessarily be a constant, but could have evolved over time.
Observationally, giants are considered rare in comparison to smaller radio galaxies.
However, because giants are presumably strongly affected by surface brightness selection, this present-day observed rarity might not translate to an intrinsic rarity.
Excitingly, by forward modelling selection effects – and in particular surface brightness selection – we can constrain the intrinsic comoving GRG number density between $z = 0$ and $z = z_\mathrm{max}$, which we denote simply by $n_\mathrm{GRG}$.

The bottom-right one-dimensional marginal of Fig. <ref> shows a strongly skewed distribution for $n_\mathrm{GRG}$, with a marginal mean $\mathbb{E}[n_\mathrm{GRG}] = 13 \pm 10\ (100\ \mathrm{Mpc})^{-3}$ and a $95\%$ probability that $n_\mathrm{GRG} > 4\ (100\ \mathrm{Mpc})^{-3}$.
These number densities are a factor of order unity higher than those of Oei et al., 2023, who inferred a marginal mean $\mathbb{E}[n_\mathrm{GRG}] = 4.6 \pm 2.4\ (100\ \mathrm{Mpc})^{-3}$ and a $90\%$ probability that $n_\mathrm{GRG} < 6.7\ (100\ \mathrm{Mpc})^{-3}$.

The joint marginal distribution of $n_\mathrm{GRG}$ and $b_{\nu,\mathrm{ref}}$ reveals a strong inverse relationship, whose origin is easy to grasp.
Models in which giants are relatively rare (i.e. with low $n_\mathrm{GRG}$) but with relatively mild surface brightness selection (i.e. with high $b_{\nu,\mathrm{ref}}$) are about as successful in reproducing the data-derived projected length–redshift histogram as models in which giants are relatively common (i.e. with high $n_\mathrm{GRG}$) but with relatively severe surface brightness selection (i.e. with low $b_{\nu,\mathrm{ref}}$).
The narrowness of the joint distribution also suggests that, if estimates of $b_{\nu,\mathrm{ref}}$ would reveal it to be $\gtrsim 10\ \mathrm{Jy\ deg^{-2}}$, it should be possible to break the degeneracy and accurately determine $n_\mathrm{GRG}$.

Recent work [Oei et al., 2023] suggests that the comoving number density of luminous, non-giant radio galaxies (LNGRGs), understood to have radio luminosities at $150\ \mathrm{MHz}$ of $l_\nu \geq 10^{24}\ \mathrm{W\ Hz^{-1}}$ and projected lengths $l_\mathrm{p} < l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$, is $n_\mathrm{LNGRG} = 12 \pm 1\ (100\ \mathrm{Mpc})^{-3}$.
Our work suggests that giants might be comparably common.
If this is indeed the case, then the widespread belief that giants form a rare population of radio galaxies must be revised.

§.§ Lobe volume-filling fraction of giant radio galaxies

The present-day volume-filling fraction of the lobes of giants in clusters and filaments of the Cosmic Web, $\mathcal{V}_\mathrm{GRG-CW}(z=0)$, is not a parameter of our model, but rather a derived quantity.
As briefly discussed in Sect. <ref>, we compute its probability distribution using the parameter hexads that we have obtained by rejection sampling from the posterior.
For each sampled hexad, we compute $\xi(l_\mathrm{p})$ using $\xi(l_\mathrm{p,GRG})$, $\Delta\xi$, and Eq. <ref>, then $f_{L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}}(l_\mathrm{p})$ using Eq. <ref>, and finally $\mathcal{V}_\mathrm{GRG-CW}(z=0)$ using $n_\mathrm{GRG}$ and Eq. <ref>.

Empirical relation between projected proper length and two-lobe proper volume for RGs from Ineson et al., 2017, Oei et al., 2022, and Oei et al., 2023.
The power law–like trend motivates Eq. <ref>.
The green fit is based on giants only, whilst the grey fit is based on all RGs.
$V_\mathrm{GRG}$ is the mean two-lobe proper volume of the shortest possible giants (i.e. giants for which $l_\mathrm{p} = l_\mathrm{p,GRG}$), while $\gamma$ is the exponent of the power law. For self-similar RG growth, $\gamma = 3$.

To arrive at Eq. <ref>, we assumed in Eq. <ref> that RG projected proper lengths and two-lobe proper volumes obey a power law relation with scatter.
To investigate the validity of this assumption, we collected the projected proper lengths and two-lobe proper volumes of all Fanaroff–Riley II RGs in the `representative' sample of Ineson et al., 2017.[We calculated projected proper lengths by summing the lobe tip distances of both lobes. For RGs for which Ineson] et al., 2017's Table 5 reports data on only one lobe, we assumed that the other lobe has an identical lobe tip distance and volume.
Among these RGs are just seven giants.
To more reliably probe the projected proper length–two-lobe proper volume relation for giants, we supplemented this sample with the giant of NGC 6185, the longest spiral galaxy–generated RG known [Oei et al., 2023], and with Alcyoneus, the longest elliptical galaxy–generated RG known [Oei et al., 2022].
Oei et al., 2023 estimated that the giant of NGC 6185 measures $l_\mathrm{p} = 2.45 \pm 0.01\ \mathrm{Mpc}$ and has a two-lobe proper volume $V = 0.35 \pm 0.03\ \mathrm{Mpc}^3$.
Similarly, Oei et al., 2022 estimated that Alcyoneus measures $l_\mathrm{p} = 4.99 \pm 0.04\ \mathrm{Mpc}$ and has a two-lobe proper volume $V = 2.5 \pm 0.3\ \mathrm{Mpc}^3$.
We show the empirical $l_\mathrm{p}$–$V$ relation for the resulting sample in Fig. <ref>.
Tentatively, we consider the assumed power law relation between projected proper length and two-lobe proper volume justified, although we warn that the sample size is small and that we have not made corrections for selection effects.
Through least-squares minimisation in log–log space, we obtained two best-fit power law relations: one for giants only (green line) and for all RGs (grey line).
In order to calculate Eq. <ref>, we adopted the parameters from the giant-based fit: $V_\mathrm{GRG} = 1.2 \cdot 10^{-2}\ \mathrm{Mpc}^3$ and $\gamma = 2.7$.

For each sampled hexad, we calculate $\mathbb{E}[L_\mathrm{p}^{\gamma}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}]$ using the law of the unconscious statistician.
The mean two-lobe volume of a giant is $\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}](z=0) = 5.1 \pm 0.3\ \cdot 10^{-2}\ \mathrm{Mpc}^3$.
As in Oei et al., 2023, we assume that clusters and filaments comprise $5\%$ of the Local Universe's volume [Forero-Romero et al., 2009]: $\mathcal{V}_\mathrm{CW}(z=0) = 5\%$.
We show the resulting posterior distribution over $\mathcal{V}_\mathrm{GRG-CW}(z = 0)$ in Fig. <ref>.

Posterior distribution for the instantaneous, present-day GRG lobe volume-filling fraction in clusters and filaments of the Cosmic Web, $\mathcal{V}_\mathrm{GRG-CW}(z = 0)$.

This probability distribution inherits its skewness from the skewed marginal of $n_\mathrm{GRG}$.
We find $\mathcal{V}_\mathrm{GRG-CW}(z = 0) = 9\substack{+15\\-4} \cdot 10^{-6}$ and a posterior mean and standard deviation of $\mathcal{V}_\mathrm{GRG-CW}(z = 0) = 1.4 \pm 1.1 \cdot 10^{-5}$.
These results appear statistically consistent with that of Oei et al., 2023: $\mathcal{V}_\mathrm{GRG-CW}(z=0) = 5\substack{+8\\-2}\cdot 10^{-6}$.

While they appear low at first sight, we speculate that these numbers are consistent with a scenario in which giants contribute significantly to cosmic magnetogenesis.
To see why, we first note that the number of giants that have ever existed might exceed those that exist now by two orders of magnitude: giants are actively powered for ${\sim}10^1$–$10^3\ \mathrm{Myr}$ <cit.>, the Universe is ${\sim}10^4\ \mathrm{Myr}$ old, and one may assume – in the absence of clear evidence to the contrary – that the comoving number density of giants $n_\mathrm{GRG}$ has remained roughly constant over cosmic time.[The dynamics of RGs are different at early epochs; reasons include a less developed large-scale structure but higher mean cosmic densities, higher AGN cold gas accretion and galaxy merger rates <cit.>, and more severe inverse Compton energy losses to the CMB <cit.>. Giants possibly representative of early epochs are discussed in e.g. Mahato] et al., 2022.
Second, a likely large fraction of all giants that existed throughout cosmic history lived when the Universe was significantly smaller.
More specifically, all giants at $z \geq z_{10} \coloneqq \sqrt[3]{10}-1 \approx 1.15$ lived when the Universe's volume was at least an order of magnitude smaller; this period covers ${\sim}38\%$ of all cosmic time.
If the mean two-lobe proper volume of giants $\mathbb{E}[V\ \vert\ L_\mathrm{p}
\geq l_\mathrm{p,GRG}]$ and the VFF of clusters and filaments $\mathcal{V}_\mathrm{CW}$ have remained roughly constant over the redshift range $[0, z_{10}]$, then by Eq. <ref> we have $\mathcal{V}_\mathrm{GRG-CW}(z = z_{10}) = 10 \cdot \mathcal{V}_\mathrm{GRG-CW}(z = 0)$.[Based on Millennium simulations, which are dark matter–only, Cautun] et al., 2014 suggest that $\mathcal{V}_\mathrm{CW}$ may instead have decreased from ${\sim}10\%$ at $z = z_{10}$ to ${\sim}5\%$ at $z = 0$. This would imply that $\mathcal{V}_\mathrm{GRG-CW}(z = z_{10})$ is a factor ${\sim}2$ smaller than we claim here.
Third, buoyant lobes might deposit magnetic fields in their wake, while diffusion, turbulence, and merger and accretion shocks from large-scale structure formation might have spread the contents of GRG lobes further through the IGM <cit.>.
The typical extent of a GRG lobe along a single spatial dimension, $\ell$, can be considered to be

\begin{align}
\ell(z) \coloneqq \sqrt[3]{\frac{\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}](z)}{2}},
\end{align}

meaning that $\ell(z = 0) = 0.295 \pm 0.006\ \mathrm{Mpc} \sim 10^{-1}\ \mathrm{Mpc}$.
If, after jet fuelling stops, lobes rise buoyantly to the edges of filaments, their total traversed path length will be ${\sim}10^0\ \mathrm{Mpc}$ <cit.>.
The columns through which GRG lobes have risen might therefore have a volume that is an order of magnitude larger than the lobes' own.
Taken together, these three effects could render the present-day VFF of magnetic fields that were once contained in the lobes of giants higher than $\mathcal{V}_\mathrm{GRG-CW}(z=0)$ by four (i.e. $2 + 1 + 1$) or more orders of magnitude.
This, in turn, suggests a significant astrophysical seeding potential.
For instance, assuming four orders of magnitude, ${\sim}10\%$ of the volume of today's Cosmic Web should have been magnetised by giants.

We finally point out that giant-induced IGM magnetic fields could have strengths consistent with observational constraints.
At the moment, the lowest magnetic field strengths measured in giant radio galaxy lobes, as inferred from images of Alcyoneus and the giant generated by NGC 6185 assuming the equipartition or minimum energy condition, are $400$–$500\ \mathrm{nG}$ [Oei et al., 2022, Oei et al., 2023].
If such field strengths would be typical, and buoyancy and diffusion lowers the density of field lines by an order of magnitude, then the typical giant-induced IGM field strength would be ${\sim}10\ \mathrm{nG}$.
This is in agreement with recent radio estimates and limits <cit.>.
We note that this argument ignores possibly significant amplification and decay mechanisms, such as turbulent amplification and decay.

§ DISCUSSION

Below, we discuss how our ML pipeline and GRG population inference compare to earlier work.

§.§ Comparison with previous machine learning search techniques

Proctor, 2016 applied an ML approach to search for GRG candidates by looking for likely pairs of (unresolved) radio lobes with the required angular length in the NRAO VLA Sky Survey <cit.>.
For this radio source component association problem, Proctor, 2016 trained an oblique classifier <cit.>, using six source finder–derived features on $51,195$ pairs of radio components, $48$ of which were verified giants.
This method proved to be useful under the assumption that giants generally appear as an isolated pair of unresolved radio blobs, which is the case for NVSS with its $45''$ resolution and $450$ $\mathrm{\mu Jy}$ $\mathrm{beam^{-1}}$ sensitivity.
Dabhade et al., 2020 visually inspected the $1,600$ GRG candidates presented by Proctor, 2016 and confirmed $151$ giants, which implies a $9\%$ precision for the GRG candidate predictions.
However, Proctor, 2016 expect that giants with resolved lobes – which rule-based source finders often incorrectly break down into multiple separate sources – require a different approach, and virtually all GRG lobes in LoTSS are resolved.[
As all giants have angular lengths $\phi \geq 1.3'$, they cover at least $13$ LoTSS $6''$ beams. This suggests that a single GRG lobe will cover multiple beams, too.]
It works in our favour that the convolutional neural network in our ML pipeline (Sect. <ref>) was specifically designed to use the morphology of the resolved, extended emission as a cue for the radio source component association.
Furthermore, as the source suggestions from our ML pipeline include optical host identifications, the candidates that we inspected not only have the required angular length but also have a host galaxy and corresponding redshift estimate assigned.
This allows us to visually inspect only those radio sources that fulfil the projected proper length $l_\mathrm{p,GRG} := 0.7$ Mpc requirement.
Overall, our ML pipeline has a precision of $47\%$ for the GRG candidates that it suggests.

§.§ Comparison with previous inference strategies

Compared to the approach of Oei et al., 2023, our approach makes better use of the redshift information available for each giant.
More specifically, we use the redshifts to make a `redshift-resolved' observed projected length histogram, while Oei et al., 2023 only compared a `redshift-collapsed' distribution of observed projected lengths to forward model predictions of $L_\mathrm{p,obs}\ \vert\ L_\mathrm{p,obs} \geq l_\mathrm{p,GRG}$.
Effectively, Oei et al., 2023 therefore used for each giant only Boolean redshift information, $\mathbb{I}(z_i < z_\mathrm{max})$: that is, a truth value indicating whether or not the giant with index $i$ resides at a redshift below the maximum considered value.

In addition, our work changed the comoving number density of giants, $n_\mathrm{GRG}$, from a derived quantity to a model parameter.
This approach acknowledges the fact that the observed number of giants, either for a specific projected length–redshift bin or for the parameter space in its entirety, scales linearly with $n_\mathrm{GRG}$ (if the selection effects remain the same).
Therefore, there is intrinsic population information contained in the observed number of giants.
However, by comparing predicted and observed probability distributions only, Oei et al., 2023 did not exploit this fact.

§.§ Future work

With the advent of large-scale, sensitive, low-frequency sky surveys such as the LoTSS, the Evolutionary Map of the Universe survey <cit.>, and the arrival of next-generation instruments such as the SKA [Dewdney et al., 2009] and the DSA-2000 <cit.> later this decade, opportunities will arise to detect many more giants than have been found hitherto.
It is therefore likely that automated approaches to giant finding and host association will become only more relevant in the future.

Regarding our own machine learning–based pipeline, there is significant room to improve both the radio component association and the host association. Visual inspection indicated a precision of $47\%$ and the empirically determined $p_{\mathrm{obs,ID}}$ in Fig. <ref> showed that even in combination with the RGZ sample, the ML pipeline recall does not surpass $70\%$.
Sensible paths to improve the radio component association within the ML pipeline architecture include switching from rectangular bounding box–based object detection (the Fast R-CNN used in this article) to pixel-based instance segmentation and using a larger convolutional backbone <cit.> or a transformer-based backbone <cit.>.
Mostert et al., 2022 conclude that a larger convolutional neural network is not effective unless one also significantly increases the quantity of high-quality training data, and in general, transformers require even more training data than convolutional neural networks <cit.>.
To that extent, adding a filtered version[For example, by identifying a handful of very active and expert volunteers and increasing the weight of their votes.] of the available LoTSS DR2 RGZ annotations [Hardcastle et al., 2023] to the training data can be considered.
Furthermore, assembling a joined data set encompassing the (labelled) survey data of other low frequency radio telescopes can be considered.
Pretraining on this data set can benefit radio galaxy component association, host identification and morphological classification tasks across the board.

Finally, there appear to be clear opportunities to make the population-based forward model presented in Sect. <ref> more accurate.
For example, at present, we have neglected photometric redshift uncertainties; however, the consequences of these uncertainties appear perfectly possible to forward model.
One such currently ignored consequence is Eddington bias: as RGs with projected lengths $l_\mathrm{p} = 0.6\ \mathrm{Mpc}$ are intrinsically more common than those with projected lengths $l_\mathrm{p} = 0.8\ \mathrm{Mpc}$, redshift error–induced projected length errors have the net effect of falsely raising the number of supposed giants with projected lengths near $l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$.
This effect could contaminate the inference of $\xi(l_\mathrm{p,GRG})$.
Somewhat more challenging, but plausibly of greater value, would be a further exploration of how surface brightness selection is effectively modelled.
A major focus of such an exploration would be to analyse the surface brightness characteristics of hitherto discovered giants.
As the masked cutouts of Fig. <ref> suggest, the machine learning–based pipeline described in this work offers the exciting potential to amass – fully automatically – surface brightness properties for thousands of giants.
The availability of such properties for a large fraction of observed giants also allows one to fit the forward model to an observed projected length–redshift–surface brightness histogram, rather than to an observed projected length–redshift histogram only.
It is highly likely that adding another dimension to the data yields tighter parameter constraints.
To make the identification probability functions of Fig. <ref> more accurate, it appears promising to have an expert visually (and exhaustively, i.e. without imposing angular length thresholds) comb through a small representative region of LoTSS DR2 in search of giants.
The resulting data set would provide a better basis for determining the identification probability functions than the RGZ–ML or Oei et al., 2023 data sets used in this work.
We note that the Boötes LOFAR Deep Field search of Simonte et al., 2022 does not appear suited for this purpose, as the increased depth of this field renders it unrepresentative of LoTSS DR2 as a whole.
Finally, the model could be expanded in an attempt to measure cosmological evolution of, for example, $n_\mathrm{GRG}$.
However, we note that adding additional parameters to the model necessitates adopting more efficient inference techniques, such as Markov chain Monte Carlo or nested sampling.
The associated numerical gain would, in part, be negated by losing the speed-up associated to the likelihood trick of Appendix <ref>.

To determine the instantaneous VFF of GRG lobes in the Cosmic Web, $\mathcal{V}_\mathrm{GRG-CW}(z)$, one needs to calculate the mean GRG two-lobe proper volume $\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_{\mathrm{p,GRG}}](z)$.
To estimate the latter quantity, we have proposed to leverage the apparent power law relation between projected length and two-lobe volume shown in Fig. <ref>.
However, the fit in this work is based on data of just nine giants.
To improve this situation, we recommend expanding the capabilities and automating the parametric Bayesian lobe volume estimation method introduced by Oei et al., 2022, Oei et al., 2023.
This method could then be applied to thousands of our ML pipeline's masked cutouts, such as the one in Fig. <ref>.
This effort would increase the number of giants on which our fit of the projected length–two-lobe volume power law is based by two to three orders of magnitude.

§ CONCLUSIONS

In this work, we concatenated an existing crowd-sourced radio–optical catalogue, a new ML pipeline to automate radio–optical catalogue creation, and a Bayesian forward model to build a next-generation giant radio galaxy discovery and characterisation machine.
Applying this setup to the LOFAR Two-metre Sky Survey, we uncovered thousands of previously unknown giants, confirmed thousands of GRG candidates, and constrained the properties of the underlying population.

* The LoTSS is an ongoing sensitive, high-resolution, low-frequency radio survey whose second data release (DR2) covers $27\%$ of the northern sky.
As the number of detected sources already ranges in the millions, it has become unfeasible (at least for small scientific teams) to conduct manual visual searches for giants, in particular for those with angular lengths close to the lower limit of $1.3'$.

* To address this challenge, we scanned all $841$ LoTSS DR2 pointings – which together cover more than 5,000 square degrees of the northern sky – with an ML pipeline that crucially includes the convolutional neural network of Mostert et al., 2022, designed for the association of radio components for highly resolved radio galaxies, and an adapted version of the automated optical host galaxy identification heuristic developed by Barkus et al., 2022.
Used as a GRG detection system, our ML pipeline has a precision of $47\%$, a significant improvement over the $9\%$ precision obtained using the previously best ML GRG detection model [Proctor, 2016, Dabhade et al., 2020].
We merged the resulting giant candidate sample with that of the RGZ citizen science campaign [Hardcastle et al., 2023], homogenised the angular lengths, and subjected the candidates to a visual quality check.
The result is a sample of more than 8,000 newly confirmed giants, of which a large fraction are considered genuine beyond reasonable doubt.
More than $10^4$ unique giants have now been identified and published.

* We expand the population-based statistical forward model of Oei et al., 2023 designed to constrain the geometric properties of giants.
In particular, by modelling the PDF of the radio galaxy projected length RV $L_\mathrm{p}$ as a curved power law, we automatically also model the PDF of the giant radio galaxy projected length RV $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ as a curved power law.
We assume that these projected length distributions do not undergo intrinsic evolution between cosmological redshifts of $z = z_\mathrm{max}$ and $z = 0$, and likewise assume an intrinsically constant comoving GRG number density throughout this redshift range.
We model surface brightness selection by assuming a lognormal lobe surface brightness distribution at the survey's central frequency $\nu_\mathrm{obs}$, which is valid for radio galaxies of intrinsic proper length $l_\mathrm{ref}$ at redshift $z = 0$.
We relate lobe surface brightness distributions for radio galaxies of other lengths and at other redshifts to this reference distribution.
In addition, we model selection caused by the imperfect ability of search methods to identify all in-principle identifiable giants.
For this purpose, we use logistic functions of projected length $l_\mathrm{p}$ and redshift $z$.

* We then sought to identify all model parameter hexads that can reproduce the projected length–redshift histogram of the joint RGZ–ML–Oei et al., 2023 LoTSS DR2 GRG sample.
Through a simple Poissonian likelihood and a uniform prior distribution, we constructed a posterior distribution over the model parameters.
By confronting the model with an observed projected length–redshift histogram, rather than with an observed projected length distribution only (as has been done in Oei et al., 2023), we obtain tighter parameter constraints.

* We find evidence in support of the claim that the projected lengths of giant radio galaxies follow a curved power-law PDF, whose tail index is equal to $\xi(l_\mathrm{p,GRG}) = -2.8 \pm 0.2$ at $l_{\mathrm{p},1} = l_\mathrm{p,GRG} \coloneqq 0.7\ \mathrm{Mpc}$ and increases by $\Delta\xi = -2.4 \pm 0.3$ (i.e. decreases by $2.4 \mp 0.3$) in the projected length interval leading up to $l_{\mathrm{p},2} = 5\ \mathrm{Mpc}$.
The predicted median lobe surface brightness at $\nu_\mathrm{obs} = 150\ \mathrm{MHz}$, $l_\mathrm{ref} = 0.7\ \mathrm{Mpc}$, and $z = 0$ is equal to $b_{\nu,\mathrm{ref}} = 30 \pm 20\ \mathrm{Jy\ deg^{-2}}$.
This surface brightness level is lower than previously thought.
Tight degeneracies resembling inverse relations exist between $b_{\nu,\mathrm{ref}}$ and the reference surface brightness dispersion measure $\sigma_\mathrm{ref}$, and between $b_{\nu,\mathrm{ref}}$ and the GRG number density $n_\mathrm{GRG}$.
The latter relation suggests that giants might be more common than previously thought.
At $n_\mathrm{GRG} = 13 \pm 10\ (100\ \mathrm{Mpc})^{-3}$, giants appear to be of an abundance comparable to that of luminous non-giant radio galaxies.
We conclude that, at any moment in time, a significant fraction of the radio galaxy population is in a GRG phase.
As an immediate consequence, the fraction of radio galaxies that end their lives as giants must be even higher.

* Finally, we generate a posterior distribution for the instantaneous volume-filling fraction of GRG lobes in clusters and filaments of the Cosmic Web, $\mathcal{V}_\mathrm{GRG-CW}(z=0)$ – a key statistic required for determining the cosmic magnetogenesis potential of giants.
We find $\mathcal{V}_\mathrm{GRG-CW}(z=0) = 1.4 \pm 1.1 \cdot 10^{-5}$.
The mean two-lobe proper volume of a giant is $\mathbb{E}[V\ \vert\ L_\mathrm{p} \geq l_{\mathrm{p,GRG}}](z=0) = 5.1 \pm 0.3 \cdot 10^{-2}\ \mathrm{Mpc}^3$.
If a GRG population similar to that in the Local Universe has existed for most of the lifetime of the Universe, and IGM mixing processes are significant, then it appears possible that magnetic fields originating from giants permeate throughout significant (${\gtrsim}10\%$) fractions of today's Cosmic Web.

Using modern automation and inference techniques (which still leave significant room for future improvements), we conducted the most detailed study of the abundance and geometrical properties of giants to date.
These cosmic colossi may provide a previously underappreciated contribution to astrophysical magnetogenesis.

The full GRG catalogue with host identifications
will soon be available on the VizieR catalogue service hosted by the Centre de Données astronomiques de Strasbourg (CDS).
M.S.S.L. Oei and R.J. van Weeren acknowledge support from the VIDI research programme with project number 639.042.729, which is financed by the Dutch Research Council (NWO).
M.S.S.L. Oei also acknowledges support from the CAS–NWO programme for radio astronomy with project number 629.001.024, which is financed by the NWO.
Finally, M.S.S.L. Oei acknowledges support from the ERC Starting Grant ClusterWeb 804208.
L. Alegre is grateful for support from the UK STFC via CDT studentship grant ST/P006809/1.
M.J. Hardcastle acknowledges support from the UK STFC [ST/V000624/1].
B. Barkus is grateful for support from the UK STFC.
We would like to thank Huib Intema for enabling the cross-institute collaboration on the Leiden Observatory computer infrastructure.
We would like to thank Frits Sweijen for coding the very useful <https://github.com/tikk3r/legacystamps>.
This research has made use of the Python  package [The Astropy Collaboration et al., 2018]; the VizieR catalogue access tool [Ochsenbein et al., 2000], CDS, Strasbourg, France (DOI: 10.26093/cds/vizier); and the `Aladin Sky Atlas' developed at CDS, Strasbourg Observatory, France [Bonnarel et al., 2000, Boch & Fernique, 2014].
LOFAR data products were provided by the LOFAR Surveys Key Science project (LSKSP; <https://lofar-surveys.org/>) and were derived from observations with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The efforts of the LSKSP have benefited from funding from the European Research Council, NOVA, NWO, CNRS-INSU, the SURF Co-operative, the UK Science and Technology Funding Council and the Jülich Supercomputing Centre.
This publication uses data generated via the Zooniverse.org platform, development of which is funded by generous support, including a Global Impact Award from Google, and by a grant from the Alfred P. Sloan Foundation.

[Ahn et al., 2012]
Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ,
203, 21

[Alegre et al., 2022]
Alegre, L., Sabater, J., Best, P., et al. 2022, , 516, 4716

[Andernach et al., 2021]
Andernach, H., Jiménez-Andrade, E. F., & Willis, A. G. 2021,
Galaxies, 9, 99

[Barkus et al., 2022]
Barkus, B., Croston, J. H., Piotrowska, J., et al. 2022, , 509, 1

[Bassani et al., 2021]
Bassani, L., Ursini, F., Malizia, A., et al. 2021, , 500, 3111

[Beckmann et al., 2024]
Beckmann, R. S., Smethurst, R. J., Simmons, B. D., et al. 2024, ,
527, 10867

[Boch & Fernique, 2014]
Boch, T. & Fernique, P. 2014, in Astronomical Society of the Pacific
Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems
XXIII, ed. N. Manset & P. Forshay, 277

[Bonnarel et al., 2000]
Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, , 143,

[Breiman, 1997]
Breiman, L. 1997, Arcing the edge, Tech. rep., Citeseer

[Carretti et al., 2023]
Carretti, E., O'Sullivan, S. P., Vacca, V., et al. 2023, , 518,

[Cautun et al., 2014]
Cautun, M., van de Weygaert, R., Jones, B. J. T., & Frenk, C. S. 2014,
, 441, 2923

[Chambers et al., 2016]
Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv
e-prints, arXiv:1612.05560

[Codis et al., 2018]
Codis, S., Jindal, A., Chisari, N. E., et al. 2018, , 481, 4753

[Condon et al., 1998]
Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, , 115,

[Connor et al., 2022]
Connor, L., Bouman, K. L., Ravi, V., & Hallinan, G. 2022, , 514,

[Cutri et al., 2021]
Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online
Data Catalog, II/328

[Dabhade et al., 2020]
Dabhade, P., Mahato, M., Bagchi, J., et al. 2020a, ,
642, A153

[Dabhade et al., 2020]
Dabhade, P., Röttgering, H. J. A., Bagchi, J., et al.
2020b, , 635, A5

[Dabhade et al., 2023]
Dabhade, P., Saikia, D. J., & Mahato, M. 2023, Journal of Astrophysics
and Astronomy, 44, 13

[de Jong et al., 2022]
de Jong, J. M. G. H. J., van Weeren, R. J., Botteon, A., et al. 2022,
, 668, A107

[de Ruiter et al., 1977]
de Ruiter, H. R., Willis, A. G., & Arp, H. C. 1977, , 28, 211

[Delhaize et al., 2021]
Delhaize, J., Heywood, I., Prescott, M., et al. 2021, , 501, 3833

[Dewdney et al., 2009]
Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W.
2009, IEEE Proceedings, 97, 1482

[Dey et al., 2019]
Dey, A., Schlegel, D. J., Lang, D., et al. 2019, , 157, 168

[Duncan, 2022]
Duncan, K. J. 2022, , 512, 3662

[Enßlin, 2003]
Enßlin, T. A. 2003, , 399, 409

[Forero-Romero et al., 2009]
Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., &
Yepes, G. 2009, , 396, 1815

[Friedman, 2001]
Friedman, J. H. 2001, Annals of statistics, 1189

[Gheller & Vazza, 2019]
Gheller, C. & Vazza, F. 2019, , 486, 981

[Gheller et al., 2018]
Gheller, C., Vazza, F., & Bonafede, A. 2018, , 480, 3749

[Girshick, 2015]
Girshick, R. 2015, in Proceedings of the IEEE international conference on
computer vision, 1440–1448

[Govoni et al., 2019]
Govoni, F., Orrù, E., Bonafede, A., et al. 2019, Science, 364, 981

[Gürkan et al., 2022]
Gürkan, G., Prandoni, I., O'Brien, A., et al. 2022, , 512,

[Hallinan et al., 2019]
Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, in Bulletin of the
American Astronomical Society, Vol. 51, 255

[Hardcastle, 2018]
Hardcastle, M. J. 2018, , 475, 2768

[Hardcastle et al., 2023]
Hardcastle, M. J., Horton, M. A., Williams, W. L., et al. 2023, ,
678, A151

[Ineson et al., 2017]
Ineson, J., Croston, J. H., Hardcastle, M. J., & Mingo, B. 2017,
, 467, 1586

[Ishwara-Chandra et al., 2020]
Ishwara-Chandra, C. H., Taylor, A. R., Green, D. A., et al. 2020,
, 497, 5383

[Lakshmanan et al., 2021]
Lakshmanan, V., Görner, M., & Gillard, R. 2021, Practical machine learning
for computer vision (" O'Reilly Media, Inc.")

[Li et al., 2023]
Li, F., Zhang, H., Xu, H., et al. 2023, in 2023 IEEE/CVF Conference on
Computer Vision and Pattern Recognition (CVPR) (Los Alamitos, CA, USA: IEEE
Computer Society), 3041–3050

[Liu et al., 2021]
Liu, Z., Lin, Y., Cao, Y., et al. 2021, in Proceedings of the IEEE/CVF
International Conference on Computer Vision, 10012–10022

[Liu et al., 2022]
Liu, Z., Mao, H., Wu, C.-Y., et al. 2022, in Proceedings of the IEEE/CVF
Conference on Computer Vision and Pattern Recognition, 11976–11986

[Lochner & Bassett, 2021]
Lochner, M. & Bassett, B. A. 2021, Astronomy and Computing, 36, 100481

[Mahato et al., 2022]
Mahato, M., Dabhade, P., Saikia, D. J., et al. 2022, , 660, A59

[Masini et al., 2021]
Masini, A., Celotti, A., Grandi, P., Moravec, E., & Williams, W. L.
2021, , 650, A51

[Mohan & Rafferty, 2015]
Mohan, N. & Rafferty, D. 2015, PyBDSF: Python Blob Detection and Source
Finder, Astrophysics Source Code Library

[Mostert et al., 2022]
Mostert, R. I. J., Duncan, K. J., Alegre, L., et al. 2022, , 668,

[Mostert et al., 2023]
Mostert, R. I. J., Morganti, R., Brienza, M., et al. 2023, , 674,

[Murthy et al., 1993]
Murthy, S. K., Kasif, S., Salzberg, S., & Beigel, R. 1993, in Proceedings of
AAAI, Vol. 93, Citeseer, 322–327

[Norris et al., 2011]
Norris, R. P., Hopkins, A. M., Afonso, J., et al. 2011, , 28, 215

[Ochsenbein et al., 2000]
Ochsenbein, F., Bauer, P., & Marcout, J. 2000, , 143, 23

[Oei et al., 2023]
Oei, M. S. S. L., van Weeren, R. J., Gast, A. R. D. J. G. I. B., et al.
2023a, , 672, A163

[Oei et al., 2022]
Oei, M. S. S. L., van Weeren, R. J., Hardcastle, M. J., et al. 2022,
, 660, A2

[Oei et al., 2023]
Oei, M. S. S. L., van Weeren, R. J., Hardcastle, M. J., et al.
2023b, , in press

[Oei et al., 2023]
Oei, M. S. S. L., van Weeren, R. J., Hardcastle, M. J., et al.
2023c, , 518, 240

[O'Leary et al., 2021]
O'Leary, J. A., Moster, B. P., Naab, T., & Somerville, R. S. 2021,
, 501, 3215

[Pedregosa et al., 2011]
Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine
Learning Research, 12, 2825

[Planck Collaboration et al., 2020]
Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, , 641,

[Proctor, 2016]
Proctor, D. D. 2016, , 224, 18

[Rice, 2006]
Rice, J. 2006, Mathematical Statistics and Data Analysis, Advanced series
(Cengage Learning)

[Richter, 1975]
Richter, G. A. 1975, Astronomische Nachrichten, 296, 65

[Ringermacher & Mead, 2009]
Ringermacher, H. I. & Mead, L. R. 2009, , 397, 164

[Schlafly et al., 2019]
Schlafly, E. F., Meisner, A. M., & Green, G. M. 2019, , 240, 30

[Shimwell et al., 2022]
Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, , 659,

[Shimwell et al., 2017]
Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017,
, 598, A104

[Simonte et al., 2022]
Simonte, M., Andernach, H., Brüggen, M., et al. 2022, , 515,

[Subramanian, 2016]
Subramanian, K. 2016, Reports on Progress in Physics, 79, 076901

[Sutherland & Saunders, 1992]
Sutherland, W. & Saunders, W. 1992, , 259, 413

[Tang et al., 2020]
Tang, H., Scaife, A. M. M., Wong, O. I., et al. 2020, , 499, 68

[The Astropy Collaboration et al., 2018]
The Astropy Collaboration, Price-Whelan, A. M. Price-Whelan, A. M.,
Sipőcz, B. M., et al. 2018, , 156, 123

[Turner & Shabala, 2015]
Turner, R. J. & Shabala, S. S. 2015, , 806, 59

[van Haarlem et al., 2013]
van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, , 556,

[Vazza et al., 2017]
Vazza, F., Brüggen, M., Gheller, C., et al. 2017, Classical and
Quantum Gravity, 34, 234001

[Vazza et al., 2021]
Vazza, F., Locatelli, N., Rajpurohit, K., et al. 2021, Galaxies, 9, 109

[Wang et al., 2023]
Wang, W., Dai, J., Chen, Z., et al. 2023, in 2023 IEEE/CVF Conference on
Computer Vision and Pattern Recognition (CVPR), 14408–14419

[Williams et al., 2019]
Williams, W. L., Hardcastle, M. J., Best, P. N., et al. 2019, ,
622, A2

[Wright et al., 2010]
Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, ,
140, 1868

[Zhang et al., 2023]
Zhang, H., Li, F., Liu, S., et al. 2023, in The Eleventh International
Conference on Learning Representations

[Zheng et al., 2017]
Zheng, H., Tegmark, M., Dillon, J. S., et al. 2017, , 464, 3486

§ CURVED POWER LAW PDF FOR $L$

In Sect. <ref>, we have started modelling the geometry of radio galaxies at the level of the projected proper length RV $L_\mathrm{p}$.
While algebraically easier – when curved power laws are considered, at least – this approach is more limited than starting the forward model at the level of the intrinsic proper length RV $L$.
In this appendix, we calculate the distribution of $L_\mathrm{p}$ upon modelling $L$ with a curved power law.
Let us assume that, for $l \geq l_\mathrm{min}$,

\begin{align}
f_L(l) \propto \left(\frac{l}{l_\mathrm{min}}\right)^{\xi(l)},
\end{align}

where $\xi(l) = al+b$.
We now use the identity that for $f(x) = \left(\frac{x}{c}\right)^{ax+b}$, one finds

\begin{align}
\frac{\mathrm{d}f(x)}{\mathrm{d}x} = \left(\frac{x}{c}\right)^{ax+b}\left(a \ln{\frac{x}{c}} + a + \frac{b}{x}\right) = f(x)\left(a \ln{\frac{x}{c}} + a + \frac{b}{x}\right).
\end{align}

\begin{align}
\frac{\mathrm{d}f_L(l)}{\mathrm{d}l} = f_L(l)\left(a \ln{\frac{l}{l_\mathrm{min}}} + a + \frac{b}{l}\right),
\end{align}

\begin{align}
\frac{\mathrm{d}f_L(l_\mathrm{p}\eta)}{\mathrm{d}l_\mathrm{p}} = f_L(l_\mathrm{p}\eta)\left(a \ln{\frac{l_\mathrm{p}\eta}{l_\mathrm{min}}} + a + \frac{b}{l_\mathrm{p}\eta}\right)\eta.
\end{align}

Therefore, finding the PDF of $L_\mathrm{p}$ requires calculating three different integrals over $\eta$:

\begin{align}
f_{L_\mathrm{p}}(l_\mathrm{p}) = &-(1+b)\int_1^\infty \sqrt{1 - \frac{1}{\eta^2}}\ f_L(l_\mathrm{p} \eta)\ \mathrm{d}\eta\nonumber\\
&-l_\mathrm{p}a\left(1+\ln{\frac{l_\mathrm{p}}{l_\mathrm{min}}}\right)\int_1^\infty \sqrt{\eta^2-1}f_L(l_\mathrm{p}\eta)\ \mathrm{d}\eta\nonumber\\
&-l_\mathrm{p}a\int_1^\infty \sqrt{\eta^2-1}f_L(l_\mathrm{p}\eta)\ln{\eta}\ \mathrm{d}\eta\ \ \ \ \text{for } l_\mathrm{p} > l_\mathrm{min}.
\end{align}

The PDF of $L_\mathrm{p}\ \vert\ L_\mathrm{p} \geq l_\mathrm{p,GRG}$ follows through Eq. <ref>.

§ LIKELIHOOD TRICK