1. Introduction
Active galactic nuclei (AGN) are now found regularly in the early Universe from searches for unobscured quasi-stellar objects (QSOs; e.g. Bañados et al. Reference Bañados2016). Indeed, the most distant AGN currently known is a QSO with a redshift of
$z=7.642$
(Wang et al. Reference Wang2021), putting it just
${\sim}670\,$
Myr after the Big Bang. The super-massive black hole (SMBH) powering this AGN is estimated to have a mass of
$1.6\times 10^9\,{\rm M}_\odot$
. Such observations naturally raise the questions of how did SMBHs form so early and/or grow so quickly? The existence of billion-solar-mass black holes in the early Universe leads to tension with current models of their formation, even with the black holes growing at their maximum Eddington rate (e.g. Volonteri Reference Volonteri2012; Smith et al. Reference Smith2017). Did these first SMBHs form from massive seeds? Did they form from small seeds but accrete and/or merge rapidly?
Many AGN are obscured, meaning that the high ultraviolet (UV)–optical luminosity from their accretion disc is hidden by dust and hence are missed by current QSO searches in the optical and near-infrared (near-IR). If obscured and radio-loud AGN, i.e. radio galaxies, can also be found at high redshifts, then a more complete census of AGN activity, and SMBH growth, in the early Universe may be obtained. Specifically, (a) such AGN do not outshine their host galaxies in the optical/near-IR, allowing their hosts to be studied in detail, and (b) the influence on the environment from their powerful jets and inverse-Compton emission (from interaction with the cosmic microwave background) can be investigated. Furthermore, if a suitably bright, radio-loud AGN could be discovered within the epoch of reionisation (EoR;
$6.5\le z\le 15$
), then the distribution of neutral hydrogen (HI) along the line of sight could be measured from absorption studies of the 21-cm line (e.g. Carilli et al. Reference Carilli, Gnedin, Furlanetto and Owen2004; Ciardi et al. Reference Ciardi2015).
Historically (1960–2000), the most distant AGN were frequently discovered from radio surveys since the powerful radio emission acts as a beacon across cosmic time (e.g. Schmidt Reference Schmidt1963; Peterson et al. Reference Peterson, Savage, Jauncey and Wright1982). These surveys found highly accreting but often obscuredAGN, where the UV–optical emission from the accretion disc is hidden by a dusty torus around the SMBH (Urry & Padovani Reference Urry and Padovani1995). The high accretion rates for obscured, type 2 AGN are only detectable via spectropolarimetry (e.g. Vernet et al. Reference Vernet2001), narrow emission lines (e.g. Nesvadba et al. Reference Nesvadba2017) or from mid-IR observations (e.g. Drouart et al. Reference Drouart2014). As discussed above, the most distant AGN known currently are QSOs found from deep and wide optical/near-IR surveys using the Lyman-break technique, whereas searches in the radio are completely independent of the Lymanbreak method.
The discovery of AGN in the distant Universe from radio surveys is undergoing a renaissance due to a new generation of deep, low-frequency radio surveys such as the TIFR GMRT Sky Survey (TGSS; 150 MHz; Intema et al. Reference Intema, Jagannathan, Mooley and Frail2017), the LOFAR Two-metre Sky Survey (LoTSS Shimwell et al. Reference Shimwell2017), and the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM) survey (70–230 MHz; Wayth et al. Reference Wayth2015). Low-frequency radio surveys provide a lever arm to aid the selection of high-redshift radio galaxies (HzRGs;
$z>5$
), which likely exist in such surveys (e.g. Wilman et al. Reference Wilman2008), from the millions of radio sources already catalogued.
Several new HzRGs have been discovered over the last few years, breaking a distance record for powerful radio galaxies that stood for nearly two decades (van Breugel et al. Reference van Breugel1999).Footnote a Saxena et al. (Reference Saxena2018a) used the 1.4-GHz Faint Images of the Radio Sky at Twenty-cm (FIRST) survey (Becker, White, & Helfand Reference Becker, White and Helfand1995) and the National Radio Astronomical Observatory Very Large Array (VLA) Sky Survey (NVSS; Condon et al. Reference Condon1998) in combination with TGSS to search for uHzRGs. Using the ultra-steep spectrum (radio spectral index
$\alpha<-1.3$
Footnote b) technique, they found a radio galaxy, TGSS J1530
$+$
1049, at
$z=5.72$
(Saxena et al. Reference Saxena2018b).
Recently, we have taken advantage of a different selection technique, exploiting the 20-band radio photometry provided by the first extragalactic data release of GLEAM (Hurley–Walker et al. Reference Hurley-Walker2017). We searched for curvature in the radio spectral energy distributions (SEDs) of candidate HzRGs, which may be caused by lower-frequency turn-overs due to synchrotron self-absorption (SSA) or free–free absorption (FFA). With this technique, and searching for sources without near-IR
$K_{\rm s}$
-band detections in the Visible and Infrared Survey Telescope for Astronomy Kilo-degree Infrared Galaxy (VIKING) survey (Edge et al. Reference Edge2013), we identified four candidate HzRGs (Drouart et al. Reference Drouart2020, hereafter D20) in the Galaxy And Mass Assembly (GAMA Driver et al. Reference Driver2011)
$9^{\rm h}$
field. Using deep
$K_{\rm s}$
-band imaging from the HAWKI instrument (Kissler-Patig et al. Reference Kissler-Patig2008) on the Very Large Telescope (VLT) to locate the host galaxy, we obtained 87–117 GHz spectra (and imaging) from the Atacama Large Millimetre Array (ALMA). From the detection of carbon monoxide (CO) lines, we found one source, GLEAM J0856
$\,{+}\,$
0223, to be a powerful radio galaxy at
$z=5.55$
(D20), more than five times as luminous at low frequencies than TGSS J1530
$\,{+}\,$
1049 at a similar redshift.
Furthermore, we found a candidate ultra-HzRG (uHzRG,
$z\,{>}\,7$
), GLEAM J0917–0012, potentially at
$z\sim 10$
based on three weak CO features (with signal-to-noise ratios, SNR, of 2.8–3.4, although a joint fit found a collective significance of
$\sim 5.4\sigma$
). Follow-up spectroscopic observations with both ALMA (targeting the same CO lines) and with the VLA (targeting lower transition lines around
$42\,$
GHz) were unable to confirm this possibility (Drouart et al. Reference Drouart2021, hereafter D21). Detailed analysis of the optical-to-far-IR SED, including the 100-GHz datum, suggested that the host of GLEAM J0917–0012 lies at either
$z\sim 3$
(optical
$z_{phot}\sim 2.2$
) or
$z\sim 8$
. We also identified a ‘companion’ source about one arcsecond to the west of GLEAM J0917–0012.
In order to elucidate the nature of GLEAM J0917–0012 we have obtained spectroscopic and photometric follow-up with the Hubble Space Telescope (HST). We report the results of these Wide Field Camera 3 (WFC3) observations here. As confirming the redshift of a galaxy in the early Universe is a complex endeavour we supplement the results of the WFC3 data with several detailed SED redshift fitting methods. The WFC3 observations and their reduction are reported in Section 2. Other data and photometry used in this work are presented in Section 3. We present the SED fitting methods in Section 4. The results of our analysis are set forth and discussed in Section 5. We conclude this paper in Section 6. In this paper, we use a flat
$\Lambda$
CDM cosmology with H
$_0$
=67.7 km s–1 Mpc–1 and
$\Omega_{\rm M}$
=0.308 (Planck Collaboration et al. 2016). All coordinates are in J2000 and all magnitudes in the AB system.
2. HST observations
GLEAM J0917–0012 was observed with HST under proposal ID 16184 using WFC3 for five orbits in late 2020. Considering the earlier indications that this radio galaxy could lie at
$z\sim 10$
(D20), the first orbit was imaging with the F105W filter (which would have been completely blue-ward of Lyman-
$\alpha$
at that redshift). However, the host of GLEAM J0917–0012 was clearly detected suggesting a lower redshift. Therefore the subsequent four orbits of grism observations were conducted with the F098M/G102 combination of filter/grism (rather than the F160W/G141 combination as originally proposed) to target the observed-frame Lyman
$-\alpha$
. The observations were carried out in two orientations of two orbits each. Unfortunately, one of the orbits was approximately colinear with the galaxy-companion direction which results in over-lapping of their spectral features. The dispersion directions are indicated in Figure 1. We processed the data using two independent software packages to ensure any significant spectral features are not an artefact of any one code.
Grey-scale F098M, F105W and
$K_s$
-band images of the host of GLEAM J0917–0012 (centre of each panel) and the companion galaxy (to the west). Overlaid on the F098M image are contours from the highest-resolution radio data (100-GHz ALMA, red) with the beam shown in the lower left. The contours are
$3,4,5,6\,\times\,10\,\mu$
Jy beam–1, the root mean square (RMS) noise of the ALMA image. The host galaxy shows a compact morphology with a size of
${\sim}0.25\,$
arcsec, and with a faint extension to the west. The two circles in the central panel indicate the
$0.7''$
radius apertures used to measure the flux densities, with the Source Extractor positions given in Table 1. The two dispersion directions of each orientation of the grism are indicated by the arrows in the third panel.
2.1. Reduction with hstaxe
2.1.1. Imaging
We first used the hstaxe v1.0.0Footnote c (Sosey & Pirzkal Reference Sosey and Pirzkal2013) software to process the data. This code is a Python and C successor to the iraf based axe code used for many years on HST data. We based our reduction on the cookbook presented in the Jupyter notebook available with HSTAXE. We used the default values in this notebook unless otherwise stated. The pipeline starts from the ‘flt’ fits files downloaded from MAST.Footnote d The AstroDrizzle routine was used to combine all the G102 grism observations to define a master pixel grid. The ‘full-depth’ F098M and F105W images were also created with AstroDrizzle from all available exposures using this master pixel grid. We used the combine_type=‘median’ and combine_nhigh=1 as recommended for more than three exposures. The pixel size of the ‘drizzled’ image matches that of the WFC3 detector (
${\sim}0.13\,$
arcsec) by default.
To detect and measure the flux densities of the host and companion, Source Extractor (Bertin & Arnouts Reference Bertin and Arnouts1996) was run on both images. We used
$0.7\,$
arcsec radius apertures to measure the fluxes and applied an aperture correction (see Section 3 for full details of the photometry). We used the detections in the F105W band (where the host galaxy has the highest SNR) to define the centres of these apertures. The only other difference from the HSTAXE cookbook Source Extractor parameters was that we used a lower DETECT_MINAREA of 4 (rather than the default of 10), and lower values of the parameters DETECT
_THRESHOLD and ANALYSIS
_THRESH (1.5 rather than the default of 3). We report the Source Extractor positions of both galaxies in Table 1.
F105W Source Extractor coordinates of the host of GLEAM J0917–0012 and the companion.
2.1.2. Grism spectra
We used the aXe Flexcube extraction method to create the 1D and 2D spectra. This method takes advantage of the mosaics at different wavelengths to aid in the flux calibration and spectral extraction. We used the following steps for both orientations separately and for all the data combined:
-
1. We used the iolprep aXe task to generate separate F098M catalogues derived from the full-depth catalogue and matched to the geometry of each F098M pre-imaging exposure (note the full-depth catalogue is created with the F105W image as a detection image).
-
2. We created a file, aXe.lis, with each row containing a G102 grism exposure name, the expected F098M catalogue names (created in the previous step) and the name of the associated F098M direct-preimage file.
-
3. We created a cubelist.lis file with each row containing the name of the full-depth images, the pivot wavelength of each filter (i.e. PHOTPLAM in nm), and the AB magnitude zero-points.
-
4. We then created the ‘flexcube’ model (with axetasks. fcubeprep), which is a combination of the mosaics and segmentation file into a master ‘FLX’ file for each full-depth image and an ‘FLX’ file for each G102 observation. These files can be used to compute the contamination of spectra and to perform the extraction.
-
5. We ran aXeprep on the data, which takes the aXe.lis file, the default grism configuration file (in this case G102.F098M.V4.32.confFootnote e) and subtracts the background. Also specified at this stage is the mfwhm parameter, which determines the extent of the area that is masked perpendicular to the trace of each object before the background is determined. We used mfwhm
$=1$
(rather than the default of 2), which is better for a lower SNR and as we have many sources with Source Extractor from the lower detection threshold (hence the field is more crowded). -
6. We ran axecore using the flexcube models, which takes the aXe.lis and configuration files. We also specify at this step the precise method to extract the spectra. In this case we turned on the titled extraction (orient
$=$
True) and slitless geometry (slitless_geom
$=$
True). We also used values of extrfwhm
$=2$
and drzfwhm
$=1$
, smaller than the defaults. See the hstaXe documentation for more details.Footnote f
-
7. We then used drzprep with optimal extraction to prepare the individual spectra from each G102 exposure for combining.
-
8. Finally, we used the axecrr task to perform an optimal extraction and drizzling of the 1D and 2D spectrum with the same aXe.lis and configuration files. This task combines all the data with default cosmic ray rejection (CRR) as provided in the configuration files. We also used values of infwhm and outfwhm equal to the values for extrfwhm and drzfwhm used with axecore (respectively). In contrast to the default pipeline, the uncertainty on each spaxel of the spectrum is determined from the RMS of the full 2D spectrum in the cross-dispersion direction then scaled to the appropriate flux by the sensitivity file: WFC3.IR.G102.1st.sens.2.fits.Footnote g
The final 1D and 2D hstaxe spectra for GLEAM J0917–0012 are presented in the top panel of Figure 2. The estimated contamination (indicated in the 1D spectra) is subtracted from both and only shown for reference. The contamination is only seen in one orientation (Or1), given it is approximately coaligned with the host galaxy/companion vector (see Figure 1). The final 1D and 2D hstaxe spectra for the companion are presented in the top panel of Figure 3, likewise with the contamination subtracted.
The hstaxe 1D and 2D (top), and grizli 1D (bottom) spectra of GLEAM J0917–0012. These are the combined spectra from both orientations. In each orientation the width of the extraction is twice the size of the source (perpendicular to the trace) of the target in the F140W image. Both reductions show the estimated contamination from other sources (indicated in red), which is subtracted from the 1D and 2D flux. This contamination occurs from the companion galaxy and only affects orientation 1. Grey shaded regions indicate where the transmission of the G102 grism drops below 10%. The spectrum is consistent with a very low or zero flux across most of the wavelength range. The hstaxe 1D spectrum shows a strong feature at
$1.15\,\mu$
m that is not seen in the grizli reduction which seems to an artefact. There is a fainter, but broader, feature at
$1.12\,\mu$
m seen in both reductions (and found with the grizli line finding routine, albeit at a SNR of
${\sim} 3$
). The grizli PDF from the template fitting (Figure 6) shows numerous maxima. We overlay the continuum (black) and line (orange) model template which is most consistent with our other SED constraints (see Section 4). This solution at
$z=8.21$
shows the Lyman-
$\alpha$
line in emission and the NV
$\lambda\lambda 1238,1242$
line in absorption. We discuss what these features in detail in Section 5.
As Figure 2, but for the companion galaxy. The grizli PDF from the template fitting (Figure 6) shows numerous maxima. We overlay the continuum (black) and line (orange) model template which is most consistent with our other SED constraints (see Section 4). This solution is
$z=1.63$
with the central line at
$0.98\,\mu$
m identified as the [OII]
$\lambda 3727$
doublet. We discuss the likelihood of this solution in conjunction with the SED fitting in Section 5.
2.2. Reduction with grizli
We also reduced the HST data using the Grism redshift & line (grizli) analysis software (Brammer Reference Brammer2019). This code is designed to offer a single tool to perform all the stages from the reduction to the spectral characterisation of current HST and future James Webb Space Telescope (JWST) Footnote h slitless spectroscopic observations. We used grizli version 1.0.dev1365 and followed the standard reduction steps, contamination modelling, and spectral template fitting for redshift determination, as are outlined in the official documentation. We now provide a summary of the steps we carried out:
-
1. Starting from the raw flt files from MAST, we first used the grizli function parse_visits, with default parameters, to associate the individual exposures into visit groups taken in the same filter and position angle. The grism and corresponding direct imaging association was also performed at this step.
-
2. We then used the preprocess function which performs the flat-fielding, pixel flagging (CRR and bad pixels masking with AstroDrizzle, and persistence masking), relative astrometric alignment, sky background subtraction, and drizzling (with AstroDrizzle) of all the direct imaging and grism visits. We only used default grizli parameters in this multi-step preprocessing and refer the reader to the official documentation for further details on these steps.
-
3. We then performed an absolute astrometric alignment of the individual exposures to the Gaia (Gaia Collaboration et al. 2016) data release two (DR2; Gaia Collaboration et al. 2018) with grizli’s fine_alignment function and default parameters.
-
4. After the fine alignment step, we created a F098M + F105W direct image mosaic with grizli’s drizzle_overlaps with a final mosaic pixel fraction of
$\texttt{pixfrac}=0.75$
, and we also masked diffraction spikes on the mosaic image with mask_IR_psf_spike for sources brighter than 17 mag AB. We then performed source detection on the mosaic image with multiband_catalog which uses the Source Extractor Python wrapper sep (Barbary 2016). We used default grizli parameters, including a detection threshold (threshold) of
$1\sigma$
above the global background RMS, a minimum source area in pixels minarea of 9, and a deblending contrast ratio deblend_cont and a number of deblending thresholds deblend_ntresh of
$0.001$
and 32, respectively. Matched-aperture photometry in the F098M and F105W bands was simultaneously performed during this step. -
5. The last step before spectral extraction is the continuum modelling of sources in the frames for contamination estimate and removal. For this step, we used grism_prep and multifit.GroupFLT as recommended in the documentation, using the mosaic catalog as the source detection catalog. The continuum modelling is an iterative process done in two steps. In the first step, a linear continuum model of all sources in the field down to the default mag_limit of 25 mag is created. Then, grizli refines the models by fitting higher-order polynomials (here
$\texttt{poly_order} = 3$
) to all sources iteratively in descending brightness, here from 18 to 24 mag, after subtracting potential contamination from neighbouring objects. grizli iterates refine_niter times over the refined modelling, with
$\texttt{refine_niter}=3$
here. -
6. We then extracted the 2D cutouts of the host of GLEAM J0917–0012 and of the companion source with extract and default parameters, and redshift-fitted the spectra with fitting.run_all_parallel in the range
$z=0.05-10$
with the default set of FSPS (Flexible Stellar Population Synthesis; Conroy et al. Reference Conroy, Gunn and White2009, Conroy & Gunn Reference Conroy and Gunn2010) and emission line templates. The advantage of the grizli redshift-fitting routine compared to traditional redshift fitters is that it fits the spectra in the native 2D frames by modelling the 1D spectral templates into the 2D grism space, taking into account source morphology and spectral smearing of pixels that are adjacent in the direct imaging.
Grizli provides the full redshift probability density function (PDF) of the fitting as well as the 1D continuum and line emission models. We show the 1D grilzi spectra and contamination estimates in the lower panels of Figures 2 and 3. Overlaid are templates selected from maxima in the PDF (Figure 6) which best correspond to the spectral line observed in each spectra and the other photometric constraints on the redshift (see Section 5).
3. Other data and photometry
3.1. Data
We compiled UV, optical, near- and mid-IR photometry from the literature in addition to the HST imaging presented here and the
$K_{\rm s}$
-band HAWKI imaging from D20 (see Table 2 for observing details of those bands). We used the two UV bands from the GALaxy Evolution eXplorer (GALEX) Medium Imaging Survey (MIS; Martin et al. Reference Martin2005), the grizy-band images from the Hyper Suprime-Cam (HSC) Subaru Strategic Program (Aihara et al. Reference Aihara2019), the JH-band imagesFootnote i from the VIKING survey, and mid-IR images from the Widefield Infrared Survey Explorer (WISE;Wright et al. Reference Wright2010).
Observing log of exposures with significant detections.
From the HSC data both the host and companion have catalogued flux densities in all five bands, although the significance of the host flux densities varies from 1.5 to 4.6
$\sigma$
(where
$\sigma$
is the RMS noise). In D21, based upon visual inspection of the HSC images presented in that work, we reported the host as a non-detection. Both source extraction and visual inspection of the images gave no hint of an obvious detection. Additional analysis of the images below (Section 3.2) finds only weak detections (
${<}3\,\sigma$
). The detection reported in the HSC data release is a ‘child’ detection of the significant detection of the companion in
$i-$
band only and is most likely erroneous (private communication, HSC team).
We used reprocessed GALEX, VIKING and WISE images in a similar fashion as described in the GAMA Panchromatic Data Release (GAMA PDR; Driver et al. Reference Driver2016). The GALEX MIS survey provided deeper images of most of the GAMA fields compared to the all-sky survey. The individual MIS frames were combined together (using swarp; Bertin et al. Reference Bertin, Bohlender, Durand and Handley2002) to form one image for each GAMA field. The resulting point spread functions (PSFs) were 4.1 arcsec in the far-UV (FUV, 1 350–1 750 Å) and
$5.2\,$
arcsec in the near-UV (NUV, 1 750–2 800 Å) channels. For the VIKING data, all individual frames with the same filter were rescaled onto a common zero-point magnitude (30) before being drizzled together with a
$0.339\,$
arcsec pixel size. Unlike the original images from GAMA PDR, these images were not convolved down to two arcsecs, but left with a mean
$0.85\,$
arcsec seeing. The WISE data comprise four bands (W1, W2, W3, W4) at 3.4, 4.6, 12 and
$22\,\mu$
m respectively with PSFs given in Table 3. These images were also reprocessed in a similar fashion to GAMA PDR (see Bellstedt et al. Reference Bellstedt2021, for full details of this version of the GAMA PDR).
Photometric properties of the UV to near-infrared and ALMA 100-GHz images of GLEAM J0917–0012, as well as the companion galaxy. We list each band along with its effective central wavelength (
$\lambda_0$
), the AB zero point (ZP), and the seeing/resolution (
$\theta_{\rm FWHM}$
- note that we provide the restoring beam parameters for the ALMA image. Where relevant, we provide the
$0.7\,$
arcsec aperture correction and uncertainty correction factor, C (see Section 3.2). We then list the aperture-corrected flux densities and uncertainties, or
$3\sigma$
limits, of the host of GLEAM J0917–0012 and its companion (see Section 3). Uncertainties include an additional 10% for the combined uncertainty in the absolute flux calibration and aperture correction.
In this work we also used the results of the far-IR-to-radio SED presented in D21. We do not repeat all those data here, but note that we have photometry across
$0.07{-}100\,$
GHz from numerous radio facilities (described in D21) in addition to
$3\sigma$
upper limits from the Herschel Space Observatory (Pilbratt et al. Reference Pilbratt2010) across 100–500
$\mu$
m (see D21). However, we do present the ALMA 100-GHz flux density in Table 3 as it plays a key role in constraining the photometric redshifts (see Section 4). We note that the ALMA image presented in Figure 1 is the most sensitive produced by giving the data natural weighting, but our best estimate (D20) of the size of the radio source at 100 GHz is
${\le}0.7\,$
arcsec from an image with Briggs (robust
$=0$
) weighting. The overall radio SED fitting in D20 and D21 did not suggest any variability in the radio and the
$1.4\,$
GHz flux densities from NVSS and FIRST are consistent within the uncertainties as are the two epochs of data from the VLA Sky Survey (Lacy et al. Reference Lacy2020). Hence, we conclude that this radio source is unlikely to be variable or to be beamed.
3.2. Photometry
We tied all 14 optical-to-mid IR (0.48–22
$\mu$
m) images to the Gaia DR2 reference frame using a catalogue obtained from the Gaia archive.Footnote j Catalogues derived from Source Extractor were cross-matched (
$\le 1\,$
arcsec) to the Gaia catalogue and the derived mean right ascension and declination offsets (all
${<}0.2\,$
arcsec) were added to the FITS sky position reference values of each image (i.e. CRVAL1 and CRVAL2 respectively). Note the offsets were less than
$0.05\,$
arcsec for HST and HSC so were not applied. As GLEAM J0917–0012 was not detected in the GALEX images we did not shift them. For the VIKING and WISE images we used the mean offsets of the four or five bands in each set, as these data have already been put onto an internally self-consistent reference frame.
Using the
$0.7\,$
arcsec radius apertures defined by the positions in Table 1, we measured the flux densities, and uncertainties, in all bandsFootnote k for the host of GLEAM J0917–0012 and the companion galaxy. The fluxes and uncertainties were directly converted from ADU to
$\mu$
Jy using the zero points and we applied an aperture correction to account for flux outside of the
$0.7\,$
arcsec aperture. The zero points and aperture corrections are reported in Table 3 along with PSF (
$\theta_{\rm FWHM}$
), the flux densities measured and their uncertainties. To determine the aperture corrections for WFC3, we took the inverse of the nearest set of values of the fractional flux enclosed in Table 7.7 of the WFC3 instrument handbookFootnote l and interpolated them linearly using the
$0.7\,$
arcsec radius and central wavelength of each filter. For the HSC, VIKING and HAWKI data we empirically measured the curve of growth for nearby non-saturated stars to determine the aperture corrections.
The Source Extractor uncertainties are derived from the local RMS multiplied by the square root of the aperture area (in units of pixels). The resampling by swarp affects the RMS through the algorithms used in the regridding process. To ensure we have accurate uncertainties, we have to apply a correction factor, C, to the uncertainties from the regridded images. On the assumption that the uncertainties should be the same for the original and regridded images, one can show that
$C=\frac{\sigma_{\rm o}}{\sigma_{\rm r}}\times f$
where
$\sigma_{\rm o}$
is the pixel rms in the original image,
$\sigma_{\rm r}$
is the rms in the resampled image, and f is the scale factor by which the one-dimensional length of the square pixels has decreased by, i.e.
$\theta_{\rm o}/0.12825$
where the
$\theta_{\rm o}$
is the one-dimensional pixel length before regridding. We add an additional 10% uncertainty, in quadrature, to all flux density measurements to account for the uncertainty in the absolute flux density scale and in the aperture corrections.
In the case of GALEX and WISE, we provide only
$3\,\sigma$
upper limits (where
$\sigma$
is the RMS noise) as GLEAM J0917–0012 and the companion are clearly not detected and would be blended if they were (and in the case of GALEX, we did not correct to the Gaia reference frame). To determine the RMS for WISE we used the RMS estimated by Source Extractor on the original image (i.e. before regridding, to avoid adding any correlated noise). These values were converted to AB magnitudes and then to
$\mu$
Jy. Determining the limit to the GALEX sensitivity was more complex than the other bands due to the low photon statistics. Indeed, the majority of pixels in the original FUV images contain no photons. Hence, we used Profound (Robotham et al. Reference Robotham2018) with specially tuned parameters to determine the RMS over an area with enough statistics and using mean not median statistics.
4. SED fitting and photometric redshifts
In addition to the spectroscopy and template fitting from grizli, we used three broadband SED fitting methods to estimate the redshifts of GLEAM J0917–0012 and the companion. These are (a) a bespoke UV-to-radio SED fitting method including the measurements and the upper limits, (b) EAZY on the optical-to-near-IR detections, and (c) RAiSERed on the radio photometry and imaging.
4.1. UV-to-radio SED fitting
Given the broad UV-to-radio wavelength coverage of GLEAM J0917–0012 and building on our analysis in D21, we have developed a bespoke tool to simultaneously fit galaxy templates plus an analytical function to the SED. In this case the analytical function is a triple power-law (TPL) to fit the synchrotron emission dominating the radio SED at
$0.07{-}\!100\,$
GHz, whereas the galaxy templates cover the UV to far-IR. Note the template and power-law can both contribute to the 100 GHz datum which motivates our approach of combining a template with an analytical model. For each class of galaxy template we determine the best fit of the combined template plus a TPL at 100 linearly spaced redshifts across
$0<z<10$
.
We employ the following templates from PÉGASE (v3; Fioc & Rocca-Volmerange Reference Fioc and Rocca-Volmerange2019, as used in D21) with the addition of some Lyman-break galaxy (LBG) ones:
-
1. An elliptical galaxy template.
-
2. A starburst template with internal dust extinction and assuming a slab geometry.
-
3. A dusty starburst template, i.e. the same starburst template but with
$\alpha_{\rm slab}\,{=}\,10$
, corresponding to a column density ten times denser.Footnote m
-
4. The 30 empirical LBG templates from Álvarez-Márquez et al. (Reference Álvarez-Márquez, Burgarella, Buat, Ilbert and Pérez-González2019) which cover a range of UV luminosities, UV continuum slopes and stellar mass.
All these templates cover from below the Lyman limit to the far-IR (
$0.01<\lambda_{\rm rest}<10\,000\,\mu$
m). For the PÉGASE templates we take 60 logarithmically spaced ages across 0–3 500 Myr, although exclude redshift/age combinations which exceed the age of the Universe. Hence, for each galaxy template (plus TPL) we calculate the highest likelihood, and hence best-fit, for all 6 000 redshift/age combinations. For the 30 LBG templates, the parameter space explored covers 3 000 redshift/template combinations.
For the TPL, we use Equation 4 from D20. The initial parameters to the TPL were those reported in Table 3 of D21, but with the low-frequency spectral break and low-frequency slope fixed. Importantly this allows the parameters in the TPL which effect the 100 GHz datum to vary as required in order to accommodate the galaxy template. Hence, for each of the 9 000 fits, the parameters varied are the template normalisation plus four out of six TPL parameters.
When fitting the data we use the same likelihood calculations as used in the MrMoose
Footnote n SED fitting code (Drouart & Falkendal Reference Drouart and Falkendal2018). This formalism has the advantage of being able to use both flux measurements, even at low or negative SNR, and flux upper limits, as we use in the lower resolution bands. The full details of this approach are explained in the appendix of Sawicki (Reference Sawicki2012), but, briefly, the minimised
$\chi^2$
is the sum of the traditional
$\chi^2$
statistic for the flux measurements and an additional term for the upper limits. This additional term accounts for the probability that the observation in a given band is drawn from a given model. As such, it is an integral from minus infinity to the user defined upper limit. Here we choose
$3\sigma$
as the upper limit although the fitting is not very sensitive to this limit. Effectively all the upper limits contribute to the logarithm of the likelihood (analogous to
$\chi^2$
if using detections only), i.e. the contribution to the likelihood progressively penalises models that lie above the
$1\sigma$
RMS noise and up to the predefined limit (here
$3\sigma$
) and favours those models that lie below the limits (see Drouart & Falkendal Reference Drouart and Falkendal2018, for more details).
We present a selection of the best-fitting SEDs to GLEAM J0917–0012 for each template class in Figure 4. We fit the models in flux density space, but the resultant SEDs are plotted in power, i.e.
$\nu\,F_\nu$
, to better illustrate the SED over almost eight decades in wavelength. The inserts show the variation of likelihood across the age/redshift parameter space (template/redshift parameter space for the LBG templates). The SEDs plotted are from representative local maxima in the likelihood distribution. We also use the same code, minus the analytical function for the synchrotron emission, to estimate the redshift of the radio-quiet companion galaxy (see Figure 5).
Observed UV-to-radio SED of GLEAM J0917–0012 including measurements (black diamonds with
$1\sigma$
errorbars) and
$3\sigma$
upper limits (triangles) overlaid with a selection of best-fitting models. Each panel corresponds to a different galaxy template class and the radio data is simultaneously fit with a synchrotron triple-power-law model (see text for full details). The insets show the distribution of likelihood values for each redshift/age combination (redshift/template combination for the LBG templates). Omitted are solutions which exceed the age of the Universe (i.e. the top right of the insets for the evolving PEGASÉ templates). The curved white line indicates when the first galaxies are thought to have formed (e.g. Laporte et al. Reference Laporte2021), around 250 Myr after the Big Bang. The symbols in the insert, with their log likelihoods labelled correspond to local maxima from which the best-fit templates are taken (plotted in the same colour). This modelling suggests an old, low-redshift (
$z\sim2$
, in magenta), or a high-redshift (
$z\sim 7.9-8.5$
, in green and blue) solution is most preferred.
As Figure 4, but for the companion galaxy. As this source is radio-quiet the templates are fit without a synchrotron model. Broadly speaking, the best solutions are at
$z< 3$
(magenta).
4.2. Optical-to-near-IR SED fitting with EAZY
Using the optical-to-near-IR detections we can implement the more traditional EAZY (version 2015-05-08) photometric redshift estimator (Brammer, van Dokkum, & Coppi Reference Brammer, van Dokkum and Coppi2008). Unlike the bespoke fitting above, EAZY cannot use the lower resolution data where we only have upper limits. Hence, we ran EAZY on the host galaxy and the companion using the data in Table 3 bar the lower-resolution UV and mid-IR data. We used the default templates provided and the redshift priors (linked to the
$K_{\rm s}$
-band magnitude in this case) although they make little difference to the best-fit redshift. We apply the intergalactic medium absorption option as we allow the estimated redshift to cover
$0.01<z<10$
. We report the marginalised redshift estimates with the prior, although the prior only has a weak effect on the final result.
4.3. Radio continuum redshift fitting with RAiSERed
Turner et al. (Reference Turner, Drouart, Seymour and Shabala2020) proposed a dynamical-model based method, called RAiSERed, for determining the redshift of a radio galaxy using exclusively radio-frequency imaging and photometry. In this code, the breadth of radio galaxies permitted by the physics of the Radio AGN in Semi-analytic Environments (RAiSE; Turner & Shabala Reference Turner and Shabala2015; Turner et al. Reference Turner, Rogers, Shabala and Krause2018a) dynamical model are compared to the properties of a given radio galaxy being investigated, yielding a PDF for its redshift. Specifically, one can compare the lobe flux density, angular size and broadband spectral shape of the dynamical model and observed source. The local environment of the radio galaxy is modelled as a prior PDF based on the local cluster mass function (Girardi & Giuricin Reference Girardi and Giuricin2000), modified for higher redshifts using semi-analytic galaxy evolution models (SAGE, Croton et al. Reference Croton2016; Raouf et al. Reference Raouf, Shabala, Croton, Khosroshahi and Bernyk2017). Hence with an increasingly dense environment at higher redshifts, degeneracies with redshift in the model are avoided.
We apply this method to GLEAM J0917–0012, taking the flux density of each lobe at 151 MHz as
$0.233\pm0.014$
Jy (half the total flux density presented in D20) whilst the angular size of each lobe is
${<}0.35$
arcsec.Footnote o The unknown lobe axis ratio is assumed to be consistent with that of Cygnus A (Turner & Shabala Reference Turner and Shabala2019) as a representative powerful radio galaxy.
The broadband spectral shape across the 0.07–20 GHz radio data may be described by two different models:
-
1. lobe-dominated
-
2. young-jet plus aged-lobe
The lobe-dominated model assumes a population of synchrotron-emitting electrons shock-accelerated uniformly in time since the commencement of jet activity, and which presently have an approximately constant magnetic field strength. The resulting spectrum is well described by the continuous injection of relativistic electrons and characterised by a power-law at low frequencies that steepens beyond an (optically thin) break frequency by
$\Delta \alpha = -0.5$
Footnote p for active sources. The initial power-law was fit in D20 and D21 as
$S_\nu \propto \nu^{-0.8\pm0.1}$
with the optically thin spectral break at
$1.6_{-0.6}^{+0.9}$
GHz. However, at any redshift, the lobe-dominated model poorly explains the steep spectral index within the MWA band (100–150 MHz) of
$S \propto \nu^{-1.0\pm0.1}$
, albeit is consistent with the low-frequency spectral index of the continuous injection model at the 2
$\sigma$
level.
For high-redshift radio sources, the much stronger cosmic microwave background radiation will rapidly deplete the energy of the synchrotron electron population in the lobe through inverse-Compton scattering. Hence, the jet contribution to the integrated luminosity will be important. Therefore, we also fit a young-jet plus aged-lobe model. Electrons in jets are continually accelerated at shocks, both internal (such as recollimation shocks) and termination shocks (i.e. hotspots), resulting in power-law spectra characteristic of the particle energy spectrum at injection (e.g. Bicknell et al. Reference Bicknell, Mukherjee, Wagner, Sutherland and Nesvadba2018). These spectra age rapidly on cessation of particle acceleration, on timescales
$<1$
Myr at GHz radio frequencies, due to the high milli-Gauss magnetic fields in the jet. Many powerful radio galaxies are intermittent (e.g. Bruni et al. Reference Bruni2019), and we therefore use only the low-frequency (
$<1.4\,$
GHz) radio data to fit the jet injection spectrum. The higher-frequency observations are therefore modelled assuming either a continuous injection or a Jaffe & Perola (Reference Jaffe and Perola1973, JP) spectrum.
Hence, the second model considers the integrated luminosity for such sources using a two-component spectrum: a young jet plus aged lobe. The young jet is assumed to have only recently accelerated synchrotron-emitting electrons, but is described by a JP continuous injection spectrum due to the extreme energy losses present at high redshift. These losses are due to a combination of the high magnetic field strength and strong inverse-Compton field which varies with redshift. Meanwhile, the aged lobe is described by a continuous injection model but with a break frequency well below the lowest observing frequency, i.e. it has a power-law break
$\Delta\alpha=-0.5$
steeper than the jet. This two-component spectrum is fitted assuming arbitrary flux density scaling for the two components and the same initial spectral index upon shock acceleration of synchrotron-emitting electrons. We constrain the initial power-law spectrum of both components as
$S_\nu \propto \nu^{-0.8\pm0.1}$
and the break frequency of the aged-lobe component as
$<100$
MHz; this model fits a spectral index at 100 MHz of
$S_\nu \propto \nu^{-1.0\pm0.1}$
, consistent with the value fitted across the MWA band. The results of this fit are shown in Figure 7.
Grizli redshift probability distribution function divided by the
$\chi^2$
, up to
$z=9$
for the GLEAM J0917–0012 and up to
$z=7$
the companion. Both distributions have many maxima occurring when one line or another is matched to the lines detected in the host and companion. In both cases we highlight (by a red vertical line) the solution which is most consistent with the SED fitting methods and is from a strong line. For example, we discount the
$z=4.14$
maxima for the companion as this would mean we have detected the [CIII]/
$1909\,\unicode{x212B}$
line which is very weak. Hence, our best solution for the host galaxy is
$z=8.21$
(Lyman-
$\alpha$
) and for the companion
$z=1.63$
([]CII]).
Results of the RAiSERed radio SED fitting. In the left panel we show the two best-fit models: a lobe model at
$z\sim 2.7$
and a young-jet plus aged-lobe model at
$z\sim 7$
. For the jet-plus-lobe model we show the two components separately, with the intermittent jet being the component dominant at higher frequencies. In the right panel we present the PDF for the redshift of the young-jet plus aged-lobe model with the smoothed distribution in black peaking for all redshifts greater than
${\sim}7$
. This fit should be considered an upper limit as we have assumed the jet component is symmetric when it may not in fact lie in the plane of the sky (see Section 5.3.3 for more details). The green dot-dashed line assumes a lower break frequency of 30 MHz, whilst the two red lines increase/decrease the density in the ambient media assumed in RAiSERed by a factor of two. Note we do not fit the low-frequency turn-over as it has no impact on the RAiSERed results.
We further consider the possibility that GLEAM J0917–0012 is a remnant radio source. The RAiSERed model fitting currently does not allow for an inactive source, so this approach cannot constrain the redshift, but only determine if this is a feasible model for the radio spectrum. We start by fitting the spectrum with the Komissarov & Gubanov (Reference Komissarov and Gubanov1994) model. This model yields a comparable power-law spectrum at low frequencies, but which steepens earlier than the lobe-dominated model for an active source such that it is
$S \propto \nu^{-0.9\pm0.1}$
at 100 MHz. The off-time is fitted as
$13\pm3$
% of the source total age under the assumption it is a remnant source. We discuss why this model is unlikely astrophysically in Section 5.3.3.
5. Results and discussion
5.1. Imaging and photometry
The host of GLEAM J0917–0012 is detected and resolved in our F105W image (Figure 1) and is detected less strongly in the F098M band (Table 3). However, even with the prioritised fitting the host of GLEAM J0917–0012 remains undetected in the HSC (shortward of the HST data) and VIKING bands. The HSC and HAWKI flux densities reported in Table 3 for the companion are consistent with those from D21 to within 10%. For the host, the
$K_{\rm s}$
-band flux density is found to be 16% lower (including the aperture correction on the flux density reported here) than that found by D21 which is most likely due to the matched aperture photometry used to measure it (although it is consistent within
$2\sigma$
).
The host is mainly compact at
$0.25''$
across, but with a faint extension to the west towards the ‘companion’ (believed to be at a different redshift; see Section 5.6). The ALMA 100-GHz continuum contours overlaid in Figure 1 clearly identify the eastern source as the host (as determined in our original analysis; see D20 and D21). This emission is potentially aligned with the faint extension to the west, but the ALMA resolution makes this inconclusive. The UV-to-radio SED fitting (Section 5.3.1) suggest that there could be both synchrotron and cold dust contributing to this emission.
The companion source is more extended (
$\sim 0.3''\times 0.4''$
) and brighter in the F105W band. It is also detected in all five HSC bands including g-band, ruling out that this source is above
$z\sim 4.2$
(using the Lyman limit). The flux densities measured here are consistent with the catalogued values from HSC to within 8%.
5.2. Spectral properties and grizli template fitting
5.2.1. GLEAM J0917–0012
The hstaxe and grizli spectra (Figure 2) both show a faint, broad emission feature at
$1.12\,\mu$
m and no discernible continuum. The grizli SED fitting routine estimates a SNR of
${\sim}3$
for this feature (see Table 6). The hstaxe extraction also shows a brighter, narrow feature at
$1.15\,\mu$
m, but this is not seen in the grizli reduction. We do see it in both orientations of the hstaxe reduction but it is likely the combined effect of cosmic rays and/or hot pixels. Furthermore, this feature occurs on the part of the G102 transmission where the throughput is decreasing very rapidly towards longer wavelengths hence we do not consider it a reliable feature.
As well as producing spectra, grizli fits templates at a range of redshifts in order to identify any features and confirm the redshift. One output from this analysis is the PDF. We show in the left panel of Figure 6 the PDF divided by
$\chi^2$
for the host. Numerous maxima occur at different redshifts, but we highlight the one we believe to be most likely at
$z=8.21$
(if the 1.12-
$\mu$
m feature is Lyman-
$\alpha$
) considering the strength of the identified line and the constraints from the SED fitting (Section 5.3). Curiously, the NV
$\lambda\lambda 1238,1242$
line is found in absorption but we do not give this much credence due to the low SNR at that wavelength. The other most likely solution is
$z=2.01$
(if the line is the [OII]
$\lambda 3727$
doublet).
5.2.2. Companion
The hstaxe and grizli spectra (Figure 3) both show a faint, broad emission feature at
$0.98\,\mu$
m and no discernible continuum. The PDF over
$\chi^2$
from the grizli template fitting (right panel of Figure 6) also shows numerous solutions. The most probable solution from a naive examination of the grizli PDF is
$z=4.14$
(which depends on the
$0.98-\mu$
m feature being the C[III] line at
$\lambda_{\rm rest}=0.1909\,\mu$
m). Given the relative weakness of this line we think this solution to be unlikely and discuss the other possibilities in Section 5.6.
5.3. SED fitting
We used three methods to estimate the redshift from the broadband photometry. The results from EAZY and RAiSERed are best summarised numerically in Table 4 where the uncertainties or limits (where available) are given as the
$68{\rm th}$
-percentile.
Constraints on the redshift of the host and companion using the EAZY and RAiSERed (which is only possible for the host). We also show for reference the results of EAZY applied to the companion in D21. The uncertainties are the
$68^{\rm th}$
-percentile confidence intervals from the resulting redshift distributions.
a From the young-jet plus aged-lobe model.
5.3.1. UV to radio SED fitting
For the host galaxy, the results of this fitting (Figure 4) for the different PEGASÉ templates (top three panels are broadly consistent. The insets in these figures show the most preferred redshift solutions and there are two consistent maxima at
$z\sim2$
and
$z\sim 8$
. The low-redshift solution is consistent with an older galaxy (with an age of several hundred Myr. The high-redshift solution is quite broad,
$7.5<z<8.5$
, and generally prefers a young age of
${\sim}10\,$
Myr. Overlaid on the data in Figure 4 are the best-fit templates at
$z=2.4, 7.9$
and
$8.5$
which are all comfortably consistent with the photometry. These results also agree well with the redshift constraints in the optical-to-IR SED fitting of D21 which found a low (
$z\sim2$
) and high (
$z>7$
) redshift solution.
As seen by the colour bar for the parameter space ‘heatmap’ inserts, while there are maxima as highlighted above, the overall variation in likelihood is not large. Hence, overall no particular redshift is strongly preferred by this method.
The LBG templates (Figure 4, bottom panel) are broadly consistent with the lower-redshift solution
$z<3$
. However, the best-fit template does not appear to be good representation of the data.
The synchrotron component is not significantly affected by the fitting of the different template classes. Only the high-frequency slope and high-frequency spectral break change with differing template contribution to the 100-GHz datum. However, the best-fitting results suggest that this far-IR datum has a significant cold dust contribution. The frequency at the low-frequency turn-over (in
$S_\nu$
space) and low-frequency slope are fixed from a fit to just the 0.07–100 GHz data. The observed frequency of the turn-over is constrained to be
$<50\,$
MHz with a 90% confidence although it appears higher in the figures as we plot in
$\nu\,S_\nu$
(rather than
$S_\nu$
).
For the companion galaxy the results of this fitting (Figure 5) for the different templates are also broadly consistent between the top three panels. These solutions generally favour a ‘low’-redshift (
$z<3$
) result compared to the potential high-redshift solution for the host galaxy. In these panels, there is a broad distribution of high likelihoods across
$0<z<3.5$
. A local maxima is highlighted in Figure 5 with the corresponding template overlaid on the data. Likewise with the LBG templates there is broad region of high likelihood. These results on their own do not further constrain the redshift of the companion.
5.3.2. EAZY
We are able to run EAZY on both the host and companion using only the optical to near-IR measurements. The results are reported in Table 4. For the companion, EAZY finds a low-redshift solution of
$z\sim 2.6$
, compatible with the EAZY result from D21 within the uncertainties. For the host galaxy of GLEAM J0917–0012, EAZY finds a result (
$z\sim 2.4$
) albeit from a higher fraction of low SNR measurements. We could not use EAZY on the host in D21 with only one detection.
5.3.3. RAiSERed
The RAiSERed fitting presents us with two different models with which to fit the radio emission of GLEAM J0917–0012 and hence estimate the redshift: a lobe-dominated and a young-jet plus aged-lobe model. We present our results in Table 5, as well as the spectral fitting results for a remnant model (no redshift constraint is available as RAiSERed assumes active sources). The lobe-dominated model constrains its redshift to be
$z = 2.74\pm0.24$
. However, this model has difficulties explaining the steep spectrum of the SED before the low-frequency turn-over. It also does not allow for the contribution of a jet to the total radio emission as is common for compact sources. The young-jet plus aged-lobe model gives a lower limit on the redshift of GLEAM J0917–0012 as
$z > 6.5$
under the assumption of a two-component spectrum. In this case, we would get
$\sim 50\%$
jet prominence at
${\sim}150\,\rm MHz$
, consistent with the high values expected from compact sources (Hardcastle, Lawrence, & Worrall Reference Hardcastle, Lawrence and Worrall1998). Figure 7 presents both the best-fit SEDs for the young-jet plus aged-lobe and lobe-dominated RAiSERed models, as well as the PDF of redshift from the young jet plus aged-lobe model.
Results of the radio SED fitting with RAiSERed and a remnant model. For each model we present the best-fit values for the low-frequency spectral index, (
$\alpha_{\rm low}$
from the electron energy distribution, i.e.
$N(E)\propto E^p$
and
$p=2\times\alpha_{\rm low}-1$
), the break-frequency (
$\nu_{\rm break}$
), and, in the case of the remnant model, the fraction (T) the age where the jet has been turned off. We also give the resultant value of the Akaike information criteria (AIC) from the fit and the constraints on the redshift for the two RAiSERed models.
a RAiSERed is only applicable for active sources, so no redshift constraint is available.
While the Akaike information criterion (AIC, Akaike Reference Akaike1974, where lower values imply a preferred fit) significantly prefers the young-jet plus aged-lobe over the lobe-dominated model (
$\Delta$
(AIC)
$>20$
), the remnant spectrum has the lowest AIC of the three models. However, following Turner & Shabala (Reference Turner and Shabala2019, their Equation (19)), we estimate that the magnetic field strength of GLEAM J0917–0012 is
$>$
10.6 nT assuming a spherical source of angular size 0.7 arcsec at
$z = 2.01$
(we assume the equipartition factor derived by Turner et al. Reference Turner, Rogers, Shabala and Krause2018a, for 3C sources); this increases to
$>$
18.5 nT assuming the axis ratio of Cygnus A. For these magnetic field strengths, the optically thin spectral break frequency of the remnant spectrum implies an active age of
$<$
0.51 and
$<$
0.23 Myr respectively (e.g. Equation (4) of Turner, Shabala, & Krause Reference Turner, Shabala and Krause2018b). The jet is unlikely to have fully escaped its host galaxy on such timescales and therefore we would expect strong FFA in our
$<$
1 GHz observations (Bicknell et al. Reference Bicknell, Mukherjee, Wagner, Sutherland and Nesvadba2018).
Another constraint from the radio emission comes from assuming that this radio galaxy is a peaked source lying on the peak rest-frame frequency versus physical size relation (O’Dea & Baum Reference O’Dea and Baum1997). This relation states that as compact radio sources expand and grow, their peak frequency shifts to lower values. The fitting of the radio SED results in an observed-frame turn-over frequency of
$<50\,$
MHz at a 90% confidence (consistent with D20 and D21). In Table 6 we report the
$1.4$
-GHz radio luminosity, the rest-frame peak frequency upper limit and hence the inferred physical size lower limit (assuming the source lies on the peak/size relation) for both the potential redshifts. Note, the uncertainty on the peak frequency dominates the uncertainty in the peak/size relation. We know that at any redshift
$>1.6$
our
$0.7\,$
arcsec radio size limit makes the projected size of this source
$<6\,$
kpc. The higher (
$z\sim 8$
) redshift solution is consistent with this limit, while the
$z\sim 2$
solution is in some tension, but potentially consistent within the range of linear sizes observed for radio galaxies with lower rest-frame peak frequencies.
Derived line and galaxy properties for two different redshift solutions for the
$\lambda=1.12\,\mu$
m line in the GLEAM J0917–0012 spectra. These redshifts correspond to known strong emission lines (unlike some solutions in the grizli PDF). For both line identifications we present the corresponding redshift, observed-frame wavelength FWHM (from the template), the line flux (from the template fitting), the rest-frame velocity FWHM, the equivalent width (EW), and luminosity of the line for each identification. We also add the stellar mass (
$M_{\rm stel}$
) from the best-fit PÉGASE templates at that redshift, the
$1.4$
-GHz radio luminosity from the synchrotron model, the (limit of the) low-frequency rest-frame turn-over (derived from the TPL fitting) and the corresponding constraint on the largest extent in the radio,
$\theta_{\rm LAS}$
, (from the 0.7 arcsec angular size) assuming this source follows the size/turn-over relation from O’Dea & Baum (Reference O’Dea and Baum1997).
5.4. Combined redshift constraints
5.4.1. GLEAM J0917–0012
Taking the observed line at
$1.12\,\mu$
m into account and assuming it is one of the stronger lines seen in galaxy spectra, then the two lines most consistent with the bespoke UV-to-radio SED fitting are the [OII]
$\lambda 3727$
doublet at
$z=2.01$
or Lyman-
$\alpha$
at
$z=8.21$
. These solutions both have high peaks in the grizli redshift PDF (Figure 6). We examine the properties of the line and the galaxy in Table 6 for these two potential redshifts. The full-width half-maximum (FWHM) of the lines comes from the best-fit grizli template so should be taken with some caution as should the inferred rest-frame velocity widths. The line fluxes are derived from the normalisation of the line component of the best-fit template. Their differing values can be explained by the larger equivalent width (EW) found for the [OII] solution which is due to the relatively higher continuum model used in that fit.
We can investigate how consistent the weak features in the grism spectra are with the best-fit templates. The Lyman-break flux for each of the PEGASÉ templates for the
$z\sim 8$
solution is just under
$10^{-14}\,$
erg s-1 cm-2. As a flux density this is equivalent to just under
$10^{-19}\,$
erg s-1 cm-2 Å-1 at
$1.12\,\mu$
m. This value is very close to the observed strength of the feature in Figure 2. The flux of the
$4\,000\,$
Å breaks seen in the templates for the
$z\sim 2$
solution is around half to one third less, hence is less consistent with the flux observed in the spectrum.
We also present in Table 6 the derived properties of this radio galaxy for either redshift option including the line luminosity, stellar mass, 1.4-GHz radio luminosity as well as constraints on the angular size derived from the rest-frame peak frequency limit (assuming that the source lies on the relation between peak frequency and physical size presented in O’Dea & Baum Reference O’Dea and Baum1997). The limit on the observed peak frequency,
$< 50\,$
MHz at a 95% confidence, comes from the radio SED fitting in D21 and is held constant in the UV-to-radio fitting (Section 4.1).
For the
$z=8.21$
solution, the inferred Lyman-
$\alpha$
luminosity is well below the knee of the
$z\sim 6.6$
luminosity function (Santos, Sobral, & Matthee Reference Santos, Sobral and Matthee2016) suggesting that the host galaxy does not have an extreme star formation rate or perhaps has moderate extinction in the UV. The stellar mass, from the PEGASÉ templates, of
$10^{10-10.7}\,{\rm M}_\odot$
is slightly higher than few known galaxies at
$8\le z\le 11$
which typically have masses of a few
$10^9\,{\rm M}_\odot$
e.g. A2744_YD4 at
$z=8.38$
and gnz11 at
$z=10.96$
Laporte et al. Reference Laporte2017; Oesch et al. Reference Oesch2016 respectively. However, at this redshift there could be a non-negligible contribution to the rest-frame UV/optical from inverse-Compton (iC) emission either from the radio jet (e.g. Ghisellini et al. Reference Ghisellini, Celotti, Tavecchio, Haardt and Sbarrato2014), or intense compact star-formation, which may lead to an over-estimation of the stellar mass. The inferred limit on the physical size of the radio emission, from the rest-frame turn-over frequency, is consistent with that from the ALMA imaging and the observed interplanetary scintillation at low frequencies radio (as discussed in D21).
For the
$z=2.01$
[OII] solution none of the rest-frame parameters are inconsistent with being at that redshift apart from the size implied by the peak frequency which is mildly inconsistent with the relation from O’Dea & Baum (Reference O’Dea and Baum1997). However, the radio luminosity,
$L_{\rm 1.4 GHz}=2.2\times 10^{27}\,$
W Hz-1, is quite high for a compact, but non-beamed, radio galaxy. Also for this solution one might expect to see [NeV]
$\lambda 3426$
, [NeVI]
$\lambda 3347$
or [MgII]
$\lambda 2799$
in the grism spectrum which we do not.
5.4.2. Companion galaxy
While the bespoke UV-to-radio SED fitting suggests only a broad redshift solution across
$0.1<z<3$
, when taken in combination with the EAZY photometric redshift presented here and the one from D21, then the companion is more likely to be
$z\sim 2$
. Hence, a probable identification of the 0.98-
$\mu$
m line is [OII]
$\lambda 3727$
at
$z=1.63$
. The different wavelengths of the weak features seen in each spectrum suggest, despite their close proximity on the sky, that even if the lower redshift solution for GLEAM J0917–0012 is correct these two sources likely lie at different redshifts.
5.5. Implications of a radio galaxy at
$z\sim 2$
Radio galaxies with a luminosity of
$L_{\rm 1.4GHz}\sim 4\times 10^{28}\,$
W Hz-1 are not unheard of at
$z\sim 2$
(e.g. De Breuck et al. Reference De Breuck2010). However, this source is remarkably compact, smaller than one arcsec, for a source of that luminosity. Most
$z\sim 2$
radio galaxies are many tens of kpc in projected extent (
$\gg 5\,$
arcsec in the observed frame). A few compact sources are comparable in luminosity, e.g. PKS B0008-421 at
$z=1.12$
with
$L_{\rm 1.4GHz}\sim 4\times 10^{28}\,$
W Hz-1 (Callingham et al. Reference Callingham2015, Reference Callingham2017). PKS B0008-421 turns over at
$\sim 1.3\,$
GHz in the rest frame and therefore is a gigahertz peaked-spectrum source (it is also observed to be
$\sim 0.2\,$
arcsec in size King Reference King1994). GLEAM J0917–0012 would turn over at the much lower rest-frame frequency of
$\sim 0.1\,$
GHz at
$z\sim 2$
; hence, from the (O’Dea & Baum Reference O’Dea and Baum1997) size-luminosity relation, we may expect it to be larger than the 0.7 arcsec we measure (see arguments in Section 5.3.3). The peak of the radio emission could be better constrained by lower-frequency radio observations at
$<70\,$
MHz, e.g. with LOFAR.
From D21, a source at this redshift would have just one CO line (the 3-2 transition at
$\sim 115\,$
GHz) observable from the spectral observations presented in D21. There is no evidence of a robust line near this frequency although it is at the upper end of the spectral band. Therefore this redshift is adjacent to a gap in our coverage at
$1.8\le z\le 2.0$
where we have no constraints on the CO line (see Figure 6d in D21). If GLEAM J0917–0012 was at a redshift of
$z=2$
, or just above, the CO luminosity of this source would be limited to
$L^\prime_{\rm CO}<$
a few
$10^9\,$
K km s-1 pc
$^2$
, otherwise it is unconstrained.
The RAiSERed modelling finds that this source would be dominated by an aged lobe if at
$z\sim 2$
rather than a young-jet plus aged-lobe model like most compact sources. This hypothesis could easily be tested by future high-resolution Very long baseline interferometric (VLBI) observations. If the radio emission is dominated by more extended synchrotron radiation then the low-frequency turn-over would be more likely due to SSA which is also testable by observations below the MWA frequency band.
5.6. Implications of a radio galaxy at
$z\sim 8$
The fitting of the 100-GHz datum suggests a likely contribution from both synchrotron emission and cold dust. If this cold dust were associated with star formation then there could be further iC emission from the compact star formation. The SED fitting in the far-IR suggests observations at
$150{-}200\,$
GHz would be very useful in confirming the presence of cold dust and star formation. From D21 a source at this redshift may have three CO lines (the 7-6, 8-7 and 9-8 transitions) that fall into the ALMA coverage, but none were detected. This would limit the CO luminosity to
$L^\prime_{\rm CO}<$
a few
$10^9\,$
K km s-1 pc
$^2$
.
This result would also confirm our original hypothesis (D20) that such a HzRG would be compact (
$<0.7\,$
arcsec) and turn over at low observed-frame frequencies. VLBI radio observations at tens of micro-arcsec would determine the true size of the radio emission. Such observations would help us determine whether the turn-over is due to SSA or FFA. Such observations would also reveal the morphology of this compact source and help test the jet verses lobe hypothesis in Section 5.3.3.
If such a luminous radio source existed in the early Universe it would permit studies of the intervening neutral hydrogen through studies of the absorbed 21-cm line shifted to low frequencies (
$154.25\,$
MHz). Therefore a sensitive spectral scan from
$154\,$
MHz to
$\sim 189\,$
MHz (
$\equiv z\sim 6.5$
, the approximate end of reionisation) would reveal the distribution of neutral hydrogen along one particular line of sight in the early Universe. Such an observation would likely require the Square Kilometre Array.
5.7. Uniqueness of radio SED
The curvature and steepness in the MWA SED used to select these sources in D20 favours sources which appear to peak below the lowest MWA frequency of 70 MHz. This curvature is a common characteristic of HzRGs. For example, taking the sample of the 158 radio galaxies known at the time to be at
$z>2$
from Miley & De Breuck (Reference Miley and De Breuck2008), we find 76 with GLEAM fluxes of
$S_{\rm 150MHz}\ge 0.5\,$
Jy. Of these, just under half (34/76) match the mildly updated curvature and steepness criteria for our full sample (Broderick et al., in preparation)—note we do not apply the compactness criteria on the Miley & De Breuck (Reference Miley and De Breuck2008) sample. Above
$z=3.5$
, most bright radio galaxies (9/11) meet our curvature and steepness criteria.
The highest redshift of these, TN J0924–2201 at
$z=5.19$
, has a curved MWA spectrum, but peaks at a higher frequency than the sources meeting our selection criteria. The only other sources similar in luminosity, i.e.
$L_{\rm 1.4GHz}>10^{26.5}\,$
W Hz-1, discovered since are GLEAM J0856
$+$
0223 at
$z=5.55$
(D21) and TGSS J1530
$\,{+}\,$
1049 at
$z=5.72$
(Saxena et al. Reference Saxena2018b), both mentioned in the introduction (also see Figure 1 from Broderick et al., in preparation).
6. Conclusions
We have presented new HST/WFC3 grism observations of the powerful radio galaxy GLEAM J0917–0012 and a companion galaxy
$\sim$
one arcsec away. The spectra of these two sources each show negligible continuum emission and a faint emission line. These single lines and the overall faintness of these two galaxies makes obtaining their redshift very difficult. Confirmation of faint, potentially very high-redshift sources requires combining numerous techniques. Hence, we supplement the spectra and the redshift estimates by grizli with three SED template fitting methods: a bespoke template-plus-synchrotron-function model to the UV-to-radio data, the UV/optical/IR photometric redshift code EAZY and the radio photometric redshift code RAiSERed. Our UV-to-radio SED fitting uniquely combines several key features including the combination of template plus analytical model fitting and proper consideration of upper limits from non-detections. While our results cannot be definitive, we present our best estimate of this source’s redshift. From combining the SED fitting as a function of redshift with the faint features in the spectra, our main conclusions are as follows:
-
For the host galaxy of GLEAM J0917–0012 the UV-to-radio SED fitting allows for two broad solutions at
$z\sim 2$
and
$z\sim 8$
. Including the
$1.12\,\mu$
m line seen in the WFC3 spectrum, the host galaxy of GLEAM 0917–0012 most likely lies at either
$z=2.01$
or
$8.21$
. The EAZY fitting supports the lower of these two solutions and the RAiSERed fitting supports the higher. The best-fitting UV-to-radio templates at
$z\sim 8$
, which have a strong Lyman break around
$\sim1.1\,\mu$
m, are more consistent with the strength of the faint feature seen in the grism spectra than the best-fitting templates at
$z\sim2$
. -
For the companion the UV to far-IR SED fitting suggests this galaxy is likely at
$z<3$
which is consistent with the results from EAZY (
$z\sim2.6\pm0.5$
). Hence, the
$0.98\,\mu$
m line in the spectrum could well be the O[II]
$\lambda 3727$
doublet at
$z=1.63$
.
However, a categorical determination of the redshift of GLEAM J0917–0012 is not yet possible. The best prospects for confirmation lie with either further deep near-IR or millimetre spectroscopy.
If at
$z=8.21$
this galaxy would be the most distant radio-loud AGN known, and well within the EoR. Therefore this source would become a high-priority target for numerous existing and future instruments. Observations with the SKA would allow spectral scans searching for direct and statistical measures of the distribution of neutral hydrogen along the line of sight. Facilities such as JWST and ground based
$30{-}40\,$
m class telescopes could be used to study the host galaxy in the near-IR as well as AGN contributions in the mid-IR.
If GLEAM J0917–0012 were to lie at
$z\sim 2$
then it would still be one of the most luminous radio-loud AGN known and be extremely compact for it rest-frame low-frequency turn-over. VLBI observations would be key to determine its size and morphology; its radio power is best explained by a lobe-dominated model, but it would be unusually compact for such a source.
Acknowledgements
We acknowledge the Noongar people as the traditional owners and custodians of Wadjak boodjar, the land on which the majority of this work was completed. We thank the referee, Gabriel Brammer, for his comments which enhanced the quality of this paper.
GN acknowledges funding support from the Natural Sciences and Engineering Research Council (NSERC) of Canada through a Discovery Grant and Discovery Accelerator Supplement, and from the Canadian Space Agency through grant 18JWST-GTO1.
This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program 16184.
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
This publication has made use of data from the VIKING survey from VISTA at the ESO Paranal Observatory, programme ID 179.A-2004. Data processing has been contributed by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh.
The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.
This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org.
This paper is based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center (ADC) at National Astronomical Observatory of Japan. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. The Subaru Telescope is honoured and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii.
This research made use of ds9, a tool for data visualisation supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEASARC) with support from the JWST Mission office at the Space Telescope Science Institute for 3D visualisation. This paper also made use of the TOPCAT software (Taylor Reference Taylor, Shopbell, Britton and Ebert2005).
JA acknowledges financial support from the Science and Technology Foundation (FCT, Portugal) through research grants PTDC/FIS-AST/29245/2017, UIDB/04434/2020 and UIDP/04434/2020.
