Open Access
Issue
A&A
Volume 668, December 2022
Article Number A27
Number of page(s) 19
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202244763
Published online 29 November 2022

© A. J. Gloudemans et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

High-redshift quasars (z >  5) are powerful tools for studying supermassive black hole (SMBH) formation and evolution and the phase transition from neutral to ionised hydrogen during the Epoch of Reionisation (EoR). In the past few decades, hundreds of quasars at z >  5 have been discovered in various optical and infrared surveys (e.g. Willott et al. 2007; Venemans et al. 2013; Bañados et al. 2016; Matsuoka et al. 2016; Wang et al. 2017; Miyazaki et al. 2018; Dey et al. 2019), which have been used to constrain SMBH growth (e.g. Onoue et al. 2019) and trace neutral hydrogen in the EoR by measuring, for example, the Lyα damping wing (e.g. Mortlock et al. 2011; Schroeder et al. 2013; Greig et al. 2017; Bañados et al. 2018).

Within this population, the subset of quasars that are also radio bright have been suggested to host more massive SMBHs at lower redshift (McLure & Jarvis 2004) and are important laboratories for studying the dominant physical mechanisms that contributed to the formation and subsequent quenching of massive galaxies, for example, gas accretion, galaxy merging, SMBH growth, and active galactic nucleus (AGN) feedback (Miley & De Breuck 2008; Saxena et al. 2017; Shao et al. 2020; Yamashita et al. 2020; Ichikawa et al. 2021; Uchiyama et al. 2022). Observations of the radio jets launched by these quasars can provide constraints on the timescales of early AGN activity by studying their synchrotron emission (Saxena et al. 2017; Shao et al. 2020). In addition, the physical mechanisms by which quasars are triggered and evolve over the Universe’s lifetime can be studied through comparisons of radio properties, such as radio sizes and spectral shapes, at high and low redshifts (Hook et al. 1998; Bañados et al. 2015; Gloudemans et al. 2021a).

The most radio-bright high-z quasars at low frequencies (≳10 mJy) can, in principle, be used to directly measure the neutral hydrogen content in the EoR by measuring the 21-cm absorption feature, because this line remains optically thin at z >  6 (Carilli et al. 2002; Mack & Wyithe 2012; Ciardi et al. 2015) and is shifted below 200 MHz at z ≳ 6. Such a measurement will have profound consequences for cosmology, since the rate and patchiness of reionisation and the contributors are still heavily debated (see e.g. Pentericci et al. 2014; Sobral et al. 2016; Naidu et al. 2020). A recent study of Bosman et al. (2022), using deep optical/NIR spectra of 30 high-z quasars, has shown that residual neutral hydrogen in the intergalactic medium (IGM) persists until z ∼ 5.3, which signifies a later end of cosmic reionisation than previously found (e.g. Becker et al. 2001; Fan et al. 2001, 2006; McGreer et al. 2013). To study the cosmic reionisation process universally via the 21-cm line, multiple sight lines are necessary, but progress has been halted by the small number of known bright high-z radio sources, with only a select number of radio-loud sources lying at z >  6 (McGreer et al. 2006; Willott et al. 2010; Belladitta et al. 2020; Bañados et al. 2021; Endsley et al. 2022b).

The new generation of low-frequency radio telescopes, such as the Low Frequency Array (LOFAR; van Haarlem et al. 2013), are able to cover large areas of the sky with sensitivities in the μJy regime, which has enabled high-z quasars to be systematically studied at low frequencies for the first time. The recently released LOFAR Two-Metre Sky Survey second data release (LoTSS-DR2; Shimwell et al. 2022) is a low-frequency radio survey of the northern sky at 120–168 MHz with a median root mean square (RMS) sensitivity of ∼83 μJy beam−1 and a resolution of 6″ covering over 5720 deg2. With this sensitivity, the survey is reaching more than an order of magnitude deeper (with ∼4× better angular resolution) than other large sky surveys at 150 MHz such as the TIFR GMRT Sky Survey (TGSS; Intema et al. 2017). Gloudemans et al. (2021a) found that 36% of the known high-z quasars are tentatively detected at > 2σ significance in LoTSS-DR2, showing the potential of using the new generation radio surveys to study this population.

In this work, we explore the potential of using LOFAR to discover new radio-loud (RL) quasars above z >  5. In principle, radio detections drastically reduce the number of stellar contaminants often found in optical and near-infrared (NIR) high-z quasar searches (e.g. Hewett et al. 2006; Findlay et al. 2012; Bañados et al. 2016) and therefore allow the exploration of new parameter spaces that were previously inaccessible due to the high contamination by stellar objects (see e.g. Bañados et al. 2015). Typically, high-z quasar candidates are selected in optical and NIR photometric surveys by probing the Lyman break and redward continuum. In this work we take a similar approach, but select only candidates with LOFAR detections, which allow us to search for quasars in a broad redshift range without requiring the full set of optical dropout bands and without drastically increasing the number of contaminants. Recently, Wagenveld et al. (2022) demonstrated the potential of finding quasars outside the classical optical colour space by discovering a quasar at z = 5.7 with an optical colour of i − z = 1.4, which is below the standard colour selection of i − z ≥ 1.5 or 2.0 (e.g. Bañados et al. 2016; Matsuoka et al. 2016), using a probabilistic candidate selection technique. This probabilistic candidate selection has also been performed by e.g. Mortlock et al. (2012), Matsuoka et al. (2016), and Nanni et al. (2022) and has been proven to be efficient in identifying high-z quasars without using traditional discrete colour cuts.

Recent wide-area high-z quasar searches have utilised the Sloan Digital Sky Survey (SDSS; York et al. 2000) and PANSTARRS-1 (PS1; Chambers et al. 2016) to select high-z quasars. However, the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys (Dey et al. 2019; hereafter referred to as ‘Legacy Surveys’) are reaching at least 1 magnitude deeper in the optical g-, r-, and z-band than PS1. By requiring a radio detection, dropout candidates between the r- and z-band can be identified without the need for an i-band detection to decrease the fraction of contaminants, therefore allowing us to search for high-z quasars in this deep large sky survey. A recent study of Wang et al. (2019) has also utilised the Legacy Surveys in combination with other optical and NIR surveys to search for new high-z quasars.

In addition to offering potential new probes of the process of reionisation, improved samples of the most RL quasars at z >  5 also enable statistical studies of the radio-loudness evolution, influence of the radio jet on the host galaxy, and quasar environment at high redshift. In this paper, we aim to constrain the relation between the radio luminosity and rest-frame UV properties such as Lyα emission. Previous studies of high-z radio galaxies (HzRGs) have found either a strong positive correlation between the Lyα and radio luminosity (Jarvis et al. 2001) or a weak correlation, consistent with no correlation (Saxena et al. 2019). However, these are based on biased samples of luminous radio galaxies with ultra-steep radio spectral indices and strong Lyα emission. Furthermore, we aim to investigate potential differences in the rest-frame UV properties of the RL and radio-quiet (RQ) high-z quasar population, since the existence of a radio-loudness dichotomy is still heavily debated (see e.g. Ivezić et al. 2002; Cirasuolo et al. 2003; White et al. 2007; Zamfir et al. 2008; Baloković et al. 2012; Beaklini et al. 2020; Gürkan et al. 2019; Macfarlane et al. 2021). Finally, the dust reddened red quasar population has been shown to have an enhanced fraction of radio-loud quasars (e.g. Fawcett et al. 2020; Rosario et al. 2020; Glikman et al. 2022). Therefore, in this study we aim to investigate whether our radio selection combined with loose optical colour criteria results in an increased fraction of red quasars.

This paper is structured as follows. In Sect. 2, we describe the high-z quasar candidate selection and the optical and radio observations on which this is based. In Sect. 3 we provide details on the conducted spectroscopic observations of our candidates, and in Sect. 4 we show their spectra (we discuss the non-quasar contaminants that our selection procedure produces in Appendix A). We present the optical and low-frequency radio properties of the newly discovered high-z quasar sample in Sect. 5. Finally, in Sect. 6 we summarise the results together with future prospects. Throughout this work, we assume a Λ-CDM cosmology with H0 = 70 km s−1 Mpc−1, ΩM = 0.3, and ΩΛ = 0.7 and use the AB magnitude system (Oke & Gunn 1983).

2. Candidate selection

Our high-z quasar candidates have been selected by a traditional dropout selection technique using the Legacy Surveys data combined with radio detections from LoTSS-DR2. In addition, we exploit the photometric redshifts (photo-zs) generated by Duncan (2022) and apply our own spectral energy distribution (SED) fitting to ensure a reliable high-z candidate sample. Here, we describe the specifics of the surveys and selection criteria used in this work to create our quasar candidate sample.

2.1. DESI legacy imaging surveys

The DESI Legacy Imaging Survey DR8 is an optical imaging survey providing g-, r-, and z-band imaging covering over 19 000 deg2 of the sky. This survey is the deepest wide-area optical survey to date and reaches median 5σ depths of 24.6, 24.0, and 23.3 mag in the g-, r-, and z-band filter, respectively. In addition to the optical photometry, the Legacy Surveys catalogues include mid-infrared photometry from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) at wavelengths of 3.4, 4.6, 12, and 22 μm from the W1, W2, W3, and W4 bands, respectively. We refer to Dey et al. (2019) for details on the photometric calibrations and measurements. The combination of depth and wide sky coverage of the Legacy Surveys together with the WISE photometry is well suited for discovering quasars at z >  5.

2.1.1. Photometric redshifts

To create our initial sample, we select optical sources from the Legacy Surveys classified as high-z galaxy candidates by the photometric redshift estimates of Duncan (2022). These photo-zs have been determined with a machine-learning approach, using sparse Gaussian processes augmented with Gaussian mixture models, together with cost-sensitive learning weight calculations. By combining these two methods, the photo-z estimates have been shown to be reliable for both the general galaxy population and rare or extreme populations at 4 <  z <  7. Adding spectroscopically confirmed quasars to the training sample and training on them separately (with a measured outlier fraction of OLF0.15 = 2.47%) makes this approach highly useful for high-z quasar searches (see Duncan 2022).

From the catalogues created by Duncan (2022), we selected optical sources with photo-z estimates of zphot >  3.5 within the LoTSS-DR2 footprint. To further ensure reliable photometry without blending of flux from neighbouring sources, the fracflux in the r- and z-band filter, which is defined by the profile-weighted fraction of the flux from other sources divided by the total flux, was set to be smaller than 0.2. In addition, we set maskbits=0 to ensure sources are not overlapping with the Legacy masked regions.

2.1.2. Colour selection criteria

After creating the initial sample, we apply optical and near-infrared (NIR) colour cuts, which are similar to those used in other high-z quasar searches (see e.g. Bañados et al. 2016; Jiang et al. 2016; Matsuoka et al. 2018). The wavelength range between the r- and z-band allows for probing the Lyman break of quasars in the redshift range of 4.9 ≲ z ≲ 6.6. Therefore, we selected sources with rz > 1.4 and required a non-detection in the g-band with a signal-to-noise ratio (S/N) below 3. In addition, to remove stellar contaminants we included a selection criterion redward of the Lyman break of −0.5 <  z − W1 <  2.5 (see Fig. 1). These selection criteria are sensitive to quasars over a broad redshift range, including the “quasar redshift gap” at z ∼ 5.3–5.7 created by conventional colour selections between the r- and i-band (see e.g. Yang et al. 2019). Furthermore, our loose colour selection allows for the possibility of finding overlooked quasars outside of the traditional colour regions. Our optical and NIR selection criteria can be summarised as follows:

(1)

(2)

(3)

thumbnail Fig. 1.

Optical and NIR colour selection criteria (given by the dashed lines) of our spectroscopically observed quasar candidates compared to literature quasars (5.0 ≤ z ≤ 7.5) and observed stellar M-, L-, and T-dwarfs.

https://docs.google.com/spreadsheets/d/1i98ft8g5mzPp2DNno0kcz4B9nzMxdpyz5UquAVhz-U8/edit?usp=sharing

2.2. LoTSS-DR2

For our radio selection, we use the LoTSS-DR2 144 MHz survey (see Shimwell et al. 2022), which contains a total of ∼4.4 million radio sources split over two fields. These are denoted as the 0hr and 13hr fields and cover sky areas of 1457 and 4178 deg2, with median rms values of 106 and 74 μJy beam−1, respectively.

The next step in our selection process was cross-matching with the > 5σ detected radio sources in the LoTSS-DR2 radio catalogues. The LoTSS-DR2 resolution is 6″ with an astrometric accuracy of ∼0.2″ at ≥20 mJy, and reaching an accuracy of 0.5″ at an S/N of 5. By cross-matching LoTSS-DR2 with the quasar catalogue from the Sloan Digital Sky Survey (SDSS; Pâris et al. 2018) at 0 ≤ z ≤ 5 and with known high-z quasars (Ross & Cross 2020), a search radius of 2″ was found to be sufficiently large to identify their radio counterparts without being dominated by false associations. Before cross-matching we have ∼1.6 million candidates in the LoTSS-DR2 region. Cross-matching the high-z quasar candidates within a radius of 2″ using the Legacy z-band and LoTSS-DR2 positions, drastically reduced our candidate sample to 310 potential high-z quasars, including 11 known quasars with z >  5 from Ross & Cross (2020). The radio-optical source separation is further discussed in Sect. 4.5. We note that we do not apply a radio compactness criterion to allow for marginally resolved radio emission from larger extended radio sources in the sample (see e.g. Miley 1971, 1968; Carilli et al. 1997; Saxena et al. 2019) and also accounting for aspects such as variable ionospheric smearing and incomplete deconvolution at low flux densities which can artificially enhance the source sizes (see Shimwell et al. 2022).

2.3. SED fitting

We carried out SED fitting on the sample of 310 potential high-z quasars using the photometric redshift code EAZY (Brammer et al. 2011), to further constrain their redshifts. Besides the Legacy Surveys photometry, we include forced photometric measurements of PS1 images in the g, r, i, z, y bands (see Chambers et al. 2016) and UKIDSS Hemisphere Survey J-band images (UHS; Dye et al. 2018) within a circular aperture (of 1.5″) at the Legacy z-band source position. The majority of our candidates (72% or 223/310) are detected in the PS1 catalogue, but including upper limits also for the non-detected sources below a 5σ significance contributes to constraining the SED shape. In the UHS J-band 33% (103/310) of our candidates are detected at ≥5σ significance in the catalogues, and 16% are not covered by UHS.

For the SED fitting, we used galaxy and quasar templates generated by Brown et al. (2014, 2019), which span a broad range of galaxy and quasar types and include additional reddening to take the full range of spectral slopes into account. From the resulting χ2 redshift distributions, we required integrated photo-z probabilities of at least 50% above z >  5 and zphot − 2 × zphot,err >  4.5 for each candidate, which yielded a total of 253 high-z candidates.

Finally, we visually inspected all remaining candidates to remove any spurious sources and sources with wrongly associated radio detections (e.g. cosmic rays and association with neighbouring optical sources), which resulted in a final sample of 142 candidates for spectroscopic follow-up observations.

All steps taken in the selection procedure are summarised in Table 1. The requirement of a radio detection drastically decreases the number of candidates from 1 573 757 to 310, which makes the final sample tractable for follow-up spectroscopy.

Table 1.

Overview of selection steps taken to create our candidate target list for spectroscopic follow-up.

3. Spectroscopic observations

To confirm our high-z candidates we have conducted spectroscopic observations using primarily the Faint Object Camera and Spectrograph (FOCAS; Kashikawa et al. 2002) on the Subaru Telescope on Mauna Kea, Hawaii and LRS2 (Chonis et al. 2016) on the Hobby-Eberly Telescope at the McDonald observatory in Texas, USA (HET; Ramsey et al. 1998; Hill et al. 2021). In addition, several of our targets were observed using the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on Keck telescope on Mauna Kea. The observational details of our newly discovered high-z quasars are summarised in Table 2. An additional 2 candidates were followed-up with the FOcal Reducer/low dispersion Spectrograph 2 (FORS2; Appenzeller & Rupprecht 1992) on the Very Large Telescope (VLT) at Paranal Observatory in Chile, but these were found to be low-redshift contaminants and therefore do not appear in Table 2.

Table 2.

Details of spectroscopic observations conducted of the newly discovered high-z quasars (sorted on RA).

The Subaru FOCAS observations were conducted using the VPH850 grating and SO58 filter (5800 − 10350 Å) with a slit width of 1″ and spectral resolution of R ∼ 600. The spectra were reduced following the standard procedure1 using Image Reduction and Analysis Facility (IRAF; Tody et al. 1986), which includes bias subtractions, flat fielding, distortion corrections, wavelength calibration and sky subtraction. Cosmic rays were removed using the algorithm L.A.COSMIC based on the Laplacian edge detection algorithm of van Dokkum (2001). The spectra were flux calibrated using standard stars HZ4, Feige34, and Feige110 observed on August 14th, November 25th, and December 27th 2021, respectively.

The HET observations reported here were obtained with the LRS2 integral field spectrograph. LRS2 comprises two double spectrographs each with a 6″ × 12″ integral field input separated by 100″ on sky: LRS2-B (3650 – 6950 Å) and LRS2-R (6450 – 10 500 Å). The integral field units have 100% fill-factor. There are two channels for each spectrograph: UV and “Orange” for LRS2-B and “Red” and “FarRed” for LRS2-R. We observed our targets with the LRS2-R spectrograph with exposure times of 1500–1800s. The HET observations were reduced using the HET LRS2 pipeline Panacea2 to perform the initial reductions including fiber extraction, wavelength calibration, astrometry, sky subtraction, and flux calibration. For the initial sky subtraction for each channel a biweight average (Beers et al. 1990) of the fibers excluding the target was used. Fringing and varying spectral resolution cause systematic residuals in the initial sky subtraction. We use a principle component analysis similar to the algorithm of Zurich Atmosphere Purge (ZAP; Soto et al. 2016) to model the sky subtraction residuals. After the final sky subtraction, we combine fiber spectra from the two channels, Red and FarRed, into a single data cube accounting for differential atmospheric refraction. We then identified the target quasar in each cube and extract our target spectra in a 1.5″ aperture. ILTJ0037+2410 has been observed with both HET LRS2 and Subaru FOCAS and can therefore be used to compare the resulting spectra. This comparison will be discussed in Sect. 4.4.1.

The Keck LRIS spectra were reduced using the data reduction package PypeIt (Prochaska et al. 2020). After basic calibrations (e.g. bias calibration, wavelength calibration), the object extraction (Horne 1986) and b-spline (Kelson 2003) sky subtraction were performed together. All individual 1D spectra were flux-calibrated with the standard star GD153.

To investigate the absolute flux scale accuracy of the reduced spectra, we convolved the spectra with the PS1 and Legacy Surveys photometric filter transmission profiles (including atmosphere) above the Lyman break and compare this to the photometric catalogue fluxes. The result is shown in Fig. 2 and the values for all the discovered quasars range between 0.58–1.78. These values were determined using the PS1 i-band if there was a photometric detection. Otherwise, the Legacy Surveys z-band values were used. The spectra are expected to result in lower fluxes than the photometry (absolute flux scale corrections > 1) due to slit losses for the Subaru/FOCAS and Keck/LRIS and losses due extraction of the quasars from the IFU for HET/LRS2. However, the quasar ILTJ2336+1842 at z = 6.6 has a correction factor of 0.58, which is likely due to the Lyman break being close to red end of the Legacy z-band and telluric lines heavily affecting that part of the spectrum. The other two quasars below 1 deviate only slightly with correction factors of 0.97 and 0.98. The deviations in absolute flux scale do not affect the redshift determination or the Lyα EW estimates of the quasars. However, they do influence the measured Lyα emission line fluxes and were therefore used in error propagation (see Sect. 4.2). Hence, we do not use the correction factors to calibrate the spectra, but do include them in the error estimation.

thumbnail Fig. 2.

Absolute flux scale correction factors between the quasar spectra and PS1 i-band and Legacy z-band photometry. A correction factor of > 1 means the spectrum continuum is lower than the photometric measurement. These corrections do not affect the redshift and Lyα EW estimates, but do decrease or increase the measured Lyα emission line flux.

4. New quasars

We have performed spectroscopic follow-up observations of 80 candidates in our sample (∼56%), which has resulted in the discovery of 20 new high-z quasars and the independent confirmation of 4 (which will be discussed below). Taking account of the 10 candidates that were not detected due to bad weather conditions, our success rate of spectroscopically confirming high-z quasars from our candidate sample was 34% (24/70). The contaminants we find are listed and discussed in Appendix A. We note that four of our quasars have been previously published in other work, but have been independently discovered in this work. These are ILTJ1419+5718 (McGreer et al. 2018), ILTJ1523+2935 (Wenzl et al. 2021), ILTJ1043+4048 and ILTJ1334+4750 (Lyke et al. 2020). Furthermore, some of our quasar candidates have been observed based on earlier catalogues and selection criteria (see Sect. 4.4 for details). Our final selection described in Sect. 2 includes 22 new quasars and 42 contaminants, yielding an identical success rate of ∼34% (22/64). This high success rate can be attributed to the radio selection by decreasing the number of stellar contaminants in the sample.

The spectra of the 24 newly discovered quasars are shown in Fig. 3. Their redshifts range from 4.9 to 6.6 (see Sect. 4.1). The discovery of quasars in this large redshift range was possible due to the wavelength gap between the Legacy survey’s r- and z-band filter. The Legacy, WISE, and LOFAR images of all new quasars are shown in Appendix B and quasar properties are summarised in Table 3. The optical spectra and measured properties are provided as a data product3.

thumbnail Fig. 3.

Observed-frame optical spectra obtained with Subaru/FOCAS, HET/LRS2, and Keck/LRIS of the newly confirmed quasars sorted on spectroscopic redshift. The total flux measurements at 144 MHz are also indicated on the left hand side of each spectrum. The detected emission lines are indicated by the grey dashed lines, which are Lyα+N V emission lines unless indicated otherwise.

Table 3.

Observed and measured physical properties of the discovered high-z quasars.

4.1. Redshifts

Determining accurate redshifts for high-z quasars at optical wavelengths is challenging, because the most prominent Lyα emission line is strongly affected by absorption. More accurate redshift estimates can be made by using atomic and molecular emission lines from the quasar, such as [CII], MgII or CO emission lines (e.g. Wang et al. 2010; Eilers et al. 2020; Schindler et al. 2020). However, the limited wavelength range of our observations does not include any of these emission lines. Therefore, in this work, we use a template fitting method of composite high-z quasar spectra (similar to Bañados et al. 2016) to constrain the redshift using the Lyman break and redward continuum. We fitted each quasar using the three composite quasar templates constructed by Bañados et al. (2016) for weak, moderate, and strong Lyα emission and the composite quasar template constructed by Selsing et al. (2016). Typically, the redshift estimates following from the χ2 minimisation of these templates agree well with maximum deviations of ∼0.05, but for our final redshift estimates we used the best fitting template. To determine errors on the spectroscopic redshift estimates, we convert the χ2 values as a function of redshift to a probability density distribution and fit the peak of the distribution using a Gaussian fit. The error on the redshift is then given by the standard deviation of the best Gaussian fit. This method yields errors on the zspec between 0.01–0.06 with a median of 0.03 (see Table 3). We emphasise that follow-up (near-)infrared or sub-millimetre observations are necessary to obtain more precise systemic redshifts for these quasars.

4.2. Emission lines

The Lyα emission line at 1216 Å is the most prominent line we detect in the optical spectra of our quasars. Additionally, for two quasars we also detect Lyβ, O I, Si IV, and C IV emission lines (see Sects. 4.4.4 and 4.4.5 for details). The strength of the Lyα emission line varies considerably between the different quasars (see Fig. 3) and also contains the contribution from the N V emission line at ∼1240 Å and Si II line at 1263 Å, with N V being the most dominant. To investigate the emission line properties as a function of redshift and radio luminosity, we estimate the Lyα+N V luminosity and equivalent width (EW) using the method of Diamond-Stanic et al. (2009). This method entails fitting a power-law of the form fλ = C × λβ to the continuum between the rest-frame wavelength regions uncontaminanted by emission lines in the range of [1285, 1295], [1315, 1325], [1340, 1375], [1425, 1470], and [1680, 1710] Å. However, note that the final wavelength range of [1680, 1710] is not covered by all our quasars above z >  5. We fix the slope of the power-law fit β to the values −1.5 and −1.7 (as done by Bañados et al. 2016) to avoid unphysical values of the continuum slope caused by noise and telluric contamination. We take into account the error on the spectroscopic redshift by fitting the power-law 10 000 times to the rest-frame spectrum determined by redshifts drawn from a normal distribution around the zspec. The Lyα+N V fluxes and EWs are determined by integrating between 1160–1290 Å above the continuum using the power-law fits (and averaging over the β = −1.5 and −1.7 solutions). For each source, the final estimated Lyα flux and EW are given by the 50th percentile on the resulting distributions and the errors by the 16th and 84th percentiles. As discussed in Sect. 3, the measured line fluxes are potentially affected by the inaccurate absolute flux scale corrections of the spectra, therefore the correction factors from Fig. 2 have been also been included in the errors on the Lyα+N V flux.

4.3. Rest-frame UV magnitudes

To determine the UV rest-frame magnitudes of our quasars we use a similar SED fitting procedure on their photometry as described in Sect. 2.3, but fix the templates to the spectroscopic redshift. From the best fits we derive the UV rest-frame magnitudes at 1450, 2500, and 4400 Å by measuring the average flux with tophat filters centred on the respective central wavelength with a width of 100 Å. To obtain the final UV magnitudes and their uncertainties we adopt a simple Monte Carlo method, where we duplicate each source 500 times and randomly adjust the fluxes by drawing from a normal distribution around each catalogue flux value with the standard deviation determined by the catalogue flux errors. The UV magnitudes are then given by the 50th percentile values, and the errors by the 16th and 84th percentile of the resulting distribution. We note that the intrinsic UV magnitudes are brighter than the observed ones due to dust reddening (e.g. Calistro Rivera et al. 2021). The dust reddening of our quasars is discussed in Sect. 5.4.

4.4. Notes on individual quasars

Here we discuss some of our newly discovered quasars that have notable properties, such as detected emission lines in addition to Lyα, and quasars that have not been selected by our main selection method as described in Sect. 2.

4.4.1. ILTJ0037+2410

The quasar ILTJ0037+2410 was observed with both by HET and Subaru, allowing a comparison of the absolute flux scale and the wavelength calibration. Both the derived Lyα+N V luminosities and the measured EWLyα+NV are in good agreement with luminosities of 10 and 10 erg s−1 and EWs of 14 and 16 for the Subaru and HET spectrum, respectively. The spectroscopic redshifts also agree very well with a difference of ∼0.01 from our fitting procedure (see Sect. 3).

4.4.2. ILTJ1013+3518

This quasar was not selected by the method described in Sect. 2, because it was omitted from the catalogues of Duncan (2022) due to its unreliable photometry caused by contamination of a neighbouring object. However, it was selected and observed prior to the creation of this final sample and confirmed to have a redshift of z = 5.03.

4.4.3. ILTJ1350+3748

This quasar has been identified as a dropout and confirmed using Keck/LRIS in 2017 (see also Banados et al. in prep, who independently discovered it).

4.4.4. ILTJ0121+2940

This quasar and ILTJ1043+4048 are the only quasars where we have detected additional emission lines. In the spectrum of ILTJ0121+2940 we have identified as Lyβ, O I, Si IV, and C IV, which can also be used to determine the spectroscopic redshift. A redshift of z = 5.24 ± 0.01 would match the Lyβ, O I, Si IV emission lines with a blueshift of C IV. However, this redshift does not match the Lyα line at ∼7700 Å, which would be shifted to ∼7585 Å for z = 5.24 and our template fitting procedure yielded a redshift of z = 5.34 ± 0.02 instead (see Sect. 4.1). Coincidentally, this wavelength region (7529–7600 Å) had to be masked out, due to a cosmic ray at that exact location. Therefore, the redshift of this quasar remains unclear and both redshift estimates are quoted in Table 3. For subsequent plots we adopt the redshift determined by the Lyβ, O I, Si IV emission lines of z = 5.24 ± 0.01. We note that potentially strong Lyα absorption cannot be explained by an IGM damping wing, since the quasar is at z <  5.5 (see e.g. Wolfe et al. 2005), nor by an absorber near the quasar, because of the presence of Lyβ emission (e.g. Bañados et al. 2019). This quasar could potentially be a broad absorption line (BAL) quasar with narrow Lyα (e.g. Ross et al. 2015) or extremely high velocity outflows (e.g. Rodríguez Hidalgo et al. 2020) causing Lyα and N V absorption.

The detection of the C IVλλ1548.2,1550.8 doublet allows for estimating the SMBH mass by measuring the full width at half maximum (FWHM) and continuum luminosity as often used to derive BH masses at high-z (see e.g. Vestergaard & Peterson 2006; Kelly et al. 2010; Park et al. 2013). To determine the continuum around the C IV line, we fit a power-law to the rest-frame spectrum using the same method as described in Sect. 4.2. Subsequently, the continuum luminosity is calculated using the power-law fit between 1500–1600 Å and the FWHM of C IV is determined by fitting a Gaussian profile to the continuum subtracted spectrum. The virial BH mass is finally calculated by

(4)

with (a, b, c)=(0.66, 0.53, 2.0) as determined by Zuo et al. (2020) using the calibrations from Vestergaard & Peterson (2006). This results in an estimated BH mass for ILTJ0121+2940 of , which is in line with previous high-z SMBH masses found by e.g. Fan et al. (2001) and Jiang et al. (2007). The errors on this BH mass estimate include the absolute flux scale uncertainties (see Fig. 2) and the redshift uncertainties, which both influence the continuum luminosity and subsequently the derived BH mass.

Note that this BH mass might be overestimated due to blueshifting of the C IV emission line, of ∼10 000 km s−1 for z = 5.24 assuming equal contribution of both components of the C IV with a rest-frame wavelength of 1549.98 Å (see Coatman et al. 2016). This blueshifting is likely caused by strong outflows and a high Eddington luminosity ratio (see e.g. Coatman et al. 2016; Marziani et al. 2019; Zuo et al. 2020).

4.4.5. ILTJ1043+4048

Similarly to ILTJ1013+3518, this quasar was omitted from the final sample due to a g-band S/N of > 5σ in the Legacy Surveys DR8 catalogue used in Duncan (2022). However, in the Legacy Surveys DR7 catalogue the quasar has a negative g-band S/N of −4.9 and was therefore previously selected for follow-up observations.

Also for this quasar we identify Lyβ, O I, Si IV, and C IV in the spectrum. In this case the observed wavelengths of the lines are consistent with the redshift determined by template fitting of z = 4.89 ± 0.03. Performing the same fitting procedure of the C IV line yields supermassive black hole mass of , which is slightly higher than expected from its absolute UV magnitude () when assuming a constant Eddington ratio and bolometric correction (Inayoshi et al. 2020). However, similar masses have been reported for other quasars around z ∼ 5, but it is at the high end of the expected SMBH mass range (see e.g. Vestergaard et al. 2008; Wang et al. 2015; Wu et al. 2015; Sbarrato et al. 2021). Observations of the Mg II line could provide a more accurate measurement of the SMBH mass of this source, due to the effects of outflows on the C IV emission line (e.g. Coatman et al. 2016; Onoue et al. 2019; Shen et al. 2019).

4.5. Optical-radio separation

Our observing campaign has been successful in selecting high-z quasars with a success rate of 34%. However, we could have further improved this by lowering the separation limit between the optical and radio positions. The separation of the quasars and contaminants as a function of the association likelihood ratio (LR; see e.g. Williams et al. 2019) are shown in Fig. 4, which clearly shows a downward trend in LR with separation, as expected. These LRs have been calculated for LoTSS-DR2 to identify the most likely optical counterparts in the Legacy Surveys for each non-complex detected radio source (see Hardcastle et al., in prep.). In short, the LR is defined as the ratio between the probability of an object being the true counterpart and being a random interloper. To calculate the LRs for optical counterparts in LoTSS-DR2 a maximum likelihood method is used (based on Williams et al. 2019), which uses a prior probability that a radio source has an optical counterpart with a certain optical colour and magnitude, the surface density of sources with specific colours and magnitudes, and the probability distribution function for the offset between the radio and optical sources. To avoid biases due to galaxy clustering, the fields around each radio source are compared to other random sky positions. Source types are not considered in the LR calculation. For more details on the LRs we refer to Williams et al. (2019) and Hardcastle et al., (in prep.).

thumbnail Fig. 4.

Association likelihood ratios versus optical-radio separations of our observed sample. The blue histograms show the distribution for the quasar sample from this work and literature combined. We demonstrate that a similar distribution (p-value = 0.5) of separations can be created by a simple Monte Carlo simulation using the positional offsets from both the Legacy and LoTSS-DR2 catalogue. From this distribution we derive the chance of a quasar having an optical-radio separation larger than 1″ to be ∼2% when assuming no intrinsic difference in the position of the radio and optical emission.

In the blue histogram in the top panel of Fig. 4 the distribution of separations is shown for the quasar sample from this work combined with the quasar literature sample. We demonstrate that a similar distribution can be created by a simple Monte Carlo simulation using the positional offsets from both the Legacy and LoTSS-DR2 catalogue. The average astrometric precision of the LOFAR sources in the LoTSS-DR2 region is σRA = 0.22″ and σDec = 0.20″ for S/N >  20 and the precision decreases gradually for lower S/N up to a maximum of 0.5″ for an S/N of ∼5 (Shimwell et al. 2022). The optical positions in the Legacy Surveys have a high (median) astrometric accuracy of ∼0.02″(Dey et al. 2019), therefore the positional uncertainties are dominated by the LOFAR astrometric offsets and by the fact that some of these have extended jet emission that can offset the radio centroid from the true core. The positional uncertainties for our quasar sample and literature sample are on average σRA = 0.40″ and σDec = 0.29″. The Monte Carlo distribution created using these positional uncertainties is statistically similar to the quasar distribution with a p-value of 0.5. From this distribution we derive the chance of a quasar having an optical-radio separation larger than 1″ to be ∼2% and larger than 1.2″ to be less than 1%. This is in agreement with all newly discovered and known quasars having optical-radio separation offsets smaller than 1.2″. Therefore taking this smaller offset (< 1.2″) as a selection criterion of our sample, the success rate increases to 42%. The increased contamination fraction with separation indicates that there might be false associations amongst the contaminants and shows that quasars at these redshifts typically do not have large offsets (> 1″) in their radio emission. For our assumed cosmology, an offset of 1″ corresponds to 5.7 ± 2.9 kpc physical size at z = 6 with the error determined by the astrometric accuracy of LOFAR. High-resolution radio observations are necessary to study the resolved jet structure (e.g. using the international LOFAR baselines; Morabito et al. 2022 or the Very Long Baseline Array (VLBA); Momjian et al. 2018).

The contaminants we find in this work are listed in Table A.1 and briefly discussed in Appendix A. Regardless of the nature of the contaminants, for future high-z quasars searches our results support the use of stricter optical-radio separation criteria to ensure maximum purity of the candidate sample. Based on the currently known radio-detected quasars, a separation of 1.2″ would be sufficient to retrieve all quasars in the LoTSS-DR2 survey.

5. New quasar sample properties

An overview of the observed properties of the newly discovered quasars is given in Table 3. To put this new sample into the context of the previously known quasar population, we examine the radio and optical properties and compare them to the known high-z quasars from literature work. For our literature sample, we primarily utilise the source list generated by Ross & Cross (2020) with low-frequency radio properties determined by Gloudemans et al. (2021a) in the LoTSS-DR2 region, and include other recently published radio-loud quasars.

5.1. Low-frequency radio properties

To measure the radio luminosities and radio-loudness of the new quasar sample from the observed LOFAR fluxes, we need to determine their radio spectral indices (using a simple power-law of Sν = να). Previous work of Gloudemans et al. (2021a) has shown that the median spectral index of the known quasar population between 144 MHz and 1.4 GHz is approximately −0.3 ± 0.1. To obtain spectral indices for our newly discov- ered quasars we cross-match our sources (within 2″) with the VLA FIRST survey at 1.4 GHz (Becker et al. 1994), which yields detections for 4 quasars. Additionally, to be able to also constrain the spectral index on marginally detected sources we perform forced photometry by taking the peak flux, which results in 1 additional ∼3σ measurement of ILTJ1334+4750.

The spectral indices derived from these detections are shown in Fig. 5, with 3σ upper limits derived for the non-detected quasars (< 3σ) in FIRST. We could not obtain a measurement for nine out of 24 quasars, which were outside of the FIRST survey footprint. We repeat this procedure using the Very Large Array Sky Survey (VLASS; Condon et al. 1998) at 2–4 GHz, which yields 3 detected sources in the VLASS catalogue and three marginally detected sources from measuring the peak flux. Since the VLASS images are known to systematically underestimate the flux densities by ≈15%, we have applied this correction to our peak fluxes. These six detected quasars in VLASS resulted in additional spectral index measurements for ILTJ1523+2935 and ILTJ2201+2338 (see Fig. 5). All measured fluxes and spectral indices are given in Table 3. From this Fig. 5 it is apparent there is a large scatter in spectral index around the median value of −0.3 ± 0.1, which is determined in Gloudemans et al. (2021a) by median stacking LOFAR and FIRST images of 93 known quasars at z >  5. Therefore, in this work we assume α = −0.3 ± 0.1 if there is no FIRST or VLASS detection. Finally, we note some of these quasars could also be peaked-spectrum sources with a turnover in the radio spectrum (see e.g. review by Urry & Padovani 1995). However, the spectral indices derived from VLASS for the other four quasars agree very well with the FIRST spectral indices (with differences of 0.04–0.15), indicating no significant spectral curvature for these sources at GHz frequencies. For the other quasars more radio observations are necessary to constrain the radio spectrum shape.

thumbnail Fig. 5.

Radio spectral index between 144 MHz (LoTSS-DR2) and 1.4 GHz (VLA FIRST) for the discovered quasars in the FIRST and VLASS survey area. Upper limits (indicated by downward-pointing triangles) are given for the quasars non-detected in FIRST. The dashed grey line shows the stacked spectral index obtained by Gloudemans et al. (2021a) using a sample of 115 known quasars at z >  5. ILTJ2201+2338 is outside of the FIRST footprint, but is detected in VLASS at 2–4 GHz, which enables a spectral index measurement. ILTJ1523+2935 is only marginally detected in FIRST (S/N ∼ 2), but is more securely detected in VLASS (S/N = 3.4).

The radio luminosities of the newly discovered quasars as a function of redshift are shown in Fig. 6. Our sample probes on average higher radio luminosities than the previously known quasars, but most of the brightest radio sources in our candidate sample were identified as low-z contaminants. However, the quasar ILTJ1601+3102 at z = 4.9 has an exceptionally high radio flux of 69 mJy. Besides the known quasar sample in LoTSS-DR2 from Gloudemans et al. (2021a), we also add other recently published radio-loud quasars: PSO J172+18 (Bañados et al. 2021), PSO J0309+27 (Belladitta et al. 2020), VIK J2318-3113 (Ighina et al. 2021), and COS-87259 (Endsley et al. 2022b). PSO J0309+27 has been detected by LOFAR outside of the current LoTSS-DR2 area with a radio flux of 52 mJy at 144 MHz (Spingola et al. in prep) and has been classified as blazar by Belladitta et al. (2020). The highest redshift radio-loud quasar COS-87259 has also been detected by LOFAR outside of the LoTSS-DR2 footprint with a 144 MHz flux of 0.48 ± 0.18 mJy (Endsley et al. 2022a). The other very high-z redshift radio-loud quasar PSO J172+18 at z = 6.8 remains non-detected by LOFAR with an RMS noise limit of ∼0.3 mJy. This non-detection by LOFAR puts a new constraint on the radio spectral index of α >  0.2 using the VLA-L band measurement at 1.52 GHz of 0.51 ± 0.02 mJy (see Bañados et al. 2021), since it was previously only non-detected by TGSS with a flux limit of 8.5 mJy. The spectral index between the VLA-L and VLA-S band was measured to be negative with . This quasar therefore shows a turnover in the spectrum. The quasar VIK J2318−3113 was the previous highest redshift radio-loud quasar, but this region is not observable by LOFAR due to the low declination and therefore a FIRST flux upper limit is given here.

thumbnail Fig. 6.

Radio luminosities at 144 MHz of newly discovered quasars in this work compared to the literature sample of known high-z quasars detected in LoTSS-DR2. Other known z >  6 radio-loud quasars PSO J172+18 (Bañados et al. 2021) and VIK J2318−3113 (Ighina et al. 2021) remain non-detected and not covered by LOFAR, respectively. PSO J0309+27 (Belladitta et al. 2020) and COS-87259 (Endsley et al. 2022b) are detected by LOFAR outside of the LoTSS-DR2 area.

The radio-loudness of quasars is traditionally defined by the ratio of the radio (5 GHz) and optical (4400Å) flux densities (radio-loudness criterion; R = f5GHz/f4400Å >  10; e.g. Kellermann et al. 1989). Therefore, to compare the radio-loudness of our new quasars with the known quasar population we estimate the radio flux at 5 GHz rest-frame by extrapolating from the observed 144 MHz flux (∼850 − 1100 MHz in rest-frame) assuming our determined radio spectral indices from Fig. 5. The resulting radio-loudness (R) compared to known quasars is shown in the left panel of Fig. 7. The radio-loudness of our quasars is in the range of 6-1080. Traditionally, quasars are defined as radio-loud if R >  10, which is satisfied by 21 out of the 24 quasars with 23 out of the 24 quasars having R >  10 within the error margin. The blazar PSO J0309+27 at z = 6.1 is still the most radio-loud quasar-like source at z >  6. We note that this result is highly dependent on the assumed radio spectral indices, which could not be determined for the full sample. For our sample of quasars a change of 0.1 in the spectral index gives rise to a difference of ∼16–19% in the radio-loudness. These errors are taken into account in Fig. 7. Follow-up radio observations at GHz frequencies are necessary to constrain the spectral indices (and potential curvature) and accurately determine the radio-loudness. This is also the case for the LOFAR detected known quasar sample, for which Gloudemans et al. (2021a) assumed a spectral index of −0.3 ± 0.1 for all non-detected quasars in the FIRST survey.

thumbnail Fig. 7.

Left: radio-loudness (R = f5 GHz/f4400 Å) as a function of spectroscopic redshift for the discovered quasars in this work and the literature quasar sample. Our quasars span a broad range in radio-loudness. Traditionally, quasars are defined as radio-loud if R >  10, which is satisfied by 21/24 of our newly discovered quasars. Right: UV rest-frame absolute magnitude at 1450 Å versus radio luminosity at 144 MHz. The new quasar sample probes fainter UV magnitudes and higher average radio luminosities than previously known quasar samples.

As shown in the right panel of Fig. 7, our quasars probe a new parameter space of fainter rest-frame UV magnitudes and on average higher radio luminosities (at 144 MHz) than previously known LOFAR detected high-z quasars. To compare, the new quasars have a median L144MHz and M1450 of 1026.02 W Hz−1 and −26.06, respectively, whereas the known quasars detected by LOFAR have a median L144MHz and M1450 of 1025.58 W Hz−1 and −26.68 respectively. Our quasars are still not radio bright enough (> 10 mJy) at high enough redshift (z >  6) for current facilities such as LOFAR to be able to measure the 21-cm absorption line. However, these measurements might be possible with similar quasars in the Southern Hemisphere with the future Square Kilometer Array (see Ciardi et al. 2015). We note that populations of fainter quasars have been found at z >  5.7 in optical surveys such as the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs; e.g. Matsuoka et al. 2016, 2022), but these generally remain non-detected by LOFAR (Gloudemans et al. 2021b).

Other recent studies in the radio regime have found that about 10% of the high-z quasars are radio-loud using the classical definition of the ratio of the radio and optical fluxes f5 GHz/f4400 Å >  10 (Bañados et al. 2015; Gloudemans et al. 2021a; Liu et al. 2021). However, in the current study we have selected quasar candidates that are relatively radio bright and therefore our fraction of radio-loud quasars is significantly higher. There is still ongoing discussion on whether the classes of radio-loud and -quiet quasars are probing different populations, since different studies find both evidence for bimodality (e.g. Ivezić et al. 2002; White et al. 2007; Zamfir et al. 2008; Beaklini et al. 2020) and against bimodality (e.g. Cirasuolo et al. 2003; Baloković et al. 2012) of the radio luminosity distribution of quasars. Recently, Gürkan et al. (2019) showed that this bimodality is not favoured by the data using a large sample of quasars (∼50 000) covered by the 150 MHz LOFAR surveys. Also, recent work of Macfarlane et al. (2021) shows that the radio luminosity distribution can be modelled using a combination of radio emission from jets (with a wide distribution of radio powers) and star formation with a smooth transition between the radio-loud and -quiet regime and without the need for bimodality.

5.2. Quasar colours

With this campaign we have selected quasars in a broad redshift range (4.9 ≲ z ≲ 6.6) and colour space (see Sect. 2.1.2), whereas previous high-z quasar searches have focused on for example 4.7 <  z <  5.4 (e.g. Wang et al. 2016), z ≳ 5.6 (e.g. Bañados et al. 2016; Jiang et al. 2016; Wang et al. 2019), or quasars in the redshift gap between z ∼ 5.3 − 5.7 (e.g. Yang et al. 2017, 2019). Despite selecting our quasar candidates using fewer and less conservative colour restrictions, our overall success rate is 34% and we have discovered 7 new quasars in the redshift gap. In other high-z quasar searches the success rate ranges between ∼5–64% depending on the selection method and redshift range (see e.g. Wang et al. 2016, 2019; Yang et al. 2017; Matsuoka et al. 2019; Nanni et al. 2022; Wagenveld et al. 2022), however this number is not always published. Our success rate can be attributed to the radio detections drastically decreasing the number of contaminants in our sample and this work therefore demonstrates that a relatively clean high-z quasar candidate sample can be created using radio detections without requiring the full set of optical and NIR filters normally employed at these redshifts. We note the WISE photometry has also contributed significantly in constraining the photometric redshifts in this work.

The optical and near-infrared colours of the new quasars do not significantly deviate from the known quasar population, as is apparent in Figs. 1 and 8. All candidates that were spectroscopically observed outside of the traditional colour space have been identified as low-z contaminants, which are discussed in Appendix A.

thumbnail Fig. 8.

Near-infrared colours of the discovered quasars in this work and contaminants compared to the literature sample of quasars and stellar dwarfs5 (for which the PS-1 z-band magnitude is used in the right panel). Despite a looser colour selection, the new quasars exhibit colours consistent with the known quasar population.

5.3. Rest-frame UV spectral properties

In this section, we compare the Lyα+N V line strengths of our RL quasars to known RQ quasars to investigate whether they may probe different populations. We compare our EW distribution to the distribution obtained by Bañados et al. (2016) for PS1 selected quasars at z >  5.6 in Fig. 9. A Kolmogorov-Smirnov (KS) test yields a p-value of 0.10 indicating no statistically significant deviation between our EW distribution and the PS1 selected distribution from Bañados et al. (2016). As stated in Bañados et al. (2016), the EW distribution determined for SDSS selected quasars at z >  3 from Diamond-Stanic et al. (2009) is likely systematically higher due to lower IGM absorption of Lyα at lower redshift. From our sample, 9 out of 24 quasars (38%) have EW values below the defined weak line quasar boundary of 15.4 Å (see Diamond-Stanic et al. 2009), which is almost 3 times higher than the fraction of weak line quasars of 14% (16/117) reported from the PS1 quasar sample of Bañados et al. (2016). This could indicate a potential difference in observed Lyα emission between radio-quiet and -loud high-z quasars, but these are still low number statistics. However, the consistency of the overall shape the EW distribution of our radio-loud quasar sample with the optically selected quasars from Bañados et al. (2016) potentially suggests there is no strong radio-loudness dichotomy between the radio-loud and radio-quiet quasar population in the sense that the optical quasar properties are similar.

thumbnail Fig. 9.

Derived Lyα+N V emission line properties of the discovered quasars. Left: normalised Lyα+N V EW distribution of the discovered quasars in this work compared to the distribution of PS1 discovered quasars at z >  5.6 by Bañados et al. (2016) and SDSS discovered quasars at z >  3 by Diamond-Stanic et al. (2009). Our EW distribution is consistent with the optically selected quasars sample of Bañados et al. (2016). Right: Lyα+N V luminosity versus the radio luminosity at 144 MHz for the quasars in this work, showing little to no correlation (see Sect. 5.3).

A previous study of high-z radio galaxies (HzRGs) by Jarvis et al. (2001) found a strong positive correlation between LLyα and L151 for their sample of 35 HzRGs at z >  1.75, however Saxena et al. (2019) only found a weak positive correlation, consistent with no correlation, for their sample consisting of 10 new HzRGs at z >  2 and literature HzRGs. To investigate if such a correlation exists for our quasar sample, we plot the radio luminosity against the Lyα luminosity in the right panel of Fig. 9. We find a weak negative correlation between the Lyα and radio luminosity with a Pearson correlation coefficient of r = −0.12. We find a very weak negative evolution of Lyα luminosity with redshift (r = −0.05), but a slightly stronger correlation of radio luminosity with redshift (r = −0.24). However, the results from this quasar sample cannot be directly compared to the radio galaxies, since the rest-frame UV and radio emission we are probing are driven by accretion that is varying on shorter timescales than radio galaxies, which generally have more extended radio-lobes (see e.g. Miley & De Breuck 2008; Saxena et al. 2019; Nyland et al. 2020).

Although we do not see any strong correlation between Lyα and radio luminosity from our new sample of faint quasars, we note that any interpretation is limited by our relatively small sample size. In addition, interpreting the emission line properties of Lyα at high redshift is complicated, since resonant scattering becomes increasingly important, and therefore the line profile depends on both the intrinsic emission and the environment of the quasar or HzRG (e.g. van Ojik et al. 1997). Future NIR spectroscopy (targeting broad emission lines such as C IV and Mg II) will enable the measurement of SMBH masses and Eddington ratios of these sources (see e.g. Onoue et al. 2019) to constrain the timescales of SMBH growth and systemic redshifts to characterise their Lyα damping wings and outflows.

5.4. Red quasars

Following Glikman et al. (2012), high-z quasars can be classified as red when the dust reddening parameter E(B − V)> 0.1. The dust reddened red quasar population has been shown to have an enhanced fraction of radio-loud quasars (e.g. Fawcett et al. 2020; Rosario et al. 2020; Glikman et al. 2022). Therefore, to investigate whether our radio selection combined with loose optical colour criteria has resulted in an increased fraction of red quasars, we examine the dust reddening properties of our sample. Previous observational studies using the Atacama Large Millimeter Array (ALMA) have found abundant dust in quasar host galaxies (see e.g. Wang et al. 2013; Venemans et al. 2016; Izumi et al. 2019), however the fraction of extinguished quasar radiation caused by dust and whether this represents an evolutionary stage in a quasar’s lifetime is still unknown. Furthermore, the FIR and millimetre emission of radio-loud quasars is likely not only caused by cold dust but also by synchrotron emission (Rojas-Ruiz et al. 2021; Khusanova et al. 2022).

To determine the E(B − V) for our quasars, we used a template fitting method on the available optical and (N)IR photometry. We fitted the composite quasar spectrum of Selsing et al. (2016) to the photometric data points redward of the Lyman break and apply both the SMC extinction law of Pei (1992) with RV = 2.93 and the starburst galaxies extinction curve of Calzetti et al. (2000) with RV = 3.1 for different values of E(B − V). The SMC dust reddening law is assumed to be similar to high-z galaxy environments due to their low metallicity (see e.g. Hopkins et al. 2004; Maiolino et al. 2004; Gallerani et al. 2010). To be able to fit the quasar spectra to the observed photometry we convolve the model spectra with the available photometric filter transmission profiles, which includes absorption by the atmosphere. Subsequently, the template scaling is calculated analytically across all bands. To determine the best fitting model a minimum χ2 is calculated for each template with E(B − V) values in the range of −0.3–1.0 (with a step size of 0.005). We also fit negative values of E(B − V) to allow for intrinsically bluer quasars than average. To compare the preference for each model we calculate the Bayesian Information Criterion (BIC), which is given by

(5)

in the case of Gaussian distributed errors, with k the number of model parameters and N the number of photometric data points. When comparing two models, a ΔBIC > 6 or ΔBIC < −6 gives a strong positive or negative evidence, respectively, for one of the models being preferred. Therefore, we compare the best fitting dust extinction models to the Selsing et al. (2016) composite quasar spectrum fit and assign a E(B − V) value of zero if this condition is not met. Furthermore, we do not include quasars with 3 or less photometric data points, to ensure a reliable fit.

[5]https://doi.org/10.5281/zenodo.4169085

The resulting E(B − V) values using the SMC extinction law are shown in Fig. 10. Only 1 quasar, ILTJ0121+2940, satisfies the criterion E(B − V)> 0.1 for both possible redshifts z = 5.24 and 5.34 as discussed in Sect. 4.4.4. However, to examine the E(B − V) values of our radio selected quasars compared to optically selected high-z quasar sample, we repeat the same fitting procedure using the z >  5 quasar sample of Ross & Cross (2020), which is also shown in Fig. 10. We divide the literature sample in radio-loud (RL) and radio-quiet (RQ) using the record of known RL quasars (see Bañados et al. 2021 and Gloudemans et al. 2021a). We note that the radio-quiet sample likely still contains a small fraction of RL quasars, which have not been identified yet. The fraction of RL quasars with E(B − V)> 0.1 (including quasars from this work) is 2.3% for our quasars compared to 1.5% for the RQ literature sample. Comparing the two distributions with a KS-test results in a p-value of 0.5, indicating no significant deviation. However, the sample size of RL quasars is still relatively small with 44 quasars and biased by observational selection. Including also the fits using the starburst galaxy extinction curve of Calzetti et al. (2000) yields higher fractions of red quasars of 9.1% and 3.2% for the RL and RQ sample, respectively. However, there is still no significant deviation between the two distributions (p-value of 0.6), therefore supporting the same conclusion.

thumbnail Fig. 10.

Distribution of E(B − V) for the discovered quasars in this work and the literature sample determined by fitting the composite quasars spectrum of Selsing et al. (2016) together with the SMC extinction law (Pei 1992). The radio-loud quasars from the literature sample are stacked on top of quasars from this work. We do not find a significant deviation between the E(B − V) distribution of the RL and RQ quasar sample. For comparison we also plot the E(B − V) Gaussian distribution found by Glikman et al. (2022) of SDSS QSOs at 0.1 <  z ≲ 3 using a similar template fitting method.

The fraction of possible red quasars we find in our sample is similar to Kato et al. (2020), who report a fraction of 2% at z >  5.6 using the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs; Matsuoka et al. 2016) sample of 93 high-z quasars. However, due to selection biases, observational limits, and definitions, red quasars are understudied at high redshift and likely biased towards lower fractions.

At lower redshift Glikman et al. (2012) have reported a red quasar fraction of ∼20% between 0.13 <  z <  3.1 using 51 radio-selected quasars and Richards et al. (2003) found a fraction of 15% between 0.3 <  z <  2.2 using a sample of 4576 Sloan Digital Sky Survey (SDSS; York et al. 2000) quasars. Recently, Glikman et al. (2022) have investigated the extinction properties of a large sample of SDSS QSOs at 0.1 <  z ≲ 3 using a similar template fitting method as was used in this work. They find that the fraction of red quasars is strongly luminosity dependent, with red quasars making up to 40% of the sources at the highest luminosities, but decreasing to a few percent at lower luminosities. To compare their findings to our sample, we plotted the E(B − V) distribution found by Glikman et al. (2022) at z <  3 in Fig. 10, which shows the distribution is similar in shape. However, we do not find the same tail of red quasars out to E(B − V)∼1.0 as Glikman et al. (2022), which is as expected since most red quasars would disappear from high-z optical and NIR surveys due to the large amount of extinction.

To investigate whether the Legacy Surveys depth allows for discovering heavily reddened quasars, we determine the maximum E(B − V) for each quasar in our sample for which it would still be visible in the Legacy z-band (reaching a 5σ depth of mz ∼ 23.3). The fitting procedure is similar to the one described above, but here we use only the W1 and W2 measurements of each quasar and increase the E(B − V) value (using the SMC extinction law) on the best fitting quasar template until the Legacy z-band depth is reached. This analysis shows that 44% of the quasars would not have been detected if they had E(B − V)> 0.25 and 93% if they had E(B − V)> 0.4, which is consistent with not finding heavily reddened quasars in a tail out to E(B − V)∼1.0 and likely contributing to our lower fractions of red quasars. Therefore, we need significantly deeper optical data to detect the dust obscured quasars that are currently too faint for ground based telescopes. The Rubin Observatory Legacy Survey of Space and Time (LSST, LSST Science Collaboration 2009) will eventually reach a 5σ depth of mz ∼ 26.1, which would yield a > 80% completeness out to E(B − V)=0.4. Furthermore, to study red quasar fractions as function of radio luminosity in a solid statistical way at high redshift we need a larger sample size of radio bright quasars, which can be discovered by combining future large area radio surveys from e.g. LOFAR and SKA with deep optical and NIR data from e.g. LSST and the forthcoming Euclid mission (Euclid Collaboration 2019).

6. Summary and future prospects

We have started a campaign to search for radio-bright high-z quasars by combining an optical colour cut technique using the Legacy Surveys with LoTSS-DR2 radio detections. The LoTSS-DR2 survey is currently reaching sensitivities within the 100 μJy regime, which enables detection of ∼30% of the high-z quasar population at low-frequencies (Gloudemans et al. 2021a). We have performed spectroscopic follow-ups for 80 of our candidates, which lead to the confirmation of 20 new quasars (and the independent confirmation of four recent confirmed high-z quasars) between 4.9 ≤ z ≤ 6.6 and a doubling of the sample of known radio-loud quasars in this epoch. The radio detections from LoTSS-DR2 decreased the contamination of stellar sources in our sample and allowed for selection of these quasars in a broad redshift range and use of the limited number of optical bands (g,  r, and z bands) from the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys, which probe at least 1 magnitude deeper than other large surveys such as PS1. Our campaign demonstrates the potential for efficiently discovering new (faint) high-z quasar populations through next generation radio continuum surveys.

Our new quasar sample probes on average fainter rest-frame UV luminosities and brighter radio luminosities than the previously known quasar population. The optical and NIR colours and Lyα EW distribution are very similar to the known quasar population, which potentially suggests that there is no strong radio-loudness dichotomy between the radio-loud and radio-quiet quasar population. However, we note that we find a higher fraction (38%) of weak line quasars compared to the 14% reported by Bañados et al. (2016). Furthermore, we find little to no correlation between the Lyα and radio luminosity of our sample with a Pearson correlation coefficient of r = −0.12. However, this analysis is limited by the small sample size and observational selection biases.

Furthermore, we estimated the red quasar fraction of the RL versus the RQ high-z quasar population by fitting the quasar composite of Selsing et al. (2016) to the photometry and applying the SMC extinction law of Pei (1992) and starburst galaxy extinction curve of Calzetti et al. (2000). Our analysis suggests no significant deviation between the E(B − V) values of the RL and RQ high-z quasar sample.

We find that the success rate of our spectroscopic follow-up campaign can be improved significantly by decreasing the positional optical-radio separation limit of our candidates; by taking the optical-radio separation criterion as 1.2″ instead of 2.0″, the 34% success rate could have been increased to 42% without loss of new or known quasars. From a simple Monte Carlo simulation we derive the probability of a quasar having an optical-radio separation larger than 1″ to be ∼2%, assuming similar positional uncertainties. The separations should therefore be carefully considered in future work. However, note that the optimal separation limit is highly dependent on the resolution, sensitivity, calibration of the radio data, and target elevation (which causes changes in the synthesised beam).

To further study the AGN and host galaxy properties of these new high-z quasars, additional observations are necessary (see e.g. Khusanova et al. 2022). The SMBH masses and broad line regions for example can be constrained by measuring the C IV and Mg II NIR lines (see e.g. Fine et al. 2010; Onoue et al. 2019) and host galaxy kinematics by measuring the sub-millimetre [C II] emission line (e.g. Neeleman et al. 2021). Furthermore, the radio spectra need to be constrained to be able to study the radio emission mechanisms of the AGN. A statistical sample of radio spectra of high-z quasars is also essential to obtain better estimates of the number counts of radio sources suitable for future 21-cm experiments such as SKA (Ciardi et al. 2015) and to further optimise strategies for discovering suitable sources.

A novel approach to finding high-z radio sources will be undertaken by the WEAVE-LOFAR (WL) survey (Smith et al. 2016). The WEAVE instrument (Dalton et al. 2014) is a multi-fiber spectrograph on the William Herschel Telescope on La Palma and is expected to have first light in the summer of 2022. The WEAVE-LOFAR survey will obtain ∼1 million spectra of radio sources from the LoTSS surveys without optical colour cuts and will therefore probe the full range of potential quasar colours, as well as build large samples. This survey has the potential of enabling the detection of dust obscured high-z radio galaxies and quasars into the EoR and discovery of suitable sources for 21-cm spectroscopy.


3

The optical spectra and quasar properties are available in electronic form at the CDS.

Acknowledgments

We are grateful for the support of Kentaro Aoki and Ichi Tanaka during our Subaru/FOCAS observations (proposal ID S21B-003). K.J.D. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 892117 (HIZRAD). Y.H. is supported by JSPS KAKENHI Grant Number 21K13953. M.J.H. and D.J.B.S. acknowledge support from the UK Science and Technology Facilities Council (STFC) under grant ST/V000624/1. P.N.B. is grateful for support from the UK STFC via grant ST/V000594/1. The Low Resolution Spectrograph 2 (LRS2) was developed and funded by the University of Texas at Austin McDonald Observatory and Department of Astronomy, and by Pennsylvania State University. We thank the Leibniz-Institut fur Astrophysik Potsdam (AIP) and the Institut fur Astrophysik Goettingen (IAG) for their contributions to the construction of the integral field units. We thank the Board of the Hobby-Eberly Telescope for the allocation of Guaranteed Time for the LRS2 instrument, which was important in enabling this investigation. This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT) under project codes LC0 015, LC2 024, LC2 038, LC3 008, LC4 008, LC4 034 and LT10 01. LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin university (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy.

References

  1. Appenzeller, I., & Rupprecht, G. 1992, The Messenger, 67, 18 [NASA ADS] [Google Scholar]
  2. Bañados, E., Venemans, B. P., Morganson, E., et al. 2015, ApJ, 804, 118 [Google Scholar]
  3. Bañados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11 [Google Scholar]
  4. Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473 [Google Scholar]
  5. Bañados, E., Rauch, M., Decarli, R., et al. 2019, ApJ, 885, 59 [CrossRef] [Google Scholar]
  6. Bañados, E., Mazzucchelli, C., Momjian, E., et al. 2021, ApJ, 909, 80 [CrossRef] [Google Scholar]
  7. Baloković, M., Smolčić, V., Ivezić, Ž., et al. 2012, ApJ, 759, 30 [CrossRef] [Google Scholar]
  8. Beaklini, P. P. B., Quadros, A. V. C., de Avellar, M. G. B., Dantas, M. L. L., & Cançado, A. L. F. 2020, MNRAS, 497, 1463 [NASA ADS] [CrossRef] [Google Scholar]
  9. Becker, R. H., White, R. L., & Helfand, D. J. 1994, in ASP Conf. Ser., 61, 165 [NASA ADS] [Google Scholar]
  10. Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850 [NASA ADS] [CrossRef] [Google Scholar]
  11. Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32 [Google Scholar]
  12. Belladitta, S., Moretti, A., Caccianiga, A., et al. 2020, A&A, 635, L7 [EDP Sciences] [Google Scholar]
  13. Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2011, ApJ, 739, 24 [NASA ADS] [CrossRef] [Google Scholar]
  15. Brown, M. J. I., Moustakas, J., Smith, J. D. T., et al. 2014, ApJS, 212, 18 [Google Scholar]
  16. Brown, M. J. I., Duncan, K. J., Landt, H., et al. 2019, MNRAS, 489, 3351 [Google Scholar]
  17. Calistro Rivera, G., Alexander, D. M., Rosario, D. J., et al. 2021, A&A, 649, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Callingham, J. R., Vedantham, H. K., Shimwell, T. W., et al. 2021, Nat. Astron., 5, 1233 [NASA ADS] [CrossRef] [Google Scholar]
  19. Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
  20. Carilli, C. L., Röttgering, H. J. A., van Ojik, R., Miley, G. K., & van Breugel, W. J. M. 1997, ApJS, 109, 1 [NASA ADS] [CrossRef] [Google Scholar]
  21. Carilli, C. L., Gnedin, N. Y., & Owen, F. 2002, ApJ, 577, 22 [NASA ADS] [CrossRef] [Google Scholar]
  22. Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints [arXiv:1612.05560] [Google Scholar]
  23. Chonis, T. S., Hill, G. J., Lee, H., et al. 2016, in SPIE Conf. Ser., 9908, 99084C [NASA ADS] [Google Scholar]
  24. Ciardi, B., Inoue, S., Mack, K., Xu, Y., & Bernardi, G. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 6 [Google Scholar]
  25. Cirasuolo, M., Celotti, A., Magliocchetti, M., & Danese, L. 2003, MNRAS, 346, 447 [NASA ADS] [CrossRef] [Google Scholar]
  26. Coatman, L., Hewett, P. C., Banerji, M., & Richards, G. T. 2016, MNRAS, 461, 647 [NASA ADS] [CrossRef] [Google Scholar]
  27. Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693 [Google Scholar]
  28. Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog: II/328 [Google Scholar]
  29. Dalton, G., Trager, S., Abrams, D. C., et al. 2014, in SPIE Conf. Ser., 9147, 91470L [Google Scholar]
  30. Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168 [Google Scholar]
  31. Diamond-Stanic, A. M., Fan, X., Brandt, W. N., et al. 2009, ApJ, 699, 782 [NASA ADS] [CrossRef] [Google Scholar]
  32. Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142 [Google Scholar]
  33. Duncan, K. J. 2022, MNRAS, 512, 3662 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dunlop, J. S. 2013, Astrophys. Space Sci. Lib., 396, 223 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dye, S., Lawrence, A., Read, M. A., et al. 2018, MNRAS, 473, 5113 [Google Scholar]
  36. Eilers, A.-C., Hennawi, J. F., Decarli, R., et al. 2020, ApJ, 900, 37 [Google Scholar]
  37. Endsley, R., Stark, D. P., Fan, X., et al. 2022a, MNRAS, 512, 4248 [CrossRef] [Google Scholar]
  38. Endsley, R., Stark, D. P., Lyu, J., et al. 2022b, MNRAS, submitted [arXiv:2206.00018] [Google Scholar]
  39. Euclid Collaboration (Barnett, R., et al.) 2019, A&A, 631, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Fan, X., Narayanan, V. K., Lupton, R. H., et al. 2001, AJ, 122, 2833 [Google Scholar]
  41. Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415 [Google Scholar]
  42. Fawcett, V. A., Alexander, D. M., Rosario, D. J., et al. 2020, MNRAS, 494, 4802 [Google Scholar]
  43. Findlay, J. R., Sutherland, W. J., Venemans, B. P., et al. 2012, MNRAS, 419, 3354 [NASA ADS] [CrossRef] [Google Scholar]
  44. Fine, S., Croom, S. M., Bland-Hawthorn, J., et al. 2010, MNRAS, 409, 591 [NASA ADS] [CrossRef] [Google Scholar]
  45. Gallerani, S., Maiolino, R., Juarez, Y., et al. 2010, A&A, 523, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Glikman, E., Urrutia, T., Lacy, M., et al. 2012, ApJ, 757, 51 [Google Scholar]
  47. Glikman, E., Lacy, M., LaMassa, S., et al. 2022, ApJ, 934, 119 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gloudemans, A. J., Duncan, K. J., Röttgering, H. J. A., et al. 2021a, A&A, 656, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Gloudemans, A. J., Duncan, K. J., Kondapally, R., et al. 2021b, A&A, 648, A7 [Google Scholar]
  50. Greig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239 [Google Scholar]
  51. Gürkan, G., Hardcastle, M. J., Best, P. N., et al. 2019, A&A, 622, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Hewett, P. C., Warren, S. J., Leggett, S. K., & Hodgkin, S. T. 2006, MNRAS, 367, 454 [Google Scholar]
  53. Hill, G. J., Lee, H., MacQueen, P. J., et al. 2021, AJ, 162, 298 [NASA ADS] [CrossRef] [Google Scholar]
  54. Hook, I. M., Shaver, P. A., & McMahon, R. G. 1998, in ASP Conf. Ser., 146, 17 [NASA ADS] [Google Scholar]
  55. Hopkins, P. F., Strauss, M. A., Hall, P. B., et al. 2004, AJ, 128, 1112 [Google Scholar]
  56. Horne, K. 1986, PASP, 98, 609 [Google Scholar]
  57. Ichikawa, K., Yamashita, T., Toba, Y., et al. 2021, ApJ, 921, 51 [NASA ADS] [CrossRef] [Google Scholar]
  58. Ighina, L., Belladitta, S., Caccianiga, A., et al. 2021, A&A, 647, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27 [NASA ADS] [CrossRef] [Google Scholar]
  60. Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Ivezić, Ž., Menou, K., Knapp, G. R., et al. 2002, AJ, 124, 2364 [CrossRef] [Google Scholar]
  62. Izumi, T., Onoue, M., Matsuoka, Y., et al. 2019, PASJ, 71, 111 [CrossRef] [Google Scholar]
  63. Jarvis, M. J., Rawlings, S., Lacy, M., et al. 2001, MNRAS, 326, 1563 [NASA ADS] [CrossRef] [Google Scholar]
  64. Jiang, L., Fan, X., Vestergaard, M., et al. 2007, AJ, 134, 1150 [NASA ADS] [CrossRef] [Google Scholar]
  65. Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222 [Google Scholar]
  66. Kashikawa, N., Aoki, K., Asai, R., et al. 2002, PASJ, 54, 819 [NASA ADS] [Google Scholar]
  67. Kato, N., Matsuoka, Y., Onoue, M., et al. 2020, PASJ, 72, 84 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B., & Green, R. 1989, AJ, 98, 1195 [Google Scholar]
  69. Kelly, B. C., Vestergaard, M., Fan, X., et al. 2010, ApJ, 719, 1315 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kelson, D. D. 2003, PASP, 115, 688 [NASA ADS] [CrossRef] [Google Scholar]
  71. Khusanova, Y., Bañados, E., Mazzucchelli, C., et al. 2022, A&A, 664, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Knapp, G. R., Leggett, S. K., Fan, X., et al. 2004, AJ, 127, 3553 [Google Scholar]
  73. Liu, Y., Wang, R., Momjian, E., et al. 2021, ApJ, 908, 124 [Google Scholar]
  74. LSST Science Collaboration (Abell, P. A., et al.) 2009, arXiv e-prints [arXiv:0912.0201] [Google Scholar]
  75. Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8 [NASA ADS] [CrossRef] [Google Scholar]
  76. Macfarlane, C., Best, P. N., Sabater, J., et al. 2021, MNRAS, 506, 5888 [NASA ADS] [CrossRef] [Google Scholar]
  77. Mack, K. J., & Wyithe, J. S. B. 2012, MNRAS, 425, 2988 [NASA ADS] [CrossRef] [Google Scholar]
  78. Maiolino, R., Oliva, E., Ghinassi, F., et al. 2004, A&A, 420, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Marziani, P., del Olmo, A., Martínez-Carballo, M. A., et al. 2019, A&A, 627, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Matsuoka, Y., Onoue, M., Kashikawa, N., et al. 2016, ApJ, 828, 26 [Google Scholar]
  81. Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2018, ApJS, 237, 5 [Google Scholar]
  82. Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2019, ApJ, 883, 183 [Google Scholar]
  83. Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2022, ApJS, 259, 18 [NASA ADS] [CrossRef] [Google Scholar]
  84. McGreer, I. D., Becker, R. H., Helfand, D. J., & White, R. L. 2006, ApJ, 652, 157 [Google Scholar]
  85. McGreer, I. D., Jiang, L., Fan, X., et al. 2013, ApJ, 768, 105 [NASA ADS] [CrossRef] [Google Scholar]
  86. McGreer, I. D., Fan, X., Jiang, L., & Cai, Z. 2018, AJ, 155, 131 [NASA ADS] [CrossRef] [Google Scholar]
  87. McLure, R. J., & Jarvis, M. J. 2004, MNRAS, 353, L45 [NASA ADS] [CrossRef] [Google Scholar]
  88. McLure, R. J., Dunlop, J. S., Cirasuolo, M., et al. 2010, MNRAS, 403, 960 [Google Scholar]
  89. Miley, G., & De Breuck, C. 2008, A&ARv, 15, 67 [Google Scholar]
  90. Miley, G. K. 1968, Nature, 218, 933 [NASA ADS] [CrossRef] [Google Scholar]
  91. Miley, G. K. 1971, MNRAS, 152, 477 [NASA ADS] [CrossRef] [Google Scholar]
  92. Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1 [NASA ADS] [Google Scholar]
  93. Momjian, E., Carilli, C. L., Bañados, E., Walter, F., & Venemans, B. P. 2018, ApJ, 861, 86 [NASA ADS] [CrossRef] [Google Scholar]
  94. Morabito, L. K., Jackson, N. J., Mooney, S., et al. 2022, A&A, 658, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616 [Google Scholar]
  96. Mortlock, D. J., Patel, M., Warren, S. J., et al. 2012, MNRAS, 419, 390 [NASA ADS] [CrossRef] [Google Scholar]
  97. Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109 [NASA ADS] [CrossRef] [Google Scholar]
  98. Nanni, R., Hennawi, J. F., Wang, F., et al. 2022, MNRAS, 515, 3224 [NASA ADS] [CrossRef] [Google Scholar]
  99. Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141 [NASA ADS] [CrossRef] [Google Scholar]
  100. Nyland, K., Dong, D. Z., Patil, P., et al. 2020, ApJ, 905, 74 [NASA ADS] [CrossRef] [Google Scholar]
  101. Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713 [NASA ADS] [CrossRef] [Google Scholar]
  102. Oke, J. B., Cohen, J. G., Carr, M., et al. 1995, PASP, 107, 375 [Google Scholar]
  103. Onoue, M., Kashikawa, N., Matsuoka, Y., et al. 2019, ApJ, 880, 77 [Google Scholar]
  104. Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, A&A, 613, A51 [Google Scholar]
  105. Park, D., Woo, J.-H., Denney, K. D., & Shin, J. 2013, ApJ, 770, 87 [NASA ADS] [CrossRef] [Google Scholar]
  106. Pei, Y. C. 1992, ApJ, 395, 130 [Google Scholar]
  107. Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, ApJ, 793, 113 [NASA ADS] [CrossRef] [Google Scholar]
  108. Prochaska, J., Hennawi, J., Westfall, K., et al. 2020, J. Open Source Software, 5, 2308 [NASA ADS] [CrossRef] [Google Scholar]
  109. Ramsey, L. W., Engel, L., Rhoads, B., et al. 1998, Am. Astron. Soc. Meeting Abstracts, 193, 10.07 [NASA ADS] [Google Scholar]
  110. Richards, G. T., Hall, P. B., Vanden Berk, D. E., et al. 2003, AJ, 126, 1131 [Google Scholar]
  111. Rodríguez Hidalgo, P., Khatri, A. M., Hall, P. B., et al. 2020, ApJ, 896, 151 [CrossRef] [Google Scholar]
  112. Rojas-Ruiz, S., Bañados, E., Neeleman, M., et al. 2021, ApJ, 920, 150 [NASA ADS] [CrossRef] [Google Scholar]
  113. Rosario, D. J., Fawcett, V. A., Klindt, L., et al. 2020, MNRAS, 494, 3061 [Google Scholar]
  114. Ross, N. P., & Cross, N. J. G. 2020, MNRAS, 494, 789 [NASA ADS] [CrossRef] [Google Scholar]
  115. Ross, N. P., Hamann, F., Zakamska, N. L., et al. 2015, MNRAS, 453, 3932 [Google Scholar]
  116. Saxena, A., Röttgering, H. J. A., & Rigby, E. E. 2017, MNRAS, 469, 4083 [Google Scholar]
  117. Saxena, A., Röttgering, H. J. A., Duncan, K. J., et al. 2019, MNRAS, 489, 5053 [Google Scholar]
  118. Sbarrato, T., Ghisellini, G., Giovannini, G., & Giroletti, M. 2021, A&A, 655, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Schindler, J.-T., Farina, E. P., Bañados, E., et al. 2020, ApJ, 905, 51 [Google Scholar]
  120. Schroeder, J., Mesinger, A., & Haiman, Z. 2013, MNRAS, 428, 3058 [Google Scholar]
  121. Selsing, J., Fynbo, J. P. U., Christensen, L., & Krogager, J. K. 2016, A&A, 585, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Shao, Y., Wagg, J., Wang, R., et al. 2020, A&A, 641, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Shen, Y., Wu, J., Jiang, L., et al. 2019, ApJ, 873, 35 [Google Scholar]
  124. Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Smith, D. J. B., Best, P. N., Duncan, K. J., et al. 2016, in SF2A-2016: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, J. Richard, L. Cambrésy, et al. (Paris: SF2A), 271 [Google Scholar]
  126. Sobral, D., Darvish, B., Matthee, J., et al. 2016, The Hosts of the Early Ionized Bubbles: the Nature and Diversity of the Most Luminous Lyman-alpha Emitters at z 6-7, HST Proposal, 24, 14699 [Google Scholar]
  127. Soto, K. T., Lilly, S. J., Bacon, R., Richard, J., & Conseil, S. 2016, MNRAS, 458, 3210 [Google Scholar]
  128. Tody, D. 1986, in SPIE Conf. Ser., 627, 7333 [NASA ADS] [Google Scholar]
  129. Uchiyama, H., Yamashita, T., Toshikawa, J., et al. 2022, ApJ, 926, 76 [NASA ADS] [CrossRef] [Google Scholar]
  130. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  131. van Dokkum, P. G. 2001, PASP, 113, 1420 [Google Scholar]
  132. van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  133. van Ojik, R., Roettgering, H. J. A., Miley, G. K., & Hunstead, R. W. 1997, A&A, 317, 358 [NASA ADS] [Google Scholar]
  134. Vedantham, H. K., Callingham, J. R., Shimwell, T. W., et al. 2020, ApJ, 903, L33 [Google Scholar]
  135. Venemans, B. P., Findlay, J. R., Sutherland, W. J., et al. 2013, ApJ, 779, 24 [Google Scholar]
  136. Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37 [Google Scholar]
  137. Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689 [Google Scholar]
  138. Vestergaard, M., Fan, X., Tremonti, C. A., Osmer, P. S., & Richards, G. T. 2008, ApJ, 674, L1 [NASA ADS] [CrossRef] [Google Scholar]
  139. Wagenveld, J. D., Saxena, A., Duncan, K. J., Röttgering, H. J. A., & Zhang, M. 2022, A&A, 660, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  140. Wang, F., Wu, X.-B., Fan, X., et al. 2015, ApJ, 807, L9 [NASA ADS] [CrossRef] [Google Scholar]
  141. Wang, F., Wu, X.-B., Fan, X., et al. 2016, ApJ, 819, 24 [NASA ADS] [CrossRef] [Google Scholar]
  142. Wang, F., Fan, X., Yang, J., et al. 2017, ApJ, 839, 27 [Google Scholar]
  143. Wang, F., Yang, J., Fan, X., et al. 2019, ApJ, 884, 30 [Google Scholar]
  144. Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699 [Google Scholar]
  145. Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44 [Google Scholar]
  146. Wenzl, L., Schindler, J.-T., Fan, X., et al. 2021, AJ, 162, 72 [NASA ADS] [CrossRef] [Google Scholar]
  147. White, R. L., Helfand, D. J., Becker, R. H., Glikman, E., & de Vries, W. 2007, ApJ, 654, 99 [Google Scholar]
  148. Williams, W. L., Hardcastle, M. J., Best, P. N., et al. 2019, A&A, 622, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  149. Willott, C. J., Delorme, P., Omont, A., et al. 2007, AJ, 134, 2435 [NASA ADS] [CrossRef] [Google Scholar]
  150. Willott, C. J., Delorme, P., Reylé, C., et al. 2010, AJ, 139, 906 [Google Scholar]
  151. Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, ARA&A, 43, 861 [NASA ADS] [CrossRef] [Google Scholar]
  152. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
  153. Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512 [Google Scholar]
  154. Yamashita, T., Nagao, T., Ikeda, H., et al. 2020, AJ, 160, 60 [NASA ADS] [CrossRef] [Google Scholar]
  155. Yang, J., Wang, F., Fan, X., et al. 2017, AJ, 153, 184 [NASA ADS] [CrossRef] [Google Scholar]
  156. Yang, J., Wang, F., Fan, X., et al. 2019, ApJ, 871, 199 [NASA ADS] [CrossRef] [Google Scholar]
  157. York, D. G., Adelman, J., Anderson, J. E. J., et al. 2000, AJ, 120, 1579 [Google Scholar]
  158. Zamfir, S., Sulentic, J. W., & Marziani, P. 2008, MNRAS, 387, 856 [NASA ADS] [CrossRef] [Google Scholar]
  159. Zuo, W., Wu, X.-B., Fan, X., et al. 2020, ApJ, 896, 40 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Contaminants

All spectroscopically observed candidates that have been identified as not being quasars are listed in Table A.1 for future reference. Based on initial visual inspection, the majority (∼60%) of the contaminants we have detected in our spectroscopic follow-up observations are M-type brown dwarf (BD) candidates. If these are indeed proven to be BDs with associated radio emission, this is an interesting discovery by itself, since the first direct low-frequency radio detections of BDs have only recently been reported by Vedantham et al. (2020) and Callingham et al. (2021). In those studies, their radio emission can partially be explained by plasma emission from the active chromospheres of the stars. However, the origin of radio emission for the in-active stars is potentially driven by star-planet interactions. The investigation of the BD nature of these sources is beyond the scope of this paper and will be discussed in separate work (Gloudemans et al. in prep.).

The other ∼40% of the contaminants have slowly rising featureless spectra and are expected to be either low-z SF galaxies with possible AGN or L- or T-type dwarfs, which are difficult to classify in the optical wavelength regime. In the case of low-z SF galaxies the red-dusty spectrum could be confused with a high-z quasar SED due to a strong Balmer break (around z∼1-2) or dust extincted continuum (see e.g. McLure et al. 2010; Donley et al. 2012; Dunlop 2013). L- and T-type dwarfs also exhibit a slowly rising spectrum in the optical, however to classify these follow-up IR observations are needed to identify for example the TiO and VO lines and methane absorption (see e.g. Knapp et al. 2004).

Table A.1.

Spectroscopically observed candidates that are confirmed not to be high-z quasars.

Appendix B: Cutouts of newly discovered quasars

thumbnail Fig. B.1.

20″ × 20″ Legacy DR8 (Dey et al. 2019), AllWISE (Cutri et al. 2021), and LOFAR (LoTSS-DR2; Shimwell et al. 2022) cutouts for the newly confirmed quasars.

thumbnail Fig. B.1.

continued.

All Tables

Table 1.

Overview of selection steps taken to create our candidate target list for spectroscopic follow-up.

Table 2.

Details of spectroscopic observations conducted of the newly discovered high-z quasars (sorted on RA).

Table 3.

Observed and measured physical properties of the discovered high-z quasars.

Table A.1.

Spectroscopically observed candidates that are confirmed not to be high-z quasars.

All Figures

thumbnail Fig. 1.

Optical and NIR colour selection criteria (given by the dashed lines) of our spectroscopically observed quasar candidates compared to literature quasars (5.0 ≤ z ≤ 7.5) and observed stellar M-, L-, and T-dwarfs.

In the text
thumbnail Fig. 2.

Absolute flux scale correction factors between the quasar spectra and PS1 i-band and Legacy z-band photometry. A correction factor of > 1 means the spectrum continuum is lower than the photometric measurement. These corrections do not affect the redshift and Lyα EW estimates, but do decrease or increase the measured Lyα emission line flux.

In the text
thumbnail Fig. 3.

Observed-frame optical spectra obtained with Subaru/FOCAS, HET/LRS2, and Keck/LRIS of the newly confirmed quasars sorted on spectroscopic redshift. The total flux measurements at 144 MHz are also indicated on the left hand side of each spectrum. The detected emission lines are indicated by the grey dashed lines, which are Lyα+N V emission lines unless indicated otherwise.

In the text
thumbnail Fig. 4.

Association likelihood ratios versus optical-radio separations of our observed sample. The blue histograms show the distribution for the quasar sample from this work and literature combined. We demonstrate that a similar distribution (p-value = 0.5) of separations can be created by a simple Monte Carlo simulation using the positional offsets from both the Legacy and LoTSS-DR2 catalogue. From this distribution we derive the chance of a quasar having an optical-radio separation larger than 1″ to be ∼2% when assuming no intrinsic difference in the position of the radio and optical emission.

In the text
thumbnail Fig. 5.

Radio spectral index between 144 MHz (LoTSS-DR2) and 1.4 GHz (VLA FIRST) for the discovered quasars in the FIRST and VLASS survey area. Upper limits (indicated by downward-pointing triangles) are given for the quasars non-detected in FIRST. The dashed grey line shows the stacked spectral index obtained by Gloudemans et al. (2021a) using a sample of 115 known quasars at z >  5. ILTJ2201+2338 is outside of the FIRST footprint, but is detected in VLASS at 2–4 GHz, which enables a spectral index measurement. ILTJ1523+2935 is only marginally detected in FIRST (S/N ∼ 2), but is more securely detected in VLASS (S/N = 3.4).

In the text
thumbnail Fig. 6.

Radio luminosities at 144 MHz of newly discovered quasars in this work compared to the literature sample of known high-z quasars detected in LoTSS-DR2. Other known z >  6 radio-loud quasars PSO J172+18 (Bañados et al. 2021) and VIK J2318−3113 (Ighina et al. 2021) remain non-detected and not covered by LOFAR, respectively. PSO J0309+27 (Belladitta et al. 2020) and COS-87259 (Endsley et al. 2022b) are detected by LOFAR outside of the LoTSS-DR2 area.

In the text
thumbnail Fig. 7.

Left: radio-loudness (R = f5 GHz/f4400 Å) as a function of spectroscopic redshift for the discovered quasars in this work and the literature quasar sample. Our quasars span a broad range in radio-loudness. Traditionally, quasars are defined as radio-loud if R >  10, which is satisfied by 21/24 of our newly discovered quasars. Right: UV rest-frame absolute magnitude at 1450 Å versus radio luminosity at 144 MHz. The new quasar sample probes fainter UV magnitudes and higher average radio luminosities than previously known quasar samples.

In the text
thumbnail Fig. 8.

Near-infrared colours of the discovered quasars in this work and contaminants compared to the literature sample of quasars and stellar dwarfs5 (for which the PS-1 z-band magnitude is used in the right panel). Despite a looser colour selection, the new quasars exhibit colours consistent with the known quasar population.

In the text
thumbnail Fig. 9.

Derived Lyα+N V emission line properties of the discovered quasars. Left: normalised Lyα+N V EW distribution of the discovered quasars in this work compared to the distribution of PS1 discovered quasars at z >  5.6 by Bañados et al. (2016) and SDSS discovered quasars at z >  3 by Diamond-Stanic et al. (2009). Our EW distribution is consistent with the optically selected quasars sample of Bañados et al. (2016). Right: Lyα+N V luminosity versus the radio luminosity at 144 MHz for the quasars in this work, showing little to no correlation (see Sect. 5.3).

In the text
thumbnail Fig. 10.

Distribution of E(B − V) for the discovered quasars in this work and the literature sample determined by fitting the composite quasars spectrum of Selsing et al. (2016) together with the SMC extinction law (Pei 1992). The radio-loud quasars from the literature sample are stacked on top of quasars from this work. We do not find a significant deviation between the E(B − V) distribution of the RL and RQ quasar sample. For comparison we also plot the E(B − V) Gaussian distribution found by Glikman et al. (2022) of SDSS QSOs at 0.1 <  z ≲ 3 using a similar template fitting method.

In the text
thumbnail Fig. B.1.

20″ × 20″ Legacy DR8 (Dey et al. 2019), AllWISE (Cutri et al. 2021), and LOFAR (LoTSS-DR2; Shimwell et al. 2022) cutouts for the newly confirmed quasars.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.