RAiSERed: radio continuum redshifts for lobed AGNs

The dynamical model-based estimate of the transverse comoving distance to an active, lobed radio source at a trial redshift $z$ can be expressed in terms of the flux density, $S_{\nu}$ (single lobe at observer-frame frequency $\nu$); angular size, $\theta$ (single lobe); axis ratio, $A$ (single lobe length divided by the lobe radius); and properties of the electron energy distribution. Particles injected into the lobe are initially described by a power law distribution of energies $N(E)=N_{0}E^{-s}$, where $N_{0}$ is a constant and $s=2\alpha_{\rm inj}+1$ for injection-time spectral index $\alpha_{\rm inj}>0.5$. The spectral index of the electron energy population steepens to $\alpha_{\rm inj}+0.5$ above the ‘optically-thin’ break frequency, $\nu_{\rm b}$, due to synchrotron radiative losses as the source ages. These parameters describing the electron population are constrained observationally by fitting the radio spectrum with the continuous injection (CI) model (Turner et al., 2018b). However, the spectral shape of remnants and jetted FR-Is are not fitted well by the functional form of the CI model (e.g. 3C28 and 3C31 in Harwood, 2017), enabling these objects to be excluded even without high-resolution imaging. Meanwhile, the flux density, source size and axis ratio are readily measured from high-resolution radio images. Importantly, the monochromatic flux density, $S_{\nu}$, used in analytic models without a full treatment of radiative loss mechanisms (cf. RAiSE; Turner & Shabala, 2015) must be measured at a frequency below the break frequency ($\nu<\nu_{\rm b}$).

Following Turner & Shabala (2019), the expected distance to a lobed FR-II radio source at some trial redshift $z$ is given by:

where $q=u_{\rm B}/u_{\rm e}\sim 0.019$ is the ratio of the energy in the magnetic field to that in the particles (i.e. equipartition factor; Turner et al., 2018b), and $\gamma\sim 300$ is the minimum Lorentz factor of the electron population after expanding from the jet terminal hotspot to the lobe (Turner & Shabala, 2019); the maximum Lorentz factor is neglected in the model as it has no significant effect on our results for $\gamma_{\rm max}\gg\gamma$. The density profile, of the total gas mass, into which the radio source expands (see Section 2.2) is approximated locally as a power law of the form $\rho(r)=\rho[r/a]^{-\beta}$; the gas density $\rho$ is defined at some galactocentric radius $a$, which we arbitrarily set as the physical equivalent of the angular size of the lobe at the trial redshift $z$. Meanwhile, the exponents $x$ and $y$ are functions of $s$, $\beta$ and the lobe magnetic field strength, $B$, defined immediately after Equation 13 in Turner & Shabala (2019). The lobe magnetic field strength at the trial redshift $z$ is estimated using the following equation derived from Equations 3, 12 and 13 of Alexander (2000):

where $\mu_{0}$ is the vacuum permeability, and $d_{{\rm M}\>\!(H_{0},\Omega_{0})}(z)$ is the transverse comoving distance at the trial redshift $z$ for the concordance cosmological model.

The constants of proportionality $f_{1}$ and $f_{3}$ in the distance equation (Equation 1) are functions of $s$, and $A$, $B$, $\beta$ and $z$, respectively defined as:

where $\sigma_{\rm T}$ is the electron scattering cross-section, $e$ and $m_{\rm e}$ are the electron charge and mass, $c$ is the speed of light, $\upsilon$ is a constant defined in Equation 5 of Turner et al. (2018b), and $\Gamma_{\rm c}=4/3$ and $\Gamma_{\rm x}=5/3$ are the adiabatic indices of the lobe plasma and external medium respectively. The time-average of the synchrotron radiative loss function, $\overline{\mathcal{Y}(t,\nu)}$, is contrained to be a constant value in the range 0.3-0.5 based on RAiSE simulations (Turner & Shabala, 2019). Meanwhile, the functions $\kappa(B,z)$ and $\sigma(B,z)$, defined in Equation 10 of Turner & Shabala (2020), are used to maintain an analytic solution across all lobe magnetic field strengths. The ratio of the lobe to expansion surface pressures is found from the numerical simulations of Kaiser & Alexander (1999) to be well modelled as (Equation 7 of Kaiser, 2000):

The key model parameters that can be constrained through observations, either for individual sources or as a population average, are summarised in Table 1.

The dynamical model-based estimate for the transverse comoving distance can be calibrated using a modest-sized sample of radio AGNs with known spectroscopic redshifts. Calibration constants are included in the distance equation (Equation 1) to provide small corrections to the most poorly constrained, albeit well-informed, model parameters. The minimum Lorentz factor and equipartition factor of the electron population in the lobe cannot realistically be constrained for individual sources and also faces moderate uncertainty in the mean value across radio sources; systematic errors in the mean values assumed for these parameters are handled using the calibration constants $b_{1}$ and $b_{2}$ respectively. The particle content in the lobe (i.e. ratio of leptons to baryons) is also modified using the $b_{1}$ calibration constant. Meanwhile, although the ratio of the lobe and expansion surface pressures is constrained by numerical simulations (Equation 4, above), model assumptions likely do not perfectly represent actual sources, leading to potential large errors in the thinnest sources (i.e. $\chi\propto A^{-1.5}$ for $\beta\sim 1$); systematic errors in the exponent on the axis ratio are handled using the calibration constant $b_{3}$. Finally, the gas density at the radius of the lobe is modelled using density profiles for observed clusters and cosmological simulations (see Section 2.2), however this modelling is subject both to small errors in the absolute scaling of the gas density and the cosmic evolution of the cluster mass function; these systematic errors are handled using calibration constants $b_{2}$ (shared with equipartition factor) and $b_{4}$ in the form $b_{2}[1+z]^{b_{4}}$. The four calibration constants are assumed to be zero by default (i.e. no correction required), however we discuss the use of these parameters in Section 4 when examining a sample of observed high-redshift sources.

The semi-analytic galaxy evolution (SAGE; Croton et al., 2016; Raouf et al., 2017) model (update of original model; Croton et al., 2006) traces the evolutionary history of baryonic matter on top of existing large scale simulations of the dark matter. In this work we use the Bolshoi model (Klypin et al., 2011) to construct the dark matter framework for SAGE. This simulation traces galaxies in a box of side length $250\rm\,Mpc$$/h$ through cosmic time with outputs at regular intervals in redshift. The Bolshoi dark matter simulation assumes a WMAP5 cosmology with $\Omega_{\rm m}=0.27$, $\Omega_{\Lambda}=0.73$, $\Omega_{\rm b}=0.0469$, $\sigma_{8}=0.82$, $h=0.70$ and $n=0.95$; the products we derive from this simulation are quite insensitive to small changes to the cosmological parameters (see Section 3.2.3). We process Bolshoi simulation outputs at 15 redshifts up to $z=5$ using SAGE to obtain mock galaxy populations across cosmic time. The RAiSERed code extrapolates empirical relationships derived for $0<z<5$ beyond this maximum redshift to provide a good approximation of the distance to ultra-high redshift sources; however, in this work we conservatively limit ourseleves to objects below $z=6$. Based on the local 200 MHz luminosity function (Franzen et al.,, in preparation), and assuming a luminosity evolution of $\sim$$(1+z)^{4\pm 0.5}$ (Seymour et al.,, in preparation), the Bolshoi simulation volume is expected to include $0.7_{-0.4}^{+0.8}$ powerful radio galaxies with luminosities $L_{200\rm\,MHz}>10^{28}\rm\,W\,Hz^{-1}$. The empirical relationships (discussed below) will therefore be derived, or more likely validated, using at least one example of even the most massive host galaxies expected for radio AGNs.

The mock galaxy populations based on the Bolshoi simulation are used to inform both the cluster mass function and the gas fraction as a function of redshift. Following Turner & Shabala (2020), we assume the low-redshift mass function of Girardi & Giuricin (2000) who find that a common Schechter function describes both galaxy groups and clusters. The Bolshoi simulations are used to extend their observations to higher redshifts by finding that the mass of the break in the Schechter function scales with redshift as approximately $(1+z)^{-3}$. The relative normalisation of the cluster mass function between redshifts is not important in this work. The cluster mass function is further weighted by the AGN duty cycle (Pope et al., 2012); high mass black holes and galaxies are known to have an enhanced probability of hosting AGNs compared to their lower mass counterparts (e.g. Sabater et al., 2019). Meanwhile, the cosmological evolution of the gas fraction as a function of halo mass is found to be well described by $f_{\rm gas}=10^{w_{0}(z)-14w_{1}}{(M_{\rm halo}/\rm M_{\odot})}^{w_{1}}$, where $M_{\rm halo}$ is the simulated dark matter halo mass taken from SAGE, and $w_{0}(z)=\text{max}\{-0.88-0.03z,-0.92+0.001z\}$ and $w_{1}=0.05$ describe the redshift dependence. These gas fractions are further scaled by a small, constant factor based on low-redshift observations (McGaugh et al., 2010; Gonzalez et al., 2013). The spread in the distribution of gas fractions is approximately constant in log-space for all halo masses, but depends on redshift as $\delta f_{\rm gas}=0.05-0.002z\,\rm dex$.

The shape of the gas density profiles is based on the Vikhlinin et al. (2006) cluster observations following the method described in Turner & Shabala (2015)¹¹1The gas density profile of the host galaxy can be neglected since core-confined radio AGNs (that are significantly impacted by their host galaxies) are excluded based on their ‘optically-thick’ low-frequency spectral turnover (i.e. due to free-free or synchrotron self-absorption).. The same shape profile can describe the environments of clusters with very different masses by scaling with the core density and virial radius of the halo. The virial radius is directly related to the mass of its halo through $r_{\rm halo}=(GM_{\rm halo}/[100H^{2}(z)])^{1/3}$, where $H(z)$ is the Hubble constant at redshift $z$, and $G$ is Newton’s gravitation constant. The Vikhlinin et al. (2006) gas density profile is integrated over the volume of the cluster within the virial radius and this unscaled mass is compared with the gas mass simulated using SAGE (based on the halo mass and gas fraction) to derive the density $\rho$ at the end of the lobe of radius $a$. The shape of the density profile is approximated locally using a power law of the form $\rho_{\rm gas}(r)=\rho\>\![r/a]^{-\beta}$ to enable analytic modelling of the radio source evolution.

The calibration constants in the RAiSERed model described in the previous sections can be constrained using observations of at least five²²2There are four calibration constants so an absolute minimum of five independent measurements are required. radio AGNs with known spectroscopic redshifts. The sum of squared differences between the radio continuum redshifts (i.e. mean of the probability density function) and the spectroscopic redshifts for the calibration sample is minimised to find the optimal values for $b_{1}$, $b_{2}$, $b_{3}$ and $b_{4}$. The optimisation is performed using the noisyopt package in python, a robust pattern search algorithm with an adaptive number of function evaluations (Mayer et al., 2016). Importantly, this algorithm does not use the function derivative to locate the minimum, an approach that would not be viable in this work given random noise is more significant that the gradient for small perturbations in the calibration constants. Despite the greatly improved computational efficiency of this algorithm over a brute force technique the computational time for a sample of 15 objects is approximately 4-6 hours on a typical laptop. The calibration of the RAiSERed model for use in large sky surveys will likely need to use supercomputer time or only a small subsample of calibrators.

Cygnus A is a low-redshift radio source (spectroscopic redshift of $z=0.056075\pm 0.000067$; Owen et al., 1997) with a double FR-II lobe morphology. Steenbrugge et al. (2010) reported core and hotspot removed measurements of the flux density in the east and west lobes of Cygnus A at six frequencies from $151\rm\,MHz$ to $15\rm\,GHz$. The size and axis ratio of the two lobes are measured from the $5\rm\,GHz$ radio images following Turner & Shabala (2019); the uncertainty in the length of each lobe is taken as the size of the synthesised beam at this frequency ($\theta_{\rm res}$; Carilli et al., 1996). The properties of the electron population in each lobe are constrained by fitting the radio spectra using the continuous injection (CI) model following the method of Turner et al. (2018b). The $151\rm\,MHz$ to $15\rm\,GHz$ flux density measurements are supplemented by $74\rm\,MHz$ observations from Cohen et al. (2007) to better constrain the electron energy injection index, $s$; the uncertainties on the Steenbrugge et al. (2010) measurements are estimated following Carilli et al. (1996). The radio continuum observations for the two lobes of Cygnus A are summarised in Table 2.

The RAiSERed model is tested for radio AGNs across redshifts $z=0$ to 6 by shifting the observed radio continuum observations of Cygnus A to higher redshifts. The flux density and lobe angular size are converted for the cosmology assumed in this work, whilst intrinsic properties of the source (axis ratio and electron energy injection index) are unchanged. The observer-frame break frequency measured for an identical source to Cygnus A (at redshift $z_{0}$) but located at redshift $z$ is found from rearranging Equation 9 of Turner & Shabala (2019),

where $B$ is the magnetic field in the lobe. The magnetic field strength is constrained using the dynamical model in Equation 2 for each lobe to provide an estimate calibrated for any systematic errors in the model, based on the obsevred properties of Cygnus A at $z=0.056$.

The radio continuum attributes of the east and west lobes of Cygnus A are shifted from $z=0.056$ to progressively higher redshifts up to $z=6$. The cosmological environments assumed in our model are not strictly valid for the Cygnus A-like sources shifted to substanially higher redshifts (i.e. the observed properties are influenced by the actual $z\sim 0.056$ environment). However, in order to maintain the correct model response to small variations in redshift we choose to use the cosmological environment appropriate to each trial redshift ³³3The findings in this section are unchanged if we assume the $z=0.056$ environment at all redshifts, although the model becomes less stable around redshifts $z=1.5-2$, coinciding with the turnover in the angular diameter distance–redshift function.. The best estimate radio continuum redshifts for the east and west lobes of these Cygnus A-like sources are shown in Figure 2 as a function of their true redshifts (i.e. shifted spectroscopic redshifts). The radio continuum redshifts are consistent with the true redshifts within the 1$\sigma$ measurement uncertainties for the majority of the simulated Cygnus A-like sources, except those at $z=0.15$-0.6 and from $z=1$-2.5 (these are consistent at the 2$\sigma$ level). We note that the absolute scaling need not be correct at this stage since our model includes calibration terms (i.e. $b_{1}$, $b_{2}$, $b_{3}$ and $b_{4}$); however, this agreement supports both our modelling of the cosmological environments and the (uncalibrated) values chosen for the less well constrained model parameters ($q$, $\gamma$, $\mathcal{Y}$).

In this section, the simulated population of Cygnus A-like sources is used to test the ability of the RAiSERed model to estimate redshifts for unresolved sources (i.e. size upper limits) and objects with poorly constrained break frequencies. The other radio continuum observables have smaller (absolute) exponents in the distance measure equation (Turner & Shabala, 2019) and thus modest uncertainties are not expected to adversely effect the redshift estimates.

The RAiSERed model is applied to the two lobes of Cygnus A, two lobes of the small subsample of 3C sources, and the combined lobes of each of the five HeRGE objects; i.e. 17 lobes are fitted with radio continuum redshifts. The redshift probability density functions for the majority of objects has a single sharp peak located approximately at the spectroscopic redshift. However, the probability density function for both lobes of 3C219 have two peaks, one at $z\approx 0.2$ and another at $z\approx 1.5$. The high-redshift peak is associated with lobe expansion rates at close to the speed of light and dissapears in the south lobe if the prior probability density function is modified to restrict velocities to be less than $0.5c$. That is, our technique can still find the correct radio continuum redshift for this source. The uncalibrated redshift distributions found using our technique for the remaining 15 lobes are plotted in Figure 7(left) as a function of their known spectroscopic redshifts. The radio continuum redshifts are strongly correlated with the spectroscopic redshifts although the high-redshift HeRGE sample has their redshifts systematically underestimated. The uncalibrated model has an average log-space error of $\delta\log(1+z^{*})=0.069\rm\,dex$ (i.e. 17% in $1+z^{*}$). Importantly, the uncalibrated radio continuum redshifts explain 70% of the variation in the known spectroscopic redshift for the high-redshift HeRGE sample, compared with at most 27% for any one of the observables in isolation.

The calibration constants are constrained (initally) using the radio continuum attributes of the 15 lobes (from Cygnus A, and the 3C and HeRGE subsamples). The fitted values are found as $b_{1}=0.76\pm 0.07$, $b_{2}=0.01\pm 0.01$, $b_{3}=-0.60\pm 0.07$ and $b_{4}=-0.33\pm 0.03$. These fits indicate that our parameter associated with the minimum Lorentz factor is a factor of 5.8 too low, the equipartition factor and gas density at redshift $z=0$ are correct, the lobe to hotspot pressure ratio has a weaker than expected dependence on the axis ratio, and that the cosmological environments predict gas densities a factor of $(1+z)^{0.33}$ too high. However, caution should be taken in directly associating these fits with their intended purpose (as stated above) as the least squares optimisation will also attempt to explain other more minor factors with these calibration constants. The redshift probability density functions for the 15 calibrator lobes are shown in Figure 7(right) as a function of their spectroscopic redshifts. The correlation between these calibrated radio continuum redshifts and their spectrocopic redshifts is now centred on the one-to-one line with an average log-space error of $\delta\log(1+z^{*})=0.040\rm\,dex$ (i.e. 9.6% in $1+z^{*}$). Combining the probability density functions of the two lobes for each low-redshift source to yield a single robust radio continuum redshift does not change these results; however, the increased uncertainties on sources with less consistent estimates improves the agreement between these robust redshifts and their spectroscopic counterparts.

We further quantify the likelihood that the calibration constants improve the accuracy of redshift measurements (compared to our uncalibrated model) for sources without known spectroscopic redshifts. This is achieved by randomly selecting six or nine objects as calibrators with the same ratio of low- and high-redshift lobes as our original sample (i.e. one-third HeRGE, two-thirds Cygnus A/3C); the remainder are specified to not have spectroscopic redshifts in the RAiSERed code. Radio continuum redshifts are calculated for these remaining nine or six lobes respectively based on the calibration constants fitted using the randomly selected calibrators. This process is repeated ten times, selecting a different set of calibrators in each iteration. The redshifts estimated for the ten randomly selected subsets of lobes not used as calibrators are plotted in Figure 8 as a function of their spectroscopic redshift.

The radio continuum redshifts estimated when using six calibrators are largely consistent with their spectroscopic counterparts and centred on the one-to-one line. However, the predicted redshift for USS 1707$+$105 is highly sensitive to the choice of calibrators, with estimates ranging from $z=3$ to 7; the axis ratio upper limit is likely the cause of the increased variability in this source. As a result of this source, the likelihood of agreement with the spectrocopic redshifts actually decreases (compared to our uncalibrated model) when calibrating the model with six objects; the average error, across both redshift and the ten random source selections, increases to $\delta\log(1+z^{*})=0.098\rm\,dex$ (i.e. 25% in $1+z^{*}$)⁴⁴4The error is 14% in $1+z^{*}$ if USS 1707$+$105 is excluded from the calculation.. By contrast, the radio continuum redshifts estimated using nine calibrators are much more stable due to the greater exploration of parameter space in the calibration. Variability in the estimated redshift of USS 1707$+$105 is now minimised between calibrator samples. Calibration of the RAiSERed model with nine sources increases the likelihood of agreement between the radio continuum and spectroscopic redshifts for all objects. The average error converges towards that found when using all 15 objects as calibrators; i.e. $\delta\log(1+z^{*})=0.058\rm\,dex$ (i.e. 14% in $1+z^{*}$). The error in our model is empirically related to the number of degrees of freedom, $\text{df}=N-4$ (for $N$ calibrators), as $\delta\log(1+z^{*})=0.139\;\!\text{df}^{\>\!-0.53}\rm\,dex$; i.e. the error is less than 10% with 14 calibrators and, extrapolating, falls to less than 5% with 40 sources. The proposed calibration technique will therefore be successful if a modest number of sources are used covering a comparable range of parameter space to the target sources lacking spectroscopic redshifts.

There are a number of metrics used for optical photometric redshifts (e.g. Tanaka et al., 2018; Schmidt et al., 2020) that can similarly be applied to test the accuracy of the redshift probability density functions as an estimator of the true spectroscopic redshift. In particular, the probability integral transformation (PIT) is the cumulative distribution function (CDF) of a redshift probability density function evaluated at its spectroscopic redshift:

The distribution of PIT values for accurate probability density functions is expected to be uniform from 0 to 1. An ensemble of overly broad probability density functions will produce an excess of PIT values around 0.5, and conversely, narrow distributions lead to an excess of the lowest and highest values. In this manner, the distribution of PIT values can be used to probe the average accuracy of the redshift probability density functions for an ensemble of radio AGNs.

The accuracy of the redshift probability density functions for the calibrated subsamples considered in the previous section is assessed in Figure 9 using the probability integral transformation. These probability density functions are updated to include the model uncertainty of the RAiSERed code. The distribution of the PIT values for the random subsamples (of six objects) calibrated using nine calibrators is evenly spread from 0 to 1, though perhaps is slightly right skewed. Meanwhile, the distribution for the random subsamples (of nine objects) calibrated using six calibrators is drawn towards the centre (i.e. 0.5). Both sets of calibrated subsamples are reasonably consistent with the expected distribution of PIT values considering the size of the sample. Importantly, there are very few sources with extreme PIT values (e.g. $<$ 0.01 or $>$ 0.99) which would indicate outliers.

The distribution of PIT values for the two sets of subsamples are assessed quantitatively using the Kolmogorov-Smirnov and Anderson-Darling statistics. The Kolmogorov-Smirnov statistic measures the maximum difference between the CDF of the distribution of PIT values for the test subsample of radio AGNs and the CDF for the expected PIT distribution (i.e. uniform from 0 to 1). The Anderson-Darling statistic is a variant of the Kolmogorov-Smirnov statistic; see Equations 6 and 8 of Schmidt et al. (2020). The Kolmogorov-Smirnov statistic for the subsample using six calibrators is $K\!S=0.19$ ($p=0.007$) and for the subsample with nine calibrators is $K\!S=0.22$ ($p=0.016$); both have $p$-values of approximately 0.01 suggesting the PIT distributions are likely not consistent with the expected uniform distribution. Similarly, the Anderson-Darling statistics for the two sets of subsamples are $A\!D=8.3$ ($p=0.001$) and $A\!D=4.1$ ($p=0.007$) for six and nine calibrators respectively. By contrast, when considering the full sample of 15 calibrators the Anderson-Darling statistic drops to $A\!D=0.42$ ($p=0.22$). The average accuracy of the RAiSERed redshift estimates therefore becomes consistent with the expected PIT distribution if a sufficient number of objects are used in the calibration.

The authors confirm that the data supporting the findings of this study are available within the article and the relevant code is included in the supplementary materials.