Radio-loud fraction of z>6 quasars

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 reionisation 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 \leq 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 \equiv f_{5\text{~GHz}}/f_{\nu,\text{4400 A}} = 2.4\pm0.5$. In addition, we detect 7 other quasars at z>6, which have previously been characterised in the literature at these frequencies. Using the full sample, we estimate the radio-loud fraction to be $3.8^{+6.2}_{-2.4}\%$, where the uncertainties are 95% 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\pm 3.9~\mu$Jy and luminosity of $\log{L_{5~\mathrm{GHz}}/(\mathrm{W~Hz}^{-1})}=24.2\pm0.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 ( > 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% of all quasars up to  ∼ 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. 1989Kellermann et al. , 2016;;Ivezić et al. 2002).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" , 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 characterised as radio-loud if it has  > 10, and as radio-quiet otherwise.Traditionally, it has been thought that AGN activity dominates the continuum emission of radio-loud ★ Email: pmk46@cam.ac.uk 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 , hence raising doubts about whether  carries any physical meaning.
While the definition of  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% at  ≈ 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% up to redshifts of  ∼ 6 (Stern et al. 2000;Ivezić et al. 2002;Bañados et al. 2015;Liu et al. 2021;Gloudemans 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 et al. 1995a) at the positions of 41 quasars at 5.5 <  < 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  > 6, and to examine the RLF at these redshifts.
In addition to investigating the SMBH and galaxy evolution, quasars at  > 6 are excellent tools for studying the Epoch of Reionization (EoR), during which the neutral hydrogen (H i) in the Intergalactic Medium (IGM) was ionised 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 et al. 2002(Carilli et al. , 2004)).Such a 21 cm absorption study would be able to probe small-scale structures in the high- 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;Trott et al. 2020;Mertens et al. 2020;HERA Collaboration et al. 2023;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  ≳ 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  > 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  ≈ 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  > 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.

The Quasar Sample
Our sample of 138 quasars at  ≥ 6 is based on the sample of the Pan-STARRS1 survey (PS1; Bañados et al. 2014Bañados et al. , 2016;;Tang et al. 2017), but also includes quasars from a number of other recent Figure 1.A sky map of the 138 quasars at redshifts  ≥ 6.We restricted the observations to declinations  > −44 • (shaded regions) for convenience with the VLA.

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 ′′ (depending on source declination and the length of the observation).At  ∼ 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 seconds.
The observing programme was designed to achieve a 5 detection threshold of 100 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 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 Redshift 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.

IMAGING AND ANALYSIS
The imaging of the VLA L-band data was performed with the stacking algorithm provided by WSClean (version 3.3, Offringa et al. 2014;Offringa & Smirnov 2017).A subset (about 50%) of the quasars was imaged independently using the w-projection algorithm (Cornwell et al. 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 ′′ , resulting in an image size of 12244×12244 pixels.We set the WSClean Briggs robust parameter1 to −0.5, resulting in an average synthesised beam FWHM of ≈ 1.8 ′′ .Parallelising 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 preimaging RFI editing.
Following the imaging, we characterised the statistical properties of the images, which ultimately allow us to determine which images exhibit reliable detections of the quasar radio continuum.The rootmean-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 ± 5 Jy beam −1 across the cutout 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 synthesised 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   > 100 Jy of approximately 10 −4 sources per square arcsec.Hence, our search area of 6 × 6 pixels (3.28 square arcsecs) contains on average 3 × 10 −4 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 (  ∝    ) for the target quasars, where we adopt a typical spectral index of   = −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(Frey et al. , 2011;;Momjian et al. 2008Momjian et al. , 2018)), which find steep radio spectra for quasars at  ∼ 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  ,5 GHz at a rest-frame frequency of 5 GHz by applying a -correction (cf.Hogg 1999;Hogg et al. 2002;Condon & Matthews 2018) where   is the luminosity distance.If a source is detected with / ⪆ 10, it can be imaged in adjacent 64 MHz VLA subbands while retaining sufficient S/N to compute an in-band spectral index (see Section A in the Appendix for a discussion on calculated inband spectral indices).For such cases, we perform the extrapolation in equation 1 using the fitted power-law rather than the assumed one.

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 which 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 ( = 6.0) and J1034-1425 ( = 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 cutouts of the L-band continuum images of the eight likely detections and J0231-2850 are shown in Figure B1 of the Appendix.
There is no evidence for extended emission in the cutout-images.
We checked this more carefully by subtracting the scaled restoring beam from the cutouts 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.

Radio Loudness
The strength of the radio emission of a quasar relative to its optical emission is characterised by the radio-loudness parameter, .There are many definitions of  in the literature (see, e.g., Hao et al. 2014, for a review).Here, we adopt the most widely used definition  =  5 GHz /  ,4400 Å (Kellermann et al. 1989), where  5 GHz and  ,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;Liu et al. 2021;Gloudemans et al. 2021), therefore allowing a fair comparison with the results therein.A quasar is generally classified as radio-loud if  > 10, and as radio-quiet otherwise.
At  = 6, the rest-frame wavelength of 4400 Å corresponds to an observed wavelength of approximately 3 m.To mitigate potential large extrapolation errors,  ,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 1 photometric data in the ALL-WISE 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  1450 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 <  < 6.5.We adopted the same error for our flux densities that are extrapolated from  1450 , with the caveat that there may be individual objects with considerably larger extrapolation errors.

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.
Figure 3 shows the distribution of the derived optical and radio luminosities along with the regions that define radio-loud and radioquiet quasars, respectively.Quasars with  < 10 are classified as radio-quiet regardless of whether or not they are detected in our Lband observations.Quasars with  > 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 characterised as radio-loud or radio-quiet.For the remaining sources, 3 are radio-loud detections and 69 are robustly characterised as radioquiet.Including literature values and excluding the radio-selected quasar J1427+3312, we have in total 65 ambiguous, 4 radio-loud, and 69 radio-quiet quasars.
Using only sources that can be robustly classified as either radioloud or radio-quiet, we find a radio-loud fraction of 4/(69 + 4) = 5.5 ± 2.7%, where the error is estimated from Poisson statistics.If we include the ambiguous quasars and assume that they are all radioquiet, we obtain a lower estimate of the RLF, 4/138 = 2.9 ± 1.4%.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  not having occurred at time  is given by the Survival function () = Pr( > ).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  () = Pr( < ) = 1 − () of the radio-loudness parameter .We estimate  () by employing the Kaplan-Meier estimator (KM, Kaplan & Meier 1958) where   is the radio-loudness parameter of a detected source and   is the number of sources that have not been detected at  >   .This can be shown to be the maximum-likelihood estimator of  () (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  () rather than ().
Our estimate of  () is shown in Figure 4 and was calculated using the Python package lifelines (Davidson-Pilon 2019).The RLF thus obtained is 3.8 +6.2 −2.4 %, where the uncertainties are 95% confidence intervals.

Does the Radio-Loud Fraction Evolve with Redshift?
Jiang et al. ( 2007) used a sample of ≳ 30, 000 optically-selected quasars at redsfhits 0 <  ≤ 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% at  = 0.5 to ≈ 8% at  ≈ 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Å.The vertical error bars indicate the errors reported in the corresponding studies.Note that uncertainties from the literature are assumed to be 1 while ours represents a 95%-confidence interval.The marker showing our estimate of the RLF is located at the median redshift of the sample, while the others are placed at the centres of the corresponding redshift ranges.Also shown is is a prediction using the fit from Jiang et al. (2007) at 2500Å absolute magnitudes of  2500 = −27.7 (grey shaded) and a dashed line representing the inverse variance weighted mean of the RLFs across redshift.See text for discussion.
Figure 5 shows a comparison between our RLF and various literature estimates of the RLF at  > 4 (Bañados et al. 2015;Yang et al. 2016;Liu et al. 2021;Gloudemans et al. 2021;Lah et al. 2023;Ivezić et al. 2002;Stern et al. 2000), along with the 1 RLF bands predicted by the fit of Jiang et al. (2007) for  2500 = −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  ≈ 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  ≳ 6 are critically needed to test the above predictions.f The detection of J0231-2850 is not reliable, as the image shows non-thermal-like artefacts.We therefore treat it as a non-detection and leave its confirmation to future work.The upper limit provided for this source is at 3 significance.

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 emission 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 ′ from the optical positions.We stack in two separate redshift bins, 6.0 <  ≤ 6.15 and 6.15 <  < 7.6 (low- and high-, hereafter, respectively), which were chosen so that both bins contain a similar number of sources ( low = 68 and  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-, hereafter), are shown as the black histograms in the upper two panels of Figure 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- and the all- 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  ≈ 6 just below the detection threshold of this survey.The apparent excess is not seen in the high- 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-, high-, and all- samples, respectively).
The lowest panel of Figure 6 show the median-stacked images of the low-, high-, and all- samples.The stacked images for the low- and all- bins have RMS noise values of ≈ 5.2 Jy and ≈ 3.9 Jy, respectively, and show detections of the stacked continuum signal at S/N > 3 significance.For the high- bin, the stacked image has Table 2. Stacking results for the  > 6 quasar sample.The columns correspond to the three redshift bins used for stacking and are defined by the minimum, maximum and median redshifts  min ,  max ,  med , respectively.The RMS noise is that of the stacked image cutouts,  ,peak is the fitted peak flux density of the median stacked images,  5 GHz is the fitted median radio luminosity at rest-frame 5 GHz, and  is the fitted median radio-loudness parameters.The uncertainties of the stacked parameters correspond to the RMS noise in the stacked images.Upper limits represent 3 uncertainty levels.

Blank Sky Mean
,peak (Jy beam −1 ) 2.6 ± 3.9 0.4 ± 4.5 1.5 ± 3.0 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  ,peak = 13.8 ± 3.9 Jy, while the median rest-frame 5 GHz radio luminosity of the quasar sample is log  5 GHz /(W Hz −1 ) = 24.2 ± 0.1.We obtain a median radio-loudness parameter of  = 0.7 ± 0.2 for the entire sample, confirming that the undetected quasar population at  ≳ 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.

SUMMARY
We report VLA 1-2 GHz continuum imaging of a sample of 138 quasars at  ≈ 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  = 6.1, and detected seven additional quasars that have previously been characterised 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 3 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 %, where the uncertainties represent 95% 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  > 6 to date.However, while the central value of our RLF is lower than earlier estimates at  ≳ 4, Note that the detected sources are outside the range of these plots but still contribute to the summary statistics.The fitted values and uncertainties of the stacked quantities are listed in Table 2.
our measurement is formally consistent with the earlier studies within ≈ 2 significance.Our results are thus consistent with an unchanging RLF at  > 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  > 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  ≥ 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.

1.3
1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 Frequency (GHz) dashed lines represent the best-fit power laws and the grey shaded regions the corresponding 95% confidence intervals.The markers with no fill represent literature values (see Table A1.) quency over a small fractional bandwidth, one usually assumes a constant spectral index.
The problem of extrapolating from a continuum measurement is two-fold.Firstly, 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 synthesised continuum image may not be the same as the central frequency of the image if the spectrum is not flat (i.e.,  0 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  0 and   .
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 Figure 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.Firstly, 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   () 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 Table A1.Spectral index fitting parameters.For each quasar the first row corresponds to the fitted parameters from this work and the second row to the values found in the literature (see table notes for the references).

APPENDIX B: VLA L-BAND IMAGES
Figure B1 shows postage stamp cutouts 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).

Figure 2 .
Figure 2. The redshift and absolute rest-frame UV magnitude ( 1450 ) distributions of the quasar sample used in this work.

Figure 3 .
Figure 3.The rest-frame radio luminosity vs. rest-frame optical luminosity at 5 GHz and 4400 Å respectively, of the  > 6 quasar sample.The blue circles with downward arrows represent 3 upper limits and the red squares represent quasars detected in the L-band.The diagonal lines indicate where the radio-loudness () parameter takes values of 1, 10, and 100 respectively.Radio detected quasars at  > 10 are classified as radio-loud and all quasars (including non-detections) at  < 10 are classified as radio-quiet (shaded region).The horizontal magenta line at log  ,5 GHz = 25.5 W Hz −1 represents an alternative separation of radio-loud and radio-quiet quasars explored inJiang et al. (2007) andBañados et al. (2015).

Figure 5 .
Figure 5. Different estimates of the RLF plotted against redshift.The horizontal error bars indicate the range of redshifts covered by each measurement.The vertical error bars indicate the errors reported in the corresponding studies.Note that uncertainties from the literature are assumed to be 1 while ours represents a 95%-confidence interval.The marker showing our estimate of the RLF is located at the median redshift of the sample, while the others are placed at the centres of the corresponding redshift ranges.Also shown is is a prediction using the fit fromJiang et al. (2007) at 2500Å absolute magnitudes of  2500 = −27.7 (grey shaded) and a dashed line representing the inverse variance weighted mean of the RLFs across redshift.See text for discussion.

Figure 6 .
Figure 6.The distributions of the flux density  ,peak (top) and 5 GHz rest frame luminosity  ,5 GHz (middle row) computed from our 1-2 GHz continuum images at the optical positions of the  > 6 quasars.The three columns correspond to the redshift ranges 6.0 <  ≤ 6.15 (left), 6.15 <  < 7.6 (centre), and 6.0 <  < 7.6 (right), respectively.The grey-shaded region shows the amplitude distribution obtained on the "blank-sky" regions located 1 ′ from the optical positions of the quasars.The grey solid vertical lines indicate the means and the blue solid lines the best-fit Gaussians of the blank-sky distributions.The inserted box whisker plots show the medians (orange vertical line), the interquartile ranges (IQR, box), the 1.5 × IQR (whiskers) and outliers (red markers).The bottom row shows the median images of the corresponding redshift ranges.The contour levels correspond to 2 and 3 with negative values shown as dashed lines.Note that the detected sources are outside the range of these plots but still contribute to the summary statistics.The fitted values and uncertainties of the stacked quantities are listed in Table2.

Table 1 .
Radio and Optical Data of the  > 6 Quasar detections at 1-2 GHz The full table including all the upper limits are provided as supplementary material to this paper.a L-band continuum observations are typically reported at 1.4 GHz.However, the exact frequency associated with the flux density in the synthesised image cannot be known without knowledge of the spectral properties of the sources in the map.The central frequency of our observations is 1.68 GHz, but for a typical source with spectral index   = −0.75,we find 1.6 GHz to be a more representative frequency of the reported flux densities (cf.Section 3).