The radio-loud fraction of quasars at z > 6 Open Access

ABSTRACT

Quasars at redshifts z > 6 are an excellent probe of the formation and evolution of supermassive black holes in the early Universe. The population of radio-luminous quasars is of particular interest, as such quasars could potentially be used to study the neutral intergalactic medium during cosmic reionization via H i 21 cm absorption studies. However, the lack of deep radio observations of z > 6 quasars leaves the population poorly constrained, and suitable candidates for an H i 21 cm absorption study have yet to be found. In this work, we present Jansky Very Large Array (VLA) 1–2 GHz radio continuum observations of 138 quasars at redshifts 6.0 ≤ z < 7.6. We detect the radio continuum emission of the z = 6.1 quasar J1034−1425, with a 1.6 GHz flux density of |$170\pm 36\, \mu$|Jy. This quasar is radio-quiet with radio-loudness, R ≡ f_5 GHz/f_{ν, 4400 Å} = 2.4 ± 0.5. In addition, we detect seven other quasars at z > 6, which have previously been characterized in the literature at these frequencies. Using the full sample, we estimate the radio-loud fraction to be |$3.8^{+6.2}_{-2.4}\ \hbox{per cent}$|⁠, where the uncertainties are 95 per cent confidence intervals. This is lower than recent estimates of the radio-loud fraction in the literature, but is still marginally consistent with no redshift evolution of the radio-loud fraction. We explore the undetected quasar population by stacking their continuum images at their optical positions and obtain a median stacked flux density of 13.8 ± 3.9 µJy and luminosity of log L_5 GHz/(W Hz⁻¹) = 24.2 ± 0.1.

1 INTRODUCTION

Quasars are among the most luminous objects in the Universe, powered by the accretion of matter onto supermassive black holes (SMBHs). Studying quasars at high redshifts (z > 6) can provide crucial information on the evolution of SMBHs, their host galaxies, and the surrounding large-scale structures in the early Universe. Approximately 10 per cent of all quasars up to z ∼ 5 are ‘radio-loud’, with strong radio synchrotron emission associated with relativistic jets launched from their active galactic nuclei (AGNs; e.g. Kellermann et al. 1989; Ivezić et al. 2002; Kellermann et al. 2016). The energy required to launch the powerful jets is thought to originate from the rotational energy of the SMBH (Blandford & Znajek 1977), the accretion disk (Blandford & Payne 1982), or a combination of both.

A useful measure of the strength of a quasar’s radio emission is given by the ‘radio loudness’ R, which is usually defined as the ratio of the rest-frame luminosities at 5 GHz and 4400 Å (Kellermann et al. 1989). A quasar is typically characterized as radio-loud if it has R > 10, and as radio-quiet otherwise. Traditionally, it has been thought that AGN activity dominates the continuum emission of radio-loud quasars, while radio-quiet emission may be due to a combination of AGN activity and host-galaxy star formation. However, using deep Low Frequency ARray observations, Gürkan et al. (2019) find that there is no clear separation between radio-loud and radio-quiet quasars i.e. no bimodality in the distribution of R, hence raising doubts about whether R carries any physical meaning.

While the definition of R is somewhat arbitrary, the fraction of radio-loud quasars at any given redshift may depend on several factors, including the SMBH accretion modes (Best & Heckman 2012), the black hole masses, and spins, etc (e.g. Rees et al. 1982; Laor 2000). Any evolution in the radio-loud fraction (RLF) would thus indicate evolution in one or more of these properties. Indeed, there has been much debate in the literature as to whether the radio properties of quasars evolve with redshift, with some studies finding that the RLF is a strong function of redshift and optical luminosity and others finding evidence for an unevolving RLF. For example, Jiang et al. (2007) find that the RLF increases with optical luminosity and decreasing redshift, predicting an RLF of about 2 per cent at z ≈ 6. This is marginally consistent with the findings of Wang et al. (2007), which, however, were obtained with a small sample size of only 13 quasars. Kratzer & Richards (2015) confirm the results of Jiang et al. (2007), but point out that the observed differences between the behaviour of the RLF and the mean radio-loudness suggest the presence of selection biases. Other studies do not find such correlations, reporting a constant RLF consistent with ≈10 per cent up to redshifts of z ∼ 6 (Stern et al. 2000; Ivezić et al. 2002; Bañados et al. 2015; Gloudemans et al. 2021; Liu et al. 2021). Studies at the highest redshifts have hitherto been limited by small sample sizes, and have hence been unable to rule out redshift evolution. Bañados et al. (2015) stacked the images from the 1.4 GHz Very Large Array FIRST survey (Becker, White & Helfand 1995a) at the positions of 41 quasars at 5.5 < z < 7.1, finding no significant detection of the stacked continuum emission. Instead, they set an upper limit of 84 µJy to the mean flux density of the undetected quasars at 1.4 GHz. Clearly, a deeper radio survey of a large quasar sample is needed to detect radio continuum emission from quasars at z > 6, and to examine the RLF at these redshifts.

In addition to investigating the SMBH and galaxy evolution, quasars at z > 6 are excellent tools for studying the Epoch of Reionization (EoR), during which the neutral hydrogen (H i) in the intergalactic medium (IGM) was ionized by the UV radiation of the first luminous sources. In particular, quasars could potentially be used to probe the neutral IGM via absorption by the forbidden H i 21 cm hyperfine transition (Carilli, Gnedin & Owen 2002; Carilli et al. 2004). Such a 21 cm absorption study would be able to probe small-scale structures in the high-z IGM without suffering from many of the systematics and wide-field effects that affect most all-sky 21 cm emission experiments (e.g. Paciga et al. 2013; Kolopanis et al. 2019; Mertens et al. 2020; Trott et al. 2020; HERA Collaboration 2022; Singh et al. 2022; Keller et al. 2023).

The detection of the cosmological H i 21 cm line in absorption requires a strong compact source (150 MHz flux density ≈10–100 mJy) at z ≳ 6.5 to achieve the necessary sensitivity with modern radio telescopes (e.g. Ciardi et al. 2013). To date, there are only four known sources at z > 6 with flux densities exceeding 10 mJy at ∼150 MHz: J1427+3312 (McGreer et al. 2006), J1429+5447 (Willott et al. 2010a; Frey et al. 2011), J2318−3113 (Decarli et al. 2018; Ighina et al. 2021), and J0309+2717 (Belladitta et al. 2020). Of these quasars, J1429+5447 has the highest redshift at z ≈ 6.2, which is marginal to probe the IGM during the EoR. Alternatively, Thyagarajan (2020) proposed to measure the H i 21 cm absorption statistically by stacking the 1-D power spectra of a large number of fainter sources (flux densities ≈1–10 mJy). The selection of a suitable observation strategy requires knowledge about the flux density distribution of the quasar population at these redshifts.

In this study, we used the L-band receivers of the Karl G. Jansky Very Large Array (VLA) to obtain deep 1–2 GHz continuum images of 138 optically confirmed quasars at z > 6, to measure their radio properties. Section 2 describes our quasar sample and the VLA observations. Section 3 describes the imaging and post-imaging data processing. Our results are presented in Section 4, and a summary of this work is provided in Section 5.

Throughout this paper, we assume a standard ΛCDM cosmology with H₀ = 67.66 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and |$\Omega _\Lambda =0.7$| as measured by Planck Collaboration et al. (2020). Magnitudes are given in the AB system.

2 DATA AND OBSERVATIONS

2.1 The quasar sample

Our sample of 138 quasars at z ≥ 6 is based on the sample of the Pan-STARRS1 survey (PS1; Bañados et al. 2014, 2016; Tang et al. 2017), but also includes quasars from a number of other recent optical and near-infrared high-redshift surveys: SDSS (Jiang et al. 2009), CFHQS (Willott et al. 2007, 2009, 2010a, b), UKIDSS (Mortlock et al. 2011), VIKING (Venemans et al. 2013), VST-ATLAS (Carnall et al. 2015), DES (Reed et al. 2015, 2017), HSC/SHELLQs (Kashikawa et al. 2015; Matsuoka et al. 2016, 2018, 2019), and DECaLS (Bañados et al. 2018). In addition, our sample contains one radio-selected quasar, J1427+3312, which was discovered by McGreer et al. (2006). Eight quasars from the aforementioned surveys were not included because we restricted the sample to declinations ≥−44°, for better observing conditions with the VLA (see Fig. 1). Except for the declination constraint, the sample included all quasars at z ≥ 6 that were known at the time of the VLA observations.

The redshift and rest-frame UV magnitude (M₁₄₅₀) distributions of the quasar sample are shown in Fig. 2. The quasar redshifts are within the range 6.0 ≤ z < 7.6, but have a relatively low median of z = 6.15. The absolute UV magnitudes range from −29.1 to −20.9, with a median of −25.6 and a mode of ≈−27.

2.2 VLA L-band observations

The 1–2 GHz continuum observations were carried out with the VLA between 2019-09-01 and 2019-09-18 (proposal 19A-056, PI: N. Thyagarajan), using the most extended (A-array) configuration, which has a maximum baseline of 35 km. The A-array was preferred due to its high angular resolution, typically ≲2 arcsec (depending on source declination and the length of the observation). At z ∼ 6, this corresponds to a projected spatial resolution of about 11 kpc.

The WIDAR correlator was used as the backend for the observations, with two 512 MHz IF bands covering the frequency range 1–2 GHz. Each IF band was divided into eight 64-MHz sub-bands, each with 64 channels. The data were recorded in full-polar mode, with a time resolution of 2 s.

The observing programme was designed to achieve a 5σ detection threshold of 100 |$\mu \mathrm{Jy\,beam}^{-1}$| at 1.4 GHz, the approximate depth required to detect prospective candidates for a statistical 21 cm absorption study (see Section 1). Each quasar target was observed for around 12 min, resulting in a total on-source time of approximately 27.6 h. The observations were grouped in 12 blocks, each starting with 5 min on a primary (flux density) calibrator and interspersed with regular 2–3 min observations of secondary (phase) calibrators. In total, the programme took up 33 h of observing time at the VLA. In each observing block, flux calibration was performed on 3C48, 3C147, or 3C286, using the flux scale of Perley & Butler (2017).

The data editing and calibration were done with standard tools available from the Common Astronomy Software Applications package (CASA version 6.5.3, McMullin et al. 2007). In addition to running the automated CASA tasks, we manually removed time and frequency ranges that are affected by strong radio frequency interference (RFI), resulting in a usable bandwidth of 500–600 MHz at a mean frequency of 1.68 GHz. Lastly, we carefully inspected the various calibration solutions to identify and remove antennas that exhibit noisy gains or other kinds of abnormality. On average, this resulted in about 25 uncorrupted antennas at any given time.

3 IMAGING AND ANALYSIS

The imaging of the VLA L-band data was performed with the w-stacking algorithm provided by WSClean (version 3.3, Offringa et al. 2014; Offringa & Smirnov 2017). A subset (about 50 per cent) of the quasars was imaged independently using the w-projection algorithm (Cornwell, Golap & Bhatnagar 2008) in CASA. The results from the two approaches were found to be in good agreement, thus validating our WSClean approach.

For each quasar, we imaged a ∼1° field of view (FoV) at a pixel resolution of 0.3 arcsec, resulting in an image size of 12244 × 12244 pixels. We set the WSClean Briggs robust parameter¹ to −0.5, resulting in an average synthesized beam FWHM of ≈1.8 arcsec. Parallelizing over 32 CPUs, the computing time is around 15 min per image. We performed self-calibration on all target fields, where an improvement of the signal-to-noise ratio was achievable. Before making the final images, we deployed a set of standard flagging routines on the residual visibilities, but at lower rejection thresholds than in the initial pre-imaging RFI editing.

Following the imaging, we characterized the statistical properties of the images, which ultimately allowed us to determine which images exhibit reliable detections of the quasar radio continuum. The root-mean-square (RMS) of the image noise was robustly estimated as the median absolute deviation (MAD) of the central 41 × 41 pixel values, multiplied by a scale factor of 1.4826 appropriate for an assumed Gaussian noise distribution (Rousseeuw & Croux 1993). Using this estimate, we find a median noise RMS of |$36\pm 5\, \mu \mathrm{Jy\,beam}^{-1}$| across the cut-out images.

We fitted a 2-D Gaussian to the central 21 × 21 pixels of each image, using a least-squares method. We fixed the widths and position angle of the Gaussian to those of the synthesized beam, but let the position vary freely within a range of ±3 pixels. The errors of the fitted amplitudes and positions were determined as described in Condon (1997). Using the fitted parameters and their associated errors, we used the following two criteria to classify a source as detected: (1) the fitted amplitude is at least three times the noise RMS and (2) the optical position of the target lies within the 2σ uncertainty contour of the fitted position. In total, there are nine sources that meet these criteria, two of which have not previously been detected at radio frequencies. However, one of these two sources is likely to be a false detection, arising from systematics in the image. This will be discussed in more detail in Section 4. Note that there is also a small probability that another source falls within the cross-matching area around the optical position of the target (source confusion). Using the source counts in Hale et al. (2023) as a guide, we find a cumulative source density at |$S_\nu > 100\,\mu$|Jy of approximately 10⁻⁴ sources per square arcsec. Hence, our search area of 6 × 6 pixels (3.28 square arcsecs) contains on average 3 × 10⁻⁴ confused sources; source confusion is thus unlikely to significantly affect our results.

Using the assumption that the sources are unresolved, we take the fitted Gaussian amplitudes as a proxy of the total source flux density. Note that the central frequency of the continuum observations is 1.68 GHz. However, without knowledge of the underlying radio spectra, it is not possible to provide an exact frequency corresponding to the measured flux density. Instead, we assume a power-law spectrum (⁠|$S_\nu \propto \nu ^{\alpha _R}$|⁠) for the target quasars, where we adopt a typical spectral index of α_R = −0.75. This value is in keeping with previous studies (Wang et al. 2007; Momjian et al. 2014; Bañados et al. 2015; Liu et al. 2021) and is in agreement with very long baseline interferometric measurements (Frey et al. 2005; Momjian, Carilli & McGreer 2008; Frey et al. 2011; Momjian et al. 2018), which find steep radio spectra for quasars at z ∼ 6. Using this assumed spectrum, we find that our measured continuum flux densities would typically correspond to spectral flux densities at an observed frequency of 1.6 GHz.

We further transform these flux densities to radio luminosity densities |$L_{\nu ,\mathrm{5\, GHz}}$| at a rest-frame frequency of 5 GHz by applying a K-correction (cf. Hogg 1999; Hogg et al. 2002; Condon & Matthews 2018)

where D_L is the luminosity distance. If a source is detected with S/N ⪆ 10, it can be imaged in adjacent 64 MHz VLA sub-bands while retaining sufficient S/N to compute an in-band spectral index (see Section  A in the Appendix for a discussion on calculated in-band spectral indices). For such cases, we perform the extrapolation in equation (1) using the fitted power law rather than the assumed one.

4 RESULTS

The flux densities, radio luminosities, and radio-loudness parameters derived from our VLA L-band detections are listed in Table 1, along with the redshifts and optical data from the literature. For quasars that are not detected in this study, we provide a table as supplementary material to this paper with 3σ upper limits on the L-band flux densities and on any other quantities derived therefrom. In total, there are nine quasars detected with ≥3σ significance, of which three are radio-loud and two do not have any previous radio detections. The two new detections, J0231−2850 (z = 6.0) and J1034−1425 (z = 6.1), are found at 3.4σ and 4.7σ significance, respectively. However, the radio image of J0231−2850 exhibits clear non-Gaussian artefacts (stripes), and the detection should hence be deemed marginal. To be conservative, we hereafter treat J0231−2850 as a non-detection. The postage stamp cut-outs of the L-band continuum images of the eight likely detections and J0231−2850 are shown in Fig. B1 of the Appendix. There is no evidence for extended emission in the cut-out images. We checked this more carefully by subtracting the scaled restoring beam from the cut-outs and inspecting the residuals around the quasar position. In these residuals, we found no sign of extended emission, thus validating our assumption that the sources are unresolved.

4.1 Radio loudness

The strength of the radio emission of a quasar relative to its optical emission is characterized by the radio-loudness parameter, R. There are many definitions of R in the literature (see, e.g. Hao et al. 2014, for a review). Here, we adopt the most widely used definition |$R=f_{5 \, \rm{GHz}}/f_{\nu ,\rm{4400 \mathring{\rm A}}}$| (Kellermann et al. 1989), where f_5 GHz and f_{ν, 4400 Å} are the flux densities at rest-frame 5 GHz and 4400 Å, respectively. This definition is also in keeping with other high-redshift studies (e.g. Bañados et al. 2015; Gloudemans et al. 2021; Liu et al. 2021) therefore allowing a fair comparison with the results therein. A quasar is generally classified as radio-loud if R > 10, and as radio-quiet otherwise.

At z = 6, the rest-frame wavelength of 4400 Å corresponds to an observed wavelength of approximately |$3\, \mu$|m. To mitigate potential large extrapolation errors, f_{ν, 4400 Å} should thus be estimated from mid-IR or near-IR observations. In our sample, there are 16 quasars with available photometric data from Spitzer’s Infrared Array Camera (IRAC; Fazio et al. 2004) at a wavelength of 3.6 µm. There are 49 additional quasars that have W1 photometric data in the ALLWISE catalogue of the Wide-Field Infrared Survey (WISE; Wright et al. 2010) at a wavelength of 3.4 µm. We use these together with a commonly assumed optical spectral index of α_ν = −0.5 (e.g. Wang et al. 2007; Bañados et al. 2015) to extrapolate the flux densities to a rest-frame wavelength of 4400 Å. For the remaining 73 sources, we extrapolate from the AB magnitudes at rest-frame 1450 Å, using the same spectral index. As pointed out by Bañados et al. (2015), some quasars might deviate considerably from the average and thus be subject to large extrapolation errors. By comparing the extrapolations from M₁₄₅₀ with those from existing IR-photometry, Bañados et al. (2015) estimate the fractional error to be 0.3 for quasars at redshifts 5.5 < z < 6.5. We adopted the same error for our flux densities that are extrapolated from M₁₄₅₀, with the caveat that there may be individual objects with considerably larger extrapolation errors.

4.1.1 Constraining the radio-loud fraction

In the following, we use the derived radio-loudness parameters to constrain the fraction of quasars that are radio-loud. In doing this, we replace our radio flux density upper limits with published results from earlier VLA observations, where available. The quasars for which such observations are available are J1148+5251 (Carilli et al. 2004), J0227−0605 (Liu et al. 2021) and J1319+0950 (Wang et al. 2011). Furthermore, we include the radio-loud quasar PSO J172.3556+18.7734 (Bañados et al. 2021), which was discovered after our VLA observations were carried out. We exclude the radio-selected source J1427+3312 from the sample to prevent a bias towards radio-loud quasars.

Fig. 3 shows the distribution of the derived optical and radio luminosities along with the regions that define radio-loud and radio-quiet quasars, respectively. Quasars with R < 10 are classified as radio-quiet, regardless of whether or not they are detected in our L-band observations. Quasars with R > 10, on the other hand, can only be conclusively classified as radio-loud if they are detected, since an upper limit in that region cannot rule out the possibility that the quasar is radio-quiet. In total, there are 65 sources that cannot be characterized as radio-loud or radio-quiet. For the remaining sources, three are radio-loud detections and 69 are robustly characterized as radio-quiet. Including literature values and excluding the radio-selected quasar J1427+3312, we have in total 65 ambiguous, four radio-loud, and 69 radio-quiet quasars.

Using only sources that can be robustly classified as either radio-loud or radio-quiet, we find a radio-loud fraction of |$4/(69+4)=5.5\pm 2.7\ \hbox{per cent}$|⁠, where the error is estimated from Poisson statistics. If we include the ambiguous quasars and assume that they are all radio-quiet, we obtain a lower estimate of the RLF, |$4/138=2.9\pm 1.4\ \hbox{per cent}$|⁠. This should be considered a lower bound to the value of RLF, as there may be some radio-loud quasars among the ambiguous cases.

A more rigorous approach to estimating the RLF makes use of survival analysis, which is normally used to analyse the expected time τ that it takes for a specific event to occur (e.g. lifespan analyses in biology). The probability of an event E not having occurred at time t is given by the Survival function S(t) = Pr(τ > t). For our purposes, we replace ‘time’ with ‘decreasing radio-loudness’, and ‘event’ with ‘detection’. In other words, we are interested in the cumulative density function F(r) = Pr(R < r) = 1 − S(r) of the radio-loudness parameter r. We estimate F(r) by employing the Kaplan–Meier estimator (KM, Kaplan & Meier 1958)

where R_i is the radio-loudness parameter of a detected source and n_i is the number of sources that have not been detected at R > R_i. This can be shown to be the maximum-likelihood estimator of F(r) (e.g. Andersen et al. 1993). The advantage of using the KM estimator is that it can deal with censored data, i.e. data where there is only partial information about the exact time (radio-loudness) at which an event (detection) occurs. In this study, we are dealing with left-censored data, or upper limits, which is why we estimate F(r) rather than S(r). The RLF is then simply Pr(R > 10) = 1 − F(r = 10).

Our estimate of F(r) is shown in Fig. 4 and was calculated using the Python package lifelines (Davidson-Pilon 2019). The RLF thus obtained is |$3.8^{+6.2}_{-2.4}\ \hbox{per cent}$|⁠, where the uncertainties are 95 per cent confidence intervals.

4.1.2 Does the radio-loud fraction evolve with redshift?

Jiang et al. (2007) used a sample of ≳30 000 optically-selected quasars at redsfhits 0 < z ≤ 5 to find that the RLF is a strong function of both luminosity and redshift: at a fixed 2500 Å absolute magnitude of −27.7, they found the RLF to decrease from ≈40 per cent at z = 0.5 to ≈8 per cent at z ≈ 3. Conversely, at any redshift, they found the RLF to be higher for more luminous quasars. We note that the median absolute magnitude of our sample at 2500 Å is −26.4, for an assumed optical-NUV spectral index of α = −0.5. However, this median magnitude is close to the lower magnitude limit at which our quasars are robustly classified as radio-loud or radio-quiet. For objects with robust classifications, the effective absolute magnitude of our sample is thus slightly brighter, ≈−27.7 at 2500 Å.

Fig. 5 shows a comparison between our RLF and various literature estimates of the RLF at z > 4 (Stern et al. 2000; Ivezić et al. 2002; Bañados et al. 2015; Yang et al. 2016; Gloudemans et al. 2021; Liu et al. 2021; Lah et al. 2023), along with the 1σ RLF bands predicted by the fit of Jiang et al. (2007) for M₂₅₀₀ = −27.7. It is clear that our estimate of the RLF is somewhat lower than, but formally consistent with, the earlier estimates of the RLF in the literature. Similar to the earlier studies, we find no significant evidence for redshift evolution in the RLF between z ≈ 4 and 6.2. However, our measurement of the RLF is also consistent with the prediction of Jiang et al. (2007) of a decline in the RLF with increasing redshift. Deeper observations of a larger sample of quasars at z ≳ 6 are critically needed to test the above predictions.

4.2 Constraining the undetected source population

In the context of radio surveys, ‘stacking’ is used to glean information on sources that lie below the flux density detection threshold of the survey, but which have been identified at other wavelengths (e.g. infrared, optical, X-ray, etc). This enables the investigation of the statistical properties of the source population below the survey detection threshold. Stacking sources in bins of different parameters (e.g. redshift, stellar mass, colour, etc) allows us to investigate how the average source properties depend on these parameters. Further, the average flux density of individually undetected sources serves as a useful guide for deeper future surveys. In this section, we investigate the evolution of the average quasar radio properties with redshift.

Several methods of stacking radio continuum emissions have been proposed in the literature. White et al. (2007) compared the stack mean with the stack median, concluding that the latter is preferable due to its lower sensitivity to the high flux density tail of the distribution, as well as higher robustness towards potential systematic effects (e.g. confusion). They further show that the median approaches the mean in the limit of a noise-dominated distribution. Other authors (e.g. Mahony et al. 2012; Condon et al. 2013; Zwart et al. 2015) point out several limitations of stacking such as the bias towards flux densities just below the detection threshold. As an alternative, they suggest fitting parametric models to the flux density distribution without the use of a single summary statistic. This ‘stack-fitting’ has been further developed and employed in Mitchell-Wynne et al. (2014), Roseboom & Best (2014), and Malefahlo et al. (2020).

As we are dealing with only 138 quasars, which is small relative to the studies mentioned above, we did not attempt to fit the flux density distribution and leave this as a possibility for future work. Instead, we stack the flux density at the positions of the optical quasars, using median-stacking. As a reference, we also stack the flux densities at ‘blank-sky’ positions located ∼1 arcmin from the optical positions. We stack in two separate redshift bins, 6.0 < z ≤ 6.15 and 6.15 < z < 7.6 (low-z and high-z, hereafter, respectively), which were chosen so that both bins contain a similar number of sources (N_low = 68 and N_high = 70). The median redshifts of the bins are 6.06 and 6.37, respectively, corresponding to a time difference of approximately 60 Myr.

For each quasar, we measured the flux density and the luminosity at both the optical quasar position and the ‘blank-sky’ position. The flux density and luminosity distributions of the quasars in the two redshift bins, and those of the complete sample (all-z, hereafter), are shown as the black histograms in the upper two panels of Fig. 6. The blank-sky distributions are shown as the grey-shaded regions; the vertical grey solid line in each panel indicates the mean of the corresponding blank-sky distribution, while the the blue solid curve shows the Gaussian fit to the blank-sky distribution. Both the low-z and the all-z bins show an apparent excess of flux density at the quasar positions, compared to the blank-sky distribution. This suggests that there may be several sources at z ≈ 6 just below the detection threshold of this survey. The apparent excess is not seen in the high-z sample. In passing, we note that a two-sided Kolmogorov–Smirnov test finds that the null hypothesis that the flux densities at both the quasar positions and the blank-sky positions are drawn from the same distribution cannot be rejected at high significance (p-values of 0.17, 0.5, and 0.15 for the low-z, high-z, and all-z samples, respectively).

The lowest panel of Fig. 6 show the median-stacked images of the low-z, high-z, and all-z samples. The stacked images for the low-z and all-z bins have RMS noise values of ≈5.2 and ≈3.9 µJy, respectively, and show detections of the stacked continuum signal at S/N > 3σ significance. For the high-z bin, the stacked image has an RMS noise of ≈5.4 µJy, and does not show any evidence for statistically-significant emission at the image centre.

Table 2 lists the flux densities, the rest-frame 5-GHz radio luminosities, and the radio-loudness parameters obtained by fitting a single Gaussian to each median-stacked image. We also list the corresponding mean values for the blank-sky stacked images. The median stacked L-band flux density is S_ν, peak = 13.8 ± 3.9 µJy, while the median rest-frame 5-GHz radio luminosity of the quasar sample is log L_5 GHz/(W Hz⁻¹) = 24.2 ± 0.1. We obtain a median radio-loudness parameter of R = 0.7 ± 0.2 for the entire sample, confirming that the undetected quasar population at z ≳ 6 is predominantly radio-quiet. As expected, the mean blank-sky flux densities are consistent with zero, underlining that the observed flux density excess at the quasar positions is likely to indeed be due to an undetected population of radio sources.

5 SUMMARY

We report VLA 1–2 GHz continuum imaging of a sample of 138 quasars at z ≈ 6.0–7.6, aiming to determine their radio properties, search for possible redshift evolution in the radio-loud fraction, and identify good targets for follow-up searches for redshifted H i 21 cm absorption. Our search revealed a new radio-quiet quasar, J1034−1425, with flux density 170 ± 36 µJy at redshift z = 6.1, and detected seven additional quasars that have previously been characterized at radio wavebands. For those sources, our results agree well with literature measurements and are consistent with no temporal variability. We find no new radio sources with 1.6 GHz flux densities ≳1 mJy amongst the sample of 138 quasars, i.e. no immediately suitable targets for follow-up H i 21 cm absorption spectroscopy.

Our survey robustly classifies three quasars as radio-loud (R >10) and 69 as radio-quiet, while the remaining 66 quasars are indeterminate as our upper limits on their L-band flux density do not rule out their being radio-loud systems. Using the Kaplan–Meier estimator, we find a radio-loud fraction of |$3.8^{+6.2}_{-2.4}\ \hbox{per cent}$|⁠, where the uncertainties represent 95 per cent confidence intervals. Due to our large sample size relative to that of Bañados et al. (2015) and Liu et al. (2021), this is the most reliable measurement of the RLF at z > 6 to date. However, while the central value of our RLF is lower than earlier estimates at z ≳ 4, our measurement is formally consistent with the earlier studies within ≈2σ significance. Our results are thus consistent with an unchanging RLF at z > 4.

By stacking our continuum images, we find that the undetected quasar population has a median flux density of 13.8 ± 3.9 µJy, a factor of ≈3 below our typical RMS noise values. Significantly deeper observations would be needed to detect the L-band continuum emission from the bulk of the quasar population at these redshifts. Studies at these redshifts are hence likely to continue to be limited to the brightest quasars in the foreseeable future. However, due to the ongoing optical search for high-redshift quasars, the number of quasar identifications at z > 6 is steadily increasing (e.g. see new discoveries by DESI; Yang et al. 2023). According to the Million Quasar Catalogue (MILLIQUAS, Flesch 2023), there are now 308 known quasars at z ≥ 6 as of August 2023, more than double the number observed in this work. This will allow future surveys to achieve much tighter constraints on the RLF, which should allow a definitive conclusion on its possible redshift evolution.

ACKNOWLEDGEMENTS

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

We thank the anonymous referee for their feedback and comments, which improved the quality of this paper. We also thank Bojan Nikolic and Chris Carilli for valuable inputs.

The software developed for this analysis uses python and the publicly-accessible and open-sourced python packages Numpy (van der Walt, Colbert & Varoquaux 2011), SciPy (Virtanen et al. 2020), Astropy (Astropy Collaboration et al. 2022), Matplotlib (Hunter 2007), and lifelines (Davidson-Pilon 2019). We furthermore acknowledge the use of CASA (McMullin et al. 2007) for data editing and calibration, WSClean (Offringa et al. 2014; Offringa & Smirnov 2017) for imaging, and AOFlagger (Offringa et al. 2010; Offringa, van de Gronde & Roerdink 2012) for RFI excision.

P.M. Keller is funded by the Institute of Astronomy and Physics Department of the University of Cambridge via the Isaac Newton Studentship. AK and NK acknowledge the Department of Atomic Energy for funding support, under project 12-R&D-TFR-5.02-0700.

DATA AVAILABILITY

The full data table of this work (Table 1) is available in machine-readable form and is provided as supplementary material to this article. The raw data from the VLA observing programme 19A-056 is publicly available from the NRAO data archive at https://data.nrao.edu/.

Footnotes

References

APPENDIX A: IN-BAND SPECTRAL INDEX FITTING

The radio-spectra of astronomical sources are typically parametrized by a power law |$S_\nu = S_{\nu _0}(\nu /\nu _0)^{\alpha _R(\nu)}$|⁠, where α_R(ν) is a frequency-dependent spectral index, and |$S_{\nu _0}$| is the flux density at frequency ν₀. When extrapolating or interpolating to another frequency over a small fractional bandwidth, one usually assumes a constant spectral index.

The problem of extrapolating from a continuum measurement is two-fold. First, the local spectral index is not necessarily known and must be either measured or given an assumed value that is representative of the source population under investigation. Secondly, the frequency corresponding to the source flux density in the synthesized continuum image may not be the same as the central frequency of the image if the spectrum is not flat (i.e. ν₀ is uncertain). For bright sources, these uncertainties can be mitigated by dividing the band into sub-bands, which are imaged separately. The extracted flux densities can then be used to fit for ν₀ and α_R.

In this work, we have three quasars that are bright enough to fit their in-band spectra: J1427+3312, J1429+5447, and J2318−3113. We image nine spectral windows, each of which has a width of approximately 64 MHz, and is relatively free of RFI contamination. The flux densities within the spectral windows are computed as described in Section 3 and are plotted in Fig. A1. For these sources, we use our fitted power-law values to obtain the results listed in Table 1.

Table A1 lists the fitted values for the flux density and the spectral index at the observed frequency of 1.4 GHz, together with values from the literature (references in the table notes). All three sources have fairly steep spectra. The literature values for J1427+3312 and J1429+5447 are consistent (within 2σ significance) with our measurements. There are, nevertheless, several effects that might cause the results to differ. First, the spectral indices from the literature were computed from two spectral measurements separated by several GHz in each case (see the table notes in Table A1). Therefore, the results could be biased at 1.4 GHz in cases where α_R(ν) varies across frequency. Secondly, the spectral index of J1429+5447 was determined from VLBI measurements (Frey et al. 2011), which might be resolving out some of the radio emission and would thus be more sensitive to the compact core, which typically has a flatter spectrum. However, the spectral index found in Frey et al. (2011) is in fact steeper than the one obtained in this work, thus ruling out this possibility. Lastly, temporal variability may also be at play with some of the sources. While the literature has not yet noted any variability for J1427+3312 and J1429+5447, there have been inconsistencies for J2318−3113 at lower frequencies (888 MHz) hinting at variability (Ighina et al. 2021, 2022). Our measurements are consistent with literature values reported at similar frequencies, so we do not find any strong evidence for variability of these sources in the L-band (see Fig. A1 and Table A1).

APPENDIX B: VLA L-BAND IMAGES

Fig. B1 shows postage stamp cut-outs of the L-band continuum images of the eight likely detections and J0231−2850, which is treated as a marginal case due to the presence of non-thermal-like artefacts in its image. These detections are discussed in more detail in the results section (Section 4).