The GLEAMing of the first supermassive black holes: III. Radio sources with ultra-faint host galaxies

Abstract

We present deep near-infrared $K_\textrm{s}$-band imaging for 35 of the 53 sources from the high-redshift ($z \gt 2$) radio galaxy candidate sample defined in Broderick et al. (2022, PASA, 39, e061). These images were obtained using the High-Acuity Widefield K-band Imager (HAWK-I) on the Very Large Telescope. Host galaxies are detected for 27 of the sources, with $K_\textrm{s} \approx 21.6$–23.0 mag (2$^{\prime\prime}$ diameter apertures; AB). The remaining eight targets are not detected to a median $3\unicode{x03C3}$ depth of $K_\textrm{s} \approx 23.3$ mag (2$^{\prime\prime}$ diameter apertures). We examine the radio and near-infrared flux densities of the 35 sources, comparing them to the known $z \gt 3$ powerful radio galaxies with 500-MHz radio luminosities $L_{500\,\textrm{MHz}} \gt 10^{27}$ W Hz$^{-1}$. By plotting 150-MHz flux density versus $K_\textrm{s}$-band flux density, we find that, similar to the sources from the literature, these new targets have large radio to near-infrared flux density ratios, but extending the distribution to fainter flux densities. Five of the eight HAWK-I deep non-detections have a median $3\unicode{x03C3}$ lower limit of $K_\textrm{s} \gtrsim 23.8$ mag (1$.\!^{\prime\prime}$5 diameter apertures); these five targets, along with a further source from Broderick et al. (2022, PASA, 39, e061) with a deep non-detection ($K_\textrm{s} \gtrsim 23.7$ mag; $3\unicode{x03C3}$; 2$^{\prime\prime}$ diameter aperture) in the Southern H-ATLAS Regions $K_\textrm{s}$-band Survey, are considered candidates to be ultra-high-redshift ($z \gt 5$) radio galaxies. The extreme radio to near-infrared flux density ratios ($\gt 10^5$) for these six sources are comparable to TN J0924$-$2201, GLEAM J0856$+$0223 and TGSS J1530$+$1049, the three known powerful radio galaxies at $z \gt 5$. For a selection of galaxy templates with different stellar masses, we show that $z \gtrsim 4.2$ is a plausible scenario for our ultra-high-redshift candidates if the stellar mass $M_\textrm{*} \gtrsim 10^{10.5}$ M$_\odot$. In general, the 35 targets studied have properties consistent with the previously known class of infrared-faint radio sources. We also discuss the prospects for finding more UHzRG candidates from wide and deep near-infrared surveys.

1. Introduction

This paper is the third in a series describing an ongoing programme to find high-redshift radio galaxies with powerful radio emission (HzRGs; redshift $z \gt 2$ and 500-MHz radio luminosity $L_{500\,\textrm{MHz}} \gt 10^{27}$ W Hz $^{-1}$ ; see review by Miley & De Breuck 2008). In particular, to find new candidates, we have cross-matched the wideband 72–231 MHz GaLactic and Extragalactic All-sky Murchison Widefield Array (MWA; Tingay et al. 2013) survey (GLEAM; Wayth et al. 2015) with the Visible and Infrared Survey Telescope for Astronomy (VISTA; Dalton et al. 2006; Emerson, McPherson, & Sutherland 2006) Kilo-degree Infrared Galaxy survey (VIKING; Edge et al. 2013). In Drouart et al. (2020), we described a pilot study of four HzRG candidates, introducing a novel selection technique that considers both radio spectral steepness and curvature within the GLEAM band. This study led to the discovery of GLEAM J0856 $+$ 0223, a powerful radio galaxy at $z=5.55$ , as well as another potentially ultra-high-redshift source that was further investigated in two companion studies (GLEAM J0917+0012; Drouart et al. 2021; Seymour et al. 2022). Then, in Broderick et al. (2022), henceforth B22, we defined a larger GLEAM–VIKING sample of 53 HzRG candidates (including GLEAM J0856 $+$ 0223 and GLEAM J0917 $+$ 0012).

Our primary aim is that the B22 sample will provide high-quality candidate powerful radio galaxies at $z \gt 5$ ; we will henceforth refer to such sources as ultra-high-redshift radio galaxies (UHzRGs). Finding UHzRGs, particularly those within the Epoch of Reionisation (EoR; $z \gtrsim 6.5$ ), would provide further vital observational constraints for theoretical explanations of the formation and rapid growth of supermassive black holes (SMBHs) in the early Universe with mass $M_\textrm{BH} \sim 10^{8-9}$ M $_\odot$ (e.g. Volonteri 2012; Johnson et al. 2013; Latif & Ferrara 2016; Smith, Bromm, & Loeb 2017; Wise et al. 2019; Smith & Bromm 2019; Inayoshi, Visbal, & Haiman 2020; Kroupa et al. 2020; Di Matteo, Angles-Alcazar, & Shankar 2023). Identifying powerful radio galaxies within the EoR would also provide bright background radio sources to facilitate searches for redshifted 21-cm hydrogen absorption by the neutral intergalactic medium, including the 21-cm forest from intervening material along the line of sight (e.g. Carilli et al. 2004; Mack & Wyithe 2012; Ciardi et al. 2015), and also allow vital constraints to be placed on active galactic nucleus (AGN) feedback (e.g. review by Hardcastle & Croston 2020) at early cosmic times.

While an increasing number of radio-loud quasars are now known at $z \gt 5$ (e.g. Bañados et al. 2018, 2021, 2023; Ighina et al. 2021, 2022, 2023, 2024a,b; Momjian et al. 2021; Gloudemans et al. 2022; Belladitta et al. 2023), all but one of these sources are not radio powerful according to our aforementioned luminosity criterion. The exception is PSO J352.4034 $-$ 15.3373 at $z=5.84$ (Bañados et al. 2018). ^a These radio-loud quasars are also not obscured (in the rest-frame ultraviolet and optical). The discovery of the obscured AGN COS-87259 at $z=6.85$ (Endsley et al. 2022, 2023b) was therefore a particularly noteworthy development. This source has $M_\textrm{BH} \gtrsim 1.5 \times 10^9$  M $_\odot$ (Endsley et al. 2023b) and is the only confirmed obscured and radio-loud (yet weakly radio emitting: $L_{500\,\textrm{MHz}} \approx 6.9 \times 10^{25}$ W Hz $^{-1}$ ) AGN currently known at $z \gt 6$ . Additionally, the source COSW-106725 is a candidate $z \sim 7.7$ obscured AGN with $M_\textrm{BH} \geq 6.4 \times 10^8$ M $_\odot$ (Lambrides et al. 2024). If the redshift is confirmed, and depending on the extrapolation of the radio spectrum below 100 MHz, this source would have $L_{500\,\textrm{MHz}} \approx 2.3 \times 10^{27}$ W Hz $^{-1}$ , that is, it would be powerful in the radio.

Gilli et al. (2022) suggested that 80–90% of SMBHs at $z = 6$ –8 are obscured, mainly by the dense interstellar medium (ISM) of the host galaxy. Moreover, Johnson & Upton Sanderbeck (2022) presented a model where sustained super-Eddington growth remains hidden due to the extreme infall rate of material, preventing any radiation from escaping. These authors demonstrated that the obscured super-Eddington accretion rates can be maintained even for non-spherical accretion (i.e. more disc-like accretion), which would allow a pathway for radio jets to emerge. Therefore, a compelling picture is beginning to emerge in which radio selection, relatively unaffected by a dense ISM in the host galaxy, may play an especially vital role in finding AGN in the early Universe. Indeed, using wide-area radio surveys such as GLEAM to find obscured AGN provides a unique way to probe this population at high redshift. We note that if the radio lobes are sufficiently compact and within the host galaxy, there could be free–free absorption by the ISM, particularly at lower frequencies (e.g. Bicknell, Dopita, & O’Dea 1997). By selection, the B22 sample will have low-frequency turnovers, which could be due to free–free absorption. Hence, it is of particular interest whether the B22 sample contains brighter radio analogues of COS-87259.

It is well known that there is a tight correlation between near-infrared K-band (2.2 $\unicode{x03BC}$ m) magnitude and redshift for powerful radio galaxies, which are among the most massive galaxies at a given redshift (the K–z relation; e.g. Lilly & Longair 1984; Eales et al. 1997; van Breugel et al. 1998; Willott et al. 2003; Rocca-Volmerange et al. 2004). This K-band light from the host galaxy traces the old stellar population at low redshift, although it gradually probes the optical and then the ultraviolet with increasing redshift. Given the K–z relation, an HzRG search can be made significantly more efficient by selecting only those sources in large radio catalogues with K-band magnitudes fainter than a given cutoff (e.g. Ker et al. 2012). Dedicated K-band follow-up of HzRG candidate samples has been commonplace in the literature (e.g. Chambers et al. 1996; van Breugel et al. 1999; Jarvis et al. 2001, 2004; De Breuck et al. 2002, 2004; Cruz et al. 2006; Bryant et al. 2009; Saxena et al. 2018, 2019).

The B22 HzRG candidate sample was partly selected on the basis of VIKING $K_\textrm{s}$ -band (2.15 $\unicode{x03BC}$ m) host galaxy non-detections: $K_\textrm{s} \gt 21.2$ mag ( $5\unicode{x03C3}$ point-source limits presented in Table 5 in B22; median seeing full width at half maximum (FWHM) 0 $.\!^{\prime\prime}$ 83; AB). As was described in B22, three of the 53 sources in our sample have deeper limits from the Southern H-ATLAS ^b Regions $K_\textrm{s}$ -band Survey (SHARKS; Dannerbauer et al. 2022). A further four sources have $K_\textrm{s}$ -band host galaxy detections from Very Large Telescope (VLT; European Southern Observatory 1998) observations using the High-Acuity Widefield K-band Imager (HAWK-I; Kissler-Patig et al. 2008), including both GLEAM J0856 $+$ 0223 and GLEAM J0917 $+$ 0012 from the Drouart et al. (2020) pilot study. In this paper, we significantly increase the number of sources in our sample with deep $K_\textrm{s}$ -band imaging by presenting the results from a HAWK-I $K_\textrm{s}$ -band observing campaign for an additional 35 targets. Such deep HAWK-I imaging allows us to further narrow down the best UHzRG candidates compared to the VIKING selection alone; for this particular observing campaign, the improvement in sensitivity is $\sim$ 1–2 mag at an equivalent detection threshold. Hence, we can significantly increase the efficiency of spectroscopic follow-up, which can be an expensive and often challenging process, particularly in terms of observing time given that (U)HzRGs have a very low surface density (e.g. Miley & De Breuck 2008).

The layout of this paper is as follows. In Section 2, we describe the HAWK-I observations as well as our method for detecting the host galaxies and measuring their magnitudes. We also outline how a compilation of radio and near-infrared flux densities was obtained for both the B22 sample and (U)HzRGs from the literature. We present our results in Section 3, including the host galaxy magnitudes as well as the $K_\textrm{s}$ -band–radio overlay plots for our best UHzRG targets. A discussion follows in Section 4, in particular focusing on (i) how our data set compares to the literature in radio–near-infrared flux density space, and (ii) plausible scenarios for our best UHzRG candidates. We then present our conclusions in Section 5. Lastly, Appendix 1 contains the remaining overlay plots as well as notes on individual sources.

In this paper, we use the following conventions. Uncertainties are given as $\pm 1\unicode{x03C3}$ . Magnitudes are given in the AB system (Oke 1974); to convert between Vega and AB magnitudes, we use the relation $K_\textrm{s, AB} = K_\textrm{s, Vega} + 1.85$ as given in, for example, Blanton & Roweis (2007). We also assume a flat, Lambda cold dark matter ( $\Lambda$ CDM) cosmology with Hubble constant $H_0=67.7$ km s $^{-1}$ Mpc $^{-1}$ , matter density parameter $\Omega_\textrm{M}=0.31$ , and vacuum density parameter $\Omega_{\Lambda}=0.69$ (Planck Collaboration et al. 2020). Additionally, $\log$ refers to the decimal logarithm (base 10).

2. Observations, data reduction and analysis

2.1 HAWK-I observing campaign and data products

Our HAWK-I observations of the B22 sample (49 targets without previous HAWK-I observations) took place from 2021 October 20 to 2022 March 30 in service mode (ESO period P108; programme ID 108.22HY.001). HAWK-I has a field of view (FoV) of $7\hbox{$.\mkern-4mu^\prime$}$ $5 \times 7\hbox{$.\mkern-4mu^\prime$}5$ and an average pixel size of 0 $.\!^{\prime\prime}$ 1066. There were some very small differences in pixel size from image to image, but only at the level of the fourth decimal place. For 33 sources, observations comprised 12 sets of $17 \times 10$ s jittered exposures, where the jitter box width was set to 30 $^{\prime\prime}$ , that is, 34 min on-source per target. Shorter observations were possible for two sources: J0129 $-$ 3109 (680 s) and J0216–3301 (510 s). To improve the seeing, we used the GRound layer Adaptive optics Assisted by Lasers (GRAAL); the angular resolution ranged from 0 $.\!^{\prime\prime}$ 33–0 $.\!^{\prime\prime}$ 78 FWHM, with a median value of 0 $.\!^{\prime\prime}$ 46. While our programme was not fully completed by the end of P108, we obtained science-quality images for 71% of the targets (35/49 sources), with most observations having a quality control grade ‘A’ and a handful ‘B’. Reduced data were mostly obtained from the automatic HAWK-I pipeline, ^c although for J0042 $-$ 3515, J0129 $-$ 3109, and J0216–3301 some manual data reduction steps were necessary.

2.2 $\textbf{K}_{\mathbf{s}}$ -band magnitude calibration

An initial magnitude calibration was provided by the automatic HAWK-I pipeline. Then, as in B22, we used the Two-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) point-source catalogue (Cutri et al. 2003) to refine the magnitude calibration of the images. We used a sample of bright stars ( $K_\textrm{s} \lt 16.9$ mag) within the HAWK-I FoV to calculate the inverse-variance-weighted mean magnitude offset to subtract from the host galaxy magnitude. These corrections were mostly small, with a median value of 0.13 mag. The one large outlier was the correction required for J0042 $-$ 3515: 3.25 mag. However, given that the standard error of the weighted mean is only 0.01 mag in this case, we have confidence that the calibrated magnitude for the host galaxy of this source is reliable. Galactic extinction was not corrected for given that the corrections are much smaller than the magnitude uncertainties for our targets.

2.3 Host galaxy magnitude measurements

We adopted a conservative approach when searching for the $K_\textrm{s}$ -band counterparts of our targets. While some sources have multiple potential host galaxy candidates, we only assigned a host when a $K_\textrm{s}$ -band source is plausibly close to the radio centroid, as one might expect for a compact radio source. Otherwise, the source was classified as having a $K_\textrm{s}$ -band non-detection. We will further discuss this strategy in Section 4.3.

Host galaxy magnitudes were determined as follows. Firstly, we used Source Extractor (sextractor; Bertin & Arnouts 1996) to efficiently identify as many host galaxy detections as possible. To increase the reliability of the results, we used a signal-to-noise (S/N) detection threshold $\geq 3$ over at least three connected pixels. However, it became apparent that some sources were not detected by sextractor using the above connected pixel criterion, yet there was visual evidence of a host galaxy (e.g. faint diffuse emission). Therefore, to ensure that the host galaxy magnitudes were determined in as consistent a manner as possible, we instead used photutils (Bradley et al. 2021) to measure a magnitude in a 2 $^{\prime\prime}$ diameter circular aperture for each source. We chose an aperture of this size as it is one of the standard choices used in the literature and many of our host galaxy detections are compact (as one might expect for a source at high redshift). In some cases, a larger circular aperture would be more appropriate to enclose all of the flux density from a given source; further information is provided in Appendix 1.

Each aperture was centred on, in order of preference, (i) the flux-weighted centre as determined by sextractor, (ii) the $K_\textrm{s}$ -band centroid as estimated by eye if the S/N was insufficient for (i), or (iii) the radio centroid when (ii) was deemed to be inaccurate due to low S/N. For (iii), we used the high-resolution Australia Telescope Compact Array (ATCA; Frater, Brooks, & Whiteoak 1992) data from B22 where available, factoring in both the angular resolution and S/N to choose the most accurate position at either 5.5 or 9 GHz. If high-resolution ATCA data were not available, we instead used the radio position from the 3-GHz VLA Sky Survey (VLASS; Lacy et al. 2020).

A host galaxy detection was regarded as significant if the integrated flux density determined with photutils within the 2 $^{\prime\prime}$ diameter aperture had S/N $\geq 3$ . The flux density statistical uncertainty was calculated by randomly placing $1\,000 \times 2$ $^{\prime\prime}$ diameter apertures within each image, avoiding sources. The $1\unicode{x03C3}$ uncertainty was then the standard deviation of these measurements. Each magnitude uncertainty was also calculated with the appropriate error propagation. In some cases (see the notes on individual sources in Appendix 1), we repeated this exercise with apertures of different sizes. For non-detections, we report a $3\unicode{x03C3}$ flux density (magnitude) upper (lower) limit. Moreover, for detections, we conservatively added a 2% flux density calibration error in quadrature with the statistical uncertainty (followed by error propagation to obtain the magnitude uncertainty).

2.4 Radio and near-infrared flux densities for the B22 sample and (U)HzRGs from the literature

In Fig. 1, we plot the 150-MHz radio flux density, $S_{150\,\rm{MHz}}$ , versus 2.2- $\unicode{x03BC}$ m near-infrared flux density, $S_{2.2\,\unicode{x03BC}\rm{m}}$ , for 40/53 sources from the B22 sample as well as (to the best of our knowledge) all of the known powerful radio galaxies at $z \gt 3$ with $L_{500\,\textrm{MHz}} \gt 10^{27}$ W Hz $^{-1}$ and $K_\textrm{s}$ -band data, plus the obscured, radio-loud AGN COS-87259 at $z=6.85$ (42 sources from the literature in total, one of which is also in the B22 sample: the UHzRG GLEAM J0856 $+$ 0223 at $z=5.55$ discovered by Drouart et al. 2020).

Figure 1. 150-MHz radio flux density (ordinate) versus 2.2- $\unicode{x03BC}$ m near-infrared flux density (bottom abscissa) for 40/53 radio sources from the B22 sample as well as for other HzRGs and UHzRGs from the literature with $3.08 \leq z \leq 6.85$ . Also plotted are four lines denoting various radio to near-infrared flux density ratios: $S_{150\,\rm{MHz}} / S_{2.2\,\unicode{x03BC}\rm{m}} =10^{3}$ – $10^{6}$ . The top abscissa shows the AB magnitudes corresponding to a selection of the 2.2- $\unicode{x03BC}$ m flux densities. Note that one of the grey ‘ $+$ ’ symbols is augmented with a triangle to denote a 2.2- $\unicode{x03BC}$ m flux density upper limit. For the sake of clarity, error bars are not plotted. Further details on the figure, including a discussion of the data presented, can be found in Sections 2 and 4. The references used for the 150-MHz flux densitie are as follows: Waldram et al. (1996), Hurley-Walker et al. (2017), Intema et al. (2017), and Endsley et al. (2022). Furthermore, the references used for the 2.2- $\unicode{x03BC}$ m flux densities and magnitudes are as follows: Eales et al. (1993), van Breugel et al. (1998), van Breugel et al. (1999), Villani & di Serego Alighieri (1999), Stern & Spinrad (1999), Jarvis et al. (2001), De Breuck et al. (2002, 2004), Brookes et al. (2006), Cruz et al. (2006), Seymour et al. (2007), Jarvis et al. (2009), Parijskij et al. (2014), Saxena et al. (2018, 2019), Drouart et al. (2020), Seymour et al. (2022), Endsley et al. (2022) and B22.

From the B22 sample, we only considered the sources that have deep HAWK-I or SHARKS $K_\textrm{s}$ -band images presented either in this paper (35 sources) or by us previously (six sources including GLEAM J0856 $+$ 0223; Drouart et al. 2020; Seymour et al. 2022; B22). The literature sources were obtained from the compilation of $3 \lt z \lt 5$ radio galaxies in Miley & De Breuck (2008, with the exception of two sources which only had optical and not near-infrared identifications), with additional sources from Jarvis et al. (2009), Saxena et al. (2019) and Yamashita et al. (2020). We also included the UHzRGs TN J0924 $-$ 2201 (z = 5.19; van Breugel et al. 1999) and TGSS J1530+1049 (z = 5.72; Saxena et al. 2018). However, we did not include the radio galaxy VLA J123642 $+$ 621331 at $z=4.42$ (Waddington et al. 1999) as it is a weaker radio source that does not meet our radio luminosity criterion; the $3\unicode{x03C3}$ upper limit at 147.5 MHz from the Tata Institute of Fundamental Research (TIFR) Giant Metrewave Radio Telescope (GMRT; Swarup 1991) Sky Survey (TGSS) Alternate Data Release 1 (Intema et al. 2017) is about 6 mJy beam $^{-1}$ .

The values of $S_{150\,\rm{MHz}}$ and $S_{2.2\,\unicode{x03BC}\rm{m}}$ per source were determined as follows.

2.4.1 Calculation of 150-MHz flux densities

Apart from two cases, the 150-MHz flux densities were obtained from GLEAM or TGSS. Where GLEAM data were available from both the Extragalactic data release (Hurley-Walker et al. 2017) as well as the deeper South Galactic Pole data release (Franzen et al. 2021), we used the data from the latter. Furthermore, where GLEAM and TGSS data were available for a particular source, we used the GLEAM flux densities only. Some of the sources in the literature sample are too far north on the sky to have been catalogued in GLEAM, but TGSS data were available instead. The sources without GLEAM or TGSS data are as follows. Firstly, the $z=3.395$ radio galaxy B2 0902 $+$ 34 is too far north to have been catalogued in GLEAM, but it was also not catalogued in TGSS despite being a bright radio source (possibly due to it being in one of the failed pointings in the survey; see Intema et al. 2017 ^d ), and so we instead used the 151-MHz flux density from the 7C survey (Waldram et al. 1996). Secondly, COS-87259 is too faint to have been catalogued in GLEAM and TGSS, and we instead used the 144-MHz flux density reported in Endsley et al. (2022) from the Low-Frequency Array (LOFAR; van Haarlem et al. 2013) Two-metre Sky Survey (LoTSS; Shimwell et al. 2017, 2019, 2022).

To enable as consistent a comparison as possible between the B22 sample and the sample from the literature, we did not use the curved fits to the broadband GLEAM data for the former presented in B22, but instead took the 143-, 151-, 158-, and 166-MHz GLEAM flux densities for each source and calculated the inverse-variance-weighted average. This step removes some of the scatter in the single-band GLEAM flux densities. At each frequency, a 2% internal calibration uncertainty (see Hurley-Walker et al. 2017) was added in quadrature to the fitting error reported in the GLEAM catalogue before the weighted average was computed.

Note that for simplicity, we did not make any corrections for the slightly different central frequencies of the data discussed above, nor for the fact that these values have small offsets from a frequency of 150 MHz. The scale offsets between the data sets used will affect the flux densities by $\lesssim10$ % (or $ \lesssim 0.04$ dex).

The inverse-variance-weighted average 150-MHz GLEAM flux densities for the 35 sources from the B22 sample with HAWK-I data presented in this paper are listed in Table 1.

Table 1. $K_\textrm{s}$ -band properties of the 35 sources observed with HAWK-I. We present the host galaxy magnitudes and corresponding flux densities determined in 2 $^{\prime\prime}$ diameter circular apertures along with uncertainties. Magnitude (flux density) lower (upper) limits are $3\unicode{x03C3}$ values. Host galaxy positions are from sextractor unless specified with a footnote. Where the $K_\textrm{s}$ -band position was sufficiently well determined, we also provide the $K_\textrm{s}$ -band–radio angular separation. Further details can be found in Sections 2.3 and 3.1. Furthermore, we list 150-MHz inverse-variance-weighted average flux densities from GLEAM; see Section 2.4.1. The uncertainties for these GLEAM flux densities include an extra 8% absolute flux density calibration uncertainty (see Hurley-Walker et al. 2017) added in quadrature.

Notes. $^\unicode{x2020}$ UHzRG candidate. $^\textrm{a}$ Radio position where the aperture was centred; see Appendix 1 for further details. $^\textrm{b}$ Separation uncertain due to an insufficiently well-constrained $K_\textrm{s}$ -band position. $^\textrm{c}$ This incipient double-lobed radio source was fitted with a single Gaussian in RACS-mid. $^\textrm{d}$ Estimated by eye. $^\textrm{e}$ Host galaxy detected in a $1.\!\!^{\prime\prime}5$ diameter aperture; see the notes on this source in Appendix 1 for further details. $^\textrm{f}$ Distance from the $K_\textrm{s}$ -band identification to the possible radio core visible at 5.5 GHz.

2.4.2 Calculation of 2.2- $\unicode{x03BC}$ m flux densities

Regarding the 2.2- $\unicode{x03BC}$ m flux densities, different strategies were used in the literature to obtain the most accurate measurement for a given source, particularly regarding the size of the aperture used. Where more than one measurement was available for the literature sample, we used the value from the largest-sized aperture considered in order to increase the likelihood that all of the flux density from the galaxy was enclosed in the aperture. In many cases, this was an 8 $^{\prime\prime}$ diameter aperture or a 64-kpc diameter aperture. Furthermore, to ensure consistency, we took the magnitude measurements from apertures larger than a 2 $^{\prime\prime}$ diameter for the sources J0007 $-$ 3040, J0034 $-$ 3112, J0048 $-$ 3540, J0301 $-$ 3132, J0326 $-$ 3013, J0909 $-$ 0154, J1032 $+$ 0339, J1040 $+$ 0150, and J1052 $-$ 0318 from the B22 sample (Appendix 1; also see Section 2.3).

In Jarvis et al. (2009), the 2.2- $\unicode{x03BC}$ m magnitude of the HzRG FIRST J163912.11 $+$ 405236.5 ( $z=4.88$ ) was estimated to be $K = 24.2$ mag (4 $^{\prime\prime}$ diameter aperture). This estimate was obtained by extrapolating from a 3.6- $\unicode{x03BC}$ m detection with an assumed colour (2.2- $\unicode{x03BC}$ m $-$ 3.6- $\unicode{x03BC}$ m = 2.25 mag). Alternatively, using the photometry reported in Jarvis et al. (2009), we followed a similar approach to the one that we used in Seymour et al. (2022) and carried out broadband spectral energy distribution (SED) fitting with the pégase code (version 3; Fioc & Rocca-Volmerange 2019) for different galaxy templates; the 2.2- $\unicode{x03BC}$ m magnitude in a 4 $^{\prime\prime}$ diameter aperture is estimated to be $K \sim 23.4$ mag (i.e. a flux density a factor of $\sim$ 2 brighter: $S_{2.2\,\unicode{x03BC}\rm{m}}\sim 1.6$ $\unicode{x03BC}$ Jy). We have used our estimated magnitude in this study; in either case, this source is the faintest at 2.2 $\unicode{x03BC}$ m from our compilation of HzRGs at $3 \lt z \lt 5$ .

The SHARKS magnitude limit for J0008 $-$ 3007 reported in B22 was converted from $5\unicode{x03C3}$ to $3\unicode{x03C3}$ (i.e. by adding 0.55 mag; $K_\textrm{s} \gtrsim 23.7$ mag and $S_{2.2\,\unicode{x03BC}\rm{m}} \lesssim 1.2$ $\unicode{x03BC}$ Jy; 2 $^{\prime\prime}$ diameter aperture). The SHARKS limit for J0007 $-$ 3040 reported in B22 is superseded by the HAWK-I data presented in this study. Regarding the SHARKS data presented in B22 for J2340 $-$ 3230, the interpretation of this source (and hence the validity of a $3\unicode{x03C3}$ limit) is uncertain (see Section 5.5 in B22), and deeper $K_\textrm{s}$ -band data are needed. We therefore did not include J2340 $-$ 3230 in our analysis, that is, the number of sources from the B22 sample with previously presented deep $K_\textrm{s}$ -band data was reduced from six to five.

For simplicity, we did not make any corrections for the slightly different central wavelengths of the image filters used in the various studies in the literature (e.g. $K_\textrm{s}$ versus K), and considered all values as measurements at 2.2 $\unicode{x03BC}$ m.

3. Results

3.1 Host galaxy magnitudes

The host galaxy magnitude measurements (2 $^{\prime\prime}$ diameter apertures) from our HAWK-I observing campaign are presented in Table 1, along with the corresponding flux density measurements, $S_{2.15\,\unicode{x03BC}\textrm{m}}$ , as well as the positions on which the apertures were centred (see Section 2.3). As can be seen in Table 1, 27 out of the 35 sources (77%) were detected in the HAWK-I images. For these detections, the host galaxy magnitudes range from $K_\textrm{s} \approx 21.6$ – $23.0$ mag. The remaining eight non-detections have $3\unicode{x03C3}$ magnitude lower limits that are generally quite uniform (median depth $K_\textrm{s} \approx 23.3$ mag). We consider five of these non-detections (J0053 $-$ 3256, J1033 $+$ 0107, J1127 $-$ 0332, J1443 $+$ 0229, and J2326 $-$ 3028; marked in Table 1 with $^\unicode{x2020}$ ) to be our best UHzRG candidates, along with J0008 $-$ 3007 from B22 (deep SHARKS limit). Notes on individual sources are given in Appendix 1.

3.2 Overlay plots

In Fig. 2, we present high-resolution radio contours overlaid on our HAWK-I images for five of the six best UHzRG candidates (see Figure 2 in B22 for the previously presented overlay plot for J0008 $-$ 3007). The remaining 30 overlay plots can be found in Fig. A1 in Appendix 1.

Figure 2. Radio contours overlaid on HAWK-I $K_\textrm{s}$ -band images for five of the six best UHzRG candidates from this study, ordered by right ascension. The VLASS (3 GHz; magenta) and ATCA (5.5 and 9 GHz; blue and red, respectively) contours were previously presented in B22, except for the new single-epoch contours from VLASS for J1033 $+$ 0107. The contours are a geometric progression in $\sqrt{2}$ , with the lowest contour at the $5\unicode{x03C3}$ level; a summary of the contour levels can be found in Table 5 in B22. For the single-epoch image from VLASS for J1033 $+$ 0107, the updated lowest contour level ( $5\unicode{x03C3}$ ) is 0.75 mJy beam $^{-1}$ . In each panel, we also show the radio synthesised beams with different hatching styles (VLASS: vertical; ATCA 5.5 GHz: forward slash; ATCA 9 GHz: backslash). The host galaxy magnitude lower limits are reported in Table 1. The overlay plots for the remaining 30 sources observed with HAWK-I can be found in Fig. A1 in Appendix 1, while the SHARKS/VLASS/ATCA overlay for J0008 $-$ 3007, our sixth UHzRG candidate, can be found in Figure 2 in B22. Furthermore, notes on individual sources can be found in Appendix 1.

Figure 3. Tracks of modelled 2.2- $\unicode{x03BC}$ m (i.e. $K_\textrm{s}$ -band) flux density from pégase for galaxies of different masses formed no earlier than $z=16$ . ‘Ell’ and ‘SB’ refer to an old, 0.5-Gyr elliptical galaxy template and a young, 30-Myr starburst galaxy template, respectively. The elliptical tracks terminate when the age implies a formation redshift above $z=16$ ; see Section 4.2 for more details. As in Fig. 1, we overlay TN J0924 $-$ 2201, GLEAM J0856 $+$ 0223, TGSS J1530 $+$ 1049 and COS-87259, as well as the HzRGs with $4 \lt z \lt 5$ from our literature sample. We also indicate the median $3\unicode{x03C3}$ limit of our deep non-detections ( $S_{2.2\,\unicode{x03BC}\rm{m}} \lesssim 1.5$ $\unicode{x03BC}$ Jy; 2 $^{\prime\prime}$ diameter apertures). Error bars ( $\pm1\unicode{x03C3}$ ) are plotted (apart from the case of the $3\unicode{x03C3}$ limit for TGSS J1530 $+$ 1049), but sometimes are smaller than the symbols. When uncertainties were not available in the literature, we assumed a 10% error. Discussion of this figure can be found in Section 4.2.

The radio data used in Figs. 2 and A1 were mostly presented in B22. Apart from four cases outlined below, we have only plotted contours that are at an angular resolution of a few arcsec or better: 3-GHz ‘quick-look’ data from the first epoch of VLASS ^e as well as 5.5- and 9-GHz ATCA data. For J1030 $+$ 0135, J1032 $+$ 0339, J1033 $+$ 0107, and J1037 $-$ 0325, we used the higher-quality, single-epoch (epoch 2.1) VLASS continuum images ^f instead, which had become available since the B22 study.

The overlays for J0007 $-$ 3040, J0034 $-$ 3112, and J1037 $-$ 0325 presented in Fig. A1 also include data from the Rapid Australian Square Kilometre Array Pathfinder (ASKAP; Johnston et al. 2007; Hotan et al. 2021) Continuum Survey at 1367.5 MHz (RACS-mid; median angular resolution 11 $\hbox{$.\!^{\prime\prime}$}$ 2× 9 $\hbox{$.\!^{\prime\prime}$}$ 3; Duchesne et al. 2023, 2024). ^g These RACS-mid data became available near the completion of our study and changed the interpretation of these three sources. Originally thought to be additional UHzRG candidates, these are larger radio sources than previously thought and therefore much more likely to lie at lower redshift. For both J0007 $-$ 3040 and J0034 $-3112$ , the high-resolution VLASS and ATCA contours trace only one of the lobes. Further details can be found in the notes on individual sources in Appendix 1.

For J2311 $-$ 3359, also presented in Fig. A1, given that this source is only faintly detected in the ATCA data presented in B22, and not detected in VLASS, we confirmed the host galaxy identification using the publicly available data (angular resolution 10 $^{\prime\prime}$ ) from a recent study by Gürkan et al. (2022), who conducted deep 887.5-MHz ASKAP observations of the GAMA-23 field from the Galaxy And Mass Assembly survey (Driver et al. 2009, 2011, 2016). ^h

$K_\textrm{s}$ -band–radio angular separations are reported in Table 1 for the 25 sources where the $K_\textrm{s}$ -band position is known to sufficient accuracy. We again used the high-resolution ATCA data where available (factoring in both the angular resolution and S/N to choose the most accurate position at either 5.5 or 9 GHz), or otherwise the VLASS position. However, for J0007 $-$ 3040, J0034 $-$ 3112, J1037 $-$ 0325, and J2311 $-$ 3359, we instead used the ASKAP data from either RACS-mid or the Gürkan et al. (2022) study. For multi-component radio sources, the geometric midpoint was calculated unless otherwise specified in Table 1. Given both the mostly compact radio morphologies of the targets and our conservative approach at assigning host galaxy identifications, it is not surprising that the offsets are small, with a median value of approximately 0 $.^{\!{\prime\prime}}$ 5.

4. Discussion

4.1 Comparison to known HzRGs

As can be seen in Fig. 1 and as was initially discussed in Section 3.1, our detection rate for the B22 sample is high in $K_\textrm{s}$ -band: 78% for the 40 sources plotted. We can see that the B22 sample overlaps with the $3 \lt z \lt 5$ powerful HzRGs, presenting an extension of the observed trend (i.e. flux density ratios $S_{150\,\rm{MHz}} / S_{2.2\,\unicode{x03BC}\rm{m}}$ $\sim$ $10^{4}$ – $10^{6}$ ) towards lower 150-MHz and 2.2- $\unicode{x03BC}$ m flux densities. Extreme flux density ratios $\gtrsim10^5$ are not surprising as they arise from a common approach in the various HzRG selection techniques applied in the literature: sufficiently bright radio flux densities to select luminous radio galaxies (most likely powered by SMBHs) combined with faint 2.2- $\unicode{x03BC}$ m magnitudes to ensure that the host galaxy is distant. By and large, these extreme ratios are only seen at $z\gt 1$ (with a notable exception being Cygnus A; e.g. see Figure 5 in Drouart et al. 2021). Fainter 150-MHz and 2.2- $\unicode{x03BC}$ m flux densities could, generally speaking, result from less massive galaxies with less massive BHs, or an increase in redshift. Which is the dominant of these two effects in the B22 sample is unclear given the lack of spectroscopic redshift information.

It should be noted that the sources presented in Fig. 1 have been selected using a variety of multi-wavelength criteria. Additionally, as one moves from lower to higher redshift, k-corrections at both radio and near-infrared wavelengths affect the flux densities. As discussed in Section 3.3 of Drouart et al. (2021), these k-corrections account for separate segments from the SED that are both approximated by power laws with negative slopes; see that study for further details.

4.2 UHzRG candidates from the B22 sample

Of particular interest in Fig. 1 are the six sources (marked in blue) previously outlined in Section 3.1 as our best UHzRG candidates. This subset of six sources has median $3\unicode{x03C3}$ limits of $K_\textrm{s} \gtrsim 23.5$ mag and $S_{2.2\,\unicode{x03BC}\rm{m}} \lesssim 1.5$ $\unicode{x03BC}$ Jy (2 $^{\prime\prime}$ diameter apertures); moreover, $S_{150\,\rm{MHz}} \gt 200$ mJy. All of these sources have extreme radio to near-infrared flux density ratios, $S_{150\,\rm{MHz}} / S_{2.2\,\unicode{x03BC}\rm{m}} \gt 10^5$ , similar to a number of the radio galaxies from the literature plotted in Fig. 1, including the UHzRGs TN J0924 $-$ 2201 ( $z=5.19$ ), GLEAM J0856 $+$ 0223 ( $z=5.55$ ), and TGSS J1530 $+$ 1049 ( $z=5.72$ ). While at least slightly fainter in $K_\textrm{s}$ -band than TN J0924 $-$ 2201 and GLEAM J0856 $+$ 0223, our UHzRG candidates have 2.2- $\unicode{x03BC}$ m flux density upper limits that are about a factor of two brighter than the upper limit for TGSS J1530 $+$ 1049: $S_{2.2\,\unicode{x03BC}\rm{m}} \lt 0.76$ $\unicode{x03BC}$ Jy and $K_\textrm{s} \gt 24.2$ mag ( $3\unicode{x03C3}$ ; $1.\!\!^{\prime\prime}5$ diameter aperture; Saxena et al. 2018, 2019). For the five sources with HAWK-I data, repeating the limit calculations in $1.\!\!^{\prime\prime}5$ diameter apertures (see the notes on individual sources in Appendix 1) so as to better match the analysis for TGSS J1530 $+$ 1049 yielded deeper median $3\unicode{x03C3}$ limits of $K_\textrm{s} \gtrsim 23.8$ mag and $S_{2.2\,\unicode{x03BC}\rm{m}} \lesssim 1.1$ $\unicode{x03BC}$ Jy. Henceforth, however, we will use the more conservative limits from the 2 $^{\prime\prime}$ diameter apertures.

Only recently has an obscured, radio-loud AGN been detected at $z \gt 6$ : COS-87259 with $z=6.85$ , previously discussed in Section 1. This source has $S_{2.2\,\unicode{x03BC}\rm{m}} = 0.67 \pm 0.09$ $\unicode{x03BC}$ Jy ( $K_\textrm{s} = 24.33 \pm 0.15$ mag; $1.\!\!^{\prime\prime}2$ diameter aperture; Endsley et al. 2022), that is, too faint to have been detected in our deepest HAWK-I observations. It was found as part of the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007), which has unparalleled multi-wavelength data, albeit over 2 deg $^{2}$ only. Assuming that the rest-frame ultraviolet and optical emission is dominated by star formation, the faint 2.2- $\unicode{x03BC}$ m flux density is consistent with a massive galaxy with stellar mass $M_\textrm{*} = (1.7 \pm 0.7) \times 10^{11}\,{\rm M}_\odot$ (Endsley et al. 2023b, although those authors explain that the mass could lie anywhere between $10^{10}$ – $10^{11}\,\textit{M}_\odot$ depending on the star formation history). The host galaxy has a very high obscured AGN infrared luminosity, $L_\textrm{IR} = (2.5 \pm 0.2) \times 10^{13}\,{\rm L}_\odot$ , which is not evident in the rest-frame ultraviolet/optical (Endsley et al. 2023b). Its stellar mass and AGN luminosity are similar to powerful radio galaxies at $2 \lt z \lt 5$ (e.g. Seymour et al. 2007; Drouart et al. 2016). While COS-87259 has a faint 144-MHz flux density ( $S_{144\,\textrm{MHz}} = 475 \pm 180$ $\unicode{x03BC}$ Jy; Endsley et al. 2022) and is therefore much less luminous in the radio than the powerful radio galaxies referenced above, it nonetheless demonstrates the existence of detectable radio emission from an obscured AGN within the EoR. Potentially, some of the prime targets from the B22 sample are the radio-luminous analogues of COS-87259, but spectroscopic redshift determinations are needed.

To obtain further evidence that our six deep non-detections are strong UHzRG candidates, in Fig. 3 we have used pégase to model 2.2- $\unicode{x03BC}$ m flux density as a function of redshift for several galaxy templates (elliptical and starburst) with different stellar masses. We followed a similar method as that described in both Drouart et al. (2021) and Seymour et al. (2022), but assumed that the galaxies formed no earlier than $z=16$ when the Universe was $\sim250\,$ Myr old. The elliptical template was assumed to be old (0.5 Gyr: the oldest it could be at $z\sim 7$ given our formation redshift and therefore the tracks terminate at $z=7.1$ ), whereas the starburst template was assumed to be young (0.03 Gyr: chosen to match the youngest galaxies seen at $7\lt z\lt 8$ ; Endsley et al. 2023a). We have overlaid the positions of the nine HzRGs with $4\lt z\lt 5$ included in Fig. 1, as well as TN J0924 $-$ 2201, GLEAM J0856 $+$ 0223, TGSS J1530 $+$ 1049, and COS-87259. In addition, we show the median $3\unicode{x03C3}$ limit of our deep non-detections. Several conclusions can be drawn from Fig. 3 and are as follows.

1. 1. If the host galaxies of our $K_\textrm{s} \gtrsim 23.5$ mag radio sources are massive like other known powerful HzRGs at $2\lt z\lt 5$ ( $M_\textrm{*} \gtrsim 10^{11}$ M $_\odot$ ; e.g. Rocca-Volmerange et al. 2004; Seymour et al. 2007; De Breuck et al. 2010; Saxena et al. 2019), then they must lie at ultra-high redshift, $z\gtrsim 5.4$ , to be so faint.

2. 2. For a lower stellar mass, $M_\textrm{*} \sim 10^{10.5}$ M $_\odot$ , the elliptical and starburst model tracks indicate that $z\gtrsim 4.2$ (elliptical) or $z\gtrsim 4.9$ (starburst) for $K_\textrm{s} \gtrsim 23.5$ mag.

3. 3. A corollary of the two preceding conclusions is that for $M_\textrm{*} \gtrsim 10^{10.5}$ M $_\odot$ and $K_\textrm{s} \gtrsim 23.5$ mag, $z \gtrsim 4.2$ .

4. 4. Another possibility is that our deep non-detections are perhaps indicative of a class of powerful radio-loud AGN in undermassive ( $M_\textrm{*} \lesssim 10^{10}\,{\rm M}_\odot$ ) host galaxies. Although not visible in the region of parameter space shown in Fig. 3, $z\gtrsim 2.4$ (elliptical) or $z\gtrsim 2.5$ (starburst) for $M_\textrm{*} \sim 10^{10}$ M $_\odot$ and $K_\textrm{s} \gtrsim 23.5$ mag. However, given the findings from previous studies of powerful HzRGs, undermassive host galaxies would be unexpected, at least at $z \lt 5$ . Furthermore, radio-loud AGN in dwarf galaxies ( $M_\textrm{*} \lesssim 10^{9.5}\,\textit{M}_\odot$ ) are at least a few orders of magnitude less radio luminous (e.g. Mezcua, Suh, & Civano 2019; Davis et al. 2022). On the other hand, at $z\gt 4$ , there is emerging evidence from the James Webb Space Telescope (JWST; Gardner et al. 2006, 2023) of a population of AGN not detected in the radio with BH to stellar mass ratios up to several dex larger than is what usually observed (e.g. Furtak et al. 2023, 2024; Goulding et al. 2023; Pacucci et al. 2023; Kokorev et al. 2023; Maiolino et al. 2024). Assuming a fiducial value of $M_\textrm{BH}$ $\gtrsim 10^9$ M $_\odot$ for a powerful radio galaxy (e.g. Nesvadba et al. 2011; Drouart et al. 2014; also see Willott et al. 1999; Kauffmann et al. 2008; Martnez-Sansigre & Rawlings 2011), then, in the case of an undermassive host galaxy scenario, our targets might have BH to stellar mass ratios exceeding the local relation (e.g. Reines & Volonteri 2015) and potentially approaching the sources analysed in the above JWST studies.

5. 5. COS-87259 lies below both of the $10^{11}$ M $_\odot$ tracks and the $10^{10.5}$ M $_\odot$ starburst track; its estimated stellar mass from broadband SED fitting is in the range $M_\textrm{*} \sim$ $10^{10}$ – $10^{11}$ M $_\odot$ (assuming no AGN contribution to the ultraviolet/optical flux density; Endsley et al. 2023b). This large range is due to the uncertainty in the star formation history of this source, which demonstrates how uncertain stellar mass estimates can be even with extensive and accurate photometry. It is worth noting that from its BH mass estimate ( $M_\textrm{BH} \gtrsim 1.5 \times 10^9$ M $_\odot$ ), COS-87259 is likely to have a relatively low stellar mass compared to local galaxies; see the discussion in Section 4 of Endsley et al. (2023b).

6. 6. A suggested stellar mass of $M_\textrm{*} \lt 10^{10.5}$ M $_\odot$ (starburst) or $M_\textrm{*} \lt 10^{11}$ M $_\odot$ (elliptical) for TGSS J1530 $+$ 1049 is also broadly consistent with the modelling carried out by Saxena et al. (2019), who determined that $M_\textrm{*} \lt 10^{10.5}$ M $_\odot$ . TN J0924 $-$ 2201 is best described by an elliptical template with $M_\textrm{*} \sim 10^{11}$ M $_\odot$ or a starburst with $M_\textrm{*} \sim 10^{10.5}\,\textrm{M}_\odot$ ; from multi-wavelength SED fitting, Seymour et al. (2007) determined the maximum stellar mass of this source to be $M_\textrm{*} = 10^{11.1}$ M $_\odot$ . Lastly, GLEAM J0856 $+$ 0223 is potentially at least slightly more massive than TN J0924 $-$ 2201 (also see discussion in Drouart et al. 2020).

7. 7. Though with some scatter, the inclusion of the $4 \lt z \lt 5$ HzRGs in Fig. 3 demonstrates an increase in 2.2- $\unicode{x03BC}$ m flux density with decreasing redshift, as expected from the K–z relation, albeit with an upper envelope described by a $\sim10^{11.5}\,\textrm{M}_\odot$ stellar mass compared to $\sim10^{12}\,\textrm{M}_\odot$ at lower redshift (e.g. Rocca-Volmerange et al. 2004). All but one of these HzRGs have 2.2- $\unicode{x03BC}$ m flux densities at least a factor of two above the median $3\unicode{x03C3}$ limit for our deep non-detections.

4.3 Alternative scenarios for the UHzRG candidates

What is the nature of the $K_\textrm{s} \gtrsim 23.5$ sources from B22 should they not be at $z \gtrsim 4$ –5? With an undermassive host galaxy scenario unprecedented at $z\lt 4$ given the powerful radio emission, these sources could instead be massive galaxies with heavy dust obscuration. As discussed in B22, none of our UHzRG candidates have detections in the mid-infrared Widefield Infrared Survey Explorer (WISE; Wright et al. 2010) AllWISE data release (Cutri et al. 2014). Additionally, no detections were found in the deeper unWISE (Schlafly, Meisner, & Green 2019) data release. However, unWISE is much shallower than our HAWK-I data; the 50% completeness limits for the 3.37- $\unicode{x03BC}$ m (W1) and 4.62- $\unicode{x03BC}$ m (W2) bands are $W1 = 20.72$ mag/ $S_{3.37\,\unicode{x03BC}\rm{m}} = 18.71$ $\unicode{x03BC}$ Jy and $W2 = 19.97$ mag/ $S_{4.62\,\unicode{x03BC}\rm{m}} = 37.33$ $\unicode{x03BC}$ Jy, respectively (Schlafly et al. 2019). JWST is opening up a new window on the study of dusty starburst galaxies at (ultra-)high redshift (e.g. Zavala et al. 2023; Barrufet et al. 2023), suggesting that we would potentially need much deeper mid-infrared data to detect our targets. Moreover, COS-87259 has 3.6- and 4.5- $\unicode{x03BC}$ m flux densities of $2.75 \pm 0.10$ and $2.77 \pm 0.12$ $\unicode{x03BC}$ Jy, respectively (Endsley et al. 2022).

Additionally, in B22, we investigated if all of the sources in the sample had been identified as pulsars in the literature, but none were. For the non-detections in the $K_\textrm{s} \gtrsim 23.5$ mag subset, a pulsar origin seems unlikely given that all of these sources are at least slightly extended in the radio (Table 2 in B22). None of the sources show evidence for variability at radio frequencies either. Furthermore, none of the sources are detected in the optical in the third data release of the SkyMapper Southern Sky Survey (SMSS; Onken et al. 2019), ⁱ nor in the 18th data release of the Sloan Digital Sky Survey SDSS; Almeida et al. 2023) for the equatorial targets. While our sources do not have X-ray data, the lack of radio variability, optical non-detections, and, additionally, the relatively steep spectral indices at GHz frequencies (B22 cf. e.g. Fender 2001 and Trushkin, Bursov, & Nizhelskij 2003) do not paint a compelling picture that we are seeing persistent radio emission from Galactic X-ray binaries (e.g. see the radio–optical parameter space presented in Stewart et al. 2018).

Lastly, our non-detections could be radio sources with one-sided lobes, that is, more extreme versions of J0007 $-$ 3040 and J0034 $-$ 3112 (as well as J2330 $-$ 3237) in Fig. A1. Therefore, the host galaxies may be located elsewhere in the HAWK-I images. For these sources, there is no evidence in RACS-mid of a second candidate lobe. While we cannot fully rule out that at least some of the $K_\textrm{s} \gtrsim 23.5$ mag sources have one-sided radio morphologies, the number density of K-band sources is $\sim$ 70 arcmin $^{-2}$ for $23 \lt K \lt 24$ (Fontana et al. 2014). Therefore, the probability of a host galaxy candidate within this magnitude range being within 2 $^{\prime\prime}$ of the radio centroid by chance is about 20%, that is, for $\sim$ 1 of our 6 UHzRG candidates with non-detections. The corresponding values for a 5 $^{\prime\prime}$ search radius are 80% and $\sim$ 5 sources. Therefore, we can see that this quickly becomes a challenging exercise to identify the true host galaxy for increasingly large, one-sided radio sources. Very long baseline interferometry (VLBI) would allow a better assessment of whether the radio emission from the UHzRG candidates is lobe-like, possibly with backflow towards a host galaxy candidate, or whether the emission is resolved into multiple lobes and also potentially a core, which would strengthen the argument for a compact radio source with an undetected host galaxy in $K_\textrm{s}$ -band. Such observations are underway and will be presented in a future study.

The above analysis of chance radio–near-infrared matches is also relevant for the host galaxy identifications presented in Section 3, although the offset from the radio centroid is typically much smaller (median 0 $.^{\!{\prime\prime}}$ 5; Table 1). We would therefore expect at most one source of the 27 with $K_\textrm{s}$ -band counterparts in Table 1 to be a chance radio–near-infrared match.

4.4 Comparison to IFRSs

Our UHzRG candidates, and indeed the majority of the other sources from the B22 sample studied in this paper, have properties consistent with the class of infrared-faint radio sources (IFRSs; e.g Norris et al. 2006, 2011; Middelberg et al. 2008, 2011; Garn & Alexander 2008; Zinn, Middelberg, & Ibar 2011; Collier et al. 2014; Maini et al. 2016). IFRSs have typically been characterised on the basis of their 1.4-GHz radio to 3.6- $\unicode{x03BC}$ m mid-infrared flux density ratios; Zinn et al. (2011) proposed standard IFRS selection criteria of (i) $S_{1.4\,\rm{GHz}} / S_{3.6\,\unicode{x03BC}\rm{m}} \gt 500$ and (ii) $S_{3.6\,\unicode{x03BC}\rm{m}} \lt 30$ $\unicode{x03BC}$ Jy. For the 40 sources from the B22 sample plotted in Fig. 1, the 1.4-GHz flux densities from the National Radio Astronomy Observatory (NRAO) VLA Sky Survey (NVSS; Condon et al. 1998) span the range $\lt $ 1.6–98.9 mJy, with a median of 18.7 mJy (Table 2 in B22). Even after factoring in a k-correction between 2.2 and 3.6 $\unicode{x03BC}$ m, which will depend on the SED, we would expect the majority of our sources to (comfortably) meet the Zinn et al. (2011) selection criteria. Indeed, this is suggested in Fig. 1, albeit with plotted flux density ratios based on the radio value at 150 MHz ( $S_{150\,\rm{MHz}} / S_{1.4\,\rm{GHz}} \sim 10$ for a radio spectral index $\alpha = -1$ ; flux density $S_{\nu} \propto \nu^{\alpha}$ ).

As mentioned above in Section 4.3, none of our $K_\textrm{s} \gtrsim 23.5$ mag UHzRG candidates were detected in AllWISE or unWISE. However, despite the difference in sensitivity between unWISE and our HAWK-I data, seven sources (with $S_{2.2\,\unicode{x03BC}\rm{m}} = 2.54$ – $8.02$ $\unicode{x03BC}$ Jy in 2 $^{\prime\prime}$ diameter apertures) from the B22 sample plotted in Fig. 1 have unWISE detections. In Table 2, we present the unWISE and NVSS flux densities for these sources, along with the values of $S_{1.4\,\rm{GHz}} / S_{3.37\,\unicode{x03BC}\rm{m}}$ . All seven sources meet the Zinn et al. (2011) IFRS selection criteria, albeit considering a slightly different mid-infrared observing wavelength. An additional handful of sources from the B22 sample may also have unWISE counterparts (e.g. J1335 $+$ 0112 and J2219 $-$ 3312, which were discussed in Section 5.5 in B22 as having potential AllWISE detections), but it was not clear whether these detections are affected by confusion from nearby sources visible at 2.2 $\unicode{x03BC}$ m in the higher-resolution HAWK-I images. A detailed analysis was beyond the scope of this paper.

Table 2. Sources in the B22 sample with unWISE detections. For each unWISE detection, we also report the 1.4 GHz flux density from NVSS and the flux density ratio $S_{1.4\,\rm{GHz}} / S_{3.37\,\unicode{x03BC}\rm{m}}$ .

Norris et al. (2011) considered the scenarios listed in Section 4.3, as well as others, to explain the nature of IFRSs. They concluded that most IFRSs are likely either radio-loud AGN at $z \gtrsim 3$ or radio-loud AGN at $1 \lt z \lt 3$ that are significantly obscured by dust. Subsequent studies have found that many IFRSs are at high redshift, up to $z=4.39$ (e.g. Collier et al. 2014; Herzog et al. 2014; Singh et al. 2017; Orenstein, Collier, & Norris 2019). This bodes well for the B22 sample.

4.5 The role of wide and deep near-infrared surveys in the study of UHzRGs

From the above discussion of Fig. 1, it is clear that the six sources from the B22 sample with the faintest $K_\textrm{s}$ -band magnitudes are the most promising UHzRG candidates and could potentially fall within the EoR. However, obtaining deep $K_\textrm{s}$ -band imaging is expensive in terms of telescope time, especially as these rare bright radio sources are well-spaced across the sky and therefore require individual observations.

Various ongoing and future near-infrared surveys that are both wide and deep offer excellent prospects for improving the efficiency of UHzRG searches. For example, SHARKS, made use of in this project (B22), will cover 300 deg $^2$ to a median $5\unicode{x03C3}$ limiting magnitude of $K_\textrm{s} = 22.7$ mag (2 $^{\prime\prime}$ diameter aperture; $S_{2.2\,\unicode{x03BC}\rm{m}} = 3.0$   $\unicode{x03BC}$ Jy; Dannerbauer et al. 2022). The $K_\textrm{s}$ -band magnitudes of TN J0924 $-$ 2201 and GLEAM J0856 $+$ 0223 (van Breugel et al. 1999; Drouart et al. 2020) are at brightness levels approximately corresponding to the $3\unicode{x03C3}$ SHARKS depth. Referring to Fig. 3, SHARKS may be able to find UHzRGs with $M_\textrm{*} \gtrsim 10^{11}$ M $_\odot$ .

The Euclid mission (Racca et al. 2016; Euclid Collaboration et al. 2024) will also significantly advance searches for UHzRGs. The Euclid Wide Survey (EWS; Laureijs et al. 2011; Euclid Collaboration et al. 2022) will cover approximately 15 000 deg $^2$ in a wide optical band (530–920 nm) and three near-infrared filters (Y-, J-, and H-bands; 0.95–2.0 $\unicode{x03BC}$ m). The near-infrared point-source depth is expected to be 24.5 mag ( $5\unicode{x03C3}$ ; radius of 50% encircled energy $\lt 0\hbox{$.\!^{\prime\prime}$}4$ ), that is, approximately 1.5 mag fainter (at an equivalent detection threshold) than the deepest HAWK-I imaging presented here. Such data will greatly speed up the discovery of UHzRG candidates going forward; as a lower limit, extrapolating from six UHzRG candidates over approximately 1 200 deg $^2$ from our study, we would expect $\gtrsim 75$ UHzRG candidates to be found in the EWS. Furthermore, while the Euclid data alone can be used to discover bright regular galaxies at $z\gt 6$ –7, having a strong radio detection as well will immediately rule out contamination by very low-mass stars or dwarf star-forming galaxies at $1\lt z\lt 2$ , which are not known to host powerful AGN. The required radio surveys are extant, with, for example, GLEAM complete and its successor (GLEAM-X; Hurley-Walker et al. 2022; Ross et al. 2024) observed. With Euclid having now launched and beginning to collect science data, further candidate UHzRGs can be identified quickly.

Looking towards the second half of this decade, the Vera C. Rubin Observatory (Ivezić et al. 2019) will have an exceptionally large FoV (9.6 deg $^{2}$ ) and survey approximately 18 000 deg $^2$ in six optical and near-infrared bands across the wavelength range 320–1050 nm. After a first pass of the survey area, the $5\unicode{x03C3}$ point-source sensitivity in Y-band will be approximately 23.0 mag (expected median zenith seeing $\approx 0\hbox{$.\!^{\prime\prime}$}7$ ), with an expected depth of 24.8 mag after 10 yr of operations. Additionally, the Nancy Grace Roman Space Telescope (Spergel et al. 2015; Akeson et al. 2019) will have an unprecedented combination of FoV and sensitivity: the Wide-Field Instrument (WFI) will have a 0.281 deg $^2$ FoV and a limiting $5\unicode{x03C3}$ point-source sensitivity in the F213 filter (2.13 $\unicode{x03BC}$ m) of 25.6 mag (1 h integration; point spread function FWHM $= 0\hbox{$.\!^{\prime\prime}$}169$ ). A combination of Roman and the upcoming Square Kilometre Array (SKA; e.g. Dewdney et al. 2009; Braun et al. 2019) would transform searches for UHzRGs, including those which are much less powerful at radio frequencies than the sources considered in this study.

5. Conclusions

In this paper, we presented deep VLT/HAWK-I $K_\textrm{s}$ -band imaging for 35 of the 53 sources from the B22 HzRG candidate sample, bringing the total number of sources from this sample with deep $K_\textrm{s}$ -band imaging to 41. Our conclusions are as follows:

1. 1. Host galaxies were detected for 27/35 sources (77%) using photutils. The corresponding $K_\textrm{s}$ -band magnitudes in 2 $^{\prime\prime}$ diameter circular apertures range from approximately 21.6–23.0 mag. The remaining eight sources were not detected in the HAWK-I images, with deep magnitude lower limits being obtained (median $3\unicode{x03C3}$ depth $K_\textrm{s} \approx 23.3$ mag; 2 $^{\prime\prime}$ diameter apertures).

2. 2. In the $S_{150\,\rm{MHz}}$ – $S_{2.2\,\unicode{x03BC}\rm{m}}$ flux density parameter space, the distribution of known powerful radio galaxies across $3\lt z\lt 6$ is extended to fainter flux densities by the inclusion of the B22 sample. Like these distant radio galaxies from the literature, the B22 sources have large to extreme radio to near-infrared flux density ratios ( $S_{150\,\rm{MHz}} / S_{2.2\,\unicode{x03BC}\rm{m}}$ $\sim$ $7.0 \times 10^3$ – $3.1 \times 10^5$ ).

3. 3. Five of the eight targets with HAWK-I non-detections, along with J0008 $-$ 3007 from B22 that was not detected in SHARKS, are of particular interest given that they are UHzRG candidates. These six sources have a median $3\unicode{x03C3}$ magnitude limit of $K_\textrm{s} \gtrsim 23.5$ mag in 2 $^{\prime\prime}$ diameter apertures; moreover, the five HAWK-I non-detections have a median $3\unicode{x03C3}$ limit of $K_\textrm{s} \gtrsim 23.8$ mag in $1.\!\!^{\prime\prime}5$ diameter apertures. All of these UHzRG candidates have extreme radio to near-infrared flux density ratios ( $S_{150\,\rm{MHz}} / S_{2.2\,\unicode{x03BC}\rm{m}} \gt 10^5$ ), similar to the three known powerful radio galaxies at $z\gt 5$ (TN J0924 $-$ 2201, GLEAM J0856 $+$ 0223, and TGSS J1530 $+$ 1049).

4. 4. If the host galaxies of these $K_\textrm{s} \gtrsim 23.5$ mag sources have stellar masses $M_\textrm{*} \gtrsim 10^{10.5}$ M $_\odot$ , then there is a distinct possibility that at least some of them lie at ultra-high redshift to be so faint, or at least at $z \gtrsim 4.2$ . They could also perhaps be powerful radio sources at $z \gtrsim 4$ hosted by undermassive ( $M_\textrm{*} \lesssim 10^{10}$ M $_\odot$ ) galaxies, that is, the radio-loud analogues of the population of AGN with large BH to stellar mass ratios being uncovered with JWST. Other possibilities are lower-redshift sources that are extremely obscured by dust or have one-sided radio emission.

5. 5. The majority of the 35 sources studied in this paper have radio and near/mid-infrared properties consistent with IFRSs, many of which have been confirmed in the literature to be at high redshift.

6. 6. Wide and deep ( $K_\textrm{s} \gtrsim 23$ –24 mag) near-infrared surveys are key to efficiently selecting the most promising UHzRG candidates.

For future work, deeper near-infrared observations are needed to identify the host galaxy for each of our UHzRG candidates. Given the tracks in Fig. 3, as well as the (very) faint near-infrared flux densities being measured with JWST for confirmed or candidate obscured AGN in massive galaxies at ultra-high redshift (e.g. Labbé et al. 2023a,b; Akins et al. 2023; Lambrides et al. 2024; Barro et al. 2024), JWST or Hubble Space Telescope (HST; e.g. Bahcall 1986) observations would be the most efficient way of securing an identification for our targets. Additionally, efforts are underway to obtain spectroscopic redshifts for the B22 sample; these data will be presented in a future paper.