# HYPERION: Broad-band X-ray-to-near-infrared emission of quasars in the first billion years of the Universe

I. Saccheo, A. Bongiorno, E. Piconcelli, L. Zappacosta, M. Bischetti, V. d'Odorico, C. Done, M. J. Temple, V. Testa, A. Tortosa, et al.

► To cite this version:

I. Saccheo, A. Bongiorno, E. Piconcelli, L. Zappacosta, M. Bischetti, et al., HYPERION: Broad-band X-ray-to-near-infrared emission of quasars in the first billion years of the Universe. Astronomy & Astrophysics - A&A, 2025, 693, 10.1051/0004-6361/202450876. insu-05101725

HAL Id: insu-05101725

https://insu.hal.science/insu-05101725v1

Submitted on 6 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# HYPERION: Broad-band X-ray-to-near-infrared emission of quasars in the first billion years of the Universe

I. Saccheo<sup>1,2,*</sup> <u>D</u>, A. Bongiorno<sup>2</sup>, E. Piconcelli<sup>2</sup> <u>D</u>, L. Zappacosta<sup>2</sup> <u>D</u>, M. Bischetti<sup>3,4</sup>, V. D'Odorico<sup>3,5,6</sup> <u>D</u>, C. Done<sup>7</sup>, M. J. Temple<sup>7</sup> <u>D</u>, V. Testa<sup>2</sup> <u>D</u>, A. Tortosa<sup>2</sup> <u>D</u>, M. Brusa<sup>8,9</sup> <u>D</u>, S. Carniani<sup>5</sup>, F. Civano<sup>10</sup>, A. Comastri<sup>9</sup> <u>D</u>, S. Cristiani<sup>3,6,11</sup> <u>D</u>, D. De Cicco<sup>12,13,14</sup>, M. Elvis<sup>15</sup>, X. Fan<sup>16</sup>, C. Feruglio<sup>3,6</sup> <u>D</u>, F. Fiore<sup>3,6</sup>, S. Gallerani<sup>5</sup>, E. Giallongo<sup>2</sup>, R. Gilli<sup>9</sup>, A. Grazian<sup>17</sup> <u>D</u>, M. Guainazzi<sup>18</sup>, F. Haardt<sup>19,20,21</sup>, R. Maiolino<sup>22,23,24</sup>, N. Menci<sup>2</sup> <u>D</u>, G. Miniutti<sup>25</sup> <u>D</u>, F. Nicastro<sup>2</sup>, M. Paolillo<sup>12,14</sup> <u>D</u>, S. Puccetti<sup>26</sup>, F. Salvestrini<sup>3</sup> <u>D</u>, R. Schneider<sup>2,27,28,29</sup>, F. Tombesi<sup>30,31,2</sup> <u>D</u>, R. Tripodi<sup>32,3,4,6</sup> <u>D</u>, R. Valiante<sup>2,28</sup>, L. Vallini<sup>9</sup> <u>D</u>, E. Vanzella<sup>9</sup>, G. Vietri<sup>33</sup> <u>D</u>, C. Vignali<sup>8,9</sup>, F. Vito<sup>9</sup> <u>D</u>, M. Volonteri<sup>34,35</sup>, and F. La Franca<sup>1,2</sup> <u>D</u>

*(Affiliations can be found after the references)*

Received 25 May 2024 / Accepted 18 October 2024

## ABSTRACT

**Aims.** We aim to characterize the X-ray-to-optical/near-infrared(NIR) broad-band emission of luminous quasars (QSOs) in the first gigayear (Gyr) of cosmic evolution in order to decipher whether or not they exhibit differences compared to the lower-z QSO population. Our goal is also to provide a reliable and uniform catalog of derivable properties for these objects (from fitting their spectral energy distribution), such as bolometric and monochromatic luminosities, Eddington ratios, dust extinction, and the strength of the hot dust emission.

**Methods.** We gathered all available photometry—from XMM-Newton proprietary data in X-rays to rest-frame NIR wavelengths—for the 18 QSOs in the HYPERION samples ($6.0 \le z \le 7.5$). For sources lacking uniform NIR coverage, we conducted observations in the *J*, *H*, and *K* bands. To increase the statistical robustness of our analysis across the UV-to-NIR region, we add 36 additional sources to our sample from the E-XQR-30 sample with $5.7 \le z \le 6.6$. We characterized the X-ray/UV emission of each QSO using average SEDs from luminous Type 1 sources and calculated bolometric and monochromatic luminosities. Finally, we constructed a mean SED extending from the X-rays to the NIR bands.

**Results.** We find that the UV-optical emission of these QSOs can be modeled with templates of $z \sim 2$ luminous QSOs. We observe that the bolometric luminosities derived while adopting some bolometric corrections at 3000 Å (BC<sub>3000 Å</sub>) largely used in the literature are slightly overestimated, by 0.13 dex, as they also include reprocessed IR emission. We estimate a revised value of BC<sub>3000 Å</sub> = 3.3, which can be used to derive $L_{bol}$ in $z \ge 6$ QSOs. We provide a subsample of 11 QSOs with rest-frame NIR photometry; these show a broad range of hot dust emission strength, with two sources exhibiting low levels of emission. Despite potential observational biases arising from nonuniform photometric coverage and selection biases, we produce an X-ray-to-NIR mean SED for QSOs at $z \ge 6$ that is a good match to templates of lower-redshift, luminous QSOs up to the UV-optical range, with a slightly enhanced contribution from hot dust in the NIR.

**Key words.** galaxies: high-redshift – quasars: general – quasars: supermassive black holes

## 1. Introduction

Thanks to the use of wide-field optical surveys such as the Sloan Digital Sky Survey (SDSS), the Canada-France High-z Quasar Survey (CFHQS), and the Panoramic Survey Telescope & Rapid Response System (Pan-STARRS), about 300 quasars (QSOs) have been discovered at $z \sim 6-7.6$ to date (e.g., Jiang et al. 2016; Bañados et al. 2016; Mazzucchelli et al. 2017; Fan et al. 2023, and references therein). These QSOs are powered by accretion onto supermassive black holes (SMBHs) with masses of $M_{SMBH} > 10^8-10^9 M_{\odot}$ and exhibit bolometric luminosities of $L_{bol} > 10^{46} \text{ erg s}^{-1}$ (e.g., Mazzucchelli et al. 2017, 2023; Shen et al. 2019).

The presence of fully grown SMBHs already at the epoch of reionization (EoR; $z \ge 5.5$) presents a challenge to theoretical models designed to help understand their formation in the early stages of cosmic evolution (e.g., Volonteri 2010; Inayoshi et al. 2020; Volonteri et al. 2021). To grow SMBHs at such high masses in a short amount of time, models usually employ two distinct pathways, namely (a) super-Eddington accretion (e.g.,

Madau et al. 2014; Volonteri et al. 2015; Pezzulli et al. 2016) or (b) massive BH seeds with $M_{BH,seed} \sim 10^3-10^4 M_{\odot}$ (e.g., Volonteri 2010; Valiante et al. 2016). Scenario (a) allows SMBH formation from lower-mass seeds, with $M_{BH,seed} \sim 10^2 M_{\odot}$, due to PopIII star remnants. Scenario (b) instead allows an Eddington-limited accretion regime.

The upcoming Euclid and Legacy Survey of Space and Time (LSST) surveys, which, according to the latest luminosity functions (e.g., Shen et al. 2020), are expected to uncover thousands of new QSOs at $z \ge 6$, will provide a complete census of these primeval sources and may help us to constrain their evolutionary history.

So far, one finding concerning these first QSOs is that their UV-to-near-infrared (NIR) broad-band emission resembles that observed in lower-redshift counterparts (e.g., Fan et al. 2004; Iwamuro et al. 2004; Reed et al. 2019). Moreover, by analyzing the UV-to-optical spectra of 50 QSOs at $z \ge 5.7$, Shen et al. (2019) found no significant major differences from those observed in the SDSS sample (Vanden Berk et al. 2001). Nevertheless, some distinctions are reported, including faster disk winds traced by CIV blueshifts (Meyer et al. 2019; Schindler et al. 2020; Yang et al. 2021) and a higher occurrence

<footnote>* Corresponding author; ivano.saccheo@uniroma3.it</footnote>
of broad absorption lines (BALs; Bischetti et al. 2022) in high-z QSOs compared to lower-redshift objects of similar luminosity. At higher energies, the monochromatic UV-to-X-ray-luminosity ratio shows no significant redshift evolution (Vito et al. 2019; Wang et al. 2021) and follows a trend akin to that seen in lower-z QSOs (e.g., Vignali et al. 2003; Just et al. 2007; Lusso et al. 2010); although Zappacosta et al. (2023) report evidence for a mild deviation, which needs to be further investigated. In contrast, in the X-ray band, Zappacosta et al. (2023) recently presented results of an X-ray spectral analysis conducted on a sample of $z \ge 6$ QSOs, namely the HYPerluminous quasars at the Epoch of Reionization (HYPERION), showing that these highly accreting sources have a significantly steeper photon index $\Gamma$ than similarly luminous $z < 6$ QSOs (see also Vito et al. 2018). This finding points toward a redshift evolution of $\Gamma$ that may be linked to a different accretion process taking place in the region closest to the SMBH.

To date, the $L_{bol}$ of high-redshift QSOs has mainly been computed through optical bolometric corrections (e.g., Reed et al. 2019; Yang et al. 2021; Mazzucchelli et al. 2023), which have been shown to be roughly constant over a wide range of luminosities and redshifts (Duras et al. 2020, but see also Trakhtenbrot & Netzer 2012). However, bolometric corrections have an intrinsic dispersion of ~0.25 dex (Duras et al. 2020) that can be accounted for by the scatter in the slopes of the power laws describing the optical continua of individual QSOs, and therefore these corrections do not necessarily provide the most reliable results for individual sources. Moreover, most measurements of the optical luminosities do not account for dust extinction, which could lead to their underestimation.

In this work, we systematically investigated the X-ray-to-NIR continuum emission of a sample of the 18 HYPERION QSOs, for which their exists exhaustive multiwavelength coverage between these bands. We also included 36 additional sources that show comparable redshift and luminosity distributions and for which there exists rest-frame UV-to-NIR photometry in order to increase the statistical significance of our analysis. Our main goal is to provide accurate measurements of accretion-disk-related properties (i.e., bolometric and monochromatic luminosities, and Eddington ratios) via fitting their spectral energy distributions (SEDs; e.g., Elvis et al. 1994), and to use these to construct a homogeneous reference catalog. In addition, we derive a mean SED for these sources ranging from the X-ray to the NIR and compare it to that measured for low-redshift QSOs (e.g., Vanden Berk et al. 2001; Krawczyk et al. 2013). Throughout this paper, we use a standard flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km/s Mpc}^{-1}$ and $\Omega_0 = 0.27$. Unless otherwise stated, uncertainties are reported at 68% confidence level and upper limits in the photometry are reported at the $3\sigma$ level.

## 2. The sample

To perform a reliable X-ray-to-NIR description of QSOs at $z > 6$, we considered the sources in the HYPERION sample (Zappacosta et al. 2023), which were recently targeted by a ~700 h XMM-Newton Heritage Programme (PI L. Zappacosta) that provided the best-quality X-ray data for these objects to date. The HYPERION sample includes the sources whose SMBHs experienced the fastest mass-growth rates during their formation. In particular, assuming an exponential, continuous growth at the Eddington rate limit, the HYPERION QSOs have been selected as the luminous ($L_{bol} > 10^{47}$ erg/s) $z > 6$ sources whose SMBH required (for its assembly) a seed black hole (BH) mass (hereafter $M_{s,Edd}$) of >1000 $M_{\odot}$ (formed at $z = 20$,

<table><caption>Table 1. The HYPERION sample.</caption><thead><tr><th>Name</th><th>Redshift</th><th>Log(<i>M</i><sub>SMBH</sub>)<br>[<i>M</i><sub>☉</sub>]</th><th>Log(<i>L</i><sub>BC</sub>)<br>[erg/s]</th><th><i>M</i><sub>s,Edd</sub><br>[<i>M</i><sub>☉</sub>]</th><th>Log(<i>F</i><sub>0.3-7 keV</sub>)<br>[erg/s cm<sup>-2</sup>]</th></tr></thead><tbody><tr><td>J1342+0928</td><td>7.541</td><td>8.90</td><td>47.19</td><td>19 120</td><td>-15.52 ± 0.06</td></tr><tr><td>J1007+2115</td><td>7.514</td><td>9.18</td><td>47.30</td><td>32 460</td><td>-15.16 ± 0.07</td></tr><tr><td>J1120+0641</td><td>7.087</td><td>9.41</td><td>47.30</td><td>18 230</td><td>-15.16 ± 0.07</td></tr><tr><td>J0038-1527</td><td>7.021</td><td>9.14</td><td>47.36</td><td>7983</td><td>&lt;-16.09</td></tr><tr><td>J0252-0503</td><td>7.0</td><td>9.15</td><td>47.12</td><td>7679</td><td>-14.12 ± 0.05</td></tr><tr><td>J0020-3653</td><td>6.834</td><td>9.24</td><td>47.16</td><td>5753</td><td>-15.33 ± 0.08</td></tr><tr><td>J0411-0907</td><td>6.824</td><td>8.80</td><td>47.31</td><td>2019</td><td>-14.49 ± 0.09</td></tr><tr><td>J0244-5008</td><td>6.724</td><td>9.08</td><td>47.19</td><td>2814</td><td>-15.33 ± 0.03</td></tr><tr><td>J231-20.8</td><td>6.587</td><td>9.50</td><td>47.31</td><td>4708</td><td>-15.32 ± 0.05</td></tr><tr><td>J036+03.0</td><td>6.533</td><td>9.49</td><td>47.33</td><td>3776</td><td>-15.63 ± 0.09</td></tr><tr><td>J0224-4711</td><td>6.526</td><td>9.36</td><td>47.53</td><td>2730</td><td>-14.23 ± 0.07</td></tr><tr><td>J011+09</td><td>6.444</td><td>9.15</td><td>47.12</td><td>1279</td><td>-15.44 ± 0.20</td></tr><tr><td>J1148+5251</td><td>6.422</td><td>9.74</td><td>47.57</td><td>4627</td><td>-14.76 ± 0.35</td></tr><tr><td>J083+11.8</td><td>6.346</td><td>9.32</td><td>47.16</td><td>1324</td><td>-14.53 ± 0.06</td></tr><tr><td>J0100+2802</td><td>6.3</td><td>10.04</td><td>48.24</td><td>5799</td><td>-14.25 ± 0.05</td></tr><tr><td>J025-33</td><td>6.294</td><td>9.57</td><td>47.39</td><td>1392</td><td>-15.62 ± 0.06</td></tr><tr><td>J0050+3445</td><td>6.245</td><td>9.68</td><td>47.29</td><td>2072</td><td>-14.50 ± 0.05</td></tr><tr><td>J029-36</td><td>6.02</td><td>9.82</td><td>47.39</td><td>1220</td><td>-15.44 ± 0.05</td></tr></tbody></table>

Notes. Redshift and $M_{SMBH}$ reported have been all estimated from the MgII line, see Zappacosta et al. (2023) for further references. X-ray fluxes from Zappacosta et al. (2023) and Tortosa et al. (2024).

Valiante et al. 2016). $M_{s,Edd}$ is a proxy for the growth-rate history experienced by SMBHs. More specifically, the HYPERION selection requires that

$$M_{s,Edd} = M_{SMBH} \times \exp(-t/t_s) \ge 1000 M_{\odot}, \quad (1)$$

where $t$ is the time elapsed since the seed formation and $t_s = 0.45\epsilon(1 - \epsilon)^{-1}\lambda_{Edd}^{-1}f_{\text{duty}}^{-1}$ is the e-folding time, where $\lambda_{Edd} = \frac{L_{\text{bol}}}{L_{\text{Edd}}} (=1)$, $f_{\text{duty}} (=1)$, and $\epsilon$ are the Eddington ratio, the duty cycle (i.e., the fraction of time during which the QSO is active), and the radiative efficiency (i.e., the fraction of accreting mass radiated away), respectively.

Table 1 presents the HYPERION sources together with some physical parameters, namely $z$, $M_{SMBH}$, $L_{BC}$ (i.e., the bolometric luminosity derived assuming a 3000 Å bolometric correction), $M_{s,Edd}$, and measured X-ray fluxes (see Zappacosta et al. 2023; Tortosa et al. 2024).

HYPERION is designed to provide the highest quality, and therefore most reliable, determination of the X-ray nuclear properties of QSOs at the EoR so far. As mentioned, the results from the analysis of first-year observations reported in Zappacosta et al. (2023) indicate an average photon index of $\Gamma \approx 2.4 \pm 0.1$, which is significantly steeper than the $\Gamma \sim 1.8-2$ typically observed in AGN at lower redshifts (e.g., Vignali et al. 2005; Piconcelli et al. 2005; Dadina 2008; Zappacosta et al. 2018, 2020). Moreover, the far-infrared and submillimeter region of the SED of these objects, as well as their dust and gas properties, are being investigated through ALMA and NOEMA data in a series of papers (Tripodi et al. 2023, 2024; Feruglio et al. 2023).

To increase the statistical significance of our analysis in the UV-to-NIR bands, we included a complementary sample from the Ultimate XSHOOTER legacy survey of quasars at $z \sim 5.8-6.6$ (D'Odorico et al. 2023, hereafter E-XQR-30), which consists of 42 bright sources ($i_{mag} < 20$) at $z \ge 5.7$. Notably, six QSOs, namely J231-20.8, J0224-4711, J029-36, J036+03.0, J025-33, and J0100+2802, belong to both HYPERION and E-XQR-30, and thus the considered E-XQR-30 subsample consists of 36 sources.
## Multiwavelength photometric data

In addition to the good-quality X-ray fluxes measured for all HYPERION QSOs by Zappacosta et al. (2023) and Tortosa et al. (2024), to construct the rest-frame X-ray-to-NIR SED of the HYPERION QSOs, we gathered all NIR and mid-infrared (MIR) photometric data available in the literature, specifically from the *z*, *Y*, *J*, *H*, and *K* filters plus the four WISE (3.4, 4.6, 12, and 22 µm, called W1, W2, W3, and W4, respectively) bands. Filters bluer than the *z* filter correspond to rest-frame wavelengths of shorter than the $Ly\alpha$ for all sources and therefore were not considered.

Only 2 sources, J1148+5251 and J0100+2802, have complete photometric coverage from *z* to W4, while for 12 QSOs we find photometry in all bands from *z* to W2; 3 sources lack data in the *H* and *K* filters, and for the remaining 4 QSOs only the *H* band is missing. For roughly half of the sources in the sample, multiple data are available in the same band for at least one filter. When possible, we gave preference to data obtained from targeted observations or data already published over data from surveys. In terms of survey catalogs, we made use of the Vista Hemisphere Survey (VHS; McMahon et al. 2021), UKIDSS-LAS (Lawrence et al. 2012), the Dark Energy Survey (DES; Abbott et al. 2021, we used the 3 arcsec aperture magnitudes), and the Pan-STARRS survey (Kaiser et al. 2010, PSF magnitudes reported). Two sources (J0411-0907 and J036+03.0) have their *Y* band covered by multiple catalogs, with their values in agreement within 0.1 mag; in this case, we decided to use the DES values.

For the WISE data, we relied on the unWISE Catalog by Schlafly et al. (2019) for W1 and W2 filters, and on the AllWISE catalog (Cutri et al. 2013) for W3 and W4. Given the wide PSF in these bands (~7.3 and 12 arcsec respectively), and the possibility of contamination from nearby projected companions, we did not use them in the fitting routine.

J1120+06 and J1148+52 have *Spitzer* IRAC and MIPS observations in 3.6 and 4.5 µm filters from Jiang et al. (2006) and 3.6, 4.5, 5.8, 8.0, and 24 µm bands from Barnett et al. (2015), respectively. For these two QSOs, as well as for the six E-XQR-30 sources with available *Spitzer* data, we used the 3.6 and 4.5 µm bands instead of W1 and W2. Moreover, for J1120+0641 we also included the 1 µm luminosity derived by Bosman et al. (2024) from its JWST/MIRI spectrum; that is, $\lambda L_{1\mu m} = 1.37 \pm 0.08 \times 10^{46}$ erg/s.

We complement the multiwavelength photometric data described above with additional proprietary NIR observations (see Appendix A) in order to achieve homogeneous and complete coverage of the emission from these sources. A detailed description of the NIR observations can be found in Appendix A, while the AB magnitudes used and their references of the HYPERION sample are available on Zenodo<fnref n="1" />. The multiwavelength data used for the E-XQR-30 sample are also hosted on Zenodo. For more details on their properties, we refer to Bischetti et al. (2022); D'Odorico et al. (2023, *J*, *H*, and *K* bands), Ross & Cross (2020, *z*, *Y* and W1 to W4 bands), and Leipski et al. (2014, *Spitzer* IRAC 3.6, 4.5, 5.8 and *Spitzer* MIPS 8.0 and 24 µm).

# 3. SED fitting

## 3.1. UV-to-NIR SED fitting

To analyze the SEDs for our QSOs, we performed SED fitting on the UV-to-NIR data using empirically derived pure AGN tem-

plates. Indeed, given the high-luminosity regime probed by the sources in our sample, we only accounted for the QSO emission, assuming a negligible contribution to the photometric points from the host galaxy (e.g., Shen et al. 2011).

Among the different templates in the literature, we used the two mean SEDs derived from samples that best match the luminosity distribution of our sources. These templates are the mean SED computed by Krawczyk et al. (2013) for their high-luminosity subsample ($\text{Log}(L_{2500\AA}) \ge 45.85$ erg/s, $0.5 \le z \le 4.8$), hereafter referred to as lum-K13, and the mean SED computed using the WISSH hyperluminous QSOs (Saccheo et al. 2023), hereafter WISSH-S23, which span a luminosity range that is even closer to that of the analyzed sample ($\text{log}(L_{bol}) > 47$ erg/s, $1.8 \le z \le 4.8$).

To account for intrinsic variations in the QSO SEDs that are not captured by the mean SED, we quantified the typical scatter of observed QSO SEDs relative to their average template. This step is needed because, even when restricting to very luminous sources, QSO SEDs show significant variations in their shapes (e.g., Richards et al. 2006; Temple et al. 2021a). Since the amount of scatter depends on the wavelength, this corresponds to assigning greater weights to the points where the SEDs show fewer variations. To calculate this scatter, we normalized the luminosity points of the QSOs originally used to compute the lum-K13 and WISSH-S23 templates using their $L_{bol}$. These normalized points were then binned into equally spaced bins on a logarithmic scale. The scatter in each bin was computed using the median absolute deviation (MAD), which is $\sigma_{SED} = \text{MAD}(\bar{y} - f)$, where *f* is the median value of the template in that bin normalized by the $L_{bol}$ of the SED, and $\bar{y}$ is the set of luminosity points within the bin.

For the lum-K13 template, we find that the scatter, $\sigma_{SED}/f$, increases as a function of wavelength, ranging from 0.06 to approximately 0.2 for $\lambda$ between 1216 Å and 1 µm, and reaching up to 0.3 in the IR at $\lambda \approx 5$ µm. In contrast, the WISSH-S23 template exhibits a more constant trend, with $\sigma_{SED}/f \approx 0.15$ between $\lambda = 1216$ Å and 3 µm, although with considerable fluctuations likely due to the small sample size.

With reference to photometric points, we accounted for the contribution of emission lines to the observed magnitudes. Emission from the BLR can contribute up to 30% of the observed luminosity (e.g., Miller et al. 2023). Therefore, when an emission line falls within a filter's transmission curve, it can lead to a significant overestimation of the primary emission of a QSO. To estimate the contribution of emission lines, we employed an approach analogous to that used by Krawczyk et al. (2013) and Saccheo et al. (2023), which involves the determination of the difference between the magnitudes obtained by convolving the composite spectrum from Vanden Berk et al. (2001) with the filters and those obtained using only the underlying continuum. The resulting magnitude difference, $\Delta m$, which varies depending on redshift and the filter considered, is then subtracted from the observed magnitudes (i.e., resulting in larger corrected magnitudes). To verify the consistency of this procedure, we repeated the analysis using the E-XQR-30 composite spectrum<fnref n="1" />. As this spectrum was derived from QSOs in our sample, it is more representative of the spectral features of the sources analyzed. However, we found only minimal differences between the two approaches ($\Delta m_{XQR_{30}} - \Delta m_{VB+01} \le 0.005$). Consequently, we chose to use the composite spectrum from Vanden Berk et al. (2001) as it extends beyond $\lambda = 3500$ Å, allowing us to estimate the contributions from the $H\alpha$ and $H\beta$ lines as well.

<footnote n="1">¹ https://github.com/XQR-30/Composites</footnote>
SED fitting on the corrected photometric points was performed for each of the two templates using the Python package EMCEE (Foreman-Mackey et al. 2019) to maximize the posterior arising from the likelihood,

$$ \mathcal{L}(K, E[B-V]) = - \sum_i \left( \frac{(y_i - K f_i 10^{-0.4A_i})^2}{\sigma_{i,tot}^2} + \log(2\pi\sigma_{i,tot}^2) \right), \quad (2) $$

where $f_i$ denotes the flux obtained by convolving the SED template with the i-th filter transmission, $A_\lambda$ is the dust-reddening law by Prevot et al. (1984), which accounts for a possible contribution from dust extinction (see also Bongiorno et al. 2012), and $\sigma_{i,tot}^2 = \sigma_{i,L}^2 + \sigma_{i,SED}^2$ is the sum of the squared uncertainties on the luminosity points and the scatter of the template SED at that wavelength, which also depends on the normalizing constant K. Data points falling at $\lambda < 1216$ Å were not taken into account in the fit. For both K and $E[B-V]$, we assumed constant, non-negative priors and derived their best-fit values as the mean of the posterior distribution, with uncertainties taken from the 16th and 84th percentiles. We note that this SED-fitting procedure assumes that each QSO has an intrinsic continuum that is the same as that of the employed template and any deviation from it is interpreted as the effect of dust extinction. As it is possible for the QSO to be intrinsically redder than the average template (see e.g., Richards et al. 2003; Krawczyk et al. 2015), we advise caution particularly with very low measured $E[B-V]$ values.

## 3.2. Results

We find that, for 73% of our sources, the lum-K13 template achieves better results –in terms of $\chi_v^2$– in describing the broad band SED of our QSOs, with a median $\chi_v^2 = 1.3$, compared to $\chi_v^2 = 1.8$ obtained with WISSH-S23. For this reason, throughout the paper, all values reported (i.e., $L_{bol}$, $E[B-V]$) are derived from the fitting obtained with the lum-K13 template. Examining the SEDs in detail (e.g., J083+11.8, shown in Fig. 1), we note that the WISSH-S23 template yields poorer results due to a flatter UV continuum for $\lambda < 2000$ Å, a feature not found in most of our EoR QSOs.

The HYPERION QSO SEDs with the final best fits for both lum-K13 and WISSH-S23 templates are reported in Fig. 2, while the figures showing the modeling with lum-K13 for the E-XQR-30 sample are available on Zenodo. Upon visually inspecting the fitted SEDs, we observe that the data points with $\lambda < 1$ µm are generally well fitted, while those at longer wavelengths exhibit larger discrepancies from the model (see also Sect. 7). This outcome is expected: as previously discussed, QSOs tend to show greater scatter around their mean value in the NIR region.

Two notable exceptions with particularly poor fits in the UV-optical region –although their $\chi_v^2$ are not the worst– are PSOJ023–02 and SDSSJ0836+0054. PSOJ023–02 has a very flat SED, which is not well modeled by reddening and may be intrinsic in nature, as it does not display the typical curling of the SED caused by dust extinction (i.e., the reddening is stronger at shorter wavelengths; see Hopkins et al. 2004). In contrast, the fit for SDSSJ0836+0054 is quite accurate for most points, except for the *Spitzer* S$_{4.5}$ data point, which is significantly more luminous than the modeled SED. While we cannot definitively explain this discrepancy, a plausible reason could be that the S$_{4.5}$ flux is overestimated. This is supported by the fact that, although

<figure><img src="image_2.png" /><figcaption>Fig. 1. SED fitting results for J083+11.8 at $z = 6.346$ using lum-K13 (blue) and WISSH-S23 (red) templates. The thickness of the lines denotes the $\pm 1\sigma$ interval. Gray points denote photometry below the $Ly\alpha$. The cyan shaded area and the dashed red lines denote the $\pm 1\sigma_{SED}$ region associated to the best fit for the lum-K13 and WISSH samples, respectively; see Sect. 3.1.</figcaption></figure>

probing almost the same wavelength range, the value reported for the W2 filter is 0.28 mag fainter.

## 3.3. X-ray and EUV SED modeling

To model the SEDs in the unobserved extreme UV (EUV; $12.4 \le \lambda/\AA < 1216$), we followed the method discussed in Lusso et al. (2012) and Shen et al. (2020), who use a power-law with a fixed slope of $\lambda L_\lambda \propto \lambda^{0.8}$ in the range between $500 \le \lambda/\AA \le 1216$, which was derived from HST observations of local AGN (Zheng et al. 1997; Telfer et al. 2002; Lusso et al. 2015), plus a power law with a free-to-vary spectral slope $\alpha$ that links $L_{500\AA}$ with the 1 keV luminosity. Thus, the full SED template takes the following form:

$$ \begin{cases} f^{\text{SED}}(\lambda) \text{ at } \lambda > 1216\AA, f^{\text{SED}}(\lambda) = \text{lum-K13 or WISSH-S23} \\ \lambda L_{\lambda} \propto \lambda^{0.8} \text{ at } 500\AA < \lambda \le 1216\AA \\ \lambda L_{\lambda} \propto \lambda^{-(\alpha+1)} \text{ at } 12.40\AA < \lambda \le 500\AA. \end{cases} \quad (3) $$

The derived $\alpha$ for HYPERION QSOs are reported in Table 2. The same is valid for the E-XQR-30 QSOs but, as there are no 1 keV luminosities, we employed the value obtained as the geometric mean of the HYPERION sample, which is $\text{Log}(\lambda L_{1\text{kev}}) = 45.01 \pm 0.08$ erg/s. The advantage of using a double power law to describe EUV emission lies in the fact that it allows us to make use of both the information available on the QSOs emission just below the $Ly\alpha$ line (from low-$z$ sources) and the high-quality X-ray data obtained for these sources at $z \ge 6$. This would not have been possible if using a fixed template or extrapolating with a single power law.

As the EUV contributes nearly 50% of the QSO emission, different modelings could lead to significantly different derived $L_{bol}$, which is the main value we want to derive. Therefore, we tested the impact of different modelings on the integrated luminosity under the EUV region. Our double-power-law parametrization yields, on average, EUV luminosities that are ~20% larger than those derived using the EUV modeling of
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Results of the SED fitting of the HYPERION QSOs using the lum-K13 (blue) and S23 (red) templates. The dotted blue and dashed red lines give the unextincted SEDs, which have been extended below the $Ly\alpha$ with a double power law as explained in Sect. 3.3; $L_{bol}$ is the integral under the dotted/dashed lines. Data at UV wavelengths shorter than $Ly\alpha$ were not considered in the fitting and are colored gray. The cyan stars indicate measured 1 keV luminosities; see Table 3. The subpanels show the residuals with respect to the K13 best-fit SED. The black circles refer to the residuals computed using only the uncertainties on the photometric points, while the red crosses represent residuals for which also the scatter on the SED template was taken into account.</figcaption></figure>

lum-K13, with differences of up to 40%. Conversely, employing a single-power-law extrapolation yields an average decrease of 8%, with a maximum deviation of 38%. Finally, by employ-

ing the far redder EUV spectral slope recently proposed by Cai & Wang (2023), we obtain EUV luminosities that are 70% lower than those obtained when using the double-power-law
<table><caption>**Table 2.** Spectral slopes of the HYPERION QSOs.</caption><thead><tr><th>Name</th><th>α<sub>EUV</sub></th><th>β<sub>UV</sub></th><th>γ<sub>opt</sub></th></tr></thead><tbody><tr><td>J1342+0928</td><td>-1.88</td><td>-0.15 ± 0.17</td><td></td></tr><tr><td>J1007+2115</td><td>-2.01</td><td>-0.74 ± 0.29</td><td></td></tr><tr><td>J1120+0641</td><td>-1.83</td><td>-1.21 ± 0.39</td><td>0.09 ± 0.11</td></tr><tr><td>J0038-1527</td><td>-2.0</td><td>-0.61 ± 0.22</td><td></td></tr><tr><td>J0252-0503</td><td>-1.84</td><td>-0.42 ± 0.25</td><td></td></tr><tr><td>J0020-3653</td><td>-1.73</td><td></td><td></td></tr><tr><td>J0411-0907</td><td>-2.04</td><td>-0.55 ± 0.15</td><td></td></tr><tr><td>J0244-5008</td><td>-1.65</td><td>-0.44 ± 0.12</td><td></td></tr><tr><td>J231-20.8</td><td>-2.07</td><td>-0.31 ± 0.25</td><td></td></tr><tr><td>J036+03.0</td><td>-2.21</td><td>-0.23 ± 0.13</td><td></td></tr><tr><td>J0224-4711</td><td>-1.89</td><td>-1.38 ± 0.14</td><td></td></tr><tr><td>J011+09</td><td>-1.81</td><td>-0.52 ± 0.11</td><td></td></tr><tr><td>J1148+5251</td><td>-1.93</td><td>-0.45 ± 0.19</td><td>-0.51 ± 0.10</td></tr><tr><td>J083+11.8</td><td>-2.28</td><td>-0.61 ± 0.14</td><td></td></tr><tr><td>J0100+2802</td><td>-1.94</td><td>-0.34 ± 0.11</td><td></td></tr><tr><td>J025-33</td><td>-2.02</td><td>-0.48 ± 0.09</td><td></td></tr><tr><td>J0050+3445</td><td>-2.12</td><td>-1.26 ± 0.22</td><td>-0.24 ± 0.16</td></tr><tr><td>J029-36</td><td>-1.9</td><td>-0.29 ± 0.16</td><td></td></tr></tbody></table>

**Notes.** EUV (12.4 ≤ λ/Å ≤ 500, not fitted), UV (1300 ≤ λ/Å ≤ 3000) and optical (3000 ≤ λ/Å ≤ 11 000) derived spectral slopes for the HYPERION sample.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Different parametrizations considered for the EUV modeling of the SEDs applied to J036.5+03.0; see Sect. 3.3.</figcaption></figure>

parametrization. Figure 3 shows these three different modelings for the case of J036+03.0, which, in terms of optical-to-X-ray ratio, can be considered as representative of the entire HYPERION sample (see Zappacosta et al. 2023).

We also tried to quantify the possible uncertainties associated with the fact that we assumed an average $L_{1\text{keV}}$ calculated from the HYPERION sample for E-XQR-30 sources. Indeed, although the proper X-ray emission –that is, ~20 Å–1 keV– only accounts for less than 1% of $L_{\text{bol}}$ in these luminous QSOs (e.g., Duras et al. 2020), the 1 keV luminosity constrains the slope of the second power law in our modeling of the unseen EUV SED, and thus may potentially affect $L_{\text{bol}}$. To quantify the impact of this assumption, we simulated a QSO with $\log(L_{2500\text{Å}}) = 47$ erg/s and we calculated how much $L_{\text{bol}}$ varies in the case where it has an X-ray luminosity of twice (half of) that predicted by the relation

by Lusso et al. (2010), finding an increase (decrease) in $L_{\text{bol}}$ of ~4%.

## 4. Bolometric luminosities, E[B-V], and bolometric corrections

The bolometric luminosities for the HYPERION and E-XQR-30 QSOs were computed by integrating their dereddened SEDs obtained using the lum-K13 template between 1 keV and 1 µm, where the choice of integration limits ensures that reprocessed radiation is not included and thus that the same contribution is not counted twice (see Marconi et al. 2004). Assuming that $L_{\text{bol}}$ were computed without dereddening the SED, we would obtain lower values for $L_{\text{bol}}$, with a difference of $\Delta \log(L_{\text{bol}}) \approx 4 \times E[B-V]$, which is nearly linear within the range of $E[B-V]$ explored by our sources. $L_{\text{bol}}$ values are provided in Table 3 (HYPERION sources) and on Zenodo (E-XQR-30 sources), together with monochromatic luminosities at various wavelengths obtained by interpolating two adjacent photometric points and correcting for the SED-fitting-derived $E[B-V]$ of the source. We also calculated $\lambda_{\text{Edd}}$ assuming the $M_{\text{SMBH}}$ derived from the MgII line, single-epoch virial-mass estimator (see Table 1 and Zappacosta et al. 2023; Mazzucchelli et al. 2023, for further references and for the associated uncertainties).

In Fig. 4, we show the luminosity distributions of both the HYPERION and E-XQR-30 samples, which range within 46.5 ≤ $\log(L_{\text{bol}}/[\text{erg/s}]) \le 48.1$; performing a Kolmogorov-Smirnov test on the $L_{\text{bol}}$ distributions of the two samples gives a p-value of ~0.05 (i.e., it does not reject the hypothesis that both samples have been drawn from the same distribution), thus strengthening our choice of using the E-XQR-30 QSOs as a complementary sample to HYPERION.

In Fig. 5 we compare $L_{\text{bol}}$ with those computed via BC and $L_{\text{BC}}$ (Mazzucchelli et al. 2023; Zappacosta et al. 2023, and references therein). We find that, when integrating the full SED, we obtain systematically lower values with $\langle\Delta(L_{\text{bol}})\rangle = 0.13 \pm 0.08$ dex on average (i.e., $L_{\text{BC}}$ are overestimated by 34%) and discrepancies of up to ~0.25 dex in a few QSOs. However, we note that $L_{\text{BC}}$ were derived using the 3000 Å bolometric correction $\text{BC}_{3000\text{Å}} = 5.15$ (Richards et al. 2006; Shen et al. 2011), which is obtained by integrating the entire SED between 10 keV and 100 µm. Therefore, including also X-ray and IR reprocessed emission leads to overestimation of $L_{\text{bol}}$ when compared to our values. This overestimation also impacts the determination of $\lambda_{\text{Edd}}$, which we find to be $0.12 \pm 0.06$ dex lower than the previous measurement on average.

Regarding the derived $E[B-V]$, the distribution is dominated by sources not affected by dust reddening, that is, with $E[B-V] = 0$; see the right panel of Fig. 4. However, we notice a slight disagreement between the HYPERION and E-XQR-30 samples; while only J0224-4711 among the HYPERION QSOs has $E[B-V] \ge 0.05$, there are 11 E-XQR-30 sources with $E[B-V]$ values exceeding this value (excluding J0224-4711). Given the relation between dust extinction and the presence of BAL features (Gallagher et al. 2007; Bischetti et al. 2022), this difference can be accounted for by the fact that, while broad absorption lines are found in 47% of E-XQR-30 sources, only J0038-1527 (Wang et al. 2018) and J231-20 (Bischetti et al. 2022) have published BAL features among HYPERION QSOs (although, to date, no BAL feature analysis has been performed as systematically and consistently as that carried out for the E-XQR-30 QSOs). A potential caveat regarding the correlation between BAL troughs and increased dust extinction is that QSOs
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Bolometric and monochromatic luminosities and $E[B-V]$ distributions for the HYPERION (blue) and E-XQR-30 (purple) samples.</figcaption></figure>

may appear reddened not due to genuine reddening but because their intrinsic fluxes are underestimated in regions of the spectrum affected by absorption features. However, while this could indeed be true for low measured $E[B-V]$ values (~0.01–0.03), the absorption features typically observed in the E-XQR-30 sample are not deep enough to influence flux measurements to the extent that they could mimic the effects of significant $E[B-V]$. Moreover, Bischetti et al. (2022) already noted that BAL QSOs in the E-XQR-30 sample exhibit redder W1–W2 colors, which are bands that do not probe wavelengths affected by absorption troughs.

As mentioned, the measurement of $L_{bol}$ for high-z QSOs usually relies on the application of a bolometric correction, especially the 3000 Å one. As the values we used were derived from lower-z sources, we investigated whether these values are appropriate for objects at the EoR or it is necessary to use different values. In the left panel of Fig. 6, we show our derived BC$_{3000\AA}$ for the HYPERION and E-XQR-30 QSOs compared to the constant bolometric correction by Richards et al. (2006), of namely BC$_{3000\AA}$ = 5.15, which has been adopted in the literature as a standard BC$_{3000\AA}$ for $z > 6$ QSOs (e.g., Wu et al. 2015; Bañados et al. 2018; Reed et al. 2019; Shen et al. 2019; Yang et al. 2021; Farina et al. 2022; Mazzucchelli et al. 2023; Fan et al. 2023; Zappacosta et al. 2023). We also show a re-evaluation of this value by Runnoe et al. (2012), of namely BC$_{3000\AA}$ = 3.11, who exclude reprocessed IR emission, which allows a fair comparison with our value. The BC$_{3000\AA}$ distribution of HYPERION and E-XQR-30 sources ranges between 2.64 and 3.94, with a mean value of $3.30 \pm 0.3$. The majority of derived BCs are higher than the recalculated BC$_{3000\AA}$ by Runnoe et al. (2012), although they are still in full agreement within the uncertainties. This small discrepancy, rather than arising from a different SED shape in our objects, can be attributed to our more careful modeling of the EUV (i.e., with the double power law) which, as mentioned above, results in higher $L_{bol}$ (and consequently higher BC) than those obtained using an average SED template. The right panel of Fig. 6 illustrates that when computing $L_{bol}$ with the lum-K13 EUV recipe, the BCs are fully consistent with the predictions by Runnoe et al. (2012).

Given our more refined treatment of the UV-X-ray SED, which gives more accurate $L_{bol}$ measurements, we recommend using the revised BC$_{3000\AA}$ value instead of 5.15 to recover a more

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** $L_{bol}$ already published (Zappacosta et al. 2023; Mazzucchelli et al. 2023) vs $L_{bol}$ computed in this work. $L_{bol}$ has been computed in the literature using the 3000 Å BC ($L_{bol}$ = 5.15 $L_{3000\AA}$). The black line gives the 1:1 ratio, while the gray shaded area delimits the 0.1 dex difference region.</figcaption></figure>

accurate $L_{bol}$ for $z > 6$ QSOs. This allows the removal of a further source of inaccuracy in the determination of $\lambda_{Edd}$, which is already affected by systematic uncertainties from the mass determination.

Additionally, we computed optical bolometric corrections at 4400 Å and 5100 Å, comparing them to relationships established for lower-z sources, as depicted in Fig. 7. For $\lambda = 4400$ Å, the bolometric corrections demonstrate a good agreement with the average value reported in Duras et al. (2020) for a large collection of AGN at $0 \le z \le 3$, resulting in a mean of $4.9 \pm 0.7$. In the case of 5100 Å BC, our derived BC$_{5100\AA}$ values are consistent with the predictions at higher luminosities extrapolated from the relationship provided by Runnoe et al. (2012). However, they almost all exceed the value reported in Krawczyk et al. (2013), although our mean value of BC$_{5100\AA}$ = $5.7 \pm 0.9$ remains in agreement within $2\sigma$ with that of these latter authors.
<figure><img src="image_7.png" /></figure>

**Fig. 6.** BC<sub>3000 Å</sub> vs L<sub>bol</sub> for HYPERION and E-XQR-30 QSOs with their weighted mean value reported as a dark-blue square. The left panel illustrates L<sub>bol</sub> computation using Eq. (3), while the right panel showcases results obtained through the lum-K13 modeling of the EUV. The constant 3000 Å BCs proposed by Richards et al. (2006) are overplotted. The original value of 5.15 is depicted as a black dot-dashed line while, the version recomputed by Runnoe et al. (2012) to exclude IR emission is represented as a red solid line.

<figure><img src="image_8.png" /></figure>

**Fig. 7.** Bolometric corrections at 4400 (left) and 5100 Å (right) with their weighted mean values reported as blue squares. Overplotted are the relationships by Duras et al. (2020), Krawczyk et al. (2013, with the associated uncertainties as dotted lines), and Runnoe et al. (2012).

**5. Predicting X-ray luminosity with a physically motivated AD model**

We made use of QSOSED (Kubota & Done 2018), a physically motivated model for disk and corona emission, to describe the nuclear emission of QSOs. Specifically, we aimed to determine whether or not these models can accurately predict the coronal emission (i.e., the 1 keV luminosities) by fitting only the accretion disk data, which include the rest-frame UV and optical photometric points. Accordingly, our analysis is limited to the 18 HYPERION QSOs. QSOSED assumes that the emitted radiation originates from three distinct regions at increasing radial distances: (i) a hot ($k_B T \ge 10$ keV) comptonizing plasma (i.e., the corona; e.g., Haardt & Maraschi 1991) at $R_{\text{ISCO}} < r \le R_{\text{hot}}$, which is responsible for the power-law-like X-ray continuum emission, (ii) a warm ($k_B T \le 5-10$ keV) comptonizing plasma in an intermediate region, $R_{\text{hot}} < r \le R_{\text{warm}}$, which contributes to the soft excess (Petrucci et al. 2018), and (iii) a standard, geometrically thin, optically thick accretion disk at $r > R_{\text{warm}}$. Furthermore, QSOSED assumes that the hot corona radiates 2% of the Eddington luminosity of a source, regardless of its L<sub>bol</sub>.

SED templates calculated by QSOSED were generated through the X-ray spectral fitting program XSPEC (Arnaud 1996) by taking into account four main free parameters, namely $M_{\text{SMBH}}$, the BH spin, the specific mass-accretion rate, $\dot{m}$, and the cosine of the inclination angle $i$, $\cos(i)$. Given the high probability of resulting in degenerate templates given the relatively low number of degrees of freedom (i.e., photometric points), we fixed the BH spin to zero, $M_{\text{SMBH}}$ to the values listed in Zappacosta et al. (2023), and the viewing angle to 30°, which is an average value among those expected for type-1 QSOs (e.g.,

Mountrichas et al. 2021). Accordingly, the only parameter left free to vary is $\log(\dot{m})$, which ranges from -1 to 0.4 in steps of 0.1. Furthermore, we did not account for dust reddening as we derived it from the empirical modeling, finding it to be negligible for the vast majority of the sources (see Table 3). In QSOSED, the luminosity from the QSOs and thus the observed flux at each wavelength are predetermined without any normalization. Given the many assumptions on our input parameters, and the systematic uncertainties on $M_{\text{SMBH}}$, (~0.55 dex; e.g., Mazzucchelli et al. 2023), we expect a limited accuracy in the SED description, which, we emphasize, is not the main result we aim to achieve with this analysis. In this model, we cannot continuously explore the parameter space, as the allowed values of $\dot{m}$ are discretized. Thus, the best fit was determined via $\chi_v^2$ minimization, while uncertainties were estimated by identifying the range of $\dot{m}$ that resulted in $\Delta\chi_v^2 \le 1$ from the best fit, corresponding to a $1\sigma$ confidence interval (e.g., Avni 1976).

The results from the fits, shown in Fig. 8 together with the derived $\dot{m}$, indeed reveal a relatively large range of $\chi_v^2$ values, varying from 1.4 for J1120+0641 to ~20, with J029-26 being a notable outlier with $\chi_v^2 = 77$ and a $\chi_v^2$ median value of 7.2. Interestingly, among the 13 QSOs for which the QSOSED template provides an acceptable fit (arbitrarily chosen to be those with $\chi_v^2 \le 15$), the 1 keV fluxes predicted by QSOSED are consistent with the measured ones within $3\sigma$ for 7 sources (54%) and within $5\sigma$ for 9 sources (70%). By examining the difference in the ratio of predicted to measured fluxes, QSOSED models deviate by a median factor of 1.6, with a maximum discrepancy factor of 6.3 in the case of J083+11.8. This outcome highlights the effectiveness of QSOSED in predicting X-ray emission based solely on the observed luminosity from the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Same as Fig. 2 but with QSOSED generated templates. Each panel presents the mass accretion and the $\chi^2$ derived by the fit.</figcaption></figure>

accretion disk even in the case of the most distant QSOs. We note that, for the three sources showing a significant difference between observed and QSOSED predicted X-ray fluxes (namely, J0411–0907, J083+11.8, and J0050+2445), the observed X-ray flux is always lower than the predicted one by factors of 2.7, 6.3, and 3.8, respectively.

The deviation of the QSOSED prediction at 1 keV from the actual data may reflect either (i) the deviations from intrinsic physical properties (e.g., mass measurement, spin assumption, source inclination, and lack of extinction) for some sources or (ii) the underlying assumptions on which this model is based, such as the fraction of the X-ray coronal emission in terms of $L_{Edd}$, and $R_{\text{hot}}$ and $R_{\text{warm}}$, the temperature of the hot and warm Comptonizing regions. An investigation into the cause of the disagreement is beyond the aim of this paper and is deferred to future work.

Finally, in light of the good predictive results demonstrated by QSOSED, we performed the fitting on E-XQR-30 QSOs, taking into account the same assumptions (i.e., $i = 30^\circ$, spin = 0),

in order to estimate their 1 keV fluxes, which are reported, converted to luminosities, on Zenodo, along with the other derived properties.

## 6. Ultraviolet and optical slopes

The optical–UV SED is usually modeled as a broken power law, with the break between the two slopes falling at $\lambda = 3000–5000$ Å (e.g., Vanden Berk et al. 2001; Temple et al. 2021a). To derive the slopes of the HYPERION and E-XQR-30 QSOs, we opted to follow Lusso & Risaliti (2016), that is, we set the break at $\lambda = 3000$ Å; although we extended the UV interval down to 1300 Å instead of 1450 Å in order to include the Y band photometric point for a greater number of QSOs, without being affected by the peak of the $\text{Ly}\alpha$ emission. Therefore, we refer to $\beta_{\text{UV}}$ as the UV slope, characterizing the accretion disk emission between $\lambda = 1300$ Å and 3000 Å, and $\gamma_{\text{opt}}$ as the optical slope, describing the emission at longer wavelengths, up to
$\lambda = 1.0$ µm, where the SED shows an inflection due to the arising contribution from hot dust. Spectral slopes are reported in luminosity density units, that is, in the form $L_v \propto \gamma^{\beta}$.

Slopes were computed independently via least-square minimization by fitting a straight line to the photometric points in the log-log space. To compute the slope, we required the QSOs to have at least three photometric points in the wavelength interval of interest; this requirement limits the number of QSOs with measured $\beta_{UV}$ and $\gamma_{opt}$ to 39 and 30, respectively. In particular, all but three sources with computed $\gamma_{opt}$ belong to the E-XQR-30 sample, as they have lower redshifts (indeed at $z \ge 6.3$, the K band moves to shorter wavelengths than 3000 Å). Figure 9 shows the computed $\beta_{UV}$ and $\gamma_{opt}$ distribution (reported in Table 2 for HYPERION QSOs and on Zenodo for E-XQR-30 ones). For comparison, we also show the contour enclosing the 39th and 86.4th percentiles (respectively 1 and $2\sigma$ under the assumption of a bivariate normal distribution) of the values derived from the Krawczyk et al. (2013) luminous subsample following an analogous methodology. Additionally, as the derived slope could vary as a function of $z$ simply because the filters are sampling different regions of the SED, we also overplot the mean $\pm 1\sigma$ values obtained by redshifting the lum-K13 template at $5.5 \ge z \ge 7.5$ and generating mock observations with UKIDSS YJHK plus W1 and W2 filters for different levels of $E[B-V]$.

The HYPERION and E-XQR-30 sources are in complete agreement with the ranges of values observed in luminous QSOs at lower $z$, and, considering the uncertainties—which are quite large due to the fact that the fitting is performed on a limited number of points—, the bulk of the sample agrees with the point obtained from the lum-K13 template without dust extinction. QSOs with a lower $\beta_{UV}$ (i.e., a flatter SED) can be explained through dust reddening and indeed they are generally those with higher $E[B-V]$ values. In particular, the three HYPERION sources with lower $\beta_{UV}$ are those with measured $E[B-V] \ge 0.03$ (see Table 3). Conversely, two E-XQR-30 QSOs exhibit bluer $\beta_{UV}$ than the rest of the population, although still within $\sim 2\sigma$ of the distribution observed in Krawczyk et al. (2013), while PSOJ023-02, located in the bottom-right corner of Fig. 9, shows a $\gamma_{opt}$ value that is in significant disagreement with the rest of the distribution. As already mentioned in Sect. 3.1, this object has a rather flat spectrum and hence the low $\gamma_{opt}$ is expected. Finally, we tested for potential correlations between $\beta_{UV}$ and both the photon index $\Gamma$ and the UV-to-X-ray ratio $\alpha_{OX}$, defined as $-0.384 \log (L_{\nu,2\text{keV}}/L_{\nu,2500\text{Å}})$ and reported in Zappacosta et al. (2023) and Tortosa et al. (2024). However, we find no evidence of correlation among these parameters. Indeed we simulated multiple realizations of the datasets by generating synthetic values normally distributed around the best-fit value and with a $\sigma$ equal to their uncertainties and, performing Spearman tests for each of them, obtained mean correlation coefficients and p-values of 0.22 and 0.4 for the $\Gamma$ versus $\beta_{UV}$ relation and -0.05 and 0.6 for $\alpha_{OX}$ versus $\beta_{UV}$.

## 7. Near-infrared modeling of the hot-dust component

It is intriguing to compare the positions of MIR points (Spitzer MIPS 24 µm, W3, and W4) to what is expected from the normalized SED templates. Previous studies (e.g., Jiang et al. 2006; Leipski et al. 2014; Bosman et al. 2024) reported a NIR emission of QSOs at the EoR resembling that of luminous objects at lower $z$, characterized by a black-body-like continuum originating from a hot-dust component close to its sublimation temperature.

<figure><img src="image_10.png" /><figcaption>Fig. 9. $\beta_{UV}$ (1300 < $\lambda$/Å < 3000) and $\gamma_{opt}$ (3000 < $\lambda$/Å < 10 000) distributions for the sample. Only the UV slope has been measured for the points along the constant line $\gamma_{opt} = -1.5$. Blue stars represent the values obtained from the lum-K13 template assuming $E[B-V]$ values of 0, 0.05, and 0.1, and assuming observations in the UKIDSS YJHK bands plus W1 and W2 for QSOs at $5.5 \ge z \ge 7.5$. The shaded light-blue areas delimit the 1 and $2\sigma$ confidence intervals for the luminous QSO subsample in Krawczyk et al. (2013).</figcaption></figure>

However, there are several indications suggesting a higher fraction of dust-poor objects at the EoR, with a lower-than-average NIR-to-optical emission ratio (Jiang et al. 2006; Leipski et al. 2014) or even dust-free sources (Jiang et al. 2010). In these cases, the NIR can be well-modeled by assuming only emission from the accretion disk, which is modeled as a power law, without the need to add emission from hot dust. The cause of the observed lack of dust in these QSOs remains uncertain. It could be attributed to a genuine deficiency of dust, possibly linked to these sources being in the early stages of the Universe (Jiang et al. 2010). Alternatively, it may be associated with a distinct torus structure, characterized by a lower covering factor, resulting in less dust being directly exposed to the primary radiation (see Lyu et al. 2017).

In total, we have 11 QSOs where the 24 µm or the W4 photometry is available, which is necessary to constrain the hot dust emission. We find that in 4 of these, the lum-K13 SED agrees with the true NIR luminosity, while 5 objects have enhanced NIR emission and 2 objects have substantially lower emission.

To better quantify the strength of the hot-dust emission, following the approach by Temple et al. (2021b), we computed $X_{HD} = L_{BB}/L_{AD}$ for our sources, which is the ratio at $\lambda = 2$ µm between the hot-dust and accretion disk components. More specifically, we modeled the QSO SEDs at $\lambda > 3000$ Å with a power law describing the emission from the accretion disk plus a single-temperature black-body, which represents the emission from hot dust located in the innermost layer of the torus, that is,

$$ \lambda L_{\lambda} = K_{1} \lambda^{\gamma_{\text{opt}}} + K_{2} \frac{1}{\lambda^{4}} \frac{1}{\exp(hc/\lambda k_{\text{B}} T) - 1} $$

Following Temple et al. (2021b), we fixed both the temperature of the black body ($T = 1280$ K) and the slope of the power
<table><caption>**Table 3.** Bolometric and monochromatic luminosities of the HYPERION QSOs.</caption><thead><tr><th>Name</th><th>L<sub>bol</sub></th><th>λL<sub>2500Å</sub></th><th>λL<sub>3000Å</sub></th><th>λL<sub>4400Å</sub></th><th>λL<sub>5100Å</sub></th><th>E[B-V]</th><th>λ<sub>Edd</sub></th><th>λL<sub>1 keV</sub></th></tr></thead><tbody><tr><td>J1342+0928</td><td>47.05 ± 0.04</td><td>46.60 ± 0.02</td><td>46.46 ± 0.02</td><td>46.20 ± 0.05</td><td>46.14 ± 0.07</td><td>0.01</td><td>1.11 ± 0.38</td><td>44.93 ± 0.09</td></tr><tr><td>J1007+2115</td><td>47.19 ± 0.07</td><td>46.75 ± 0.03</td><td>46.67 ± 0.02</td><td>46.50 ± 0.03</td><td>46.44 ± 0.04</td><td>0.03</td><td>0.81 ± 0.16</td><td>44.87 ± 0.18</td></tr><tr><td>J1120+0641</td><td>47.22 ± 0.05</td><td>46.84 ± 0.05</td><td>46.76 ± 0.05</td><td>46.59 ± 0.01</td><td>46.49 ± 0.01</td><td>0.05</td><td>0.50 ± 0.28</td><td>45.19 ± 0.07</td></tr><tr><td>J0038-1527</td><td>47.25 ± 0.06</td><td>46.87 ± 0.03</td><td>46.78 ± 0.03</td><td>46.53 ± 0.02</td><td>46.42 ± 0.03</td><td>0.03</td><td>1.00 ± 0.23</td><td>44.96 ± 0.16</td></tr><tr><td>J0252-0503</td><td>47.00 ± 0.06</td><td>46.59 ± 0.04</td><td>46.48 ± 0.04</td><td>46.24 ± 0.04</td><td>46.15 ± 0.06</td><td>0.01</td><td>0.56 ± 0.08</td><td>44.95 ± 0.19</td></tr><tr><td>J0020-3653</td><td>47.10 ± 0.06</td><td>46.69 ± 0.05</td><td>46.66 ± 0.04</td><td>46.46 ± 0.02</td><td>46.35 ± 0.03</td><td>0.02</td><td>0.58 ± 0.13</td><td>45.21 ± 0.09</td></tr><tr><td>J0411-0907</td><td>47.17 ± 0.04</td><td>46.74 ± 0.02</td><td>46.66 ± 0.02</td><td>46.47 ± 0.02</td><td>46.43 ± 0.03</td><td>0.01</td><td>1.84 ± 0.24</td><td>44.82 ± 0.08</td></tr><tr><td>J0244-5008</td><td>47.07 ± 0.04</td><td>46.58 ± 0.03</td><td>46.54 ± 0.02</td><td>46.41 ± 0.01</td><td>46.37 ± 0.02</td><td>0.01</td><td>0.77 ± 0.27</td><td>45.31 ± 0.05</td></tr><tr><td>J231.6-20.8</td><td>47.13 ± 0.05</td><td>46.70 ± 0.04</td><td>46.59 ± 0.07</td><td>46.29 ± 0.06</td><td>46.21 ± 0.10</td><td>0.01</td><td>0.34 ± 0.15</td><td>44.73 ± 0.13</td></tr><tr><td>J036.5+03.0</td><td>47.23 ± 0.04</td><td>46.79 ± 0.03</td><td>46.72 ± 0.03</td><td>46.50 ± 0.02</td><td>46.40 ± 0.03</td><td>0.0</td><td>0.44 ± 0.13</td><td>44.61 ± 0.14</td></tr><tr><td>J0224-4711</td><td>47.47 ± 0.05</td><td>47.04 ± 0.02</td><td>47.04 ± 0.02</td><td>46.82 ± 0.01</td><td>46.75 ± 0.01</td><td>0.07</td><td>0.75 ± 0.19</td><td>45.35 ± 0.06</td></tr><tr><td>J011+09</td><td>46.87 ± 0.04</td><td>46.39 ± 0.02</td><td>46.42 ± 0.02</td><td>46.06 ± 0.04</td><td>45.92 ± 0.05</td><td>0.01</td><td>0.42 ± 0.04</td><td>44.86 ± 0.11</td></tr><tr><td>J1148+5251</td><td>47.48 ± 0.04</td><td>46.95 ± 0.02</td><td>46.91 ± 0.02</td><td>46.81 ± 0.01</td><td>46.76 ± 0.01</td><td>0.02</td><td>0.44 ± 0.05</td><td>45.30 ± 0.09</td></tr><tr><td>J083.8+11.8</td><td>47.13 ± 0.04</td><td>46.71 ± 0.03</td><td>46.64 ± 0.03</td><td>46.48 ± 0.02</td><td>46.38 ± 0.03</td><td>0.01</td><td>0.51 ± 0.16</td><td>44.41 ± 0.07</td></tr><tr><td>J0100+2802</td><td>47.99 ± 0.03</td><td>47.54 ± 0.03</td><td>47.57 ± 0.06</td><td>47.37 ± 0.01</td><td>47.32 ± 0.01</td><td>0.0</td><td>0.70 ± 0.11</td><td>45.79 ± 0.03</td></tr><tr><td>J025-33</td><td>47.38 ± 0.04</td><td>46.93 ± 0.01</td><td>46.89 ± 0.01</td><td>46.69 ± 0.01</td><td>46.64 ± 0.01</td><td>0.01</td><td>0.52 ± 0.04</td><td>45.06 ± 0.14</td></tr><tr><td>J0050+3445</td><td>47.24 ± 0.06</td><td>46.84 ± 0.02</td><td>46.78 ± 0.03</td><td>46.55 ± 0.01</td><td>46.51 ± 0.01</td><td>0.04</td><td>0.29 ± 0.06</td><td>44.76 ± 0.09</td></tr><tr><td>J029-36</td><td>47.17 ± 0.06</td><td>46.66 ± 0.03</td><td>46.60 ± 0.02</td><td>46.48 ± 0.01</td><td>46.42 ± 0.01</td><td>0.02</td><td>0.17 ± 0.03</td><td>45.03 ± 0.09</td></tr></tbody></table>

**Notes.** Logarithm of bolometric and monochromatic luminosities of the HYPERION sample in units of [erg/s] derived from the lum-K13 best fit. Monochromatic luminosities are corrected for the $E[B-V]$. The Eddington ratios $\lambda_{Edd}$ were computed assuming the $M_{BH}$ reported in Table 1. $L_{1\text{keV}}$ was derived from the observed 0.3–7 keV fluxes by adopting the photon index $\Gamma$ reported in Zappacosta et al. (2023) and Tortosa et al. (2024).

law ($\gamma_{opt} = -0.16$). We kept T fixed, as the only point useful to constrain the black-body emission is the Spitzer 24 µm one (or, alternatively, W4), while we fixed $\gamma_{opt}$ to compare our results against those by Temple et al. (2021b), who showed that the derived $X_{HD}$ values strongly depend on the assumed $\gamma_{opt}$. However, we note that the assumed slope is remarkably close to the one we obtain in Sect. 8 for the mean SED and can therefore be considered as an average value; although it is slightly different from the mean of the individual $\gamma_{opt}$ computed in Sect. 6, which is $\langle\gamma_{opt}\rangle = -0.23$. In addition, by visually inspecting the fitting results, we find that the optical emission of the 11 analyzed sources is generally well described by this power law. As done for the SED fitting routine, the constants of normalization $K_1$ and $K_2$ were derived by minimizing the likelihood.

The results of the fitting are shown in Fig. 10. For each of the analyzed QSOs, it is necessary to add a hot-dust component, and therefore there are no QSOs without a dust component in our sample. We find the median value of $X_{HD}$ to be 2.4, with a MAD of 1.4, which is in agreement with the 2.5 median value reported in Temple et al. (2021b); moreover, a K-S test ($p$-value = 0.19) indicates no evidence of a different underlying $X_{HD}$ distribution compared to that derived by Temple et al. (2021b).

Specifically, two sources, namely J0100+2802 and SDSSJ0836+0054, exhibit $X_{HD} \le 1$. While only J0100+2802 can be considered –within $1\sigma$– as properly dust poor according to the criterion set by Jun & Im (2013), which requires $X_{HD} \le 0.15$ (see Temple et al. 2021b), both sources clearly show much lower dust emission than expected. Notably, SDSSJ0836+005 was previously identified as a QSO with dust deficiency by Leipski et al. (2014). On the other hand, three QSOs, ULASJ0148+0600, SDSSJ0842+1218 and SDSSJ0927+20, are found to have $X_{HD} \ge 4.7$.

Temple et al. (2021b) reported a positive correlation between $X_{HD}$ and the velocity of the CIV emission line relative to the systemic MgII-derived redshift (hereafter, CIV blueshifts), suggesting objects with faster winds have stronger hot-dust emission. Figure 11 shows the addition of the 11 HYPER-

ION and E-XQR-30 sources to the original plot of these latter authors. Our sources do not show any clear CIV-$X_{HD}$ correlation. However, this result is primarily due to two factors. First, the very small sample size poses challenges in determining the correlation. As a test, we randomly selected 11 sources from the data presented in Temple et al. (2021b), and only in 25% of the cases did we recover a positive correlation with a $p$-value < 0.05. Second, and more importantly, our sources exhibit significantly faster winds compared to the QSOs in Temple et al. (2021b), as expected for such luminous sources (Fiore et al. 2017; Vietri et al. 2018; Meyer et al. 2019; Timlin et al. 2020; Schindler et al. 2020). Consequently, a majority of our objects lie outside the contour enclosing 86.4% (i.e., $2\sigma$) of the Temple et al. (2021b) QSOs.

## 8. The mean spectral energy distribution of luminous, $z > 6$ QSOs

We generated a mean SED based on the combined HYPERION and E-XQR-30 samples to obtain a more comprehensive view of the broadband emission of luminous high-redshift QSOs. As our sample consists of QSOs with a relatively narrow redshift range, each observed band probes similar rest-frame wavelengths across all QSOs. Therefore, instead of constructing a continuous SED template, we determined the mean luminosity value for each band. With this approach, we do not need to extrapolate the SED to wavelengths not covered by certain QSOs. We also grouped Spitzer IRAC 3.6 and W1, Spitzer IRAC 4.5 and W2, Spitzer MIPS 24 and W4 as the same band; the J1120+0641 1 µm luminosity point was included with Spitzer MIPS 8 µm. To minimize the scatter due to the small sample size of our analyzed QSOs, we normalized their SEDs at 3500 Å. The normalization is particularly relevant for the W3 and W4 bands, where most sources in our sample have not been detected. Indeed, as brighter QSOs are more likely to be detected, using non-normalized data would result in a mean NIR emission that
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** NIR SEDs for the 11 QSOs with MIR photometry. The blue line gives the accretion disk emission modeled as a power law with a spectral index of -0.16, while the gold line describes the emission from hot dust modeled as a black body with T = 1280 K. The red lines show the sum of the two components.</figcaption></figure>

is not representative of the true mean properties of the sources. Although there is no physical reason to choose a specific normalization wavelength, we tested various wavelengths within the range of 2000 Å < λ < 7000 Å and found that normalizing at 3500 Å resulted in the lowest scatter across all bands. Therefore, we adopted this wavelength as our normalization point.

Figure 12 shows the resulting mean SED, with the uncertainties given as the standard deviation divided by the square root of the number of QSOs in each band, along with a comparison with the mean SEDs derived in Krawczyk et al. (2013) for both their luminous and whole QSO sample (dot-dashed blue line and dotted purple line, respectively) and WISSH-S23 (dashed green line).

As visible in Fig. 12, we do not find any significant deviations from the mean SEDs of luminous QSOs at lower z (i.e., lum-K13 and WISSH-S23), while the average SED by Krawczyk et al. (2013), labeled “all-K13”, has a flatter UV-optical slope and, depending on the chosen normalization wavelength, predicts larger emission at λ ~ 6000–9000 Å as in Fig. 12, or a lower UV bump. Therefore, we do not find any evolution of the SED with redshift, at least for the 1000 Å–1 µm wavelength interval. While there might be a selection effect, given that all our sources were identified by targeting objects with typical colors of lower-z QSOs (e.g., Reed et al. 2019), and therefore may not be universally applicable to the entire population of high-z AGN, this result strengthens our choice to employ templates of luminous QSOs for computing bolometric luminosities.

At wavelengths above 1 µm, the distribution of rest-frame NIR points shows a considerable spread; there is an indication of even stronger NIR emission than in the IR-selected WISSH

QSOs, although this needs to be further investigated, as the currently limited data do not allow robust conclusions to be made. Moreover, the W3 and W4 points are likely biased toward NIR-bright sources, as the limiting magnitudes in these bands are much shallower than those of W1 and W2.

We fitted the mean SED with a broken power law jointed at λ = 3000 Å, plus a blackbody with a fixed temperature of T = 1280 K, as described in Sect. 7. For all parameters, we assumed flat priors. We find βUV = -0.63 ± 0.05, which is in agreement with both the value reported in Telfer et al. (2002, i.e., -0.69, from the analysis of 184 QSOs with HST spectra) and that presented by Lusso et al. (2015, i.e., -0.61, derived from 53 luminous sources at redshift z ~ 2.4). This value of βUV is instead softer than both what is found in the composite spectrum by Selsing et al. (2016) derived from 102 QSOs at z = 1–2.1 and also that found by Temple et al. (2021a) from a subsample of bright SDSS QSOs (18.6 < i < 19.1), that is, βUV = -0.30 and -0.349, respectively. The power law modeling the mean SED redwards of 3000 Å is found to be steeper with a slope of γopt = -0.14 ± 0.04. Interpreting γ as the mean optical slope, we confirm that high-luminosity QSOs have a steeper optical continuum (e.g., Richards et al. 2006; Krawczyk et al. 2013) even at high redshift with respect to the bulk of the population at lower luminosities (i.e., γopt = -0.46; Vanden Berk et al. 2001). Interestingly, computing the hot dust-to-AD ratio with this modeling gives XHD = 3.28 ± 0.5, which is higher than both the median value found by Temple et al. (2021b) and that calculated as the mean of the individual QSOs.

The continuous SED template obtained by the broken power-law plus black-body modeling is shown in Fig. 13. To connect the template with the mean 1 kev luminosity, we extended the
<figure><img src="image_12.png" /></figure>

**Fig. 11.** *X*<sub>HD</sub> vs CIV blueshift adapted from Fig. 3 in Temple et al. (2021b) with the inclusion of HYPERION and XQR-30 QSOs. E-XQR-30 blueshifts are taken from Mazzucchelli et al. (2023), while HYPERION values are computed as the mean of all the blueshifts reported in the literature for each QSO; see Tortosa et al. (2024). Errors on *X*<sub>HD</sub> are given as the 16th and 84th percentiles. Cyan stars indicate QSOs with Spitzer MIPS 24 µm fluxes, while red stars have W4 photometry. The histograms show the distributions of CIV blueshift and *X*<sub>HD</sub> in Temple et al. (2021b, shown here in gray, and normalized by a factor of 500) and for the HYPERION and E-XQR-30 QSOs (red). The *p*-values obtained by performing a K-S test on the two distributions are reported next to the histograms.

SED in the EUV region using the same recipe as that reported in Eq. (3).

We also computed the mean SEDs of several subsamples obtained by splitting the sample based on the median of several physical properties (*M*<sub>SMBH</sub>, CIV blueshift, *L*<sub>bol</sub>) or on the criterion adopted to assemble the HYPERION sample (i.e., using *M*<sub>s,Edd</sub> = 1000 *M*<sub>☉</sub> as threshold) but found no significant difference either between them or with respect to the overall average SED (see Fig. 14). This result holds true even when comparing the mean values of *β*<sub>UV</sub> and *γ*<sub>opt</sub> of each subsample, as shown in Fig. 15. Indeed, each pair of values is in agreement with that of its complementary subsample within 2σ.

**9. Summary and conclusions**

In this work, we characterize the X-ray-to-NIR broad-band emission of luminous QSOs at the EoR with the main goal being to provide the first systematic investigation of the SED of these sources. Our main results can be summarized as follows:

- We produced X-ray-to-NIR SEDs for the 18 QSOs belonging to the HYPERION sample (Zappacosta et al. 2023; Tortosa et al. 2024). Exploiting the unprecedented quality of X-ray data from the HYPERION QSOs, we employed a double power law to characterize the X-ray to UV SED (Lusso et al. 2012, see Eq. 3); given the absence of data in that wavelength range, this modeling ensures the most reliable representation of the EUV region achievable at present. The UV-to-NIR SED of the QSOs was instead modeled

<figure><img src="image_13.png" /></figure>

**Fig. 12.** Photometry points for the HYPERION and E-XQR-30 samples, normalized at $\lambda = 3500$ Å. The red diamonds show the geometric mean value for each band with the associated $1\sigma$ as an error bar. Overplotted are also the mean SEDs by Krawczyk et al. (2013) (whole and luminous QSOs samples, purple dotted and blue dot-dashed lines, respectively) and the WISSH hyperluminous mean SED (dashed green line). The vertical dotted lines indicate the position of the $Ly\alpha$ at 1216 Å and the normalization wavelength.

<figure><img src="image_14.png" /></figure>

**Fig. 13.** Continuous mean SED obtained by fitting the mean points (red diamonds) with a broken power law ($\lambda_{break} = 3000$ Å) plus a single temperature ($T = 1280$ K) black body. The EUV SED ($\lambda < 1216$ Å) was reconstructed as discussed in Sect. 3.1. The light blue and gold shaded area describe the best-fit broken power law and black body, respectively. The gray diamond was excluded from the fit as it is also the average of points falling below the $Ly\alpha$. The black dotted lines indicate the $Ly\alpha$ and the normalization wavelength at 3500 Å.

using templates derived from luminous QSOs at lower z, namely lum-K13 (Krawczyk et al. 2013) and WISSH-S23 (Saccheo et al. 2023). We find that, for all QSOs, the optical–UV region is well described by the lum-K13 SED (see Fig. 2), while we obtained, on average, poorer results when employing the WISSH-S23 SED; although these latter span a luminosity range closer to the HYPERION sources.

- We increase the statistical significance of our analysis in the UV-to-NIR region by including an additional 36 QSOs drawn from the E-XQR-30 sample (D'Odorico et al. 2023), which share similar properties with HYPERION sources in terms of redshift and luminosity. We confirm that the lum-K13 template also provides excellent results in describing the SED of these sources.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Mean SEDs derived by splitting the sample according to the median value of several physical properties. Each panel reports the threshold value used to split the sample. Filled markers refer to QSOs above the threshold, while empty ones indicate QSOs below the threshold.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Mean $\gamma_{opt}$ vs mean $\beta_{UV}$ computed for each subsample. Errors are reported as the standard deviation of the values divided by the square root of the number of sources in the subsample. The yellow triangle shows the mean value for the full sample. As in Fig. 14, filled markers refer the subsample above the threshold, while empty ones refer to the one below.</figcaption></figure>

- By integrating the SEDs, we homogeneously computed bolometric luminosities for all 54 objects. Their values range from $46.8 \le \log(L_{bol}) \le 48.1$, as shown in Fig. 4. We find that our derived $L_{bol}$ are systematically lower than the ones already reported in the literature, with an average discrepancy factor of 0.13 dex.

- We calculated bolometric corrections at 3000, 4400, and 5100 Å for QSOs with redshifts of $z \ge 5.5$ and check their consistency with those derived from lower $z$ sources. Given the similarity of the SEDs of our QSOs to those of their lower $z$ counterparts, the obtained bolometric corrections are broadly consistent with literature findings. However, we observe that the values for $\lambda = 3000$ Å and 5100 Å are systematically larger than the bolometric correction val-

ues largely used in the literature. As an improved estimate for these sources, we show that better results are obtained when using $BC_{3000\AA} = 3.3$ and $BC_{5100\AA} = 5.7$ to derive $L_{bol}$.

- Taking advantage of the X-ray coverage of HYPERION, we verified the accuracy of theoretical models (Kubota & Done 2018) in predicting coronal emission from that of the accretion disk for these high-$z$ sources. Despite the many crude assumptions we make regarding the parameters, we surprisingly find that models and observations deviate by a median factor of just 1.6.

- We investigated the hot-dust emission in 11 QSOs with available MIR data by modeling their rest-frame NIR using a combination of a power law and a black body (see Sect. 7 and Fig. 10). Our analysis reveals that, for all these sources, the inclusion of a dust component is necessary, although to a varying extent. We quantified the hot-dust emission strength by computing the hot-dust to accretion-disk emission ratio $X_{HD}$ (Temple et al. 2021b). The overall $X_{HD}$ distribution is consistent with that found by Temple et al. (2021b, cf. our Fig. 11), but two QSOs exhibit a notably low $X_{HD}$, with one source, J0100+2802, being within $1\sigma$ of the threshold outlined in Jun & Im (2013) for classification as dust poor. On the other hand, three sources have a $X_{HD}$ that is approximately twice the average value observed in luminous QSOs.

- Finally, we derived a mean SED for these high-$z$ QSOs extending from the X-ray to NIR (i.e., from 1 keV to about 3 µm); see Figs. 12 and 13. This average SED is in excellent agreement with the mean SEDs of luminous QSOs, demonstrating that, for the analyzed sample, we do not find any redshift evolution of the shape of the SED, at least in the optical-UV region. As typically reported for luminous QSOs at lower-$z$, we find that optical emission at 3000 Å $< \lambda < 1$ µm is well described by a power law with a steeper slope ($F_{\nu} \propto \nu^{-0.14}$) than observed in the bulk of the low-luminosity AGN population; that is, $F_{\nu} \propto \nu^{-0.46}$ (Vanden Berk et al. 2001).

As a natural extension of this work, we aim to (i) increase the sample size of QSOs at $z > 6-7$ with good X-ray data and (ii) include rest-frame photometric points at $\lambda > 1$ µm thanks to JWST/MIRI observations in order to obtain a uniform NIR coverage. Improving these aspects will be fundamental for highlighting potential redshift-dependent properties, which we have so far been unable to constrain in a significant way. To this end, a recently accepted *XMM-Newton* Large Program (PI Zappacosta) will allow the investigation of X-ray emission also in sources with low $M_{s, Edd}$. In the NIR, it will instead be critical to increase the number of sources with photometric coverage to confirm the tentative indication from the computed mean SED that, on average, EoR sources exhibit enhanced emission from hot dust compared to average templates. JWST/MIRI observations will play a crucial role in adequately sampling the SED in the 1 to 2.5 µm region, providing a systematic characterization of hot-dust emission in these early quasars.

## Data availability

All the photometric data used in this work are available on Zenodo at https://zenodo.org/records/14181275. The same repository also contains the derived bolometric and monochromatic luminosities, as well as the spectral slopes, for both the HYPERION and E-XQR-30 samples.

*Acknowledgements.* We are grateful to the anonymous referee for their useful comments and suggestions which helped us to improve the paper. The authors acknowledge financial support from the Bando Ricerca Fondamentale
INAF 2022 Large Grant “Toward an holistic view of the Titans: multi-band observations of z > 6 QSOs powered by greedy supermassive black holes”. AB, MB, MB, SC, VD, FF, CF, SG, VT, NM, LZ acknowledge support from the European Union – Next Generation EU, PRIN/MUR 2022 2022TKPB2P – BIG-z. MB acknowledges support from INAF project 1.05.12.04.01 – MINI-GRANTS di RSN1 “Mini-feedback” and from UniTs under FVG LR 2/2011 project D55-microgrants23 “Hyper-gal”. DD acknowledges PON R&I 2021, CUP E65F21002880003. FT acknowledges funding from the European Union – Next Generation EU, PRIN/MUR 2022 2022K9N5B4. GM acknowledges financial support by grant PID2020-115325GB-C31 funded by MICIN/AEI/10.13039/501100011033. This work was supported by STFC grant ST/X001075/1.

## References

Abbott, T. M. C., Adamów, M., Aguena, M., et al. 2021, ApJS, 255, 20

Arnaud, K. A. 1996, ASP Conf. Ser., 101, 17

Avni, Y. 1976, ApJ, 210, 642

Bañados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11

Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473

Barnett, R., Warren, S. J., Banerji, M., et al. 2015, A&A, 575, A31

Bischetti, M., Feruglio, C., D’Odorico, V., et al. 2012, Nature, 605, 244

Bongiorno, A., Merloni, A., Brusa, M., et al. 2017, MNRAS, 427, 3103

Bosman, S. E. I., Álvarez-Márquez, J., Colina, L., et al. 2024, Nat. Astron., 8, 1054

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, https://doi.org/10.5281/zenodo.7419741

Cai, Z.-Y., & Wang, J.-X. 2023, Nat. Astron., 7, 1506

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, The IRSA 2MASS All-Sky Point Source Catalog

Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products

Dadina, M. 2008, A&A, 485, 417

D’Odorico, V., Bañados, E., Becker, G. D., et al. 2023, MNRAS, 523, 1399

Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73

Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1

Fan, X., Hennawi, J. F., Richards, G. T., et al. 2004, AJ, 128, 515

Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373

Farina, E. P., Schindler, J.-T., Walter, F., et al. 2022, ApJ, 941, 106

Feruglio, C., Maio, U., Tripodi, R., et al. 2023, ArXiv e-prints [arXiv:2304.09129]

Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143

Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, J. Open Source Softw., 4, 1864

Gallagher, S. C., Richards, G. T., Lacy, M., et al. 2007, ApJ, 661, 30

Haardt, F., & Maraschi, L. 1991, ApJ, 380, L51

Hopkins, P. F., Strauss, M. A., Hall, P. B., et al. 2004, AJ, 128, 1112

Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27

Iwamura, F., Kimura, M., Eto, S., et al. 2004, ApJ, 614, 69

Jiang, L., Fan, X., Hines, D. C., et al. 2006, AJ, 132, 2127

Jiang, L., Fan, X., Brandt, W. N., et al. 2010, Nature, 464, 380

Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222

Jun, H. D., & Im, M. 2013, ApJ, 779, 104

Just, D. W., Brandt, W. N., Shemmer, O., et al. 2007, ApJ, 665, 1004

Kaiser, N., Burgett, W., Chambers, K., et al. 2010, SPIE Conf. Ser., 7733, 77330E

Krawczyk, C. M., Richards, G. T., Mehta, S. S., et al. 2013, ApJS, 206, 4

Krawczyk, C. M., Richards, G. T., Gallagher, S. C., et al. 2015, AJ, 149, 203

Kubota, A., & Done, C. 2018, MNRAS, 480, 1247

Lawrence, A., Warren, S. J., Almaini, O., et al. 2012, VizieR Online Data Catalog: II/314

Leipski, C., Meisenheimer, K., Walter, F., et al. 2014, ApJ, 785, 154

Lusso, E., & Risaliti, G. 2016, ApJ, 819, 154

Lusso, E., Comastri, A., Vignali, C., et al. 2010, A&A, 512, A34

Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623

Lusso, E., Worseck, G., Hennawi, J. F., et al. 2015, MNRAS, 449, 4204

Lyu, J., Rieke, G. H., & Shi, Y. 2017, ApJ, 835, 257

Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38

Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169

Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, ApJ, 849, 91

Mazzucchelli, C., Bischetti, M., D’Odorico, V., et al. 2023, A&A, 676, A71

McMahon, R. G., Banerji, M., Gonzalez, E., et al. 2021, VizieR Online Data Catalog: II/367

Meyer, R. A., Bosman, S. E. I., & Ellis, R. S. 2019, MNRAS, 487, 3305

Miller, J. A., Cackett, E. M., Goad, M. R., et al. 2023, ApJ, 953, 137

Mountrichas, G., Buat, V., Georgantopoulos, I., et al. 2021, A&A, 653, A70

Petrucci, P. O., Ursini, F., De Rosa, A., et al. 2018, A&A, 611, A59

Pezzulli, E., Valiante, R., & Schneider, R. 2016, MNRAS, 458, 3047

Piconcelli, E., Jimenez-Bailón, E., Guainazzi, M., et al. 2005, A&A, 432, 15

Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389

Reed, S. L., Banerji, M., Becker, G. D., et al. 2019, MNRAS, 487, 1874

Richards, G. T., Hall, P. B., Vanden Berk, D. E., et al. 2003, AJ, 126, 1131

Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470

Ross, N. P., & Cross, N. J. G. 2020, MNRAS, 494, 789

Runnoe, J. C., Brotherton, M. S., & Shang, Z. 2012, MNRAS, 422, 478

Saccheo, I., Bongiorno, A., Piconcelli, E., et al. 2023, A&A, 671, A34

Schindler, J.-T., Farina, E. P., Bañados, E., et al. 2020, ApJ, 905, 51

Schlafly, E. F., Meisner, A. M., & Green, G. M. 2019, ApJS, 240, 30

Selsing, J., Fynbo, J. P. U., Christensen, L., & Krogager, J. K. 2016, A&A, 585, A87

Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45

Shen, Y., Wu, J., Jiang, L., et al. 2019, ApJ, 873, 35

Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252

Telfer, R. C., Zheng, W., Kriss, G. A., & Davidsen, A. F. 2002, ApJ, 565, 773

Temple, M. J., Hewett, P. C., & Banerji, M. 2021a, MNRAS, 508, 737

Temple, M. J., Banerji, M., Hewett, P. C., Rankine, A. L., & Richards, G. T. 2021b, MNRAS, 501, 3061

Timlin, J. D., Brandt, W. N., Qi, Q., et al. 2020, MNRAS, 492, 719

Tortosa, A., Zappacosta, L., Piconcelli, E., et al. 2024, A&A, 691, A235

Trakhtenbrot, B., & Netzer, H. 2012, MNRAS, 427, 3081

Tripodi, R., Feruglio, C., Kemper, F., et al. 2023, ApJ, 946, L45

Tripodi, R., Scholtz, J., Maiolino, R., et al. 2024, A&A, 682, A54

Valiante, R., Schneider, R., Volonteri, M., & Omukai, K. 2016, MNRAS, 457, 3356

Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549

Vietri, G., Piconcelli, E., Bischetti, M., et al. 2018, A&A, 617, A81

Vignali, C., Brandt, W. N., & Schneider, D. P. 2003, AJ, 125, 433

Vignali, C., Brandt, W. N., Schneider, D. P., & Kaspi, S. 2005, AJ, 129, 2519

Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Vito, F., Brandt, W. N., Bauer, F. E., et al. 2019, A&A, 630, A118

Volonteri, M. 2010, A&ARV, 18, 279

Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148

Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nat. Rev. Phys., 3, 732

Wang, F., Yang, J., Fan, X., et al. 2018, ApJ, 869, L9

Wang, F., Fan, X., Yang, J., et al. 2021, ApJ, 908, 53

Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512

Yang, J., Wang, F., Fan, X., et al. 2021, ApJ, 923, 262

Zappacosta, L., Comastri, A., Civano, F., et al. 2018, ApJ, 854, 33

Zappacosta, L., Piconcelli, E., Giustini, M., et al. 2020, A&A, 635, L5

Zappacosta, L., Piconcelli, E., Fiore, F., et al. 2023, A&A, 678, A201

Zheng, W., Kriss, G. A., Telfer, R. C., Grimes, J. P., & Davidsen, A. F. 1997, ApJ, 475, 469

<footnote n="1">¹ Dipartimento di Matematica e Fisica, Università Roma Tre, Via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="2">² INAF – Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monte Porzio Catone, Italy</footnote>

<footnote n="3">³ INAF – Osservatorio Astronomico di Trieste, Via G. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="4">⁴ Dipartimento di Fisica, Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="5">⁵ Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="6">⁶ IFPU – Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34151 Trieste, Italy</footnote>

<footnote n="7">⁷ Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="8">⁸ Dipartimento di Fisica e Astronomia ‘Augusto Righi’, Università degli Studi di Bologna, Via P. Gobetti, 93/2, 40129 Bologna, Italy</footnote>

<footnote n="9">⁹ INAF – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti, 93/3, I-40129 Bologna, Italy</footnote>

<footnote n="10">¹⁰ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="11">¹¹ INFN – National Institute for Nuclear Physics, Via Valerio 2, I-34127 Trieste, Italy</footnote>

<footnote n="12">¹² Department of Physics, University of Napoli ‘Federico II’, Via Cinthia 9, 80126 Napoli, Italy</footnote>
13 Millennium Institute of Astrophysics (MAS), Nuncio Monseñor Sotero Sanz 100, Providencia, Santiago, Chile

14 INAF – Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

15 Center for Astrophysics – Harvard & Smithsonian, Cambridge, MA 02138, USA

16 Steward Observatory, University of Arizona, Tucson, Arizona, USA

17 INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy

18 European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

19 DiSAT, Università degli Studi dell'Insubria, Via Valleggio 11, I-22100 Como, Italy

20 INFN – Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

21 INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807 Merate, Italy

22 Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., Cambridge CB3 0HE, UK

23 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

24 Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK

25 Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, ESAC Campus, 28692 Villanueva de la Cañada, Spain

26 ASI – Agenzia Spaziale Italiana, Via del Politecnico snc, I-00133 Roma, Italy

27 Dipartimento di Fisica, Università di Roma La Sapienza, Piazzale Aldo Moro 2, I-00185 Roma, Italy

28 INFN – Sezione Roma1, Dipartimento di Fisica, Università di Roma La Sapienza, Piazzale Aldo Moro 2, I-00185 Roma, Italy

29 Sapienza School for Advanced Studies, Viale Regina Elena 291, I-00161 Roma, Italy

30 Physics Department, Tor Vergata University of Rome, Via della Ricerca Scientifica 1, 00133 Rome, Italy

31 INFN – Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy

32 Department of Astronomy, University of Maryland, College Park, MD 20742, USA

33 University of Ljubljana, Department of Mathematics and Physics, Jadranska ulica 19, SI-1000 Ljubljana, Slovenia

34 INAF – Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, Via A. Corti 12, 20133 Milano, Italy

35 Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 Bis bd Arago, 75014 Paris, France
<table><caption>**Table A.1.** Diary of NIR imaging observations for the HYPERION sample.</caption><thead><tr><th>Name</th><th>Instrument</th><th>Filters</th><th>Exp. time</th><th>Date</th></tr></thead><tbody><tr><td>J1342+0928</td><td>NICS</td><td>J, H, K</td><td>62, 60, 123</td><td>25/05</td></tr><tr><td>J1120+0641</td><td>NICS</td><td>J, H, K</td><td>7, 15, 10</td><td>18/04</td></tr><tr><td rowspan="2">J0038-1527</td><td>SOFI</td><td>J, H, K</td><td>12, 6, 6</td><td>03/11</td></tr><tr><td>NICS</td><td>J, H, K</td><td>6, 7, 9</td><td>02/11</td></tr><tr><td rowspan="2">J0252-0503</td><td>SOFI</td><td>J, H, K</td><td>7, 10, 15</td><td>03/11</td></tr><tr><td>NICS</td><td>Y, J, H, K</td><td>10, 15, 18, 22</td><td>02/11</td></tr><tr><td>J0020-36</td><td>SOFI</td><td>J, H, K</td><td>9, 15, 15</td><td>03/11</td></tr><tr><td rowspan="2">J0411-0907</td><td>SOFI</td><td>J, H, K</td><td>18, 21, 30</td><td>03/11</td></tr><tr><td>NICS</td><td>Y, J, H, K</td><td>6, 11, 11, 12</td><td>02/11</td></tr><tr><td>J0244-5008</td><td>SOFI</td><td>J, H, K</td><td>21, 10, 10</td><td>03/11</td></tr><tr><td>J231-20.8</td><td>NICS</td><td>Y, J</td><td>10, 23</td><td>24/04</td></tr><tr><td>J036+03.0</td><td>NICS</td><td>Y, J, H, K</td><td>5, 15, 20, 31</td><td>16/09</td></tr><tr><td>J0224-4711</td><td>SOFI</td><td>J, H, K</td><td>6, 6, 6</td><td>03/11</td></tr><tr><td rowspan="3">J011+09</td><td>NICS</td><td>Y, J, H</td><td>44, 50, 40</td><td>04/09</td></tr><tr><td>NICS</td><td>K</td><td>60</td><td>16/09</td></tr><tr><td>NICS</td><td>Y, J, H, K</td><td>21, 44, 50, 60</td><td>02/11</td></tr><tr><td>J083+11.8</td><td>SOFI</td><td>J, H, K</td><td>18, 21, 27</td><td>03/11</td></tr><tr><td rowspan="2">J0050+3445</td><td>NICS</td><td>Y, J, H, K</td><td>15, 20, 30, 30</td><td>16/09</td></tr><tr><td>NICS</td><td>Y, J, H, K</td><td>6, 10, 11, 16</td><td>02/11</td></tr><tr><td>J029-36</td><td>SOFI</td><td>J, H, K</td><td>18, 18, 18</td><td>03/11</td></tr></tbody></table>

**Notes.** Exposure times are in minutes. All observations were conducted in 2022.

both TNG and NTT, we found that SOFI images provided more accurate (in terms of S/N ratio) but still compatible data and thus we report the SOFI measurements. In 3 cases (J1120+06, J231-20 and J036+03) observations were significantly affected by poor sky conditions (i.e. seeing > 1.7 arcsec); we discarded those observations.

## Appendix A: Proprietary NIR observations and data reduction

New NIR (Y, J, H and K bands) observations have been obtained for 15 sources with NICS at the INAF funded telescope TNG during programmes AOT45 and AOT46 (PI L. Zappacosta and F. La Franca respectively) and with SOFI at ESO/NTT, run 110.244M.001, (PI I. Saccheo). Details about the observations and the net exposure time per band are reported in Table A.1. Data were collected adopting a dithering strategy and acquiring at least 5 images per filter and were reduced with the following steps: we first applied a median stacking technique with a clipping algorithm to produce an empty sky-frame. This frame was then subtracted from all the raw frames to obtain a set of images with a median background level equal to zero. This step also removes bias and dark current levels. We then divided the sky-subtracted frames by the normalized flat-field image, obtained by twilight sky frames (for TNG) or already provided by the telescope staff (SpecialDomeFlat, see SOFI User manual). The flat-fielded images were then astrometrically registered and co-added obtaining the final frames for each filter. We performed photometry extraction using the python package photutils (Bradley et al. 2022) and utilized an 8 pixel radius (e.g., ~ 2 arcsec) to perform aperture photometry since we deal with point-like sources in not crowded fields. Errors on fluxes were computed as is done in IRAF package DAOPHOT, i.e. $\sigma_f^2 = f/epadu + A_f\sigma_{bkg}^2(1 + A_f/A_{bkg})$ where f is the flux measured in the aperture area $A_f$, $\sigma_{bkg}$ is the standard deviation of the background computed in a circular annulus of area $A_{bkg}$ and epadu is the telescope gain expressed in electrons per count. Instrument magnitudes were then calibrated by using nearby standard stars observed immediately before and after the QSO observations (J, H and K bands for NTT and Y band at TNG) or by using stars in the field also observed in the 2MASS survey (J, H and K bands at TNG Cutri et al. 2003). For sources observed at