2.1. Input catalogues

2.2. Selection criteria for HzRG candidates

2.2.2. GLEAM spectral properties: Steepness and curvature

Using a similar approach to our pilot study, we then fitted a model to the GLEAM flux density data for each of the remaining sources. In the pilot project, a second-order polynomial was fitted in $\log\!{(S_{\nu}})$ – $\log\!{(\nu)}$ space:

(1) \begin{equation}\log\!(S_{\nu}) = \alpha\log\!\left(\frac{\nu}{\nu_{0}}\right) + \beta\log^2\!\left(\frac{\nu}{\nu_{0}}\right) + \log\!(S_{0}),\end{equation}

where $\alpha$ is the spectral index at reference frequency $\nu_{0} = 151$ MHz, i.e. at the centre of the GLEAM band, $\beta$ the curvature term, and $S_{0}$ the flux density at $\nu_{0}$ . In this study, however, the fitting was done in linear space to preserve the Gaussian characteristics of the flux density errors, which is especially important for the fainter sources that we considered (see below). Equation (1) is then equivalent to

(2) \begin{equation}S_{\nu} = S_{0}\!\left(\frac{\nu}{\nu_{0}}\right)^{\alpha} \times 10^{\beta\log^2\!\left(\frac{\nu}{\nu_{0}}\right)}.\end{equation}

For each sub-band flux density, the uncertainty was calculated by combining the fitting uncertainty from the catalogue and the internal GLEAM flux density calibration uncertainty, the latter being 2% for the targets of interest in this study (Hurley-Walker et al. Reference Hurley-Walker2017; Franzen et al. Reference Franzen, Hurley-Walker, White, Hancock, Seymour, Kapińska, Staveley-Smith and Wayth2021). Correlations between the sub-band flux densities (e.g. see Hurley-Walker et al. Reference Hurley-Walker2017) were not modelled; each sub-band flux density was assumed to be an independent measurement. To first order, this is not expected to affect the accuracy of the spectral steepness/curvature selection technique. We did, however, take the correlations into account when modelling the broadband radio spectra (Section 5.4).

The next step was to isolate those sources with significantly curved spectra (step 6 in Table 1). The scientific rationale for this step follows the same argument presented in D20, that is many well-studied lower-redshift radio galaxies have observed-frame radio spectra that begin to flatten or turn over at low frequencies, and by ‘shifting’ these sources to larger distances (higher redshifts) we can predict the optimal region within the $\alpha$ – $\beta$ parameter space to search for HzRGs. To carry out step 6, we used fitting criteria of (i) a reduced chi-squared goodness-of-fit statistic $\chi^2_{\rm red} < 2$ (probability of obtaining a more extreme $\chi^2_{\rm red}$ by chance is $p \approx 0.008$ ) and (ii) $\lvert\beta\rvert > \sigma_{\beta}$ , where $\sigma_{\beta}$ is the uncertainty for $\beta$ . The latter criterion generally filters out those sources that are better fitted with a single power law, as $\beta$ will be close to zero for these cases. However, for three sources in our sample, J0053 $-$ 3256, J1037 $-$ 0325 and J2311 $-$ 3359, $\chi^2_{\rm red}$ is slightly larger for a curved fit compared with a single-power-law fit (i.e. the first part of Equation (2): $S_{\nu} = S_{0}(\nu/\nu_{0})^{\alpha}$ ) across the GLEAM band: $\Delta \chi^2_{\rm red}$ is in the range $2.6 \times 10^{-4}$ – 0.036. These sources do not fully meet all of our selection criteria, and we provide further details in Section 2.3. More generally, a comparison of the reduced chi-squared values for curved and single-power-law fits suggests that we may have overfitted about 6% of the sources classified as having curved GLEAM spectra.

In addition, we used a fitted flux density cutoff $S_0 \geq 40$ mJy, i.e. an order of magnitude fainter than in the D20 pilot project. Such a cutoff enables the discovery of less luminous sources such as J1530 $+$ 1049, in addition to powerful radio galaxies such as J0856 $+$ 0223 and J0924 $-$ 2201 (and possibly J0917 $-$ 0012). We also removed the $S_{151} < 1.0$ Jy upper flux density cutoff that was used in the pilot project. For example, some radio galaxies with $z > 4$ have $S_{150} > 1$ Jy (see e.g. Table 4 in Saxena et al. Reference Saxena2018b).

The distribution of the remaining 3327 sources in the $\alpha$ – $\beta$ parameter space is shown in Figure 1. Using the same argument as in the pilot study for the tracks that sources follow in this parameter space as they are progressively redshifted (Figure 1 in D20; additional examples shown in Figure 1), our next selection criterion (step 7 in Table 1) was $\alpha \leq -0.7 \cap \beta \leq -0.2$ . The first part of this expression is almost identical to the spectral steepness criterion used in D20 ( $\alpha < -0.7$ ), while the second part was relaxed from $-1.0 < \beta < -0.4$ in the pilot project to a wider range, given the potential trajectories of the tracks mentioned above. The total number of sources that remained after this step was 1187.

Figure 1. The $\alpha$ – $\beta$ parameter space for all isolated and compact GLEAM sources in the VIKING survey area that have curved low-frequency radio spectra with $S_{0} \geq 40$ mJy (see Section 2.2.2). Error bars on individual data points are not shown for the sake of clarity; however, in the bottom-right corner of the panel, we indicate the median uncertainties in $\alpha$ ( $\pm 0.08$ ) and $\beta$ ( $\pm 0.55$ ) for the sources in our sample. The paucity of sources near $\beta=0$ is due to our $\lvert\beta\rvert > \sigma_{\beta}$ selection criterion. HzRG candidates were mostly selected in the region of the plot delineated by the blue solid lines ( $\alpha \leq -0.7 \cap \beta \leq -0.2$ ). The seven sources that do not fully meet our selection criteria are marked (see discussion in Section 2.3). We do not show the full range of $\alpha$ and $\beta$ values in the plot; indeed, extreme values often correspond to much poorer fits. One source in our sample, J1141 $-$ 0158, is not visible in the figure as it has a very large curvature term ( $\beta = -6.5 \pm 2.7$ ); in this case, the fit is likely significantly affected by the S/N in the lower part of the GLEAM band (see Figure 2). J0856 $+$ 0223 and J0917 $-$ 0012 from the pilot study, also in our larger sample, are marked separately. Furthermore, we show the HzRGs J0924 $-$ 2201 ( $\alpha = -0.26 \pm 0.04$ ; $\beta = -1.42 \pm 0.27$ ) and J1530 $+$ 1049 ( $\alpha = -1.31 \pm 0.28$ ; $\beta = -4.8 \pm 2.2$ ), both of which are detected in GLEAM in other parts of the sky not covered by VIKING. We can see that J0924 $-$ 2201 would be missed by our selection criteria; this is because the source turns over in the middle of the GLEAM band (Figure 11 in Callingham et al. Reference Callingham2017; also see Section 6.2 and Figure 3 in this paper). J1530 $+$ 1049 has a similar USS spectral index in the GLEAM band compared with its spectrum between TGSS and NVSS/FIRST ( $\alpha^{1400}_{147.5} = -1.4 \pm 0.1$ ; Saxena et al. Reference Saxena2018a,b), and $\lvert\beta\rvert$ , while large, also has a large uncertainty due to low S/N measurements at 92, 204, and 212 MHz. We analyse the broadband spectrum of J1530 $+$ 1049 in Section 6.2 and Figure 3. Lastly, we plot five tracks corresponding to the predicted observed-frame values of $\alpha$ and $\beta$ of a model source at a given redshift with particular rest-frame spectral indices and break frequencies. Along each track, redshift increases from right to left, and the star markers indicate $z=$ 2, 4, 6 and 8. The tracks shown are not an exhaustive range of possibilities, but instead illustrate how the track trajectories can change depending on the underlying rest-frame spectrum of the source. The track model parameters correspond to the smoothly varying double- and triple-power-law fits that we use in this paper to model the broadband spectra of our targets (Equations (5) and (6)). Tracks 1, 2, 3, 4 and 5 move into our selection region at redshifts of 9.3 (not shown in the panel), 4.5, 2.7, 3.4 and 6.0, respectively. The track trajectories would continue to the upper left if the redshift is increased. If the (lower) rest-frame break frequency is below/above 500 MHz, then a particular track trajectory is still the same as plotted, but the position corresponding to a given redshift is further along the track to the left/right.

Although we applied a number of selection criteria above to restrict the list of sources to those that are compact and isolated, source blending remains a potential issue that must be considered carefully given the relatively low angular resolution of the GLEAM data. This is particularly relevant regarding the reliability of the $\alpha$ and $\beta$ measurements. We discuss this potential issue further in Section 5.4.3.

2.3. The GLEAM–VIKING HzRG candidate sample

After applying the criteria described in Section 2.2.3, we were left with a sample of 53 sources within the VIKING survey region (step 10 in Table 1). The sample comprises 51 new HzRG candidates and both J0856 $+$ 0223 and J0917 $-$ 0012 from our pilot study (Section 1), which satisfy the refined selection criteria used in this paper.Footnote ^e Note that the sample includes six sources that do not fully meet the $\alpha$ and $\beta$ selection criteria: J0034 $-$ 3112, J0053 $-$ 3256, J0129 $-$ 3109, J1037 $-$ 0325, J1246 $-$ 0017 and J2311 $-$ 3359. For all of these sources bar one, this is because (i) they marginally fall outside of our selection region in $\alpha$ – $\beta$ parameter space (Figure 1), yet within the $1\sigma$ uncertainties they are consistent with the selection criteria, or (ii) $\lvert\beta\rvert$ is marginally smaller than $\sigma_{\beta}$ . Our original approach for the $\alpha$ – $\beta$ fitting had been to use Equation (1), and due to the propagation of the flux density uncertainties in $\log\!{(S_{\nu}})$ – $\log\!{(\nu)}$ space, these sources initially fully satisfied our selection criteria and had been followed up with the ATCA (Section 4). The remaining source, J2311 $-$ 3359 from the GAMA-23 field, was removed in step 2 in Table 1, as it is too faint to be detected in NVSS; moreover, it does not satisfy our $\lvert\beta\rvert > \sigma_{\beta}$ criterion, and ASKAP early science data suggests an LAS of $5.\!^{\prime\prime}6$ (see Section 5.5). J2311 $-$ 3359 was identified as an HzRG candidate of particular interest given its USS nature ( $\alpha = -1.58 \pm 0.19$ in GLEAM) and followed up with the ATCA before the selection criteria for this project were fully defined.

In addition to these six sources, our sample includes the source J1335 $+$ 0112. This source meets all of our selection criteria except that it has an AllWISE identification (which we discuss later in Section 5.5). Unfortunately, this source was followed up with the ATCA before our $2^{\prime\prime}$ AllWISE criteria was introduced (it had been $1^{\prime\prime}$ ). As we outline in Section 5.5, J1335 $+$ 0112 may still be a high-redshift target, and therefore we decided that there was sufficient scientific interest to include the source in our sample.

The sample is presented in Table 2, including the $\alpha$ and $\beta$ values for each source and the fitted 151-MHz GLEAM flux density. The fitted 151-MHz flux densities span the range 47.6–2439 mJy, with a median of 260.6 mJy. We have marked the locations of the sources in our sample in the $\alpha$ – $\beta$ parameter space in Figure 1. The median spectral index of the sample at 151 MHz is $-0.96$ , and only three sources would be traditionally classified as USS with $\alpha \leq -1.3$ : J0007 $-$ 3040, J2311 $-$ 3359 (discussed above) and J2314 $-$ 3517.

4.1. Projects CX437 and C3377

An observing log for our ATCA observations is presented in Table 3. Observations were carried out simultaneously at 5.5 and 9 GHz using the Compact Array Broadband Backend (CABB; Wilson et al. Reference Wilson2011), with a nominal bandwidth of 2.048 GHz at each frequency comprising $2048 \times 1$ -MHz channels. We used a general observational strategy of target snapshots interleaved with scans of phase calibrators. When more than one target was observed, we ensured that the individual snapshots were sufficiently well spread in hour angle to give sufficient (u, v) coverage for imaging. The standard primary calibrator PKS B1934 $-$ 638 was observed in each run, with the flux density scale reported in Reynolds (Reference Reynolds1994).

Table 3. ATCA observing log for our 5.5- and 9-GHz observations. On 2020 December 2/3, the effective frequency of the upper band was 8.8 rather than 9 GHz. Noise levels are those measured near the target; for the C3377 data, we give the median noise levels at 5.5 and 8.8/9 GHz. Furthermore, for the C3377 data, the reported angular resolutions were determined by taking the separate medians of the major and minor axis FWHMs as well as the BPAs. PKS B1934 $-$ 638 was the primary calibrator in all of the observing runs. Further details on the ATCA observations can be found in Section 4.

In project CX437, we used available Director’s Discretionary Time to observe J2311 $-$ 3359. The array was in the 750C configuration. 20-min target scans were interleaved with 2-min phase calibrator observations. Note that these observations were set up differently to what we describe below for C3377; this was due to the fact that the expected very faint flux densities for J2311 $-$ 3359 at 5.5 and 9 GHz required a significantly longer on-source integration time to increase the likelihood of a detection.

In project C3377, we observed the remaining 48 sources. Observations were carried out with the 6A array configuration for the SGP targets and with the hybrid H168 array configuration for the EQU sources. The H168 configuration was needed for the EQU sources to ensure adequate (u, v) coverage. The individual target scans were mostly 3 min in duration, although with some variations. Phase calibrators were observed every $\sim$ 10–25 min and the scans were either 1 or 1.5 min in duration.

4.2. Data reduction and imaging

The data were reduced and imaged using standard procedures in miriad (Sault, Teuben, & Wright Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995). Radio-frequency interference (RFI) was particularly problematic in our observing run on 2020 December 2/3; hence, a significant amount of primary calibrator data in the upper part of the 9-GHz band had to be flagged. To make the subsequent flux density calibration of the phase calibrator and target data as straightforward as possible, we also flagged the corresponding channels in these data, resulting in a nominal sensitivity penalty of about 10–15% and a shift in the effective frequency from 9 to 8.8 GHz.

When the S/N was sufficient, we used several iterations of imaging and phase-only self-calibration, the latter based on the clean component models. We used multifrequency deconvolution (Sault & Wieringa Reference Sault and Wieringa1994) and the robust weighting parameter (Briggs Reference Briggs1995) was set to 0.5. Note that we did not include baselines to antenna 6 in the imaging and self-calibration step for the data from project CX437 as well as the H168 observations from C3377. This was to remove the large gap in the (u, v) coverage, with a nominal sensitivity penalty of about 20% and poorer angular resolution. However, high-resolution data were already available from FIRST and VLASS for these targets. The maximum antenna spacings for our ATCA data sets were then 192, 750 and 5939 m for the H168, 750C and 6A observations, respectively. A primary beam correction was applied to all cleaned maps, although this did not make a significant difference to the integrated flux densities as all of our targets were at or very close to the pointing centre. Angular resolutions and noise levels can be found in Table 3.

After imaging the data, we then used pybdsf (Mohan & Rafferty Reference Mohan and Rafferty2015) for source finding and integrated flux density measurements. Each flux density uncertainty was determined by combining the fitting uncertainty from pybdsf, the internal calibration uncertainty and the flux density scale uncertainty in quadrature. We estimate that the internal calibration uncertainty is at the $\sim$ 5% level at both 5.5 and 8.8/9 GHz; the flux density scale uncertainty is estimated to be an additional $\sim$ 5% using the information available in Reynolds (Reference Reynolds1994) and Perley & Butler (Reference Perley and Butler2017).

The 5.5- and 8.8/9-GHz flux densities are given in Table 2. In a handful of cases noted in this table as well as in Section 5.1, due to sufficiently low S/N it was necessary to estimate the flux densities from an additional analysis of the images in question, rather than using pybdsf. Overlay plots showing the high-resolution ATCA radio contours from the 6A configuration observations and the medium-resolution contours from the 750C configuration observation of J2311 $-$ 3359 are presented in Section 5.3 (Figure 2).

5.1. Overview of available radio data

Our radio data from GLEAM and the ATCA were supplemented by flux density measurements from other radio surveys (Tables 2 and 4). We used the following catalogues: the 74-MHz VLA Low-Frequency Sky Survey Redux (VLSSr; Cohen et al. Reference Cohen, Lane, Cotton, Kassim, Lazio, Perley, Condon and Erickson2007; Lane et al. Reference Lane, Cotton, van Velzen, Clarke, Kassim, Helmboldt, Lazio and Cohen2014), TGSS at $147.5$ MHz, 325-MHz GMRT survey data of H-ATLAS/GAMA fields from Mauch et al. (Reference Mauch, Klöckner, Rawlings, Jarvis, Hardcastle, Obreschkow, Saikia and Thompson2013), the 365-MHz Texas Survey (TXS; Douglas et al. Reference Douglas, Bash, Bozyan, Torrence and Wolfe1996), the 408-MHz Molonglo Reference Catalogue (MRC; Large et al. Reference Large, Mills, Little, Crawford and Sutton1981; Large, Cram, & Burgess Reference Large, Cram and Burgess1991), the 843-MHz Sydney University Molonglo Sky Survey (SUMSS; Bock, Large, & Sadler Reference Bock, Large and Sadler1999; Mauch et al. Reference Mauch, Murphy, Buttery, Curran, Hunstead, Piestrzynski, Robertson and Sadler2003; Murphy et al. Reference Murphy, Mauch, Green, Hunstead, Piestrzynska, Kels and Sztajer2007), RACS at 887.5 MHz, FIRST at 1.4 GHz, NVSS at 1.4 GHz and VLASS at 3 GHz.Footnote ^g Following De Breuck et al. (Reference De Breuck, van Breugel, Röttgering and Miley2000), we only considered those sources that are well modelled (‘ $+++$ ’ flag) in the TXS catalogue. SUMSS data are only available for SGP sources and FIRST for EQU sources. Furthermore, we used the VLASS catalogue derived from first-epoch data presented and analysed in Gordon et al. (Reference Gordon2020, Reference Gordon2021), applying a multiplicative flux density scale correction of 1/0.87 to the catalogued flux densities and flux density fitting uncertainties (see discussion and analysis in Gordon et al. Reference Gordon2021). For the VLASS data, we also assumed a conservative 10% calibration uncertainty, adding this value in quadrature to the corrected fitting uncertainties from the catalogue to obtain the final uncertainties listed in Table 2.

Table 4. Additional flux density measurements for sources in our sample. We have rescaled the values in the literature so as to place them on the Baars et al. (Reference Baars, Genzel, Pauliny-Toth and Witzel1977) flux density scale; the multiplicative correction factors that we used are stated. VLSSr uncertainties are accurate to two significant figures. See Sections 5.1 and 5.4 for further details.

Notes. ${}^{\rm a}$ The flux density calibration uncertainty was also increased by 5% following the recommendation presented in Section 5.3 in Lane et al. (Reference Lane, Cotton, van Velzen, Clarke, Kassim, Helmboldt, Lazio and Cohen2014). References: Lane et al. (Reference Lane, Cotton, van Velzen, Clarke, Kassim, Helmboldt, Lazio and Cohen2014, VLSSr), Mauch et al. (Reference Mauch, Klöckner, Rawlings, Jarvis, Hardcastle, Obreschkow, Saikia and Thompson2013, GMRT), Douglas et al. (Reference Douglas, Bash, Bozyan, Torrence and Wolfe1996, TXS) and Large et al. (Reference Large, Mills, Little, Crawford and Sutton1981, Reference Large, Cram and Burgess1991, MRC).

The radio data presented in Tables 2 and 4 are from surveys that are not always tied to the same flux density scale. We return to this topic in Section 5.4, as it is an important consideration for the modelling of the broadband radio spectra of our sample. However, for the comparison between the 147.5-MHz TGSS and 151-MHz GLEAM flux densities described below in Section 5.2, we note that we have used rescaled TGSS flux densities from Hurley-Walker (Reference Hurley-Walker2017); these measurements are reported in Table 2. The average multiplicative correction factor applied for the sources in our sample is 0.97, with minimum and maximum values of 0.79 and 1.20, respectively. The rescaling moves the TGSS flux densities onto the flux density scale of Baars et al. (Reference Baars, Genzel, Pauliny-Toth and Witzel1977) that was used to calibrate GLEAM, rather than the Scaife & Heald (Reference Scaife and Heald2012) scale used by Intema et al. (Reference Intema, Jagannathan, Mooley and Frail2017) for the main TGSS catalogue. It also corrects for position-dependent flux density scale variations in TGSS.

There are three sources in the sample with one or more non-detections at the $3\sigma$ level: J0034 $-$ 3112 (843 MHz), J1521 $-$ 0104 (8.8 GHz), and J2311 $-$ 3359 (843 MHz, 1.4 GHz, and 3 GHz). For J0034 $-$ 3112 and J1521 $-$ 0104, we made use of a technique often used when studying the light curves of radio transients (e.g. Swinbank et al. Reference Swinbank2015): rather than including upper limits in the analysis described below, we instead measured the flux density from a forced point-source fit at the target location using imfit in miriad. This ensured the consistency of the characteristics of the data points, rather than a combination of detections and upper limits. The flux densities from these fits are $S_{843} = 9.0 \pm 3.3$ mJy for J0034 $-$ 3112 and $S_{8800} = 0.14 \pm 0.06$ mJy for J1521 $-$ 0104. Both of these flux densities are consistent with the formal $3\sigma$ upper limits ( ${<}10$ and ${<}0.2\ \mathrm{mJy\, beam}^{-1}$ , respectively). For J2311 $-$ 3359, the SUMSS and NVSS upper limits are not significantly constraining, and we therefore excluded them from our analysis. However, for the VLASS data point, we also measured the flux density from a forced fit: $S_{3000} = 0.18 \pm 0.13$ mJy. This measurement is consistent with the formal $3\sigma$ upper limit ( ${<}0.38\ \mathrm{mJy\, beam}^{-1}$ ) to within the $2\sigma$ error on the forced-fit value.

For J0856 $+$ 0223 and J0917 $-$ 0012, these sources were analysed in detail in D20, Drouart et al. (Reference Drouart2021) and Seymour et al. (Reference Seymour2022) over a wider frequency range than is considered here. We report radio properties for these two sources in Tables 2 and 4, but only over the same frequency range as for the other sources in the sample.

5.2. Comparison of radio flux densities from data sets matched or closely spaced in observing frequency

To assess the reliability of the radio flux densities, we compared several of our data sets: FIRST and NVSS, SUMSS and RACS, and GLEAM and TGSS. Moreover, as our sample was selected using the criterion LAS ${\leq}5^{\prime\prime}$ , there was the possibility that one or more very compact components in a given source had resulted in significant variability between the epochs of the surveys used in the three comparisons listed above. Variability resulting from refractive interstellar scintillation (Shapirovskaya Reference Shapirovskaya1978; Rickett, Coles, & Bourgois Reference Rickett, Coles and Bourgois1984; Rickett Reference Rickett1986, Reference Rickett1990) has been observed often at frequencies over the range that we are considering here (e.g. Hunstead Reference Hunstead1972; Gaensler & Hunstead Reference Gaensler and Hunstead2000; Bannister et al. Reference Bannister, Murphy, Gaensler, Hunstead and Chatterjee2011; Ofek & Frail Reference Ofek and Frail2011; Bell et al. Reference Bell2019; Hajela et al. Reference Hajela, Mooley, Intema and Frail2019; Ross et al. Reference Ross2021; Murphy et al. Reference Murphy2021; but see Koay et al. Reference Koay2012 for the case of high-redshift AGN). Other possibilities include intrinsic variability from a compact core and/or jet component (e.g. Mooley et al. Reference Mooley2016; Nyland et al. Reference Nyland2020; Ross et al. Reference Ross2021; Wołowska et al. Reference Wołowska2021) or a combination of both refractive scintillation and intrinsic variability, the former occurring on shorter timescales (e.g. Rickett, Lazio, & Ghigo Reference Rickett, Lazio and Ghigo2006; Bhandari et al. Reference Bhandari2018; Sarbadhicary et al. Reference Sarbadhicary2021).

5.3. Overlay plots

In Figure 2, we present $K_{\rm s}$ -band/radio overlay plots for the sources in our sample. We show contours from our highest-resolution radio data sets, i.e. FIRST, VLASS (angular resolution ${\approx} 2 .\!^{\prime\prime} 5$ ) and the ATCA 5.5- and 9-GHz data from the 6A array configuration observations. The deepest $K_{\rm s}$ -band image available for a given source has been used (i.e. from VIKING, SHARKS or HAWK-I). Note, however, that we do not include J0856 $+$ 0023 and J0917 $-$ 0012, which have been analysed extensively elsewhere (D20; Drouart et al. Reference Drouart2021; Seymour et al. Reference Seymour2022). For the overlay plots, a summary of the $K_{\rm s}$ -band host galaxy magnitudes or lower limits and the lowest radio contour levels used can be found in Table 5.

The host galaxies of J0133 $-$ 3056 and J0842 $-$ 0157 are detected in the HAWK-I $K_{\rm s}$ -band images. Otherwise, any possible host galaxy detections in the VIKING or SHARKS images (e.g. for J1329 $+$ 0133) are tentative at best and below a $5\sigma$ brightest pixel value, which indeed is why the sources are included in our sample (Section 2.2.3).

Table 5. Summary of the $K_{\rm s}$ -band images and lowest radio contour levels used in the overlay plots in Figure 2. The reported radio contour levels are usually $5\sigma$ , but in a handful of cases they are either $3\sigma$ or $4\sigma$ . We also report the $K_{\rm s}$ -band host galaxy magnitudes (limits at the $5\sigma$ level); see Section 3 for further details.

Notes. ${}^{\rm a}$ Not corrected using the 1/0.87 scaling factor discussed in Section 5.1.

${}^{\rm b} 3\sigma$ contour.

${}^{\rm c} 4\sigma$ contour.

Given our LAS ${\leq}5 ^{\prime\prime}$ selection criterion and the angular resolutions of the VLASS and ATCA data (which are similar for the SGP sources), it is not particularly surprising that all but one of our sources have radio morphologies that can be classified as one of the following: single component, incipient double or resolved double. The one exception is the triple source J0309 $-$ 3526. Offsets between the VLASS and ATCA centroids are within the astrometric uncertainties of the VLASS quick-look images (up to ${\approx} 1^{\prime\prime}$ ; Lacy et al. Reference Lacy2019). Phase errors remaining in the VLASS maps can lead to erroneous extension visible in the overlay plots (e.g. for J0301 $-$ 3132 and J0326 $-$ 3013), and care must be taken interpreting the radio morphology at 3 GHz (see Lacy et al. Reference Lacy2019 for further details).

Further information can be found in the notes on individual sources in Section 5.5.

5.4. Radio spectral modelling

The radio data for each source in our sample spans a maximum frequency range of 74 MHz–9 GHz, with flux density measurements at up to 29 different frequencies. We were therefore able to explore the broadband radio spectral properties of our sample, which could then be compared with the GLEAM-only $\alpha$ – $\beta$ fitting that was used as part of the sample selection (Section 2.2.2). For each source for which we present an overlay plot in Figure 2, we have also plotted the corresponding observed-frame broadband radio spectrum in the right-hand column of this figure. Overlaid on each spectrum is the preferred model: either a single or double power law. Radio spectra of particular interest are discussed in the notes on individual sources in Section 5.5, followed by further discussion in Sections 6.1 and 6.2.

For J0856 $+$ 0223 and J0917 $-$ 0012, radio spectra over a wider frequency range were presented in D20, Drouart et al. (Reference Drouart2021) and Seymour et al. (Reference Seymour2022). We do not show spectra for these sources in Figure 2, but instead analyse these sources later in this paper in Section 6.2 and Figure 3.

We now describe the steps taken to construct the broadband radio spectra and carry out the spectral modelling.

5.4.2. Modelling the data

As in D20, we fitted a single, double and triple power law to each spectrum. We made use of the emcee (Foreman-Mackey et al. Reference Foreman-Mackey2013a, b) and george (Ambikasaran et al. Reference Ambikasaran, Foreman-Mackey, Greengard, Hogg and O’Neil2015; Foreman-Mackey Reference Foreman-Mackey2015) modules in python, so as to carry out a Markov Chain Monte Carlo analysis of the parameter space. A single-power-law fit was defined as follows:

(4) \begin{equation}S_{\nu} = N\!\left(\frac{\nu}{\nu_{0}}\right)^{\alpha},\end{equation}

where N is a constant, $\alpha$ the spectral index across the observed-frequency range, and $\nu_{0}$ the reference frequency, which we chose to be 1000 MHz. Thus, N is the fitted flux density at 1000 MHz for the single-power-law fit. For the smoothly varying double-power-law fit,

(5) \begin{equation}S_{\nu} = N\!\left(\frac{\nu}{\nu_{0}}\right)^{\alpha_{\rm l}} \times \left[1 + \left(\frac{\nu}{\nu_{\rm b}}\right)^{\lvert\alpha_{\rm h} - \alpha_{\rm l}\rvert}\right]^{\mathrm{sgn}{(\alpha_{\rm h} - \alpha_{\rm l})}},\end{equation}

where $\alpha_{\rm l}$ and $\alpha_{\rm h}$ are the spectral indices either side of the break frequency $\nu_{\rm b}$ and $\mathrm{sgn}{}$ the signum function. Lastly, the smoothly varying triple-power-law fit was of the form

(6) \begin{align}S_{\nu} =\, & N\!\left(\frac{\nu}{\nu_{0}}\right)^{\alpha_{\rm l}} \times \left[1 + \left(\frac{\nu}{\nu_{\rm bl}}\right)^{\lvert\alpha_{\rm m} - \alpha_{\rm l}\rvert}\right]^{\mathrm{sgn}{(\alpha_{\rm m} - \alpha_{\rm l})}} \nonumber\\[4pt]& \times \left[1 + \left(\frac{\nu}{\nu_{\rm bh}}\right)^{\lvert\alpha_{\rm h} - \alpha_{\rm m}\rvert}\right]^{\mathrm{sgn}{(\alpha_{\rm h} - \alpha_{\rm m})}},\end{align}

where $\alpha_{\rm l}$ and $\alpha_{\rm m}$ are the spectral indices either side of the lower break frequency $\nu_{\rm bl}$ and similarly for the spectral indices $\alpha_{\rm m}$ and $\alpha_{\rm h}$ as well as the higher break frequency $\nu_{\rm bh}$ .

Table 6. Non-informative priors for each of the parameters in Equations (4)–(6). Section 5.4 describes how we modelled the broadband radio spectra.

The data were fitted with input units of MHz and mJy. We used non-informative priors for each of the parameters in Equations (4)–(6); the priors are listed in Table 6. The priors were chosen such that we assumed a triple-power-law fit with a low-frequency turnover (that could result from synchrotron self-absorption and/or free–free absorption). For both the break frequency in the double-power-law fit and the higher break frequency in the triple-power-law fit, the prior was sufficiently general such that we could model spectral steepening at higher frequencies due to one or more energy loss mechanisms, or high-frequency spectral flattening due to a radio core component beginning to dominate over the lobe emission.

When fitting the data, the uncertainty for each GLEAM flux density was calculated by combining the fitting and absolute calibration uncertainties, the latter being 8% for all of the sources in our sample (Hurley-Walker et al. Reference Hurley-Walker2017; Franzen et al. Reference Franzen, Hurley-Walker, White, Hancock, Seymour, Kapińska, Staveley-Smith and Wayth2021). Furthermore, it was necessary to take into account the correlations that exist between the GLEAM sub-bands, so as to avoid erroneous fits (see discussion in Hurley-Walker et al. Reference Hurley-Walker2017). To do so, we used a blocked Matérn covariance function for the GLEAM flux densities (Rasmussen & Williams Reference Rasmussen and Williams2006).

The preferred model for each radio spectrum was determined using the sample-size-corrected Akaike Information Criterion (AICc; Akaike Reference Akaike1974; Burnham & Anderson Reference Burnham and Anderson2002). In particular,

(7) \begin{equation}{\rm AICc} = \chi^2 + 2k + \frac{ 2k(k+1)}{n - k - 1},\end{equation}

where n and k are the number of data points and free parameters, respectively, and $\chi^2$ is the standard goodness-of-fit statistic. For all possible realisations from the emcee fitting, we determined the minimum value of AICc, $\min$ (AICc), for each of the possible three models. We then selected the preferred model by using the standard convention of examining the difference

(8) \begin{align}\Delta{\rm AICc} & = \min\!({\rm AICc})_{i}\nonumber\\&\quad - \min\{\min(\rm{AICc})_{\rm SPL}, \min\!({\rm AICc})_{\rm DPL}, \min\!({\rm AICc})_{\rm TPL}\},\end{align}

where the subscripts SPL, DPL and TPL refer to the single-power-law, double-power-law, and triple-power-law fits, respectively, and the ith model can be one of the three possibilities. Our preferred model first satisfied the condition $0 \leq \Delta{\rm AICc} < 4$ and secondly was the model with the fewest free parameters satisfying this condition. Note that this can mean, for example, that a double-power-law fit has the lowest AICc value, but the single power law was selected as the preferred fit because $\Delta{\rm AICc}$ is sufficiently small enough.

The fitted model parameters are presented in Table 7. For 34 sources, the preferred model is a double power law, with the remaining 17 sources described by a single power law. The triple power law is not the preferred model for any of the sources (but see Drouart et al. Reference Drouart2021 and Seymour et al. Reference Seymour2022, where the radio emission from J0917 $-$ 0012 was modelled with a triple power law).

Table 7. Fitted parameters for the broadband radio spectra. We give the preferred model type: SPL for a single power law and DPL for a double power law. The 16th, 50th and 84th percentiles are reported for the marginalised parameter distributions. We also list the parameter values corresponding to the best fit with the smallest value of AICc as well as the reduced chi-squared goodness-of-fit statistic for this model. Further information on the modelling can be found in Section 5.4.

5.5. Notes on individual sources

In this section, we discuss the overlay plots and/or radio spectra of a number of sources in the sample.

J0007 $-$ 3040: This source is of particular interest, with $K_{\rm s} > 23.1$ from SHARKS, a relatively compact radio morphology with LAS $=3.\!^{\prime\prime}0$ in both the ATCA and VLASS images, and a double-power-law spectrum with best-fitting spectral indices of $\alpha_{\rm l} = -1.385$ and $\alpha_{\rm h} = -2.75$ . While $\alpha_{\rm l}$ indicates a USS spectrum at frequencies below $\nu_{\rm b} = 1.97$ GHz, $\alpha_{\rm h}$ is exceptionally steep and relatively well constrained. This source therefore appears to be a promising HzRG candidate; possible scenarios for explaining the properties of the radio spectrum are discussed in Section 6.1. Another possibility is that this source may be an as of yet undetected pulsar, but this seems less likely given the LAS and the fact that there are hints of an incipient double radio morphology for this source.

We also note that we checked whether emission potentially being resolved out in the ATCA maps could explain at least some of the observed spectral curvature. Adjusting the robust weighting parameter at 5.5 GHz to $-$ 2 (i.e. close to uniform weighting) gave an image with a very similar resolution to the 9-GHz map generated with a robust weighting parameter of 0.5. The measured 5.5-GHz flux density was not significantly different to the value reported in Table 2. Additionally, a 5.5-GHz map generated with a robust weighting parameter of 2 (i.e. close to natural weighting) also gave a very similar flux density to the value in Table 2.

J0133 $-$ 3056: This source is an incipient double in VLASS with LAS $= 4.\!^{\prime\prime} 3$ . The likely host galaxy is seen in the HAWK-I $K_{\rm s}$ -band image close to the peak of the radio emission, with magnitude $K_{\rm s} = 22.23 \pm 0.08$ .

J0309 $-$ 3526: The ATCA 5.5-GHz and VLASS data both suggest a multi-component morphology, perhaps a triple source, although as can be seen in the corresponding panel in Figure 2, the morphologies are not fully consistent between the two frequencies. There is clear extension at 5.5 GHz in the direction of the synthesised beam, almost orthogonal to what appears to be the main axis of the radio emission. Given the VLASS morphology and the fact that we could not phase self-calibrate the 5.5-GHz data due to low S/N, it seems most plausible that the extension at 5.5 GHz in the direction of the synthesised beam is spurious. The spectral index between 5.5 and 9 GHz is extremely steep: $\alpha^{9000}_{5500} = -2.8 \pm 1.1$ . There is evidence that the accuracy of $\alpha^{9000}_{5500}$ is affected by diffuse emission being resolved out; this can be seen, for example, by decreasing the robust weighting parameter from 0.5 to $-$ 2 (i.e. uniform weighting in the latter case) and reimaging the 5.5-GHz data. Higher S/N data at mid/high frequencies in array configurations with sufficient low-surface-brightness sensitivity are needed for this source.

J0842 $-$ 0157: The FIRST and VLASS morphologies are very compact: LAS $=1.\!^{\prime\prime}7$ and $1.\!^{\prime\prime}1$ , respectively. We take the HAWK-I $K_{\rm s}$ -band source closest to the radio centroid as the host galaxy identification; the magnitude of the host galaxy is $K_{\rm s} = 22.96 \pm 0.12$ .

J0909 $-$ 0154, J1141 $-$ 0158 and J1443 $+$ 0229: The best-fitting model for J0909 $-$ 0154 has significantly flattened at the bottom end of the GLEAM band and would be expected to begin to turn over at frequencies $\lesssim 70$ MHz. J1141 $-$ 0158 and J1443 $+$ 0229 are the only sources in the sample where the best fit peaks and turns over (within the GLEAM band at 205 and 91 MHz, respectively); higher S/N data at the lower GLEAM frequencies are needed to confirm the turnover in each case, however.

J1112 $+$ 0056: There is a second NVSS source $51.\!^{\prime\prime} 3$ from the GLEAM position (beyond the region shown in the panel for this source in Figure 2) that has a 1.4-GHz flux density that is about 40% of the NVSS flux density of the HzRG candidate (NVSS J111212 $+$ 005519 with $S_{1400} = 6.7 \pm 0.5$ mJy; cf. $S_{1400} = 17.1 \pm 0.7$ mJy for the HzRG candidate). Both of these sources have very similar two-point spectral indices $\alpha^{1400}_{887.5}$ ; if this spectral similarity is also the case at low frequencies, the accuracy of the GLEAM $\alpha$ – $\beta$ fitting should not be affected significantly. The GLEAM/TGSS flux density ratio in Table 2 is $1.27 \pm 0.15$ ; this tentatively suggests that there could be an excess in GLEAM due to source blending, albeit not statistically significant at e.g. the $3\sigma$ level.

J1125 $-$ 0342 and J1317 $+$ 0339: These are the USS-selected sources TN J1125 $-$ 0342 and TN J1317 $+$ 0339, respectively, from De Breuck et al. (Reference De Breuck, van Breugel, Röttgering and Miley2000). Their redshifts remain unknown and the VIKING non-detections are the deepest constraints on their $K_{\rm s}$ -band magnitudes: $>21.4$ (J1125 $-$ 0342) and $> 21.3$ (J1317 $+$ 0339).

For J1125 $-$ 0342, two sources separated by $8.\!^{\prime\prime} 3$ are visible in the VLASS Epoch 1 image. The 3-GHz flux density of the southern source is ${\approx} 1.4$ mJy, a factor of $\sim$ 20 fainter than the much brighter source to the north. J1125 $-$ 0342 could therefore be a very asymmetric double, or the fainter source to the south is unrelated. In the case of the former scenario, the angular size would then suggest that the source is too extended to be at a very high redshift. For the purpose of analysis in this paper, we have assumed the latter scenario. Deeper $K_{\rm s}$ -band imaging and higher-resolution radio imaging at e.g. 5.5 and 9 GHz are needed.

For J1317 $+$ 0339, it can be seen in Figure 2 that the 365-MHz flux density point is a clear outlier; we did not include this data point when fitting the radio spectrum. The $\alpha_{\rm l}$ value in Table 7 implies that the spectrum is not as steep as suggested by the two-point spectral index between the TXS and NVSS surveys. This is the case for J1125 $-$ 0342 as well, although not to the same extent as for J1317 $+$ 0339. We included the 365-MHz data point when fitting the radio spectrum of J1125 $-$ 0342.

J1329 $+$ 0133: There is a hint of a host galaxy identification in the VIKING $K_{\rm s}$ -band image, but the brightest pixel value is only at the $3.5\sigma$ level (which we do not consider to be a secure detection).

J1335 $+$ 0112: As previously discussed in Section 2.3, this source has a detection in AllWISE: the W1-band magnitude is $20.033 \pm 0.130$ (Cutri et al. Reference Cutri2014). The source is not detected in any of the other three longer-wavelength AllWISE bands. The $K_{\rm s}$ -band magnitude is $> 21.2$ (Table 5); the ${\gtrsim} 1.2$ mag break between 2.15 and $3.37\, \unicode{x03BC} \mathrm{m}$ might potentially indicate a redshifted 4000 Å break (i.e. $z \gtrsim 4.4$ ). The AllWISE W2-band magnitude is $> 19.064$ ( $5\sigma$ ). Additional follow-up and analysis is needed.

J1340 $+$ 0009: There is a hint of a host galaxy identification in the VIKING H-band image, but the brightest pixel value is only at the $4.2\sigma$ level (which we do not consider to be a secure detection).

J1351 $-$ 0209: This is the brightest radio source in our sample (by a factor of about three at 151 MHz) with $S_{151} = 2.44$ Jy. It is also catalogued as the Parkes source PKS B1349 $-$ 019 (Wright & Otrupcek Reference Wright and Otrupcek1990). The PKS 2.7- and 5-GHz flux densities from Wright & Otrupcek (Reference Wright and Otrupcek1990) are 130 and 50 mJy, respectively, consistent with the VLASS and ATCA 5.5-GHz data. Given this consistency at very similar frequencies as well as the fact that uncertainties are not reported for these Parkes measurements, we chose not to use the Parkes flux densities in the radio spectrum modelling.

Additionally, Downes et al. (Reference Downes, Peacock, Savage and Carrie1986) reported 1.5- and 4.9-GHz flux densities from their high-resolution VLA data (250 and 60 mJy, respectively); the former is consistent with the NVSS flux density and the latter with both the Parkes 5-GHz and ATCA 5.5-GHz measurements. Again, there is no significant advantage in including these data in our radio spectrum modelling, particularly as it is unclear what the flux density uncertainties are in these cases as well. The radio morphology in the 4.9-GHz VLA map shows two components; the angular extent ( ${\sim} 2.\!^{\prime\prime}2$ ; position angle ${\sim}{-}45^\circ$ ) is reasonably consistent with the FIRST and VLASS LAS measurement ( $2.\!^{\prime\prime}8$ ).

Dunlop et al. (Reference Dunlop, Peacock, Savage, Lilly, Heasley and Simon1989) did not detect the host galaxy in optical imaging; the reported B- and R-band magnitude limits are ${\gtrsim} 24$ . The source is listed in Dunlop & Peacock (Reference Dunlop and Peacock1990) as a candidate high-redshift object, but with $K \approx 19.75$ ; this is most likely a misidentification or erroneous catalogue entry given that our overlay plot shows no evidence of a $K_{\rm s}$ -band identification to a much greater depth ( $ K_{\rm s} > 21.2$ ).

J1521 $-$ 0104: There is a hint of extension in RACS about $20^{\prime\prime}$ to the north-west (beyond the region shown in the panel for this source in Figure 2), albeit at a marginal level ( ${\sim} 3.9\sigma$ ). The extension is coincident with the near-infrared source 2MASS J15215337 $-$ 0104034, where $K_{\rm s} = 14.449 \pm 0.033$ (Cutri et al. Reference Cutri2003). This 2MASS source is $22.\!^{\prime\prime}2$ from our HzRG candidate. The possible radio/near-infrared association might then suggest that the HzRG candidate shown in Figure 2 is part of a larger, head–tail source. On the other hand, there is no evidence of similar extension in our other radio data sets. Given that the extended 887.5-MHz radio emission and in turn the radio/near-infrared association is tentative, we regard the source shown in Figure 2 as an HzRG candidate with the requisite compact radio morphology. A deeper $K_{\rm s}$ -band image would allow us to determine if a host galaxy is associated with this radio source.

Thyagarajan et al. (Reference Thyagarajan, Helfand, White and Becker2011) classified J1521 $-$ 0104 as variable at 1.4 GHz based on analysis of the three separate snapshot observations taken for this source as part of FIRST. These authors reported a minimum variability timescale of 8 d. The catalogued FIRST and NVSS flux densities in Table 2 are consistent on a longer timescale (2.1 yr), however.Footnote ⁱ Short-timescale variability from a scintillating radio core would rule out that J1521 $-$ 0104 is a component in a larger radio source. Alternatively, the variability may have been due to a scintillating lobe hotspot; this alone would not provide conclusive evidence of the true angular extent of this radio source.

J2219 $-$ 3312: While this source passed step 8 in Table 1, there is an AllWISE source $1.\!^{\prime\prime}2$ from the ATCA position. The ATCA centroid is slightly further to the south-south-east than in TGSS (offset $1.\!^{\prime\prime}3$ ), and there is also evidence in Figure 2 of a small offset between the ATCA and VLASS centroids. The AllWISE source has W1 and W2 magnitudes of $20.089 \pm 0.161$ and $20.188 \pm 0.347$ , respectively (Cutri et al. Reference Cutri2014), but is not detected in the longer-wavelength AllWISE bands. The $K_{\rm s}$ -band limit is $> 22.3$ (Table 5); the ${\gtrsim} 2.2$ mag break between $K_{\rm s}$ -band and W1 is larger than for J1335 $+$ 0112 discussed above. We regard the AllWISE source as a tentative host galaxy identification; this association needs to be confirmed with follow-up work, particularly deep $K_{\rm s}$ -band imaging.

J2311 $-$ 3359: This source, not selected from our $\alpha$ – $\beta$ criteria (as discussed in Section 2.3), has an extremely steep spectrum well described by a single power law with $\alpha = -1.842$ . Therefore, J2311 $-$ 3359 is faint in our higher-frequency images. The source is unresolved in the ATCA data, but the low S/N does not allow an accurate LAS determination. This source is also unresolved in RACS, but we measured an LAS of $5.\!^{\prime\prime}6$ in an 887.5-MHz ASKAP early science image of the GAMA-23 field (see Seymour et al. Reference Seymour2020). Such an LAS value falls outside of our LAS ${\leq}5^{\prime\prime}$ criterion, but needs to be confirmed with additional radio data. Note in Figure 2 that the VLASS contours are offset from the ATCA contours, but the former are from a line-like artefact in the map.

J2330 $-$ 3237: We interpret this source as an asymmetric double with LAS $4.\!^{\prime\prime}2$ , where the components have significantly different flux densities (the western lobe being much fainter). Using the available high-resolution radio data, the flux densities of the brighter eastern lobe are $6.3 \pm 0.7$ , $3.74 \pm 0.27$ and $2.11 \pm 0.16$ mJy at 3, 5.5 and 9 GHz, respectively. Similarly, the values for the significantly fainter component to the west are ${<}0.45$ , $0.22 \pm 0.06$ and ${<}0.078$ mJy, respectively (upper limits at the $3\sigma$ level). Given that the western lobe is detected at 5.5 GHz only, we verified that this source was not spuriously created as a result of our phase-only self-calibration step. The two-point spectral indices are $\alpha^{5500}_{3000} = -0.86 \pm 0.22$ and $\alpha^{9000}_{5500} = -1.16 \pm 0.21$ for the eastern lobe; similarly these values are $\alpha^{5500}_{3000} \gtrsim -1.2 $ and $\alpha^{9000}_{5500} \lesssim -2.1$ for the western lobe. The significant spectral steepening of the western lobe at higher frequencies could at least be partly due to flux density possibly being resolved out in the higher-resolution 9-GHz map. Further evidence in favour of a double-lobed morphology is the hint of a $K_{\rm s}$ -band identification between the two components (brightest pixel value $= 3.7\sigma$ ); there is similar marginal evidence in H-band (brightest pixel value $= 4.4\sigma$ ).

J2340 $-$ 3230: Two radio sources are visible, with the north-eastern source having a possible marginal (brightest pixel value $= 4.6\sigma$ ) $K_{\rm s}$ -band detection in SHARKS. The radio flux densities of this source are $\sim$ 0.75, $0.79 \pm 0.11$ and $0.58 \pm 0.07$ mJy at 3, 5.5 and 9 GHz, respectively. Therefore, there is tentative evidence that this source is turning over at GHz frequencies. For the purpose of analysis in this paper, we have assumed that the HzRG candidate is the south-western source (which is e.g. an order of magnitude brighter at 5.5 GHz) and that the north-eastern source is an unrelated source that is nearby in projection. Alternatively, this could be a lower-redshift source with a one-sided jet (assuming that one of the two sources is the radio core), where the LAS is $6.\!^{\prime\prime}2$ . A deeper $K_{\rm s}$ -band image is needed to identify a possible near-infrared counterpart coincident with the south-western source; there is tentative evidence of a SHARKS $K_{\rm s}$ -band detection (brightest pixel value $= 3.8\sigma$ ).