XQz5: A New Ultraluminous z$\sim$5 Quasar Legacy Sample

TripleSpec4.1 (Wilson et al., 2004) is a fixed-assembly cross-dispersed long-slit near-infrared imaging spectrograph with slit dimensions of 1.1” by 28” mounted on SOAR, a 4.1m aperture telescope situated in Cerro Pachón, Chile. The instrument covers a simultaneous wavelength range from 0.94 to 2.47 $\mu$m, with spectral resolution R$\sim$3500.

We observed 26 targets in the classical observing mode over 9 nights between 2021-08-22 and 2022-05-15, under the programme IDs 2021B-0036 and 2022A-389756. Using the standard ABBA observing strategy, the total integration time for each target was calculated to achieve a signal-to-noise ratio (SNR) over 20 around the Mg ii emission-line. For targets in the H-band magnitude range $15.8<\textit{H}<17.0$, the integration times ranged from 2880s to 12480s between the brightest and faintest targets. We carried out the data reduction and extraction of the infrared spectra using an adaptation of the Spextool IDL package (v4.1) (Cushing et al., 2004), specifically designed for TripleSpec4.1. The post-extraction flux calibration and telluric correction were performed with the xtellcorr IDL package (Vacca et al., 2003) using standard stars observed immediately before and after each target, which are well matched in airmass. The resulting mean $\rm{SNR}_{\rm{2600}}\simeq 15.8$ per 40 km s^-1 pixel measured from the median SNR in the 100Å window between 2550–2650Å in the continuum region near the Mg ii line.

VLT’s X-shooter instrument (Vernet et al., 2011) is a multi-wavelength, medium resolution spectrograph with three spectroscopic arms: UVB (300–560nm), VIS (550–1020nm), and NIR (1020–2480nm). We proposed observations for 40 target quasars under the programme IDs 108.22H9.001, 109.23D1.001, and 109.23D1.002. Observations were completed for 32 targets in “service mode” between 2021-10-20 and 2022-08-12, 13 of which were also observed with TripleSpec4.1. We also include one target, J2335–5901, which was observed under programme ID 0104.A-0410(A). Under service mode, we defined Observation Blocks (OBs) containing the instrument setup and requested weather conditions.

For the instrument setup, the adopted slit widths are 1.6” in the UVB, 1.5” in VIS, and 1.2” in NIR. Each slit width is the widest slit available for science observations to maximise the signal while providing a nominal resolving power of 3200, 5000, and 4300, respectively. The integration time for each target is calculated to achieve $\rm{SNR}>20$ around the Mg ii line in the NIR arm after resampling the extracted spectra into 50 km s^-1 velocity bins. Of the executed OBs, 71% were within observing condition specifications with an additional 5% almost within specifications, resulting in repeat observations of two targets (J091656–251146 and J121921–360933). Where the executed OBs did not meet the desired observing specifications, we manually check the data quality for each individual frame and find that despite the lower quality, all of the data could still be used for further analysis. The total VLT time investment for these new targets is 22 hours. We add two additional targets from the 0104.A-0410(A) programme which invested 0.8 hours of on-target exposure time for J215728.21–360215.1 at $z=4.692$ and 1.6 hours for J233505–590103 at $z=4.53$, resulting in high quality spectra (SNR${}_{2600}\sim$40/pixel).

We further supplement our sample with additional X-shooter data from the ESO Science Archive Facility of 36 targets, originating from the programmes 084.A-0574(A), 084.A-0780(B), 087.A-0125(A), 094.A-0793(A), 098.A-0111(A). We selected these additional targets from the Million Quasars (Milliquas) v7.7 Catalogue²²2https://quasars.org/milliquas.htm (Flesch, 2021), most of which are bright quasars in the redshift range $5.0<z<5.3$ concentrated at more northern declinations as seen in Figure 1. The targets from 084.A-0574(A), 087.A-0125(A), and 094.A-0793(A) were observed with slit dimensions such that the effective resolution is R = 8900, 5600 for the VIS and NIR arms, respectively. Targets from 098.A-0111(A) and 0100.A-0243(A) were observed with higher resolution NIR arm, where R = 8900, 8100, and BR 1202–0725 from 084.A-0780(B) was observed with R = 11400, 8100.

Due to the high sky-background level, extraction of the NIR spectra can be non-trivial. Although ESO provides a pipeline for X-shooter data reduction (Freudling et al., 2013), we find that the post-extraction spectra often exhibited significant sky-subtraction residuals. Consequently, we opted to use the semi-automated reduction procedure implemented in the PypeIt package (Prochaska et al., 2020), which contains code infrastructure for X-shooter data reduction. The initial flux calibration is derived from a standard star typically observed within a day of the observation. However, the standard is not often matched to the airmass of our targets, so we later rescale the spectra to the observed photometric magnitudes as described in Section 2.4. Furthermore, the telluric correction algorithm in PypeIt utilises a large telluric grid to find the best-fitting atmospheric model in place of deriving the correction from a telluric standard.

After the extraction of the 1D spectra of all targets observed with X-shooter, the resulting mean $\rm{SNR}_{1450}\simeq\rm{SNR}_{2600}\simeq 16.8$ per pixel measured from the median SNR between 1400–1500Å and 2550–2650Å.

The Wide Field Spectrograph (WiFeS) instrument (Dopita et al., 2007; Dopita et al., 2010) installed on the fully automated Australian National University 2.3m telescope is a double-beam integral field spectrograph with a wavelength coverage of 330–930 nm. For this sample, WiFeS observations are conducted on the 13 targets without optical coverage, using R = 3000 gratings. The data reduction is performed with the standard Python-based pipeline, PyWiFeS (Childress et al., 2014) and extraction of the 1D spectra from the IFU data cubes was performed with PySpecExtract³³3https://github.com/nishamrutha/PySpecExtract. The resulting mean $\rm{SNR}_{1450}=6.1$ per 50 km s^-1 pixel. We re-calibrate the flux normalisation using optical photometry with the same method applied for X-shooter data as described in the folllowing Section 2.4.

In Table 1, we list properties of all targets and extracted spectra in our sample. The columns from left to right are the target name, coordinates, redshift derived from the peak of the Mg ii emission-line (the details of which is presented in Section 3), programme ID, Milky Way extinction colour excess (see Section 2.4), magnitudes in the i, z, H-bands and signal-to-noise ratios in 100Å windows centered at 1450 and 2600Å. The optical i and z-band photometry are obtained from SMSS DR4 (Onken et al., prep), if available, or Pan-STARRS DR1 (Chambers et al., 2016). The near-infrared H-band photometry is from one of various surveys: VISTA Hemisphere Survey (VHS; McMahon et al., 2013) DR5/6, VISTA Kilo-degree Infrared Galaxy Survey (VIKING; Edge et al., 2013) DR2, UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al., 2007) DR9, or Two Micron All-Sky Survey (2MASS; Skrutskie et al., 2006). In some cases, the H-band magnitude is interpolated (using the relations from O22) from the VHS J and K-bands or the SMSS/Pan-STARRS z and ALLWISE/CatWISE2020 W1 bands (Wright et al., 2010; Cutri et al., 2013; Marocco et al., 2021). The survey sources of the optical and near-infrared photometry are indicated in the supplementary table. The full table, extended with associated uncertainties, is available online as supplementary material.

Our combined sample consists of 83 ultraluminous quasars based on new observations from SOAR/TripleSpec4.1 and VLT/X-shooter alongside additional X-shooter data from the ESO Science Archive Facility. These quasars are found in the redshift range $4.5<z<5.3$ with H-band magnitudes $15.44<\textit{H}<18.83$ and their sky distribution in equatorial coordinates is presented in Figure 1. Here, we describe the post-processing steps taken to clean the spectra and prepare them for continuum and emission-line modelling.

1. 1.

   We transform the spectrum to the rest-frame by estimating an initial redshift from the peak of the observed Mg ii data. We then apply a per-pixel SNR floor of 5 and mask all data below.

2. 2.

   We apply a sigma-clipping mask across box-widths of 40 pixels, using a 3$\sigma$ threshold to remove narrow absorption features and noise above 3$\sigma$. Narrow absorption features, if left in the spectrum, can affect the modelling of the quasar continuum or intrinsic line profiles of broad emission lines.

3. 3.

   The spectra are resampled from their original post-extraction resolution (13.5/40 km s^-1 for X-shooter NIR/TripleSpec4.1 and <11.6/50 km s^-1 for X-shooter VIS/WiFeS) into a common wavelength grid from rest-frame 1050–3600Å using a flux-conserving algorithm, SPECTRES (Carnall, 2017), with a velocity resolution of 50 km s^-1 for each bin. Rebinning and consolidating the flux in 50 km s^-1 bins increases the SNR of the spectral data, most notably for the X-shooter data, and standardises the output from different spectroscopic instruments for continuum and emission-line modelling.

4. 4.

   As the initial X-shooter flux calibration is based on a standard that is not reliably observed in the same night or under similar conditions, we independently recalibrate the X-shooter spectra to photometry. We integrate the observed spectroscopic data over the bandpass transmission curves obtained from the SVO Filter Profile Service⁴⁴4http://svo2.cab.inta-csic.es/theory/fps/ (Rodrigo & Solano, 2020) and scale the flux to match the observed magnitudes. For the X-shooter VIS arm, we use the inverse variance to weight the flux calibration between the i-band and z-band. We then scale the X-shooter NIR arm according to the calibration of the VIS arm, enforcing continuity between arms. Targets observed with TripleSpec4.1 are already calibrated to an airmass-matched standard observed immediately before or after the quasar, so we take its original calibration to be most reliable. If there is TripleSpec4.1 near-infrared spectroscopic data available for an object, then we use the TripleSpec4.1 spectrum to scale all other spectroscopic data. If we calibrated TripleSpec4.1 spectra to H-band photometry taken at a variety of different epochs, we find that the median absolute correction for the TripleSpec4.1 sample to be 15% of its original flux calibration with an additional 15% RMS deviation. In some cases, prior to the rescaling and splicing of the optical and near-infrared spectra, a discontinuity can be observed from the combined effects of chromatic slit losses and quasar variability. We note that the 15% median absolute calibration discrepancy would result in $<0.1$ dex uncertainty in the estimated continuum luminosity.

5. 5.

   We correct for Milky Way extinction by adopting $R_{\rm{v}}=3.1$ with the Fitzpatrick (1999) extinction curve and colour excess from the Schlegel, Finkbeiner & Davis (SFD; Schlegel et al., 1998) extinction map. We then apply a suggested 14% recalibration to the colour excess, $E(B-V)=0.86\times E(B-V)_{\rm{SFD}}$, informed by fits to the blue tip of the stellar locus (Schlafly et al., 2010). We assume no further dust extinction from the quasar host galaxy, following the common practice for similar studies of quasar line properties, including Shen et al. (2011) – which defines our mass calibration (see Section 3.2.1) – and others (e.g., Rakshit et al., 2020; Wu & Shen, 2022; Lai et al., 2023b). Furthermore, luminous quasars are intrinsically less likely to exhibit notable dust extinction, with the majority of quasars consistent with null colour excess (e.g., Krawczyk et al., 2015). Due to the ultraluminous nature of the quasars in this sample, we do not account for the starlight contribution in the spectra, which is insignificant compared to the AGN.

After post-processing of our sample, we find the median $\rm{SNR}=27.2$ per 50 km $s^{-1}$ velocity bin measured from the median SNR between 2700–2750Å and 2850–2900Å. In Figure 2, we show our sample, XQz5, compared to other high-quality X-shooter programmes, XQ-100 and Enlarged XQR-30 Plus (E-XQR-30+), which are concentrated on lower and higher redshifts, respectively. The E-XQR-30+ sample is composed of the core 30 quasars from XQR-30 and supplemented with high-redshift literature quasars with NIR spectra as described in D’Odorico et al. (2023); in addition for this paper, we have extended E-XQR-30+ with 112 quasars listed in the Fan et al. (2022) review paper. We note that the gravitationally lensed object QSOJ0439$+$1634 has been removed. We also show quasars from T11, which is a lower luminosity sample at a similar redshift range, and we plot three notable quasars: J0529–4351 (Onken et al., 2023), J2157–3602 (Wolf et al., 2018b; Onken et al., 2020), and J0100$+$2802 (Wu et al., 2015), the three most luminous quasars in the respective redshift ranges of each sample, of which we include J2157–3602 in our sample. Other luminous quasar samples such as ELQS (Schindler et al., 2017) and QUBRICS (Cristiani et al., 2023) also reach up to $z\sim 5$, but they have not been observed with a comparable broad spectral coverage as the X-shooter samples, which include observed-frame near-infrared.

We provide all of the observed-frame telluric-corrected spectra in the supplementary material. The spectra are provided in the native resolution and do not include the post-processing steps described in this section, except for the flux re-calibration. Raw data files are available in their respective data archive facilities, the ESO Science Archive Facility⁵⁵5http://archive.eso.org/cms.html and NOIRLab Astro Data Archive⁶⁶6https://astroarchive.noirlab.edu/.

Before we model the Mg ii line, we first subtract a model of the quasar continuum. Our quasar continuum model is composed of two components: a power law and Fe ii template, which is altogether called the pseudo-continuum. The pseudo-continuum is then uniquely defined by five free parameters and is fit simultaneously to the spectroscopic data contained within selected rest-frame windows: 1973–1983Å, 2060–2340Å, 2600–2740Å, and 2840–3100Å, chosen to minimise contribution from strong emission-lines while also covering wavelength regions in close proximity to the Mg ii line with a wide diversity of expected Fe ii flux.

The Mg ii feature is sensitive to the Fe ii continuum emission, which extends underneath the emission-line, and it has previously been established that the choice of Fe ii template can bias the result. Notably, adopting the Vestergaard & Wilkes (2001) template, which is often used to calibrate the Mg ii-based single-epoch virial black hole mass estimate (e.g., Vestergaard & Osmer, 2009; Shen et al., 2011), results in larger line widths than with other models of the Fe ii emission (Schindler et al., 2020; Lai et al., 2023b). Thus, we rely on four empirical and semi-empirical Fe ii templates: Vestergaard & Wilkes (2001, VW01), Tsuzuki et al. (2006, T06), Bruhweiler & Verner (2008, BV08), and Mejía-Restrepo et al. (2016, M16) to minimise bias induced by any single template. The templates are renormalised and broadened with a Gaussian kernel to match the features of each observed spectrum.

A third potential component to the pseudo-continuum, the Balmer continuum, is sometimes introduced in similar studies of the Mg ii emission-line in AGN (e.g., Wang et al., 2009). However, the flux of the Balmer continuum decreases continuously on the blue side of the Balmer edge located at 3645Å, and is relatively weak near the Mg ii line (Vestergaard & Wilkes, 2001; Dietrich et al., 2002). We also note that the Balmer continuum flux can be degenerate with the other pseudo-continuum components. Thus, despite an improved global continuum fit, the isolation of the Mg ii emission-line profile is often independent of a reasonably constrained Balmer continuum model (e.g., Shen & Liu, 2012; Lai et al., 2023b). Therefore, we do not include the Balmer continuum in our pseudo-continuum model, although we note that the exclusion of the Balmer continuum can change the flux of the power-law and Fe ii components.

Even though we do not observe evidence of a strong Balmer flux contribution across the Balmer edge in the luminous quasar spectra of XQ-100 (Lai et al., 2023b), we estimate the systematic uncertainty related to the exclusion of the Balmer contribution to the continuum model. A common model of the Balmer continuum, adopted in other studies (e.g. Dietrich et al., 2003; Kurk et al., 2007; De Rosa et al., 2011; Shin et al., 2019), normalises the flux at the Balmer edge to 30% the flux of the power-law component at 3675Å, where the Fe ii emission is absent. We assume a lower normalisation of 10% that is also commonly used in the literature (Calderone et al., 2017; Rakshit et al., 2020). Following these studies, we adopt the same electron temperature ($T_{\rm{e}}$ = 15,000 K) and optical depth ($\tau_{\lambda}$ = 1), finding that the Balmer flux contribution at 3000Å would be between 7–8% of the power-law flux, depending on the power-law slope, where we have used the 16^th and 84^th percentiles of all measured slopes from our sample, $\gamma=-1.36^{+0.40}_{-0.33}$. Thus, the power-law flux at 3000Å could be overestimated by up to 10% depending on how the Balmer flux is proportioned between the power-law and Fe ii components. We later quantify on how this could affect our results.

Broad emission-line profiles can be found with a wide range of properties and morphological complexities, which render emission-line models described by a single analytical function unsuitable, particularly with a Doppler broadened doublet emission-line like Mg ii. As such, quasar spectral modelling studies typically utilise a multiple Gaussian approach to fit each individual emission feature (e.g., Greene & Ho, 2005; Shen et al., 2011; Rakshit et al., 2020; Wu & Shen, 2022). Accordingly, we model the Mg ii broad emission-line with three Gaussian components and use an additional component to fit the narrow line, if present. The maximum width of the narrow component is set at 1000 km s^-1.

After subtracting the pseudo-continuum model, we fit the Mg ii line model within rest-frame 2730–2870Å. We then determine the Mg ii emission-line properties from the total line profile composed of the three broad components. The final mean line properties and measurement uncertainties are determined from the emergent emission-line models assuming the four different Fe ii templates. We further add in quadrature the additional uncertainty measured by fitting 50 synthetic spectra, which are generated by randomly redistributing the flux in each fixed-velocity bin according to its Gaussian uncertainties.

In Table 2, we describe each of the line properties measured from the quasar continuum and Mg ii broad emission-line model. For each target, we update the measured redshift based on the wavelength of the peak flux in the reconstructed Mg ii line profile. From the combined broad emission-line profile, we measure the Mg ii FWHM and line dispersion, Sigma, which is the second moment of the line profile. The blueshift of the line is measured from the flux-bisecting wavelength, as

where $\lambda_{\rm{ref}}$ is the vacuum wavelength of the Mg ii doublet which is set to 2799.117Å, and $\lambda_{\rm{med}}$ is the wavelength bisecting the total continuum-subtracted emission-line flux. We also measure the line rest-frame equivalent width (EW) and integrated luminosity. The monochromatic luminosity at 3000Å is measured from the power-law continuum model. For the 14 targets with more than one spectrum, the final line properties are measured from the weighted average of multiple observations. In our sample, we find one target, J051509–431854, with clear Mg ii absorption, which is a common characteristic of low-ionisation broad absorption line (LoBAL) quasars. This is consistent with the expectation that quasar with BAL outflows are found in less than 20% of quasars (e.g. Gibson et al., 2009) and LoBAL outflows are even rarer, consisting of 10–15% of BAL quasars. However, there is some evidence that the BAL fraction increases with redshift (Bischetti et al., 2023). A more detailed analysis of the BAL quasar fraction will be in an separate forthcoming study.

In Figure 3, we show two targets which represent a range of line model susceptibility to the Fe ii template. The target on the left panel (J013539–212628) is highly sensitive to the Fe ii model, with a Mg ii FWHM that varies from 2910–5400 km s^-1. In contrast, the target on the right panel (J040733–281031) is minimally affected as a result of its weak Fe ii emission. Its Mg ii FWHM range is from 2420–2540 km s^-1.

To measure the bolometric luminosity of a quasar exactly, it is necessary to observe the quasar at all wavelengths from all possible viewing angles and have an understanding of the source component decomposition for the observed luminosity at different wavelengths. Because that is not feasible, it is common among studies of quasars to apply a first-order correction to the observed monochromatic luminosity in the form, $\rm L_{\rm{bol}}=k_{\lambda}\rm L(\lambda)$, where $k_{\lambda}$ is the bolometric correction factor and $\rm L(\lambda)$ represents the monochromatic luminosity, $\lambda\rm L_{\lambda}$ at rest wavelength $\lambda$.

The bolometric correction factor, $k_{\lambda}$, is calibrated for a mean quasar spectral energy distribution (SED), integrated across a particular wavelength range. In this study, we use $k_{3000}=5.15$, a first-order correction from Richards et al. (2006) that is calibrated to the mean SED of 259 SDSS Type 1 quasars, integrated from 100 $\mu$m to 10 keV and assumed to be emitting isotropically. However, as the Type 1 quasars used to construct the mean SED are likely to be biased towards face-on orientations, the resulting luminosity will be overestimated. Thus, we adopt the model presented in Lai et al. (2023a) and assume that the maximum opening angle of the obscuring torus for our sample, which defines the boundary between a Type 1 and Type 2 quasar classification, is 65^∘. We then use a thin disc model (Li et al., 2005) to estimate the correction between the observed and bolometric luminosities, assuming the quasars are randomly distributed in observed orientations and black hole spins. Using this method, the mean correction factor for anisotropy is calculated to be a 25% luminosity suppression to the observed luminosity, which is consistent with the correction in Runnoe et al. (2012), derived from similar arguments. Thus, for this study, we adopt $\rm L_{\rm{bol}}=0.75\times 5.15\times\rm L(3000)$ and we use this bolometric correction consistently for the comparison with literature samples, unless otherwise stated. Due to the spectral diversity of quasars, the systematic uncertainty from deriving the assuming a mean SED can be as high as 50% (Richards et al., 2006), and so a systematic error of $\Delta\log{\rm L_{\rm{bol}}}=0.3$ is assumed. We further note that the bolometric correction is expected to be mass or luminosity-dependent (e.g. Runnoe et al., 2012; Netzer, 2019). For an ultraluminous 24 billion solar mass black hole, the bolometric correction factor at 3000Å could be as low as $k_{3000}=1.62$ (Lai et al., 2023a), which is less than half of the bolometric luminosity derived with this method. An additional 0.05 dex adjustment may also result from a 10% overestimation of continuum flux by not explicitly modelling the Balmer component as cautioned in Section 3.1 and calibration errors may contribute up to 0.1 dex, but these corrections are not as significant as the assumed systematic error. We emphasise the possibility that the brightest and most massive black holes in this sample could have overestimated bolometric luminosities, but we adopt a consistent correction for all targets to facilitate comparison.

Using $k_{3000}=5.15$ with a 0.75 anisotropy correction factor, the estimated bolometric luminosity range occupied by the XQz5 sample is then $\log{(\rm{L_{bol}}/\rm{erg\,s^{-1}})}=46.5-48.34$ with Eddington ratios spanning nearly an order of magnitude from 0.19–1.8.

In this section, we place our black hole mass, luminosity, and Eddington ratio estimates in context with the literature. First, we compare the measurements obtained from the SDSS quasars in our sample with the results from the C iv-based virial estimator from SDSS DR16Q (Wu & Shen, 2022). We also compare our estimates of bolometric luminosity with those from O22⁹⁹9O22 does not include black hole mass estimates, which derived its $L_{\rm{3000}}$ luminosity estimates by matching each target’s photometric information to a composite spectrum of bright quasars from $z=1-2$ (Selsing et al., 2016). These comparisons are presented in Figure 4, from which we observe that the bolometric luminosity estimates between XQz5 and the literature are in tight agreement, despite the fact that the DR16Q estimates come from a different wavelength regime as $\rm L_{\rm{bol}}=3.81\,\rm L_{\rm{1350}}$. The literature luminosity of one outlier (not shown), identified as SDSS J001714.67–100055.4, is 0.7 dex lower than our measured luminosity and we found its continuum model in DR16Q to have been affected by a broad absorption feature. As for the measured black hole masses, the C iv-based virial mass estimator applied to SDSS data significantly overestimates the $M_{\rm{BH}}$ for our sample on average. When the 1690–1705Å continuum fitting window falls near the red end of the wavelength coverage affected by enhanced sky background and molecular oxygen skylines, there is no anchor point for the continuum model on the red side of C iv. For this overlapping sample, we find that the cubic polynomial continuum model component adopted in (Wu & Shen, 2022) significantly underestimates the continuum underneath the C iv line when fitted only to the blue side of C iv, which results in inflated C iv FWHMs and overestimated black hole masses.

Before comparing the XQz5 to other quasar samples, we first caution that systematic biases can arise from heterogeneous subsample selection criteria, as even well-defined subsamples can be biased if the sample completeness is unknown (Shen & Kelly, 2012; Kelly & Shen, 2013; Wu et al., 2022). Furthermore, we have defined the bolometric luminosity estimate using a static correction to the UV continuum luminosity. Implicitly, this assumes a consistent overall mean SED (Richards et al., 2006), which could again be a source of bias between heterogeneously selected samples. However, as XQz5 is a well-defined sample with high completeness at the bright end, we make these comparisons in order to examine relative distributions of other quasar samples compared to our sample. We include our 0.75 anisotropy correction factor consistently in our comparisons to literature estimates of bolometric luminosity.

We first construct a comparison to the T11 sample in Figure 5, using their Mg ii measurements that we have rescaled to the Shen et al. (2011) calibration. Both XQz5 and the flux-limited T11 sample constructed from SDSS DR6 are biased towards high Eddington ratios with similar Eddington ratio distributions. However, the XQz5 targets are measured to be both more luminous and have higher mass black holes than T11. Both samples occupy the high mass and bright end of the subsample of SDSS DR16Q with well-constrained continuum luminosities and Mg ii-based mass measurements between z = 0.35 – 2.3. In total, we find 11 targets in XQz5 with mildly super-Eddington luminosity.

Using SDSS DR16Q quasars with Mg ii measurements as the parent sample, we construct a luminosity-matched sample by collecting two quasars with the closest matching luminosity to each quasar in XQz5, resulting in a subsample composed of 164 quasars in the redshift range $0.6<z<2.0$. The median $\rm L_{\rm{bol}}$ difference of a selected SDSS DR16Q quasar and its associated XQz5 quasar is $\sim 0.19$ dex. Despite this, due to the scarcity of high luminosity SDSS quasars, the exact luminosity distribution of XQz5 could not be reproduced. In the top panels of Figure 6, we compare the luminosity, black hole mass, and Eddington ratio distribution between XQz5, the luminosity-matched DR16Q subsample, and the WISSH sample (Bischetti et al., 2017; Vietri et al., 2018), which is a high-luminosity SDSS and WISE-selected set of quasars across a wide range of redshifts ($1.8<z<4.7$) designed to reflect the most luminous known AGN in the universe. The black hole parameters of the WISSH sample are typically derived from the C iv emission line and $\rm L(1350)$, which are inter-calibrated to the results from the Mg ii and H$\beta$ line but contribute additional scatter due to their weak correlation with other virial mass estimators (e.g., Baskin & Laor, 2005; Shen & Liu, 2012). Nevertheless, we find in Figure 6, that both of the luminosity-matched lower redshift samples host black holes with higher mass than in XQz5. Statistical analysis of the XQz5 and lower-redshift mass distributions with the Kolmogorov-Smirnov (KS) test yields p-values less than 1%, indicating a statistically significant difference in the parent population. This also implies that ultraluminous quasars observed at lower redshift are generally less active with lower Eddington ratios than their counterparts at higher redshift, which is in concordance with recent literature (e.g., Yang et al., 2021; Farina et al., 2022). The mean Eddington ratio difference between XQz5 and the two lower-redshift samples is $\Delta\lambda_{\rm{Edd}}\sim 0.3$

When comparing XQz5 with high-quality brightness-selected X-shooter quasar samples (XQ-100 and E-XQR-30+) in the bottom panels of Figure 6, we find all samples to be distributed towards high Eddington ratios, suggesting that all three samples are composed of rapidly accreting massive quasars. The mean(median) Eddington ratio of XQz5 is $\sim 0.59(0.58)$, which is the highest mean(median) value of the three samples. The E-XQR-30+ sample has a long low-luminosity tail which consists entirely of literature quasars compiled in the Fan et al. (2022) database. Unlike XQz5 and XQ-100, these quasars are not preferentially selected from the bright quasar population. By excluding the low luminosity tail of E-XQR-30+ such that the distribution better reflects the brightest known quasars at $5.8<z<6.6$, we find that the three Eddington ratio distributions are not statistically different as measured by the KS test.

In this comparison, we find XQz5 to be the brightest X-shooter quasar sample by proportion and by number of quasars brighter than $\log\left(L_{\rm{bol}}\rm{[ergs^{-1}]}\right)=47.3$. As both XQz5 and XQR-30 are selected from the brightest known quasars at the time each sample was constructed, they are highly complete at the bright end of the quasar luminosity function. In particular, the completeness of the parent sample suggests that it is unlikely that a sample of brighter targets would exist within the redshift range of XQz5.

From XQR-30 at $z\sim 6.1$ to the later cosmic epoch of XQz5 at $z\sim 5.0$, we observe a significant increase in luminosity with decreasing redshift and a milder boost in black hole mass over a period of $\sim 260$ Myr. This reflects an underlying population of quasars that are on a rising branch stage of their cosmic evolution, as their SMBHs have been feeding continuously from an earlier seed. A forthcoming paper will model and examine the exponential mass growth of quasars at this redshift range.