Open Access
Issue
A&A
Volume 689, September 2024
Article Number A220
Number of page(s) 29
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202349054
Published online 13 September 2024

© The Authors 2024

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

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

1. Introduction

In the past few decades, many advances have been made in unveiling the properties of the host galaxies of the first quasi-stellar-objects (QSOs) from theory and observations. The Atacama Large Millimeter/sub-millimeter Array (ALMA), along with the Northern Extended Millimeter Array (NOEMA), the Very Large Array (VLA), and Herschel, have probed the cold gas and dust of the QSO host galaxies. The dust continuum was detected in many z ∼ 6 QSOs with estimated far-infrared (FIR) luminosities of LFIR = 1011−13 L and dust masses of about Mdust = 107−9 M (Decarli et al. 2018; Carniani et al. 2019; Shao et al. 2019).

The rest-frame FIR continuum emission originates from dust heated by the ultraviolet (UV) radiation from young and massive stars (Decarli et al. 2018; Venemans et al. 2020; Neeleman et al. 2021) and the active galactic nucleus (AGN) radiation (Schneider et al. 2015; Di Mascia et al. 2021; Walter et al. 2022). The latter contribution is usually neglected when modeling the FIR spectral energy distribution (SED) of z ∼ 6 QSOs, although the AGN heating can contribute 30 − 70% of the FIR luminosity (Schneider et al. 2015; Duras et al. 2017). Moreover, dust masses are often determined with huge uncertainties and only rely on single-frequency continuum detections. However, if multifrequency ALMA and/or NOEMA observations are available in the rest-frame FIR, which probe both the peak and the Rayleigh Jeans tail of the dust SED, the dust temperature and mass can be constrained with statistical uncertainties < 10 − 20% (e.g., Carniani et al. 2019; Tripodi et al. 2022, 2023a; Witstok et al. 2023). This results in a high accuracy in the determination of the star formation rate (SFR).

Accurate estimates of the dust masses would also allow us to derive the molecular gas mass of the host galaxy through the gas-to-dust ratio (GDR). Although it is possible to directly probe the molecular reservoirs of the QSO host galaxies using the rotational transitions of the CO (e.g., Vallini et al. 2018; Madden et al. 2020), very few high-z QSOs are observed in CO because this emission line is typically faint at high z. The GDR has indeed often been assumed in order to compute the gas mass, resulting in a high degree of uncertainty in its estimate. Studies of z ∼ 2.4 − 4.7 hyperluminous QSOs showed that the GDR spans a broad range of values, 100−300, with an average GDR ∼180 (Bischetti et al. 2021). This is consistent with the results found for submillimeter (submm) galaxies to z ∼ 3 − 5 with a GDR ∼ 150 − 250 (e.g., Saintonge et al. 2013; Miettinen et al. 2017). Recent studies suggested that the GDR for SMGs may be lower than 100 (Birkin et al. 2021; Liao et al. 2024). In low-z galaxies, a GDR ∼ 100 is typically observed (Draine et al. 2007; Leroy et al. 2011), implying that the GDR increases with redshift. However, the combination of gas-mass estimates from CO with the dust mass from finely sampled FIR SEDs results in an accurate determination of the GDR, and this allows us to study its evolution with redshift.

Another central question is to understand the efficiency of large halos in forming stars during the first cosmic billion year. In this context, luminous QSOs are efficient signposts of large halos at high redshift (e.g., Wang et al. 2024), which allows us to investigate the complex physics at play in the building of massive galaxies. The efficiency of a galaxy in forming stars is an intricate mechanism that depends on the efficiency with which cold gas is formed from baryons and on the efficiency with which the cold gas is converted into stars. The former is influenced by the capability of gas cooling and the capability of gas removal or gas heating by feedback. The gas star formation efficiency (SFE) has been investigated for decades in both the local and high-redshift universe. It is parameterized by the Kennicutt–Schmidt (KS) relation (see e.g., Kennicutt 1998; Genzel et al. 2010; Speagle et al. 2014; Calabrò et al. 2024). Because spatially resolved observations and of robust estimates of the SFR are still lacking, the KS relation is still poorly studied in high-z QSOs. Many questions about the SFE in these objects are therefore still unanswered.

Targeting the QSO host galaxies at these redshifts provides a unique opportunity to characterize the formation and concurrent build-up of SMBHs and their host galaxies, and also the physical properties of the ISM in the first billion year of the Universe (e.g., Decarli et al. 2018; Venemans et al. 2020; Neeleman et al. 2021). The local correlation between massive BHs and bulges is indeed tight, suggesting that the same process that assembled galaxy bulges may cause most of the growth of massive BHs. QSOs at z ≳ 6 appear to lie above the local MBH − M* correlation, and thus, the BH growth seems to precede that of its host galaxy (Ding et al. 2023; Stone et al. 2024; Yue et al. 2024). However, before the advent of the James Webb Space Telescope (JWST), it was difficult to obtain reliable and accurate estimates of the stellar mass in QSO hosts at high z, even with deep HST observations (Mechtley et al. 2012; Marshall et al. 2020). Therefore, instead of the stellar mass, the dynamical mass is used at high z, exploring the MBH − Mdyn relation. For instance, Feruglio et al. (2018) and Pensabene et al. (2020) studied the MBH − Mdyn relation in luminous high-z QSO hosts and reported that this relation evolves with redshift and that high-z QSOs lie above the local relations. This implies that the SMBHs formed significantly faster than their hosts. In this context, Izumi et al. (2018, 2019) studied the MBH − Mdyn using a sample of seven z > 6 low-luminosity quasars. They found that while the luminous quasars typically lie above the local relation (Kormendy & Ho 2013) with BHs overmassive compared to local AGNs, the discrepancy becomes less evident for low-luminosity quasars (Maiolino et al. 2023). In order to avoid being biased by the properties of higher luminosity sources, it is essential to study the whole population of QSOs and galaxies at high z, including the low-luminosity sources. This is now possible with data from JWST (see e.g. Santini et al. 2023; Harikane et al. 2023), which enable us to determine the stellar mass especially in low-luminosity sources based on a sharper PSF (especially in the short-wavelength bands) and longer-wavelength coverage out to 5 micron compared to HST, which is less affected by dust attenuation. Recently, Harikane et al. (2023) analyzed a first statistical sample of faint type 1 AGNs at z > 4 identified by JWST/NIRSpec deep spectroscopy. They compared their AGNs to the AGNs at z ∼ 0 (Reines & Volonteri 2015) and to high-z QSOs. Their AGNs have similar BH masses, but systematically lower stellar masses than the local AGNs. Moreover, Maiolino et al. (2023, 2024) showed even more extreme MBH − M* ratios in their sample of 12 AGNs at z = 4 − 7 that belong to the JADES survey. Similar results were also obtained in previous studies with a smaller number of AGNs (Übler et al. 2023). This indicates that the BH grows faster than its host galaxy at high redshift. The fast BH growth was also suggested by previous studies at z ∼ 2 (e.g., Zhang et al. 2023a). These overmassive (compared to their host stellar masses) BHs are indeed predicted in some theoretical models simulations (e.g., Toyouchi et al. 2021; Trinca et al. 2022; Inayoshi et al. 2022; Hu et al. 2022; Zhang et al. 2023b; Pacucci et al. 2023; Pacucci & Loeb 2024). In summary, recent observations showed that high-z QSOs may lie above the local MBH − M* correlation, and they are thus likely to follow the BH-dominance growth or BH-dominance evolutionary path (i.e., the green line in Fig. 3 of Volonteri 2012).

The significance of our work becomes evident in the context of these prior studies, as we aim to explore the interconnected growth of supermassive black holes and their host galaxies. Additionally, we seek to understand whether the high-redshift hosts of quasars serve as the progenitors of the massive galaxies observed in the local Universe. These two queries are pivotal to the study of galaxy evolution and are tightly interconnected. They can only be addressed through a reliable and precise understanding of the properties of QSOs and their massive hosts.

In this work, we present the analysis of the cold-dust SEDs of four QSOs at z > 6, PSO J036.5078+03.0498 (hereafter J036+03, Venemans et al. 2015) at z = 6.5405, VDESJ022426.54−471129.4 (hereafter J0224−4711, Reed et al. 2017) at z = 6.5222, PSO J231.6576−20.8335 (hereafter J231−20, Mazzucchelli et al. 2017) at z = 6.587, and SDSS J205406.49−000514.8 (hereafter J2054−0005, Reed et al. 2017) at z = 6.0391, based on new ALMA band 8 (B8) and/or band 9 (B9) observations, archival ALMA observations from band 3 (B3) to band 6 (B6), and a NOEMA observation at 3 mm. These observations allowed us to retrieve reliable and accurate estimates of the dust parameters and the SFR of the QSO host galaxies in our sample. Furthermore, from the analysis of archival observations, we were also able to detect the CO(7−6) and [CI] emission lines for J0224−4711 and the CO(6−5) emission line for J2054−0005, from which we were able to estimate the gas mass for these objects.

In order to investigate the dust properties, SFR, GDR, and gas SFE in a statistically sound sample of QSOs at z > 6, we considered another six high-z QSOs whose dust properties and SFR were already derived with high accuracy. The sample is presented in the following section.

The paper is organized as follows. In Sect. 2, we describe our sample of QSOs at z > 6. In Sect. 3, we describe the observations used in this work. In Sect. 4, we present the results for the continuum and line emissions. In Sect. 5, we analyze the cold-dust SED of the QSOs in our sample, and we derive the gas mass for QSOs J0224−4711, J1319+0950, and J2054−0005. In Sect. 6, we contextualize our findings regarding the dust and gas in high-z QSO, and we discuss the observed scenario for the evolutionary paths of these objects. In Sect. 7, we summarize the content of this work.

Throughout the paper, we adopt the ΛCDM cosmology from Planck Collaboration VI (2020): H0 = 67.4 km s−1 Mpc−1, Ωm = 0.315, and ΩΛ = 0.685. Thus, the angular scale is 5.66 kpc/arcsec at z = 6.3.

2. Sample

We selected from the currently known ∼300 QSOs at z > 6 (see, e.g., Fan et al. 2023; Jiang et al. 2016; Mortlock et al. 2011; Wang et al. 2019a; Inayoshi et al. 2020) all the QSOs at z > 6 for which we were able to investigate and compare the evolutionary path of the SMBHs and their host galaxies, and we linked it to feedback processes. This implied that we should be able derive or retrieve accurate estimates of the SMBH properties, dust properties, SFR and gas mass, and to detect outflowing emission. In other words, all these QSOs have high-quality rest-frame optical-UV spectra and ALMA and/or NOEMA observations targeting the continuum emission from low to high frequencies (in the range of observed frequency ∼100−600 GHz), that is, probing both the Rayleigh-Jeans and peak region of the cold-dust SED, and targeting the CO and [CII] line emissions (see, e.g., Zappacosta et al. 2023; D’Odorico et al. 2023; Wang et al. 2013; Neeleman et al. 2019; Witstok et al. 2023; Decarli et al. 2022, and other refs throughout this work). Moreover, we selected QSOs with Lbol > 1047 erg s−1. Consequently, our findings are not directly applicable to the entire population of AGNs at z > 6; instead, they specifically pertain to high-luminosity sources. Nevertheless, this luminosity bias aligns with our objectives, as we are particularly interested in investigating the role of AGN feedback in influencing the SMBH-galaxy evolution, and notably, high-luminosity QSOs exhibit compelling evidence of powerful outflows (see, e.g., Bischetti et al. 2022; Shao et al. 2022; Tripodi et al. 2024; Salak et al. 2024).

Our final sample consisted of ten QSOs at z > 6. Together with the four QSOs presented in the previous section (J036+03, J0224−4711, J231−20, and J2054−0005), we considered six other QSOs: SDSS J010013.02+280225.8 (hereafter J0100+2802), SDSS J231038.88+185519.7 (hereafter J2310+1855), J1319+0950, ULAS J134208.10+092838.35 (hereafter J1342+0928), SDSS J114816.64+525150.3 (hereafter J1148+5251), and PSO J183.1124+05.0926 (hereafter J183+05).

We divided the sample into two subsamples of six and four QSOs each, hereafter called HYPERION QSOs and z > 6 QSOs, respectively (see Table 1).

Table 1.

General properties of the QSOs in our sample.

The first subsample is composed of the six QSOs belonging to the sample called HYPerluminous QSOs at the Epoch of ReionizatION (HYPERION). HYPERION comprises the titans of z > 6 QSOs Zappacosta et al. (2023), that is, those whose SMBH had the fastest mass growth history. In particular, these QSOs were selected so that the SMBHs powering them required to form at least a 1000 M BH seed (at z = 20) under the hypothesis of a continuous exponential accretion at the Eddington rate. These SMBHs likely assembled from the largest BH seeds, or alternatively, they experienced peculiar, possibly super-Eddington, mass accretion histories. Of the ∼300 QSOs known at the epoch of reionization (EoR), the HYPERION QSOs comprise 18 QSOs with a mean redshift z ∼ 6.7, and average log(Lbol/erg/s)∼47.3, and a BH mass in the range 109 − 1010 M. HYPERION is based on a 2.4 Ms XMM-Newton Multi-Year Heritage Programme (PI: Zappacosta) to provide for the first time for such a large sample of z > 6 QSOs a uniform high quality X-ray spectral characterization for a detailed investigation of the nuclear properties and the accretion and ejection processes that are tied to the fast build-up experienced by their SMBH. The first results suggest a genuine redshift evolution of their X-ray spectral slopes, which appear to be steeper than reported in z < 6 QSOs with a similar luminosity and accretion rate. This supports a different regime for the X-ray nuclear properties of the first quasars that might be linked to the presence of fast disk-driven winds (Zappacosta et al. 2023). While the nuclear properties of HYPERION QSOs are well constrained, as are the dynamical masses of most QSOs based on archival [CII] observations, their dust properties and SFR are still mostly unconstrained.

The z > 6 QSO subsample comprises the QSOs that did not satisfy the criteria for belonging to the HYPERION sample (see Sect. 2 of Zappacosta et al. 2023), that is, they did not experience a rapid SMBH mass growth.

We divided our final sample based on the HYPERION survey selection criteria in order to assess whether the BH accretion history has strong implications for the properties of the host galaxy.

Seven QSOs in our sample1 are also part of the XQR-30 ESO Large Program (D’Odorico et al. 2023). Therefore, they have high S/N X-shooter spectra, from which Mazzucchelli et al. (2023) derived BH masses using the C IV and Mg II emission lines. The Mg II BH masses derived for J0100+2802, J036+03, J0224−4711, and J231−20 agree well with those in Zappacosta et al. (2023).

The properties of the ten QSOs in our final sample are summarized in Table 1, and a more detailed presentation of each object can be found in Appendix A.

3. Observations

All ALMA observations targeting QSOs J036+03, J0224−4711, J231−20, J1319+0950, J2054−0005 were calibrated and imaged as outlined in the following.

The visibility calibration of the observations was executed by the ALMA science archive. The imaging was performed through the Common Astronomy Software Applications (CASA; McMullin et al. 2007), version 5.1.1-5. We applied tclean using natural weighting and a 3σ cleaning threshold. We imaged the continuum using the multifrequency synthesis (MFS) mode in all line-free channels. To image the line emissions, we used the CASA task uvcontsub to fit the continuum visibilities in the line-free channels with a first-order polynomial for QSO J0224−4711, since the continuum showed a non-negligible slope, and zeroth-order polynomial for the other QSOs. We then obtained continuum-subtracted cubes to be used in our analyses. We produced line maps using the MFS mode in the channels enclosing the emission lines (see Table C.2).

We analyzed all archival observations available for J036+03 and J0224−4711. We also used the results on the continuum emission at ∼107 GHz of J036+03, obtained from a NOEMA observation (project ID: S17CD).

Because many ALMA observations are available for J2054−0005, especially in B6, we considered for each band those whose continuum sensitivity and resolution were suited for our analysis of the SED. That is, when multiple observations per band were available, we considered the observation with the highest continuum sensitivity and with a resolution that allowed us to spatially resolve the source and/or that was as close as possible to the angular resolution of the B9 observation, in order to ensure a reliable and consistent analysis of the cold-dust SED (for more details, see Sect. 4.1).

For J231−20, we used the results presented in Pensabene et al. (2021), who analyzed all the ALMA observations from B3 to B6 available for this QSO, and we independently analyzed a new B8 ALMA observation targeting this QSO (ID: 2021.2.00064.S). Since J231−20 was found to have a close companion (at a distance < 10 kpc), we discuss the analysis of the B8 observation in more detail in Sect. 4.1.

Information about the project ID, synthesized beam, lines detected, line channels, and root mean square (rms) noise of the continuum map and of the continuum-subtracted cube for each observation analyzed in this work is reported in Tables C.1 and C.2.

4. Data analysis

In the following, we report the results for our proprietary observations targeting the continuum (Sect. 4.1) and/or line emissions (Sect. 4.2) in ALMA bands 3, 8, and/or 9 of QSOs J036+03, J0224−4711, J1319+0950, and J2054−0005. As an exception, we discuss the archival observation targeting the continuum emission in B8 of J231−20 because it was found to have a close companion in lower-frequency observations. Details of the analysis of archival observations of the continuum emissions of QSOs J036+03, J0224−4711, and J2054−0005 can be found in Appendix B. The results obtained from the analysis of archival observations are reported in Table C.1.

4.1. Continuum emission

4.1.1. QSO J036+03

The top left and right panels of Fig. 1 show the continuum emission at 404.9 GHz and 670.9 GHz of J036+03, respectively. Neither emission is spatially resolved, and therefore, we considered the peak flux as the total flux of the source, that is, 6.63 ± 0.39 mJy/beam at 404.9 GHz and 5.60 ± 0.69 mJy/beam at 670.9 GHz. The emission at 670.9 GHz presents a secondary peak at RA, Dec = [ − 1.5, −1.5], which likely is an artifact due to the low resolution of the B9 observation. This interpretation is supported by the absence of line and continuum emission at the same spatial position in all lower-frequency observations and by the presence of another dimmer peak that is located symmetrically with respect to the QSO position (RA, Dec = [1.0, −1.5]).

thumbnail Fig. 1.

Dust continuum maps. Top left panel: 404.9 GHz dust continuum map of QSO J036+03 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.6 mJy/beam). The clean beam (2.10 × 1.43 arcsec2, PA = 37.40°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Top right panel: 670.9 GHz dust continuum map of QSO J036+03 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.5 mJy/beam). The clean beam (1.12 × 0.97 arcsec2, PA = 21.11°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Bottom left panel: 405.2 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.88 mJy/beam). The clean beam (3.87 × 2.79 arcsec2, PA = −76.45°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Bottom right panel: 670.9 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, 8, and 11σ, σ = 1.3 mJy/beam). The clean beam (1.08 × 0.95 arcsec2, PA = 14.05°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak.

4.1.2. QSO J0224−4711

The continuum emission of J0224−4711 at 405.2 GHz and 670.9 GHz is shown in the bottom left and right panels of Fig. 1, respectively. Neither emission is spatially resolved, and therefore, we consider the peak flux as total flux of the source, which is 8.73 ± 0.38 mJy/beam at 405.2 GHz and 19.9 ± 0.96 mJy/beam at 670.9 GHz.

4.1.3. QSO J2054−0005

We show the continuum map at 92.3 GHz underlying the CO(6−5) emission line in the bottom row of Fig. 2. Performing a 2D Gaussian fit, we found that the peak flux is 0.066 ± 0.004 mJy/beam, the integrated flux is 0.082  ±  0.009 mJy, and the emission is spatially resolved with a size of (0.27 ± 0.07)  ×  (0.07 ± 0.09) arcsec2.

thumbnail Fig. 2.

Dust continuum and CO(6−5) emission line maps of QSOs J1319+0950 and J2054−0005. Top left panel: 103 GHz dust continuum map of QSO J1319+0950 (levels −3, −2, 2, 3, 5, 10, and 20σ, σ = 0.005 mJy/beam). The clean beam (0.30 × 0.30 arcsec2, PA = −78.56°) is indicated in the lower left corner of the diagram. Top right panel: CO(6−5) map of QSO J1319+0950 (levels −3, −2, 2, 3, 5, 10, and 15σ, σ = 0.017 mJy/beam). The clean beam (0.32 × 0.31 arcsec2, PA = −78.56°) is indicated in the lower left corner of the diagram. Bottom left panel: 92 GHz dust continuum map of QSO J2054−0005 (levels −3, −2, 2, 3, 5, 7 and 10σ, σ = 0.006 mJy/beam). The clean beam (0.42 × 0.32 arcsec2, PA = −61.3°) is indicated in the lower left corner of the diagram. Bottom right panel: CO(6−5) map of QSO J2054−0005 (levels −3, −2, 2, 3, 5, 7 and 10σ, σ = 0.025 mJy/beam). The clean beam (0.39 × 0.30 arcsec2, PA = −61.4°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak for each source.

4.1.4. QSO J231–20

J231−20 was found to have a close companion that is detected in multiple emission lines and in the continuum emission from ∼100 GHz up to ∼250 GHz at a distance of ∼2 arcsec from the QSO (Pensabene et al. 2021; Decarli et al. 2017; Neeleman et al. 2019). Pensabene et al. (2021) were able to disentangle the continuum emission of the QSO from that of the companion using ∼1 arcsec resolution ALMA observations. Unfortunately, the low resolution of the new band 8 observation (beam of 4.3 × 2.9 arcsec2) did not allow us to distinguish the two emissions (see Fig. B.1). Performing a 2D fit with CASA on the continuum map, we found a peak flux of 8.43 ± 0.39 mJy/beam, and the source was unresolved. The B8 flux is indeed contaminated by the continuum emission of the companion, and this can therefore bias the SED fitting and the determination of the dust properties. We corrected for the companion contribution using the results found by Pensabene et al. (2021). In their SED modeling, the flux of the companion at ∼400 GHz is more than 0.6 dex lower than that of the QSO. Conservatively considering the highest-temperature model for the companion, that is, that the flux of the companion is 0.6 dex lower than that of the QSO, we can correct our flux in band 8. We obtain 6.74 ± 0.31. We considered this flux as a detection in the SED fitting, but it might also be seen as a lower limit to the flux of the QSO.

4.2. Emission lines

4.2.1. CO(7–6) and [CI] emission lines in J0224–4711

We used the continuum-subtracted data cube in B3 to study the CO(7−6) and [CI] emission lines of J0224−4711. The central and right panels of Fig. 3 present the maps of CO(7−6) and [CI], respectively, imaged using the MFS mode in the channels specified in Table C.2. The CO(7−6) and [CI] emission are detected with a statistical significance of ∼8σ and ∼3σ, respectively, and neither is spatially resolved. Performing a 2D Gaussian fit, we obtained a peak flux of 0.36 ± 0.02 mJy/beam for CO(7−6) and of 0.19 ± 0.03 mJy/beam for [CI].

thumbnail Fig. 3.

Dust continuum and emission line maps of QSO J0224−4711. Left panel: 95.33 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, and 5σ, σ = 0.016 mJy/beam). The clean beam (3.95 × 2.32 arcsec2, PA = −89.43°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Central panel: CO(7−6) emission line map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.042 mJy/beam). The clean beam (3.49 × 2.05 arcsec2, PA = −89.06°) is indicated in the lower left corner of the diagram. The cross indicates the position of the B3 continuum peak. Right panel: [CI] emission line map of QSO J0224−4711 (levels −3, −2, 1, 2, and 3σ, σ = 0.054 mJy/beam). The clean beam (3.49 × 2.05 arcsec2, PA = −89.06°) is indicated in the lower left corner of the diagram. The cross indicates the position of the B3 continuum peak.

Fig. 4 shows the moment-0, −1, and −2 maps of the CO(7−6) emission, obtained by applying a 2.5σ threshold to the continuum-subtracted cube in the line channels, and the spectrum extracted from the region with S/N > 2 in the CO(7−6) map. The moment-0 map shows a velocity gradient oriented east to west with Δv = 100 km s−1, and the moment-2 map shows a range of the velocity dispersion2 between 20 and 120 km s−1. The CO(7−6) line profile peaks at a frequency of 107.239 GHz, corresponding to z = 6.5220 ± 0.0002, consistent with previous determinations (Reed et al. 2017) and with the redshift derived from MgII emission line (D’Odorico et al. 2023). From the fit with a single Gaussian, the FWHM of the line is 307 ± 25 km s−1 and the integrated flux is 0.34 ± 0.05 Jy km s−1, which corresponds to a luminosity of LCO(7 − 6) = (1.6 ± 0.2)×108L and L CO ( 7 6 ) = ( 9.6 ± 1.4 ) × 10 9 $ L^\prime_\mathrm{CO(7{-}6)} = (9.6\pm 1.4)\times 10^9 $ K km s−1 pc2 (following Eqs. 1 and 3 of Solomon & Vanden Bout 2005). The [CI] line is slightly blueshifted, peaking at 107.610 GHz, which corresponds to z = 6.5211 ± 0.0004. Performing a single-Gaussian fit, we obtained that the FWHM of the line is 165 ± 29 km s−1 and the integrated flux is 0.108 ± 0.036 Jy km s−1, which corresponds to a luminosity3 of L[CI] = (5.1 ± 1.7)×107 L and L [ CI ] = ( 3.0 ± 1.0 ) × 10 9 $ L^\prime_{\mathrm{[CI]}} = (3.0\pm 1.0)\times 10^9 $ K km s−1 pc2. All line properties are reported in Table 2.

thumbnail Fig. 4.

Moment maps of the CO(7−6) emission line and spectrum of CO(7−6) and [CI] emission lines of J0224−4711. From left to right: Integrated flux, mean velocity map, velocity dispersion map, and continuum-subtracted spectra of CO(7−6) and [CI]. The clean beam is plotted in the lower left corner of the moment maps. The cross indicates the peak position of the Band 3 continuum emission. The spectrum was extracted from the region included within ≥2σ in the CO(7−6) map.

Table 2.

Line properties of QSOs J0224−4711, J1319+0950, and J2054−0005.

4.2.2. CO(6–5) emission line in J1319+0950

In order to analyze the CO(6−5) emission line of J1319+0950, we rebinned the continuum-subtracted data cube in B3 to 50 km s−1-pagination.

The top row of Fig. 2 presents the map of the CO(6−5) emission line of J1319+0950 and its underlying continuum, imaged using the MFS mode in the channels specified in Table C.2. The CO(6−5) is detected with a statistical significance of ∼15σ, and performing a 2D Gaussian fit, we obtained a peak flux of 0.303 ± 0.014 mJy/beam and an integrated flux of 0.570 ± 0.038 mJy. The emission is spatially resolved with a size of (0.35 ± 0.03)×(0.23 ± 0.03) arcsec2.

The top row of Fig. 5 shows the moment-0, −1, and −2 maps of the CO(6−5) emission of J1319+0950, obtained by applying 3σ threshold to the continuum-subtracted cube in the line channels, and the spectrum extracted from the region with S/N > 2 in the CO(6−5) map. The moment-0 shows a velocity gradient oriented northeast to southwest with Δv = 400 km s−1, and the moment-2 map shows a range of the velocity dispersion between 20 and 160 km s−1. The CO(6−5) line profile peaks at a frequency of 96.939 GHz, corresponding to z = 6.1331 ± 0.0004. Performing a single-Gaussian fit to the line spectral profile, we obtained an FWHM of the line of 529 ± 41 km s−1 and an integrated flux of 0.662 ± 0.094 Jy km s−1, which corresponds to a luminosity of LCO(6 − 5) = (2.4 ± 0.3)×108 L and L CO ( 6 5 ) = ( 2.3 ± 0.3 ) × 10 10 $ L^\prime_\mathrm{CO(6{-}5)} = (2.3 \pm 0.3)\times 10^{10} $ K km s−1 pc2.

thumbnail Fig. 5.

Moment maps of the CO(6−5) emission lines and spectra of J1319+0950 (top row), J2054−0005 (bottom row). From left to right: integrated flux, mean velocity map, and velocity dispersion map, continuum-subtracted spectrum of CO(6−5). The clean beam is plotted in the lower left corner of the moment maps. The cross indicates the peak position of CO(6−5). The spectra were extracted from the region included within ≥2σ in the CO(6−5) map of each source.

4.2.3. CO(6–5) emission line in J2054–0005

In order to analyze the CO(6−5) emission line of J2054−0005, we rebinned to 50 km s−1 the continuum-subtracted data cube in B3.

The bottom row of Fig. 2 presents the map of the CO(6−5) emission line of J2054−0005 and its underlying continuum, imaged using the MFS mode in the channels specified in Table C.2. The CO(6−5) is detected with a statistical significance of ∼10σ. Performing a 2D Gaussian fit, we obtained a peak flux of 0.349 ± 0.036 mJy/beam, an integrated flux of 0.496 ± 0.080 mJy, and the emission is spatially resolved with a size of (0.30 ± 0.09)×(0.08 ± 0.16) arcsec2.

The bottom row of Fig. 5 shows the moment-0, −1, and −2 maps of the CO(6−5) emission of J2054−0005, obtained by applying 3σ threshold to the continuum-subtracted cube in the line channels, and the spectrum extracted from the region with S/N > 2 in the CO(6−5) map. The moment-0 shows a velocity gradient oriented southeast to northwest with Δv = 120 km s−1, and the moment-2 map shows a range of the velocity dispersion between 20 and 80 km s−1. The CO(6−5) line profile peaks at a frequency of 98.233 GHz, corresponding to z = 6.0391 ± 0.0002. Performing a single-Gaussian fit to the line spectral profile, we obtained a FWHM of the line of 229 ± 20 km s−1 and an integrated flux of 0.288 ± 0.047 Jy km s−1, which corresponds to a luminosity of LCO(6 − 5) = (1.0 ± 0.2)×108L and L CO ( 6 5 ) = ( 9.8 ± 1.6 ) × 10 9 $ L^\prime_\mathrm{CO(6{-}5)} = (9.8\pm 1.6)\times 10^9 $ K km s−1 pc2.

5. Analysis

5.1. Dust properties and star formation rate

In order to ensure a self-consistent analysis of the cold-dust SEDs of J036+03, J0224−4711, J231−20 and J2054−0005, we performed an SED fitting of the flux densities reported in Table C.1 considering the tapered fluxes for the higher-resolution observations if needed (see Appendix B for details of the tapering procedure). The observations of J231−20 in Pensabene et al. (2021) did not need additional tapering because their resolution was well matched and low enough to account for the fainter and more extended emission. Moreover, for J231−20, we considered the flux corrected for the contribution of the companion to the QSO emission as explained in Sect. 4.

We modeled the dust continuum in an optically thick regime with a modified black body (MBB) function given by

S ν obs obs = S ν / ( 1 + z ) obs = Ω ( 1 + z ) 3 [ B ν ( T dust ( z ) ) B ν ( T CMB ( z ) ) ] ( 1 e τ ν ) , $$ \begin{aligned} {S\!}_{\nu _{\rm obs}}^\mathrm{obs} = {S\!}_{\nu /(1+z)}^\mathrm{obs} = \frac{\Omega }{(1+z)^3}[B_{\nu }(T_{\rm dust}(z))-B_{\nu }(T_{\rm CMB}(z))](1-e^{-\tau _{\nu }}), \end{aligned} $$(1)

where Ω = ( 1 + z ) 4 A gal D L 2 $ \Omega = (1+z)^4A_{\mathrm{gal}}D_{\mathrm{L}}^{-2} $ is the solid angle as a function of the surface area of the galaxy, Agal, and of the luminosity distance of the galaxy, DL. The dust optical depth is

τ ν = M dust A gal k ν = M dust A gal k 0 ( ν ν 0 ) β , $$ \begin{aligned} \tau _{\nu }= \frac{M_{\rm dust}}{A_{\rm gal}}k_\nu = \frac{M_{\rm dust}}{A_{\rm gal}}k_0\biggl (\frac{\nu }{\nu _0}\biggr )^{\beta }, \end{aligned} $$(2)

where β is the emissivity index, kν is the opacity, k0 = 0.45 cm2 g−1 is the mass absorption coefficient, and ν0 = 250 GHz. The latter two terms define the opacity model adopted here (Beelen et al. 2006; Carniani et al. 2019). In this model, kν carries huge systematics because the actual composition of the dust is currently unknown. Variations in the opacity model can, in principle, lower dust masses by a factor of ∼3 or increase them by a factor of ∼1.5. The solid angle for each source was estimated using the continuum mean size from resolved observations (see Table C.1). The effect of the CMB on the dust temperature is given by

T dust ( z ) = ( ( T dust ) 4 + β + T 0 4 + β [ ( 1 + z ) 4 + β 1 ] ) 1 4 + β , $$ \begin{aligned} T_{\rm dust}(z)=((T_{\rm dust})^{4+\beta }+T_0^{4+\beta }[(1+z)^{4+\beta }-1])^{\frac{1}{4+\beta }}, \end{aligned} $$(3)

with T0 = 2.73 K. We also considered the contribution of the CMB emission given by Bν(TCMB(z) = T0(1 + z)) (da Cunha et al. 2013).

Therefore, the MBB model has three fitting parameters: the dust temperature (Tdust), the dust mass (Mdust), and β. We explored the 3D parameter space using a Markov chain Monte Carlo (MCMC) algorithm implemented in the EMCEE package (Foreman-Mackey et al. 2013). We assumed uniform priors for the fitting parameters: 10 K < Tdust < 300 K, 105 M < Mdust < 109 M, and 1.0 < β < 3.0. The best-fitting Tdust, Mdust, and β, obtained from a MCMC with 40 chains, 6000 trials and a burn-in phase of ∼100 for each QSO, are reported in Table C.3. The errors on the best-fitting parameters were computed considering the 16th and 84th percentiles of the posterior distribution of each parameter. Overall, we found Tdust = 50 − 80 K, Mdust ∼ 108 M, and β = 1.6 − 2.5. Fig. 6 shows the observed SEDs with the best-fitting function and the posterior distributions of Tdust, Mdust, and β for J036+03, J0224−4711, and J231−20, and Fig. 7 shows this for J2054−0005.

thumbnail Fig. 6.

Observed SEDs and best-fitting results. Left panels: Observed SEDs of QSOs J036+03, J0224-4711, and J231−20 (top, central, and bottom row). Our new ALMA B8 and B9 data are shown as yellow stars, and other archival observations are plotted as cyan diamonds. The best-fitting curve is shown as the solid blue line. Right panels: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid orange lines indicate the best-fitting value for each parameter, and the dashed lines mark the 16th and 84th percentiles for each parameter. For J231−20, the empty diamond is the flux in B8, not corrected for the presence of the companion, and the dashed line is the best-fitting SED considering this flux (for details, see Sect. 5.1).

thumbnail Fig. 7.

Same as Fig. 6 for QSO J2054−0005.

We estimated the total infrared (TIR) luminosity for the best-fit model for each source by integrating from 8 to 1000 μm rest-frame, and we derived the SFR considering that SFR ∝ LTIR. The proportionality factor between SFR and LTIR depends on the initial mass function (IMF) chosen. We adopted a Kroupa IMF (Kroupa & Weidner 2003). The values of LTIR and SFR derived for our QSO sample are reported in Table C.3. An Salpeter or Chabrier IMF (Salpeter 1955; Chabrier 2003) would imply an SFR that is higher by a factor of 1.16 or lower by factor of 0.67, respectively. When comparing with results of the SFR from the literature, we rescaled the SFR according to the Kroupa IMF (see also Kennicutt & Evans 2012):

SFR ( M yr 1 ) = 1.496 × 10 10 L IR 8 1000 μ m ( L ) . $$ \begin{aligned} \mathrm{SFR}\ (\mathrm{M_\odot \,yr^{-1}}) = 1.496 \times 10^{-10}\ L_{\rm IR\ 8{-}1000\,\upmu \mathrm{m}}\ (\mathrm{L}_\odot ). \end{aligned} $$(4)

As a word of caution, we recall that we adopted a lower limit for the flux in B8 of J231−20, and thus, Tdust and SFR can be underestimated in principle. By performing another fit considering the B8 flux not corrected for the contribution of the companion (which is an upper limit to the flux of the QSO), we indeed derived Tdust = 61 K and SFR = 1913 M yr−1, which is a factor of 2 higher than that obtained considering the lower limit in B8, while the dust mass and emissivity index agree well with the previous values. The best-fitting curve is displayed as a dashed line in the bottom left panel of Fig. 6, along with the uncorrected B8 flux as an empty diamond. Hereafter, we conservatively use the results of the SED fitting with the B8 flux corrected for the companion emission, considering that the SFR can vary of a factor of 2 at most.

We adopted a contribution of the AGN emission to the dust heating of 50% considering that all our sources are hyperluminous QSOs at z ≳ 6, likely sharing similar properties with the hyperluminous QSOs at z = 2 − 4 in Duras et al. (2017) and with J1148+5251 in Schneider et al. (2015), which also belongs to the HYPERION sample. This assumption yields SFR = 466 ± 75 M yr−1 for J036+03, an SFR = 2485 768 + 1682 M yr 1 $ 2485^{+1682}_{-768}\,\mathrm{M_{\odot}\,yr^{-1}} $ for J0224−4711, an SFR = 496 ± 118 for J231−20, and an SFR = 730 ± 75 for J2054−0005. The role of AGN in heating the dust is further discussed in Sect. 6.1.6.

We wished to perform a statistically sound investigation of the dust properties in a sample of z > 6 QSOs and therefore explored the archives and the literature in our search for other suitable candidates for our analysis of cold-dust SED, that is, with ALMA and/or NOEMA observations of the continuum emission probing both the Rayleigh-Jeans regime and the peak region of the SED. To our knowledge, only four other QSOs at z > 6 currently have multiple-frequency observations sampling the Rayleigh-Jeans part and, barely, the peak of the cold-dust SED: QSOs J1319+0950, J1148+5251, J1342+0928, and J183+05. Their SEDs were analyzed in Carniani et al. (2019) (J1319+0950 and J1148+5251), in Novak et al. (2019), Witstok et al. (2023) (J1342+0928), and in Decarli et al. (2023) (J183+05). Since we adopted the same method as Carniani et al. (2019), we used the results for the dust properties and SFRs of J1319+0950 and J1148+5251 found in Carniani’s work, which are reported here in Table C.3. For consistency, we modeled the observed SEDs of the other two QSOs using the procedure described in the previous paragraph because the method and/or the opacity models adopted in Novak et al. (2019), Witstok et al. (2023), Decarli et al. (2023) were different from ours. The best-fitting SEDs are shown in Fig. 8, and the results can be found in Table C.3. Within the uncertainties, our results agree well overall with those in the previous papers. The discrepancies in the best-fitting values arise from the different regime assumed (optically thin in Novak et al. 2019, and thick in Witstok et al. 2023; Decarli et al. 2023) and in the different opacity models.

thumbnail Fig. 8.

Observed SEDs and best-fitting results. Left panels: Observed SEDs of QSOs J183+05, J1342+0928 (top and bottom row) as cyan diamonds. The best-fitting curve is shown as a solid blue line. Right panels: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid orange lines indicate the best-fitting value for each parameter, and the dashed lines mark the 16th and 84th percentiles for each parameter.

It is worth noting that the uncertainties on Tdust (SFR) increase to ∼40% (> 80% for SFR) when the peak of the cold-dust SED is not probed: This is the case of J1319+0950 and J1342+0928, which lack ALMA observations in bands 8 and 9, and of J0224−4711, for which even band 9 is unable to reach the peak of the SED given that this QSO has an extremely bright dust emission. In the latter case, an ALMA B10 observation is required. The strongest effect is seen in J1342+0928, for which Tdust is determined with an uncertainty of ∼40%, leaving the SFR basically unconstrained. To investigate the impact of B8–B9 observations on the uncertainties of Tdust, we again fit the SEDs for all the QSOs in our sample, excluding the data available at λrest > 100 μm for each source, that is, corresponding to ALMA B8 and B9. We found that in this case, the uncertainties on Tdust become ∼70% on average (in some cases, up to 90%). This highlights the importance of high-frequency observations (ALMA bands 8, 9, and 10) in providing precise and reliable estimates of the dust properties and SFR.

In Fig. 9 we show the cold-dust SEDs analyzed in this work. All SEDs were normalized at λ = 850 μm rest-frame in order to facilitate the comparison. It is worth noting the difference in slope and position of the peak, which are clearly related to the variation in the dust properties, namely Tdust, Mdust, and β.

thumbnail Fig. 9.

Compilation of cold-dust SEDs of the six QSOs analyzed in this work, normalized at 850 μm.

Additionally, we developed an observed mean cold-dust SED by averaging (mean) the observational data available for all the QSOs in our sample in six wavelength bins4. They are shown as dashed gray lines in Fig. 10, along with the resulting averaged fluxes as violet squares. The error bars are the standard errors on the mean associated with each flux. We fit this mean SED following the procedure described above, and we obtained T dust = 54 4 + 8 $ T_{\mathrm{dust}}=54^{+8}_{-4} $ K, Mdust = (4.6 ± 3.0)×108 M, and β = 1.6 ± 0.55. Fig. 10 shows the best-fitting function and the posterior distribution of each fitting parameter.

thumbnail Fig. 10.

Results for the observed mean cold-dust SED. Left panel: Observed mean fluxes and best-fitting function. The dashed gray lines mark the six bins we used to develop the mean SED. Right panel: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid green lines indicate the best-fitting value for each parameter, while the dashed lines mark the 16th and 84th percentiles for each parameter.

5.2. Molecular gas mass

In this section, we infer the molecular gas mass of QSOs J0224−4711, J1319+0950, and J2054−0005 from molecular tracers such as the CO(7−6) and CO(6−5) emission lines. The molecular gas mass estimates in high-z QSO host galaxies rely on CO observations, and especially on intermediate (Jup = 5−7) transitions (Venemans et al. 2017; Yang et al. 2019; Decarli et al. 2022), which are found to be at the peak of the CO spectral line energy distribution (CO SLED, Li et al. 2020; Feruglio et al. 2023). In principle, low-J transitions should be preferred, as they are less sensitive to uncertainties on the CO excitation. However, these are quite challenging to detect because of their intrinsic faintness.

The molecular gas mass was therefore derived as

M H 2 , CO = α CO r J ( J 1 ) 1 L CO ( J ( J 1 ) ) , $$ \begin{aligned} M_{\rm H2, CO} = \alpha _{\rm CO} r^{-1}_{J-(J-1)} L_{\mathrm{CO}(J-(J-1))}, \end{aligned} $$(5)

where αCO = 0.8 M (K km s−1 pc2)−1 is the CO-to-H2 conversion factor typical for ULIRG and QSOs (Carilli & Walter 2013), and rJ − (J − 1) is the CO(J − (J − 1))/CO(1−0) luminosity ratio. Fig. 11 shows the CO SLEDs normalized to the Jup = 6 transition of J036+03 (Decarli et al. 2022) and J2054−4711, limited to the (6−5) and (7−6) transitions (for CO(7−6) of J2054−0005 see Decarli et al. 2022), compared with other z > 6 QSOs, and QSO APM0879+5255 at z = 3.9 and QSO Cloverleaf at z = 2.56 (Li et al. 2020; Feruglio et al. 2023), for which the CO SLED is well constrained down to Jup = 1. We note that the CO SLED for QSO at z > 6 presents a flattening at the CO(6−5) and (7−6) transitions on average, and that all CO SLEDs have similar slopes from CO(6−5) to CO(2−1) transitions, while at higher J, the CO SLED are very different.

thumbnail Fig. 11.

CO SLED of J036+03 and J2054−0005 compared with those of other QSOs at z > 6 and at lower redshift. The CO SLEDs for J036+03 and J2054−0005 are shown as pink circles and cyan markers, respectively. For J2054−0005, the star refers to the CO(6−5) line studied in this work and the circle to the CO(7−6) in Decarli et al. (2022). The results for J036+03 are taken from Decarli et al. (2022). All other sources are displayed as shadowed squares: J1007+2115 at z = 7.5419 in red (Feruglio et al. 2023); J0439+1634 at z = 6.511 in gray (Yang et al. 2019); J1148+5251 at z = 6.419 in brown (Riechers et al. 2009; Gallerani et al. 2014); J0100+2802 at z = 6.327 in purple (Wang et al. 2019b); J2310+1855 at z = 6.003 in blue (Li et al. 2020); APM08279+5255 at z = 3.911 in green (Papadopoulos et al. 2001; Weiß et al. 2007); and Cloverleaf at z = 2.560 in orange (Bradford et al. 2009; Uzgil et al. 2016).

So far, APM0879+5255 and Cloverleaf are the only objects for which the CO(1−0) transition is detected and the ratio of the transition CO(2−1) and (1−0) is consistent between the two sources. The upper limit on the CO(1−0) of J1148+5251 also indicates a similar CO(2−1)-to-CO(1−0) ratio. Therefore, in order to estimate the molecular gas mass for J1319+0950, and J2054−0005 from CO(6−5), we assumed r65 = CO(6−5)/CO(1−0) = 1.23 ± 0.66, considering the average between the maximum and minimum values of r65 in APM0879+5255 and Cloverleaf. We then obtained MH2 = (1.5 ± 0.2)×1010 for J1319+0950, and MH2 = (6.4 ± 1.0)×109 M for J2054−0005. Regarding J0224−4711, we assumed that CO(7−6) ∼ CO(6−5), given the averaged CO SLED of high-z QSO (see Fig. 11), and r76 = 0.76 ± 0.41, considering the average between APM0879+5255 and Cloverleaf. This gave MH2 = (1.0 ± 0.1)×1010 M for J0224−4711.

For the purposes of this work, we also investigated the properties of QSO J036+03, and thus, we computed its molecular gas mass from CO(6−5) consistently with the rest of the sample. Decarli et al. (2023) reported the detection of the CO(6−5) emission line of J036+03, finding L CO ( 6 5 ) = ( 1.27 ± 0.11 ) × 10 10 $ L^\prime_\mathrm{CO(6{-}5)}=(1.27\pm 0.11)\times 10^{10} $ K km s−1 pc2. This yielded MH2 = (7.5 ± 0.64)×109 M following our assumptions.

Our analysis revealed that the gas masses in our sample of QSOs, estimated with the smallest statistical uncertainties, are ∼1010 M on average. However, it is worth noting that our assumptions for the CO(6−5)-to-CO(1−0) (or CO(7−6)-to-CO(1−0)) line ratio introduce significant systematic uncertainties to our estimates of ≳50%. Our results agree well on average with the gas masses found in Neeleman et al. (2021) derived from converting the [CII] underlying continuum flux into a dust mass and then assuming a constant GDR (=100). In contrast, the gas mass estimated from converting the [CII] luminosity directly into a molecular mass using the conversion of Zanella et al. (2018) mostly is higher by one order of magnitude than ours. Specifically, molecular gas masses for QSOs J1319+0950 and J2054−0005 were derived in Neeleman et al. (2021), who reported higher values than ours using the GDR and the Zanella conversion factor. In the former case, the GDR assumed in Neeleman et al. (2021) (GDR = 100) may not be valid for every QSO (see the results and discussion of the GDR in our sample in Sect. 6.1.5). In the latter, the Zanella conversion factor, which is calibrated for z ∼ 2 main-sequence galaxies, may not be applicable for high-z QSO hosts. This was also reported in Tripodi et al. (2022) for the case of QSO J2310+1855 at z ∼ 6.

The molecular gas masses and line properties are summarized in Table 2. The uncertainties reported on the gas mass are only statistical. Systematics uncertainties are induced by the choice of αCO (a factor of 0.2−0.3 dex) and by the scaling from high-J CO transition to J = 1 using the CO excitation ladder (a factor of 20−30%).

6. Discussion

In this section, we first discuss our findings concerning the dust and gas within our high-z QSO sample by conducting a comparative analysis with the properties observed in lower-z sources. Finally, we investigate the evolutionary paths of the ten QSOs in our sample, exploiting the results for the properties of their host galaxies derived in Sects. 5.1 and 5.2.

6.1. Properties of the QSO host galaxies

Our final sample comprises ten QSOs at 6 ≲ z < 7.5, out of which six belong to the HYPERION sample. In the previous sections, we analyzed the cold-dust SED in a homogeneous way for all quasars, and the results are summarized in Table C.3. In the following, we briefly describe the compilation of galaxy and AGN host samples we used for the comparison.

For the local Universe, we considered the JCMT dust and gas the In Nearby Galaxies Legacy Exploration (JINGLE) survey, the Herschel (U)LIRG Survey (HERUS), and a sample of QSOs selected from the Palomar-Green (PG) survey (as also done in Witstok et al. 2023). The JINGLE sample is composed of 192 nearby (0.01 < z < 0.05) galaxies, for which Lamperti et al. (2019) studied their cold-dust SED using photometric data in the 22−850 μm range from Herschel, applying a hierarchical Bayesian fitting approach. We used the cold-dust properties derived from their two MBB models, which also take into account the warm dust component in the optically thin regime. These results agreed well with those from a single MBB model. The HERUS sample comprises 43 z < 0.3 (ultra)-luminous infrared galaxies or (U)LIRGs observed with the Infrared Astronomical Satellite (IRAS) and Herschel (Sanders et al. 2003; Clements et al. 2018). Clements et al. (2018) derived the dust properties for this sample assuming an MBB function in the optically thin regime. Witstok et al. (2023) recalculated the dust properties of the JINGLE and HERUS datasets, employing a method analogous to ours (see Sect. 5.1 and Witstok et al. 2022), except for the variation in the chosen opacity model. Their results agreed with those presented in Lamperti et al. (2019) and Clements et al. (2018) within the uncertainties. Therefore, for the sake of simplicity, we used the results from the original papers. The dust properties of the PG sample, consisting of 85 nearby (z < 0.5) QSOs, were obtained by Petric et al. (2015) through modeling the photometry taken by Herschel by means of an MBB function with β fixed to 2.0 and the dust model of Draine et al. (2007)5.

At higher redshift, we selected a sample of 81 2 < z < 7 strongly gravitationally lensed, dusty star-forming galaxies identified by the South Pole Telescope (SPT). Reuter et al. (2020) analyzed the cold-dust SEDs of the objects in this sample employing an MBB function in the optically thick regime with β fixed to 2. We also included the dust properties inferred for 16 QSOs belonging to the WISE-SDSS Selected Hyper-luminous (WISSH) sample at 2 < z < 5 (Duras et al. 2017). Their SEDs were analyzed accounting for the cold dust, dusty torus, and warm dust emission components. In particular, the cold-dust emission was modeled by an MBB function in the optically thin regime with β = 1.6. Finally, at 4 < z < 7.5, we considered a sample composed of three submm galaxies (SMGs) and eight star-forming (SF) galaxies from Witstok et al. (2023). As stated before, they adopted a model for the SED fitting analogous to ours. For the QSO host galaxy samples, we note that the PG QSOs have a lower bolometric luminosity Lbol compared to our z ≳ 6 sample (Lbol, PG ∼ 1044 − 47 erg s−1), whereas the WISSH quasars are the most luminous, with Lbol, WISSH ≳ 1047.5 erg s−1 (see, e.g., Vietri et al. 2018).

6.1.1. Dust temperature

In the past few years, the trend of Tdust as a function of time has been the object of multiple studies (e.g., Faisst et al. 2017; Schreiber et al. 2018; Sommovigo et al. 2020; Bakx et al. 2021; Viero et al. 2022). However, the behavior of the Tdust − z relation, which increases or plateaus at z > 4, has remained unclear. The top panel of Fig. 12 presents the redshift distribution of the dust temperatures in our sample compared to the low-z samples described above. Additionally, we display the average Tdust derived for REBELS6 galaxies by Sommovigo et al. (2022) as a light blue triangle, and Tdust found for three REBELS galaxies at z > 6.7 by Algera et al. (2023) as light blue diamonds. As a word of caution, the dust temperatures for REBELS galaxies derived in Sommovigo et al. (2022) carry large uncertainties since they are based on a combination of a single photometric point in band 6 and the [CII] line emission. Hence, in this context, we exclusively present the mean value. Algera et al. (2023) analyzed band 8 and 6 observations of REBELS-127, REBELS-25, and REBELS-38 considering both optically thin and thick cases. We present the results for the optically thick case, fixing β = 1.58. Because we relied on only two photometric points per galaxy, the inferred Tdust are still uncertain (ΔTdust = 30 − 60%); the results obtained in the thin and thick regime agree within the large uncertainties. For clarity, we show the mean values of the temperature distribution for the JINGLE and HERUS samples with their corresponding standard deviation9. Interestingly, there is no significant difference in Tdust between QSOs and normal SF galaxies at fixed redshift, also considering that our sample is biased toward high luminosity and therefore possibly higher dust temperatures.

thumbnail Fig. 12.

Dust temperature as a function of redshift. Our sample is shown as stars (red for Hyperion QSOs, green otherwise). Top panel: Comparison of our results with samples of local QSOs (PG, blue colormap with contours, Petric et al. 2015), local star-forming galaxies (JINGLE, mean value as pink dot, Lamperti et al. 2019), local ULIRGs (HERUS, mean value as orange dot, Clements et al. 2018), 2 < z < 7 star-forming galaxies (SPT, pink colormap with contours Reuter et al. 2020), 2 < z < 5 QSOs (WISSH, yellow squares, Duras et al. 2017), 4 < z < 7.5 SMGs and star-forming galaxies (light blue dots, Witstok et al. 2023), three z > 6.7 REBELS galaxies (light blue diamonds, Algera et al. 2023), an average-REBELS galaxy (light blue triangle, Sommovigo et al. 2022), and two galaxies with very low dust temperatures at z > 6 (black crosses, Witstok et al. 2022; Harikane et al. 2020). The observed trends inferred by Viero et al. (2022) and Schreiber et al. (2018) are shown as an orange and green lines and shaded regions, respectively. Our best-fitting curve is shown as a purple line with a shaded region. The theoretical relation found in Sommovigo et al. (2022) is shown as a dashed gray line. The solid black line marks the CMB temperature level. Bottom panel: Focus on the comparison to the samples of local PG QSOs (blue color map with contours) and 2 < z < 5 WISSH QSOs. Our best-fitting curve considering QSOs alone is shown as a purple line with a shaded region.

We observe an increasing trend of Tdust with redshift that is naturally expected from a theoretical perspective as a result of decreasing gas depletion times, as seen in Sommovigo et al. (2022). Sommovigo et al. theoretically derived that T dust t dep 1 / 6 ( 1 + z ) 5 / 2 ( 4 + β ) $ T_{\mathrm{dust}}\propto t_{\mathrm{dep}}^{-1/6}\propto (1+z)^{5/2(4+\beta)} $, where β is the dust emissivity index. This Tdust − z relation is shown in Fig. 12 (dash-dotted gray line) for β = 2.03 and assuming optical depth and metallicity (in solar unit) both equal unity, as presented in Sommovigo et al. (2022), implying Tdust ∝ (1 + z)0.42. This theoretical relation is able to reproduce the trend observed in many SF galaxies (e.g., considering REBELS and ALPINE galaxies), and in Schreiber et al. (2018) from the stacking of star-forming galaxies at 0.5 < z < 4 in the deep CANDELS fields. The stacked SEDs were fit with a library of template SEDs generated with β = 1.5. Schreiber et al. (2018) found a linear trend with redshift, and we also show an extrapolated curve at z > 4 as dashed green line.

On the other hand, the population of SPT SF galaxies, QSOs ,and SF with Tdust > 60 K shows a steeper increase in Tdust with z, which is well captured by the observed relation inferred by Viero et al. (2022). They employed stacked maps in the FIR/submm of 111227 galaxies at 0 < z < 10 from the COSMOS-MOS2020/FARMER catalog to derive dust temperatures at different redshifts. Their observed relation for Tdust(z), which is a second-order polynomial in z, agrees well with the observed trend of the WISSH QSOs, with the sources at z > 5 that shows the most extreme temperatures (∼70 − 80 K), and with the lower limit of Bakx et al. (2020) and estimate of Behrens et al. (2018) based on data from Laporte et al. (2017) (see Viero et al. 2022). However, it fails to describe the bulk of the PG QSOs and HERUS ULIRGs and the sources at z > 5 with temperatures below 60 K (Faisst et al. 2020; Béthermin et al. 2020; Hashimoto et al. 2019; Sommovigo et al. 2022, see Viero et al. 2022). As a word of caution, these latter low-temperature estimates are very uncertain because they were derived from ALMA observations that did not bracket the peak of the SED, and they are therefore not adequate for modeling hotter dust components.

In order to find a general relation that can be applied to SF galaxies, SMGs, and QSO host galaxies at high z, we fit the observed data (excluding the stacked results) adopting the parameterization of Tdust(z) that was theoretically found in Sommovigo et al. (2022), therefore using a power law of the form

T dust ( z ) = A × ( 1 + z ) B , $$ \begin{aligned} T_{\rm dust}(z) = A \times (1+z)^B, \end{aligned} $$(6)

where A and B are free fitting parameters. The latter is linked with β, since β = −4 + 5/2B (considering the derivation of Sommovigo et al. 2022). We found A = 19 ± 2 K and B = 0.7 ± 0.1, implying a value of β that is unphysical. This underlines that Tdust does not depend uniquely on redshift, especially when considering different galaxy populations. Sommovigo et al. (2022) pointed out that at fixed redshift, the scatter in Tdust derives from variations either in optical depth or in metallicity. This is a likely scenario for different galaxy populations. Our best-fitting function, shown as a purple line and shaded region in the top panel of Fig. 12, is slightly flatter than Viero’s at z > 1, but still agrees with it pretty well in both the low- and high-z regimes. It captures both the flattening and the increasing trends at z > 4, that is, it describes the populations with both higher and lower Tdust.

Taking advantage of our robust estimates for Tdust in our QSO sample, we also performed a fit only to the observed QSO dust temperatures at different redshift in order to investigate the Tdust − z relation in QSOs for the first time. We found A = 29 ± 2 K and B = 0.35 ± 0.04. This trend is flatter than Sommovigo’s, but still agrees well with observed SF galaxies. The majority of the sources are within the uncertainties, including those from Schreiber et al. (2018). Three sources with low Tdust belonging to the sample of Witstok et al. (2023) (SPT0311−58W, SPT0311−58E, and A1689−zD1) agree within 2σ with our relation, along with another three sources that are well known for their low dust temperatures (J0217−0208, COS−3018555981, and REBELS-25, see Witstok et al. 2022; Algera et al. 2023; Harikane et al. 2020). As a word of caution, the dust temperature estimated for J0217−0208 and COS−3018555981 has large uncertainties because the data available for the two sources are poor (one or two detections and an upper limit) and because of the assumption of β = 1.5. In particular, for COS−3018555981, we show the upper limit on Tdust derived in Witstok et al. (2022). Even if these two sources were excluded from our fit of the Tdust − z relation, our results would not change. As discussed in Sommovigo et al. (2022), the scatter seen in Tdust at fixed redshift can be explained by variations in column density and metallicity within sources for optically thin and thick galaxies, respectively. This best-fitting relation is shown as a purple line with shaded region in the bottom panel of Fig. 12.

6.1.2. Dust emissivity

We now focus on the investigation of the dust emissivity index, β, which is physically related to the microscopic properties of dust. Fig. 13 shows the redshift distribution of β for our sample and the comparison samples. As was also found in Witstok et al. (2023), β does not evolve with redshift, with a mean value of ∼1.6 ± 0.2 (see also Beelen et al. 2006), indicating that the average dust properties do not change drastically. Interestingly, other studies of SMGs found a value of β that was closer to 2 on average (e.g., da Cunha et al. 2021; Cooper et al. 2022; Bendo et al. 2023; McKay et al. 2023; Liao et al. 2024). There is no difference (on average) in this case either between the QSOs belonging to the HYPERION sample and the others. As a word of caution, we note that high-z samples are biased toward high luminosity, and this can introduce a bias toward lower values of β (see also Witstok et al. 2023). β is significantly higher than 2 in only two exceptions, all of which lie in the high-z regime: J0100+2802, J036+03 from our sample, and GN10 from Witstok et al. (2023).

thumbnail Fig. 13.

Dust emissivity index, β, as a function of redshift. The symbols and colors are the same as in Fig. 12. For the PG, WISSH, and SPT samples, the symbols are empty to indicate that the value of β in these cases was fixed and not derived from SED analyses.

It is not straightforward to assess the physical reasons of these high β value. In principle, β depends on the physical properties and chemical composition of the grains, and possibly on environment and temperature. There are cases in which β can be higher than 2 (see, e.g., Valiante et al. 2011; Galliano et al. 2018). Spatially resolved studies in low-z galaxies showed a spread of β within a single galaxy, probably due to the temperature mixing, the different properties of the grain populations, or a combination thereof (Galliano et al. 2018). Liao et al. (2024) discussed how to interpret β in terms of dust properties and suggested that the grain composition plays a significant role in determining β. Simulations conducted by Hirashita et al. (2014) indicated that a β value of 2 can result from emission by either graphite or silicate grains. However, the correlation between the exact grain composition and β is unclear. In our case, a physical explanation for the high value of β would require studies of the dust grain properties and/or a detailed analysis of the temperature mixing, which is beyond the scope of this work.

6.1.3. Dust mass and star formation ratio

In Fig. 14, we compare the SFR versus dust mass of our sample with results from the PG and WISSH sample of QSOs. In our sample, five QSOs are the first at z > 6 for which the SFR was derived with the smallest statistical uncertainty based on high-frequency observations. As already mentioned in Sect. 5.1, the uncertainty of ∼40% on Tdust of J1342+0928 strongly affects the determination of the SFR, which has a large uncertainty. Therefore, this QSO is marked with a green circle in the plot. A correlation between the dust content and the star formation activity is evident in all the three samples, which are at different redshifts, with some scatter. This correlation is thought to be a consequence of the Schmidt–Kennicutt relation, linking the SFR to the dust content through the GDR (see, e.g., Santini et al. 2014). Overall, the SFR is higher in the z > 2 samples (by ∼2 orders of magnitudes on average), supporting the well-known concept that high-z QSOs are hosted in highly star-forming galaxies. The dust masses are higher at high z by ≲2 orders of magnitudes on average as well.

thumbnail Fig. 14.

Star formation rate as a function of dust mass. The symbols and colors are the same as in Fig. 12. The SFRs for the WISSH and our sample are corrected for a factor of 50% to take the AGN contribution to the dust heating into account.

Observations of galaxies at low and high z (Dunne et al. 2011; Driver et al. 2018; Pozzi et al. 2020; Beeston et al. 2018) found that the amount of dust in galaxies has decreased by a factor of ∼2 − 3 during the last ∼8 Gyr, and this was recently reproduced by Parente et al. (2023) using the semi-analytic model (SAM) L-GALAXIES2020. Considering an average value, we also observe a mild increase in the dust mass from our HYPERION sample ( M dust , HYP mean = 1.8 × 10 8 M $ M_{\mathrm{dust, HYP}}^{\mathrm{mean}}=1.8 \times 10^8\,\mathrm{M}_\odot $) to the WISSH sample at lower-z ( M dust , WISSH mean = 3 × 10 8 M $ M_{\mathrm{dust, WISSH}}^{\mathrm{mean}}=3 \times 10^8\,\mathrm{M}_\odot $), and followed by a drop of 2 orders of magnitude at z ∼ 0. Interestingly, here, we observe a difference within our sample: Four HYPERION QSOs have lower dust mass than the other z > 6 QSOs, up to more than about one order of magnitude and a factor of 2 on average. However, this is a preliminary result because two-thirds of the sample still has to be analyzed, and we will design future ALMA and NOEMA observations to complete the overview of the sample.

Even though we focused on QSOs, we recall that the relation between SFR and dust mass has also been studied in normal galaxies (Santini et al. 2014; da Cunha et al. 2010; Hjorth et al. 2014) and in SMGs at z = 1 − 4 (e.g., Dudzevičiūtė et al. 2020). Santini et al. (2014) used Herschel observations to estimate the dust mass of a large sample of galaxies extracted from the GOOD-S, GOOD-N, and COSMOS fields, and they performed a stacked analysis on a grid of redshifts, stellar masses, and SFR. Similarly to us, they found correlations between SFR and dust mass at different redshifts, from the local Universe out to z = 2.5. Their analysis revealed no clear evolution of the dust mass with redshift at a given SFR and stellar mass, indicating that galaxies with similar properties (in terms of SFR and stellar mass) do not show a significant difference in terms of dust content across the cosmic epochs, out to z = 2.5.

It is challenging to assess whether or not there is an evolution of Mdust with z at fixed SFR for the samples of QSOs we considered because the statistics is still poor at high z, especially for the low Mdust regime. Moreover, we are not able to explore the relation involving the stellar mass in our sample because M* estimates are not yet available for our high-z sample. Future JWST campaigns devoted to the investigation of the stellar content in QSOs and galaxies at z ≳ 6 will certainly allow us to perform these studies in detail (as done in, e.g., Harikane et al. 2023; Santini et al. 2023).

6.1.4. Star formation efficiency

The efficiency with which gas is converted into stars may vary by some orders of magnitude depending on the QSO and galaxy properties, such as luminosity, mass, the presence of a disk, bulges or mergers, and on cosmic time. At fixed redshift, the observed strong variation in the gas star formation efficiency (SFE = SFR/MH2) has usually been attributed to the existence of two distinct star formation laws: a secular mode in main-sequence (MS) galaxies, and a star-bursting (SB) mode characterized by short gas depletion timescales (τdep = 1/SFE ≲ 100 Myr; e.g., Solomon & Vanden Bout 2005; Sargent et al. 2014; Daddi et al. 2010; Scoville et al. 2023; Genzel et al. 2010; Vallini et al. 2024). The Kennicutt–Schmidt (KS) relation between the ΣMH2 and ΣSFR (Kennicutt 1998) is the most common method for studying the SFE in galaxies. However, we lack spatially resolved observations for all the QSOs in our sample. As an alternative, we discuss the integrated KS relation, considering MH2 and SFR.

The wide range of SFEs, and thus τdep, that is spanned by different galaxies at different cosmic times is shown in Fig. 15. We compared our results (red and green stars) with the 2 < z < 5 WISSH QSO sample, 1 < z < 1.5 QSOs (Frias Castillo et al. 2024), a sample of narrow-line Seyfert 1 galaxies at z < 0.1 (Salomé et al. 2023), 23 PG QSOs at z < 0.1 (Shangguan et al. 2020a,b), 0.025 < z < 0.05 xCOLDGASS galaxies (Saintonge et al. 2011, 2017), 0.002 < z < 0.09 LIRGs and ULIRGs with CO(1-0) detections from the literature (Tacconi et al. 2018), and AGNs at z < 0.05 selected from the Swift-BAT all-sky catalog (Rosario et al. 2018; Koss et al. 2021). As a reference, we also plot the gas MS computed by Sargent et al. (2014) considering a sample of galaxies at 0 < z < 2 (solid black line). While discussing this comparison, we recall the possible sources of systematic uncertainties (see also Sects. 5.2, 6.1.6): the assumption of the αCO, the scaling from the CO ladder for the CO(J → J − 1) transition with J > 1, and the assumption of Tdust for the cold-dust SED fitting and of the AGN contribution to the dust heating (for AGN-QSOs only).

thumbnail Fig. 15.

SFR as a function of the molecular gas mass, MH2. The symbols for our sample are the same as in Fig. 12. We compare our results with the WISSH sample (Bischetti et al. 2021), with 1 < z < 1.5 QSOs (pink dots), and with z ∼ 0 sources such as galaxies, ULIRGs, QSOs, AGNs, and narrow-line Seyfert 1 galaxies (symbols and colors in the legend; Saintonge et al. 2011, 2017; Salomé et al. 2023; Shangguan et al. 2020a,b; Rosario et al. 2018; Koss et al. 2021). The gray lines represent fixed values of the gas depletion time (i.e., the inverse of SFE), which is reported at the top of the line. The solid black line is the galaxy main sequence derived by Sargent et al. (2014) considering massive galaxies (M* > 1010 M) up to z ∼ 2.

Our sample exhibits τdep ∼ 10−2 Gyr, which is smaller by 1 − 2 orders of magnitude than those of the comparison samples up to z ∼ 2. As pointed out in Sect. 6.1.1, this may be connect to the high dust temperatures found in our QSOs. Vallini et al. (2024) modeled the variation in the dust temperature in terms of gas depletion timescales, optical depth, and metallicity. Based on the values of Tdust in Table C.3, the implied τdep of the Vallini et al. (2024) models agrees with those found from the observed SFR and MH2 (see Table C.3).

Moreover, considering the high-z MS derived by Scoville et al. (2023), our QSO host galaxies are found to be star-bursting sources. The star-bursting nature of our sample is also supported by the fact that τdep is lower than that found considering the relation of Tacconi et al. (2020) for τdep(z) in SB galaxies extrapolated up to z ∼ 6 (see also Vallini et al. 2024). According to the results of Scoville et al. (2023), for SB galaxies such as ours, 70% of the increased SFR relative to the MS is due to the elevated SFEs and not to the increased gas masses at early epochs.

6.1.5. Gas-to-dust ratio

We computed the GDR for the QSOs in our sample using the dust and gas masses reported in Table C.3, which were either estimated in this work or were taken from previous works in the literature (all references are reported in the table). In Fig. 16, we present the redshift distribution of the GDR in our sample compared with the WISSH sample (Bischetti et al. 2021), a sample of 2 < z < 5 star-forming galaxies hosting a heavily obscured AGN in the Chandra Deep Field-South (magenta dots, D’Amato et al. 2020), and a sample of z > 5.5 QSOs (blue squares Calura et al. 2014). D’Amato et al. (2020) derived the dust mass by modeling the SEDs with an MBB in the optically thin regime assuming β = 2.0, and Calura et al. (2014) also adopted a MBB in the optically thin regime, but assuming both Tdust = 47 K and β = 1.6 because these two works both mostly relied on a single data point at 250 GHz. Therefore, these latter estimates for Mdust are highly uncertain and, for instance, Mdust (GDR) would be lower (higher) by ∼2.4 times if a lower temperature, Tdust = 33 K, were assumed. For the comparison samples in Fig. 16, the vertical black line represents an average uncertainty on the GDR of about 0.3 dex, which accounts for the uncertainties on the cold-dust SED modeling. For our sample, the dust masses are in constrast constrained with uncertainties smaller than 30% on average. Systematic uncertainties are also significant when the gas mass is derived from CO transitions higher than J = 1 − 0 because this requires assuming a scaling between the luminosity of the considered CO transition and of CO(1−0) (see Sect. 5.2). This correction depends on the CO excitation ladder, which may vary from one source to the next depending on the ISM conditions. Given the large systematics in the gas mass determination (0.2−0.3 dex for αCO and 20−30% for r65 or r76), we did not include error bars in the GDR plot. We stress, however, that the dust masses in our sample are derived with the smallest statistical uncertainties. This has never been achieved before in a sample of QSOs at z ≳ 6. In low-z galaxies, a GDR ∼ 100 is typically observed (e.g., Draine et al. 2007; Leroy et al. 2011), while studies of massive star-forming galaxies and SMGs out to z ∼ 3 − 5 found a GDR that might increase with z, with a typical GDR ∼ 120 − 250 at z ∼ 2 − 4 (e.g., Saintonge et al. 2013; Miettinen et al. 2017). In some cases, the GDR was also found to be < 100, as in recent studies of SMGs (Birkin et al. 2021; Liao et al. 2024). Overall, the WISSH sample and our sample of QSOs show GDRs above 100 on average (this value is also commonly assumed when deriving the dust mass). In particular, our HYPERION QSOs show the highest GDRs on average (GDRs > 100) of the sources at z > 6. Two QSOs, J1319+0950 and J231−20, exhibit particularly low (< 50) GDRs that are comparable with other QSOs from Calura et al. (2014). The tail of low GDRs in our sample, that is, GDR < 50, can mostly be attributed to the different dust mass values because the gas masses are ∼ 1010 M for all the sources (see Sect. 5.2). J1319+0950 and J231−20 both have the highest dust masses in our sample, Mdust ∼ (5 − 6)×108 M.

thumbnail Fig. 16.

Redshift distribution of the GDR. The symbols and colors for the WISSH and our sample are the same as in Fig. 12. We compare our results with the WISSH sample (Bischetti et al. 2021), a sample of 2 < z < 5 star-forming galaxies hosting a heavily obscured AGN in the Chandra Deep Field-South (magenta dots, D’Amato et al. 2020), and a sample of z > 5.5 QSOs (blue squares Calura et al. 2014). The two vertical lines at the bottom right side are the systematic uncertainties induced by the choice of αCO and dust mass estimation (black line, ∼0.3 dex) and r65 or r76 in computing the gas mass (brown line, ∼0.1 dex).

6.1.6. Caveat: Effect of active galactic nuclei on the dust heating

As mentioned in Sect. 5.1, when the SFR is derived, the presence of an AGN at the center of the host galaxy likely plays a role in heating the surrounding dust. In theory, the galaxy is characterized by a distribution of dust temperatures that rise toward the center. Both observations and simulations have shown that in the innermost region of the host galaxy (∼a few hundred parsec), Tdust can be as high as some hundred K, and this is mainly due to the presence of the AGN as a heating source (Walter et al. 2022; Di Mascia et al. 2021, 2023). The resolution of our observations prevents us from performing a spatially resolved study of the dust emission, and therefore, from isolating the warmer central dust component. We therefore modeled the dust emission with a single temperature BB, that is, we mixed the dust emission that is mainly heated by the stellar distribution with the emission that is mainly heated by the central AGN. This implies that the SFR derived above from the best-fit MBB might be contaminated by the contribution of the AGN to the heating of the dust in the center.

To overcome this problem, there are two different possibilities that can also be used simultaneously. First, observations with a resolution of ∼100 − 300 pc can be used to disentangle the warmer dust component, which mainly resides in the center of the colder and more extended component. Then, it is possible to fit the warm and cold dust component separately, either using two MBBs or with an MBB for the cold dust emission and a radiative transfer modeling for the warm dust heated by the AGN. This approach has been adopted at high-z only by Tsukui et al. (2023) for a QSO host galaxy at z = 4.4. They found a warm-dust component with Tdust, warm = 87 K and an AGN contribution to the dust heating of ∼60%.

The availability of very high-resolution data is still quite rare at high z, and this approach has indeed not yet been applied to any QSO at z > 4.4. As an alternative to the spatially resolved study of the SED, radiative transfer models can be used to determine the AGN contribution to the dust heating. Schneider et al. (2015) used a radiative transfer code to follow the transfer of radiation from the central source and from stellar sources through the dusty environment of the host galaxy of QSO SDSS J1148+5251 at z = 6.4. For the stellar sources in the host galaxy, they adopted the SED computed with the PÉGASE population synthesis model using as input the star formation histories, age, and metallicities of the stellar populations predicted by GAMETE/QSOdust. To account for the emission of the dust in the host at longer wavelength, they considered two heating sources, (1) the stellar component and (2) the central AGN. This latter can account for up to ∼70% of the dust emission, supporting the idea that AGN have a strong impact on heating the surrounding dust. The exact amount of the contribution of AGN to the dust heating depends on the prescription adopted to model the central source and the torus. Schneider et al. (2015) tested alternative but still reasonable models and found that the AGN contribution can vary from 30% to 70%. Duras et al. (2017) also investigated the effect of the AGN to the FIR emission in the WISSH quasars at 1.8 < z < 4.6, using the same approach as in Schneider et al. (2015). In particular, considering the least and most luminous QSOs in their sample, they found that the AGN contribution to the FIR fluxes is 43% and 60%, respectively, pointing toward a mild trend with luminosity. Therefore, they assumed an average contribution of ∼50% to the total FIR luminosity, which also accounts for the uncertainties in the radiative transfer model and is in line with the average result found in Schneider et al. (2015).

Considering that the dynamic range of bolometric luminosity, redshift, and properties of the host galaxy of the QSOs in our sample is quite small and that these properties are remarkably like those of J1148+5251, the impact of the AGN on the dust heating may be similar. Therefore, we considered an AGN contribution of ∼50% to the dust heating. A precise estimate of the AGN impact on the FIR emission would require either a spatially resolved study of the SED or radiative transfer modeling. The former is not yet feasible because of the lack of high-resolution data, and the latter is beyond the scope of this project, but will be the goal of future investigations.

6.2. Evolutionary paths

As mentioned in Sect. 1, the BH masses of luminous QSOs at z ≳ 6 can be extremely high for their cosmic epoch, of about ≳109 M. The temporal window within which they could have accreted their mass spans approximately 650−750 Myr (considering a possible seeding era at z = 15 − 22). These masses are not lower than those of hyperluminous QSOs at lower redshift, meaning that the BH growth had to be a fast process and that the process had to stop with a similarly high efficiency after the rapid build-up. QSOs at z ≳ 6 indeed appear to lie above the local MBH − Mdyn correlation, and thus the BH growth seems to precede that of its host galaxy (Volonteri 2012; Pensabene et al. 2020; Fan et al. 2023).

Focusing on the BH dominance evolutionary scenario, which seems to be the most likely formation path of local massive galaxies, we can distinguish two large regimes: the first regime is characterized by an intense and predominant growth of the BH, and the second regime is marked by an intense and predominant growth of the host galaxy10. Additionally, a transition phase can be imagined in which the growth of the two components is in almost perfect balance (a symbiotic growth). Already at 6 < z < 8, SMBHs have reached masses similar to those observed in the most massive local galaxies. This implies that the SMBH growth has to slow down. At these redshifts, we may therefore be able to witness QSOs in their transition phase toward the galaxy regime or in the galaxy regime, for instance, moving toward the local relation in the MBH − Mdyn plane.

In the following, we quantitatively characterize the three different regimes outlined above (BH dominated, galaxy dominated, and symbiotic). For this purpose, we defined the growth efficiency of the galaxy as Ggal = SFR/Mgal, where Mgal = Mdyn − MBH is the mass of the galaxy, and the SFR was corrected for the QSO contribution. This term is a lower limit to the specific SFR of the galaxy defined as SFR/M*, which might be a better probe of the galaxy growth in principle. We did not use the specific SFR because the stellar mass is not yet available for most high-redshift QSOs. We derived the BH growth efficiency as GBH =(1−ϵ)BH/MBH, where ϵ is the radiative efficiency, and BH is the BH accretion rate that depends on the bolometric luminosity of the source11. We assumed ϵ = 0.1 (e.g., Marconi et al. 2004), and we used the BH mass derived from the MgII emission line. Comparing these two terms, we distinguished among (1) GBH > Ggal or black-hole-dominance regime, (2) GBH = Ggal or symbiotic growth, and (3) GBH < Ggal or galaxy-dominance regime. The proportionality factor between GBH and Ggal that allowed us to distinguish among these different regimes can be straightforwardly translated into an angle, that is, an angle > 45° corresponds to the BH-dominance regime, an angle ∼45° to the symbiotic growth, and an angle < 45° to the galaxy regime. The properties derived for our sample of QSOs are summarized in Table 3. In particular, the SFRs were derived in this work, the BH masses and bolometric luminosities for BH (Zappacosta et al. 2023; Mazzucchelli et al. 2023), and the dynamical masses (Neeleman et al. 2021; Riechers et al. 2009; Tripodi et al. 2022) were taken from the literature as reported in the last column of Table 3.

Table 3.

Properties of QSOs in our sample.

Before we discuss our results, it is important to briefly illustrate the possible caveats. As reported in Table 3, both MBH and Mdyn carry large systematic uncertainties12 (0.1−0.3 dex). Additionally, even though the SFR was derived with high statistical accuracy and corrected for a reasonable AGN contribution (50%; see Sect. 6.1.6), this contribution may in principle vary from ∼30% to 70%. These uncertainties (a factor of ∼2 for MBH and Mdyn, and a factor of ∼1.5 for the SFR) affect the inferred evolutionary scenario because GBH may vary by a factor of 2 and Ggal by a factor of 3 at most. Keeping in mind these systematics, we discuss our results below, based on the best-fitting values for MBH, Mdyn, and SFR (reported in Table 3).

We evaluated the evolutionary state of the QSOs in our sample. Fig. 17 presents our results, where the HYPERION QSOs in our sample are shown as a star with green contours, and the other QSOs in our sample as a star with red contours. The stars are color-coded as a function of Ggal. The slope of the arrows associated with each star corresponds to the angle reported in Table 3, that is, to the proportionality factor of the GBH − Ggal relation. Given the large systematics quoted above, we did not draw the uncertainty of the arrows for clarity. Overall, the growth of the QSOs in our sample is mainly dominated by the galaxy, or it is symbiotic (see also Table 3). In particular, the color-code of the stars shows that the closer a QSO to the local relation, the slower the growth of the galaxy on average.

thumbnail Fig. 17.

BH mass vs. dynamical mass for our sample (stars with green contours for HYPERION QSOs, and red contours otherwise), WISSH QSOs at z ∼ 2 − 4 (gray diamonds; Bischetti et al. 2021), and luminous z ∼ 4 − 7 QSOs (gray dots and gray squares; Venemans et al. 2016, 2017; Willott et al. 2013, 2015, 2017; Kimball et al. 2015; Trakhtenbrot et al. 2017; Feruglio et al. 2018; Mortlock et al. 2011; De Rosa et al. 2014; Kashikawa et al. 2015; Neeleman et al. 2021). The dashed black line (and shaded area) is the local relation from Kormendy & Ho (2013). The light green dots are the remaining HYPERION QSOs for which we were not yet able to perform a detailed study of the dust properties due to a lack of observations in the submm regime. The stars are color-coded based on the value of Ggal. The thin red arrows indicate upper limits on the dynamical mass.

In Tripodi et al. (2022, 2023b), we investigated the evolutionary path of J2310+1855 and J0100+2802 in detail. For J0100+2802, we found GBH > Ggal, suggesting that the BH still dominates the process of BH-galaxy growth in this QSO at z = 6.32713. On the other hand, in QSO J2310+1855 at z ∼ 6, AGN feedback might be slowing down the accretion onto the SMBH, while the host galaxy grows fast (Tripodi et al. 2022; Bischetti et al. 2022). One of the likely causes of the slow-down of the SMBH accretion is radiatively driven AGN winds that impact on the accreting matter, providing enough momentum to stop further accretion, and which can further propagate outward on the scale of the host galaxy. In J2310+1855, the SMBH accretion may be limited by the ionized wind traced by a CIV broad absorption line (BAL) system (Bischetti et al. 2022). J2310+18655 also shows evidence of a [CII] outflow approximately located in the central kiloparsec, with an outflow mass Mout = 3.5 × 108 M. Additionally, Shao et al. (2022), Butler et al. (2023) detected molecular outflows traced by OH and OH+. QSO J231−20, which shows a galaxy growth similar to J2310+1855, is also a BAL QSO (Bischetti et al. 2022). Additionally, evidence of high-velocity [CII] emission wings and/or OH absorption wings indicating powerful and fast outflows was found in J2054−0005 (Salak et al. 2024), J1148+5251 (Maiolino et al. 2012; Cicone et al. 2015), and tentatively in J1319+0950 (Herrera-Camus et al. 2019). For all of these QSOs, the angles of the GBH − Ggal relation are ≲45°. Therefore, the onset of the symbiotic or galaxy dominated regime may be linked to a phase of strong feedback that hampers the BH growth. Three QSOs still experience a BH-dominated growth: J0100+2802, J1342+0928, and J183+05. J0100+2802 and J1342+0928, both belonging to the HYPERION sample, have the highest and lowest BH mass in our sample, respectively. Given that the locally observed SMBH masses do not exceed 1011 M and that the BH of J0100+2802 is the most massive observed at z > 6 with already 1010 M, a substantial BH growth is an unlikely prospect. BH-dominated growth therefore is indeed a peculiar scenario for J0100+2802. Conversely, it is likely that J1342+0928 is in the process of strong BH growth given its low BH mass. Interestingly, J183+05 is closest to the local relation. Both in J0100+2802 and J183+05, high-velocity wings were detected in [CII] and OH, respectively (Tripodi et al. 2024; Butler et al. 2023). The outflows here might also be related to the fact that both QSOs may be approaching a phase of symbiotic growth. For the uncertainties on MBH, Mdyn, and SFR quoted above, both J0100+2802 and J183+05 would fall in the symbiotic-growth regime on average. The other QSOs in our sample would keep the same classification as reported in Table 3.

To summarize, our study suggests that QSOs at z ≳ 6 experience a phase of intense galaxy growth. This may be connected to the emergence of strong outflows that are able to regulate the BH growth. On the MBH − Mdyn plane, high-z QSOs appear to converge toward the massive end of the local relation. This suggests that they are viable and plausible candidate progenitors of the massive galaxies found in the local Universe.

7. Summary

We exploited newly acquired and archival ALMA observations to perform an analysis of the cold-dust SED in a sample of ten QSOs at z ≳ 6, deriving the dust properties and SFR with the smallest statistical uncertainty (< 10 − 25%). We also investigated the molecular gas in our sample in order to estimate the GDR. The final goal was to discuss the evolutionary path of the QSOs in our sample by assessing whether high-z QSOs can be considered the progenitors of local massive galaxies. We divided our sample into two subsamples depending on whether they belonged to the HYPERION sample. The HYPERION QSOs experienced an intense and rapid SMBH growth, and we therefore aimed to examine whether they showed similar or distinct properties in terms of host galaxies compared to other QSOs at the same redshift. We summarize our main findings below.

  • The analysis of the CO(6−5) and CO(7−6) emission lines in select QSOs provides insights into their molecular gas masses, which average around 1010 M, which is consistent with typical values for high-redshift QSOs. Our findings support the picture in which high-z QSO host galaxies have large gas reservoirs that constitute the fuel for star formation. The combination of this information with the SFR and dust masses estimated from the analysis of the SEDs can provide crucial insights into the galaxy assembly and evolution.

  • Proprietary and archival ALMA observations in bands 8 and 9 enabled precise constraints on the dust properties and SFR of four QSOs at z > 6 for the first time. Taking advantage of these new results, we developed a mean cold-dust SED by combining all the observations available at 50 μm < λ < 500 μm for the ten QSOs in our sample. This offered a comprehensive view of the dust properties in the QSO hosts.

  • Our investigation of the redshift distribution of dust temperatures revealed a large scatter in dust temperature between QSOs and typical star-forming galaxies at fixed redshift, but indicated a general trend of increasing Tdust with redshift, as theoretically expected. Considering the whole population of galaxies at 0 < z < 7, our best-fitting Tdust − z relation is of the form Tdust ∝ (1 + z)0.7 ± 0.1, which is steeper than expected from theory (Tdust ∝ (1 + z)0.4 (see Sommovigo et al. 2022). This implies that the variation in Tdust in different sources has non-negligible dependences on other physical properties, such as optical depth and metallicity. When we focused on QSOs alone, we recovered the theoretical results with Tdust ∝ (1 + z)0.35 ± 0.04.

  • Investigating the variation in the dust emissivity index β with redshift, we found that it is approximately constant with z (β ∼ 1.6), indicating that all sources share similar dust properties. Two QSOs, J0100+2802 and J036+03, show β > 2, suggesting that they may have peculiar dust properties. The analysis of the dust extinction curves in the rest-frame UV-optical could provide complementary information on the properties of dust in these QSOs.

  • All the QSOs in our sample are highly star forming, with an SFR ∼ 200 − 2000 M yr−1. Compared with local QSOs, the SFR is higher at z > 2 by ∼2 orders of magnitudes on average, ranging from few hundred to thousands of M yr−1. This supports the well-known concept that high-z QSOs are hosted in highly star-forming galaxies. The dust masses are higher at high z as well, by ≲2 order of magnitudes on average. As a preliminary result, the dust masses of the HYPERION QSOs in our sample are lower on average (about a factor of 2) than the other QSOs in our sample.

  • The observed high SFR in our sample yields high SFEs, and thus, very low gas depletion timescales (τdep ∼ 10−2 Gyr). The latter is connected to the observed high dust temperatures and indicates that the nature of our sample’s host galaxies is indeed star bursting. In order to derive firm conclusions on the SFE in high-z QSOs, we need a systematic study of the SFE in a larger sample that includes QSOs at lower luminosities.

  • We found a large scatter of the GDR in our sample, from 30 to 250. The lowest measured GDRs are due to massive reservoirs of dust, Mdust ∼ 5 − 6 × 108 M, which pose challenges to theoretical modeling of dust formation. Interestingly, the HYPERION QSOs show the highest GDRs in our sample owing to their lower dust masses, Mdust ∼ 2 − 9 × 107 M, whereas their H2 gas reservoirs are in line with those previously found in QSOs at the same z (MH2 ∼ 1010 M).

  • We were able to investigate and discuss the evolutionary path of our sample of ten QSOs with an accurate determination of the dust properties and SFR. Our study suggests that QSOs at z ≳ 6 experience a phase of rapid galaxy growth. This may be connected to the emergence of strong outflows that are able to regulate the BH growth. BALs were indeed detected in J2310+1855 and J231−20 (Bischetti et al. 2022), and evidence of powerful and fast outflows was found in J2310+1855 (Tripodi et al. 2022; Shao et al. 2022; Butler et al. 2023), J2054−0005 (Salak et al. 2024), J1148+5251 (Maiolino et al. 2012; Cicone et al. 2015), and tentatively in J1319+0950 (Herrera-Camus et al. 2019). On the MBH − Mdyn plane, high-z QSOs appear to be converging toward the massive end of the local relation. This makes high-z QSOs viable and plausible candidate progenitors of massive galaxies found in the local Universe. Interestingly, the average pathway pursued by high-z QSOs to end up as local massive galaxies involves an intense BH growth, which is supported by the upward offset from the local MBH − Mdyn relation, followed by a substantial growth of the galaxy. This is in contrast with the picture of the formation of massive local galaxies via symbiotic growth. Our scenario is further supported by the evidence of a stellar bulge in J2310+1855 (Tripodi et al. 2023a), which indicates that the structure of QSOs at z ∼ 6 is surprisingly similar to that typical of local massive galaxies.


1

These QSOs are J0100+2802, J036+03, J0224−4711, J231−20, J2310+1855, J1319+0950, and J183+05.

2

The maximum value of the velocity dispersion toward the nucleus is usually affected by beam smearing (Davies et al. 2011).

3

All luminosities were computed following Eqs. 1 and 3 of Solomon & Vanden Bout (2005).

4

The bins were defined considering the different ALMA bands and trying to have a similar number of data in each bin.

5

They also derived the dust properties using the dust model of Compiègne et al. (2011) with β = 1.91 finding that the dust masses are systematically larger by about 20%−40% when adopting the model of Draine et al. (2007).

6

‘Reionization Era Bright Emission Line Survey’ (REBELS; PI: Bouwens) is an ALMA large program targeting 40 of the brightest known galaxies at z > 6.5 (Bouwens et al. 2022).

7

The Band 8 continuum of REBELS-12 was undetected, therefore they derived an upper-limit.

8

Fixing β introduces an uncertainty of ∼25% in the estimate of Tdust.

9

For both samples, the mean value of the distribution corresponds to the median.

10

As a note of caution, we must emphasize that during the BH- (galaxy-) dominated regime, the galaxy (BH) is still growing, but more slowly and/or less efficiently than the BH (galaxy).

11

BH =Lbol/(ϵ c2), where c is the speed of light.

12

The errors reported on Mdyn take into account the systematics due to the circular velocity and the radius of the galaxy used to derive Mdyn (see also the discussion in Neeleman et al. 2021). For the majority of the QSOs in our sample, Mdyn was derived in Neeleman et al. (2021) in a consistent way.

13

The slope of the arrow in Tripodi et al. (2023b) is higher than reported in this work. This is because the SFR was computed using the Chabrier IMF rather then the Kroupa IMF (as in this work), which yielded an SFR that was lower by a factor of ∼0.7 than in this work, and therefore, to a lower Ggal and a higher slope of the arrow.

Acknowledgments

We thank the referee for providing insightful comments that helped the clarity and robustness of this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01633.S, 2015.1.00399.S, 2018.1.01790.S, 2021.2.00151.S, 2017.1.01472.S, 2018.1.01188.S, 2021.1.00934.S, 2021.2.00064.S, 2018.1.01289.S, 2019.1.00672.S, 2017.1.01195.S, 2016.1.01063.S. ALMA is a partnership of ESO (representing its member states), NFS (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The project leading to this publication has received support from ORP, that is funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 101004719 [ORP]. RT acknowledges financial support from the University of Trieste. RT, CF, FF, MB, EP acknowledge support from PRIN MIUR project “Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”, contract #2017PH3WAT. C.-C.C. acknowledges support from the National Science and Technology Council of Taiwan (NSTC 111-2112M-001-045-MY3), as well as Academia Sinica through the Career Development Award (AS-CDA-112-M02). EP, LZ and CF acknowledge financial support from the Bando Ricerca Fondamentale INAF 2022 Large Grant “Toward a holistic view of the Titans: multi-band observations of z > 6 QSOs powered by greedy supermassive black-holes”. RM acknowledges ERC Advanced Grant 695671 QUENCH, and support from the UK Science and Technology Facilities Council (STFC). RM also acknowledges funding from a research professorship from the Royal Society. SC acknowledges support from the European Union (ERC, WINGS,101040227). Facilities: ALMA. Software: astropy (Astropy Collaboration 2022), Matplotlib (Hunter 2007), SciPy (Gommers et al. 2023), CASA (v5.1.1-5, CASA Team 2022).

References

  1. Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2023, MNRAS, 518, 6142 [Google Scholar]
  2. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473 [Google Scholar]
  4. Bakx, T. J. L. C., Tamura, Y., Hashimoto, T., et al. 2020, MNRAS, 493, 4294 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bakx, T. J. L. C., Sommovigo, L., Carniani, S., et al. 2021, MNRAS, 508, L58 [NASA ADS] [CrossRef] [Google Scholar]
  6. Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694 [NASA ADS] [CrossRef] [Google Scholar]
  7. Beeston, R. A., Wright, A. H., Maddox, S., et al. 2018, MNRAS, 479, 1077 [NASA ADS] [Google Scholar]
  8. Behrens, C., Pallottini, A., Ferrara, A., Gallerani, S., & Vallini, L. 2018, MNRAS, 477, 552 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bendo, G. J., Urquhart, S. A., Serjeant, S., et al. 2023, MNRAS, 522, 2995 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bertoldi, F., Carilli, C. L., Cox, P., et al. 2003, A&A, 406, L55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2 [Google Scholar]
  12. Birkin, J. E., Weiss, A., Wardlow, J. L., et al. 2021, MNRAS, 501, 3926 [Google Scholar]
  13. Bischetti, M., Feruglio, C., Piconcelli, E., et al. 2021, A&A, 645, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bischetti, M., Feruglio, C., D’Odorico, V., et al. 2022, Nature, 605, 244 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112 [NASA ADS] [CrossRef] [Google Scholar]
  17. Butler, K. M., van der Werf, P. P., Topkaras, T., et al. 2023, ApJ, 944, 134 [NASA ADS] [CrossRef] [Google Scholar]
  18. Calabrò, A., Pentericci, L., Santini, P., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449768 [Google Scholar]
  19. Calura, F., Gilli, R., Vignali, C., et al. 2014, MNRAS, 438, 2765 [NASA ADS] [CrossRef] [Google Scholar]
  20. Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105 [NASA ADS] [CrossRef] [Google Scholar]
  21. Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939 [NASA ADS] [Google Scholar]
  22. CASA Team (Bean, B., et al.) 2022, PASP, 134, 114501 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  24. Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv e-prints [arXiv:1612.05560] [Google Scholar]
  25. Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Clements, D. L., Pearson, C., Farrah, D., et al. 2018, MNRAS, 475, 2097 [NASA ADS] [CrossRef] [Google Scholar]
  27. Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103 [Google Scholar]
  28. Cooper, O. R., Casey, C. M., Zavala, J. A., et al. 2022, ApJ, 930, 32 [NASA ADS] [CrossRef] [Google Scholar]
  29. da Cunha, E., Eminian, C., Charlot, S., & Blaizot, J. 2010, MNRAS, 403, 1894 [Google Scholar]
  30. da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13 [Google Scholar]
  31. da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, ApJ, 919, 30 [NASA ADS] [CrossRef] [Google Scholar]
  32. Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJ, 714, L118 [NASA ADS] [CrossRef] [Google Scholar]
  33. D’Amato, Q., Gilli, R., Prandoni, I., et al. 2020, A&A, 641, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Davies, R., Förster Schreiber, N. M., Cresci, G., et al. 2011, ApJ, 741, 69 [NASA ADS] [CrossRef] [Google Scholar]
  35. Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Nature, 545, 457 [Google Scholar]
  36. Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97 [Google Scholar]
  37. Decarli, R., Pensabene, A., Venemans, B., et al. 2022, A&A, 662, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Decarli, R., Pensabene, A., Diaz-Santos, T., et al. 2023, A&A, 673, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. De Rosa, G., Venemans, B. P., Decarli, R., et al. 2014, ApJ, 790, 145 [Google Scholar]
  40. Di Mascia, F., Gallerani, S., Behrens, C., et al. 2021, MNRAS, 503, 2349 [NASA ADS] [CrossRef] [Google Scholar]
  41. Di Mascia, F., Carniani, S., Gallerani, S., et al. 2023, MNRAS, 518, 3667 [Google Scholar]
  42. Ding, X., Onoue, M., Silverman, J. D., et al. 2023, Nature, 621, 51 [NASA ADS] [CrossRef] [Google Scholar]
  43. D’Odorico, V., Feruglio, C., Ferrara, A., et al. 2018, ApJ, 863, L29 [Google Scholar]
  44. D’Odorico, V., Bañados, E., Becker, G. D., et al. 2023, MNRAS, 523, 1399 [CrossRef] [Google Scholar]
  45. Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866 [NASA ADS] [CrossRef] [Google Scholar]
  46. Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891 [Google Scholar]
  47. Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828 [Google Scholar]
  48. Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510 [NASA ADS] [CrossRef] [Google Scholar]
  49. Duras, F., Bongiorno, A., Piconcelli, E., et al. 2017, A&A, 604, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Eilers, A.-C., Simcoe, R. A., Yue, M., et al. 2023, ApJ, 950, 68 [CrossRef] [Google Scholar]
  51. Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21 [Google Scholar]
  52. Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020, MNRAS, 498, 4192 [NASA ADS] [CrossRef] [Google Scholar]
  53. Fan, X., Strauss, M. A., Schneider, D. P., et al. 2003, AJ, 125, 1649 [Google Scholar]
  54. Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373 [NASA ADS] [CrossRef] [Google Scholar]
  55. Feruglio, C., Fiore, F., Carniani, S., et al. 2018, A&A, 619, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Feruglio, C., Maio, U., Tripodi, R., et al. 2023, ApJ, 954, L10 [NASA ADS] [CrossRef] [Google Scholar]
  57. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  58. Frias Castillo, M., Rybak, M., Hodge, J., et al. 2024, A&A, 683, A211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Gallerani, S., Ferrara, A., Neri, R., & Maiolino, R. 2014, MNRAS, 445, 2848 [Google Scholar]
  60. Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673 [Google Scholar]
  61. Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091 [NASA ADS] [CrossRef] [Google Scholar]
  62. Gommers, R., Virtanen, P., Burovski, E., et al. 2023, https://doi.org/10.5281/zenodo.7655153 [Google Scholar]
  63. Greiner, J., Bolmer, J., Yates, R. M., et al. 2021, A&A, 654, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93 [Google Scholar]
  65. Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hashimoto, T., Inoue, A. K., Tamura, Y., et al. 2019, PASJ, 71, 109 [NASA ADS] [CrossRef] [Google Scholar]
  67. Herrera-Camus, R., Tacconi, L., Genzel, R., et al. 2019, ApJ, 871, 37 [NASA ADS] [CrossRef] [Google Scholar]
  68. Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2020, A&A, 633, L4 [Google Scholar]
  69. Hines, D. C., Krause, O., Rieke, G. H., et al. 2006, ApJ, 641, L85 [NASA ADS] [CrossRef] [Google Scholar]
  70. Hirashita, H., Ferrara, A., Dayal, P., & Ouchi, M. 2014, MNRAS, 443, 1704 [NASA ADS] [CrossRef] [Google Scholar]
  71. Hjorth, J., Gall, C., & Michałowski, M. J. 2014, ApJ, 782, L23 [NASA ADS] [CrossRef] [Google Scholar]
  72. Hu, H., Inayoshi, K., Haiman, Z., et al. 2022, ApJ, 935, 140 [NASA ADS] [CrossRef] [Google Scholar]
  73. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  74. Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27 [NASA ADS] [CrossRef] [Google Scholar]
  75. Inayoshi, K., Nakatani, R., Toyouchi, D., et al. 2022, ApJ, 927, 237 [NASA ADS] [CrossRef] [Google Scholar]
  76. Iwamuro, F., Kimura, M., Eto, S., et al. 2004, ApJ, 614, 69 [NASA ADS] [CrossRef] [Google Scholar]
  77. Izumi, T., Onoue, M., Shirakata, H., et al. 2018, PASJ, 70, 36 [NASA ADS] [Google Scholar]
  78. Izumi, T., Onoue, M., Matsuoka, Y., et al. 2019, PASJ, 71, 111 [CrossRef] [Google Scholar]
  79. Jiang, L., Fan, X., Hines, D. C., et al. 2006, AJ, 132, 2127 [NASA ADS] [CrossRef] [Google Scholar]
  80. Jiang, L., Fan, X., Annis, J., et al. 2008, AJ, 135, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  81. Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222 [Google Scholar]
  82. Kashikawa, N., Ishizaki, Y., Willott, C. J., et al. 2015, ApJ, 798, 28 [Google Scholar]
  83. Kennicutt, R. C. 1998, ARA&A, 36, 189 [NASA ADS] [CrossRef] [Google Scholar]
  84. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531 [NASA ADS] [CrossRef] [Google Scholar]
  85. Kimball, A. E., Lacy, M., Lonsdale, C. J., & Macquart, J. P. 2015, MNRAS, 452, 88 [NASA ADS] [CrossRef] [Google Scholar]
  86. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
  87. Koss, M. J., Strittmatter, B., Lamperti, I., et al. 2021, ApJS, 252, 29 [CrossRef] [Google Scholar]
  88. Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076 [NASA ADS] [CrossRef] [Google Scholar]
  89. Lamperti, I., Saintonge, A., De Looze, I., et al. 2019, MNRAS, 489, 4389 [Google Scholar]
  90. Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21 [CrossRef] [Google Scholar]
  91. Leipski, C., Meisenheimer, K., Walter, F., et al. 2013, ApJ, 772, 103 [NASA ADS] [CrossRef] [Google Scholar]
  92. Leipski, C., Meisenheimer, K., Walter, F., et al. 2014, ApJ, 785, 154 [Google Scholar]
  93. Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, ApJ, 737, 12 [Google Scholar]
  94. Li, Y., Hopkins, P. F., Hernquist, L., et al. 2008, ApJ, 678, 41 [NASA ADS] [CrossRef] [Google Scholar]
  95. Li, J., Wang, R., Riechers, D., et al. 2020, ApJ, 889, 162 [NASA ADS] [CrossRef] [Google Scholar]
  96. Liao, C.-L., Chen, C.-C., Wang, W.-H., et al. 2024, ApJ, 961, 226 [NASA ADS] [CrossRef] [Google Scholar]
  97. Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66 [Google Scholar]
  99. Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023, A&A, submitted [Google Scholar]
  100. Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59 [NASA ADS] [CrossRef] [Google Scholar]
  101. Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169 [Google Scholar]
  102. Marshall, M. A., Mechtley, M., Windhorst, R. A., et al. 2020, ApJ, 900, 21 [NASA ADS] [CrossRef] [Google Scholar]
  103. Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, ApJ, 849, 91 [Google Scholar]
  104. Mazzucchelli, C., Bischetti, M., D’Odorico, V., et al. 2023, A&A, 676, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. McKay, S. J., Barger, A. J., Cowie, L. L., Bauer, F. E., & Rosenthal, M. J. N. 2023, ApJ, 951, 48 [CrossRef] [Google Scholar]
  106. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127 [Google Scholar]
  107. Mechtley, M., Windhorst, R. A., Ryan, R. E., et al. 2012, ApJ, 756, L38 [NASA ADS] [CrossRef] [Google Scholar]
  108. Miettinen, O., Delvecchio, I., Smolčić, V., et al. 2017, A&A, 606, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Mortlock, D. J., Patel, M., Warren, S. J., et al. 2008, ArXiv e-prints [arXiv:0810.3859] [Google Scholar]
  110. Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616 [Google Scholar]
  111. Neeleman, M., Bañados, E., Walter, F., et al. 2019, ApJ, 882, 10 [Google Scholar]
  112. Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141 [NASA ADS] [CrossRef] [Google Scholar]
  113. Novak, M., Bañados, E., Decarli, R., et al. 2019, ApJ, 881, 63 [NASA ADS] [CrossRef] [Google Scholar]
  114. Pacucci, F., & Loeb, A. 2024, ApJ, 964, 154 [NASA ADS] [CrossRef] [Google Scholar]
  115. Pacucci, F., Nguyen, B., Carniani, S., Maiolino, R., & Fan, X. 2023, ApJ, 957, L3 [CrossRef] [Google Scholar]
  116. Papadopoulos, P., Ivison, R., Carilli, C., & Lewis, G. 2001, Nature, 409, 58 [NASA ADS] [CrossRef] [Google Scholar]
  117. Parente, M., Ragone-Figueroa, C., Granato, G. L., & Lapi, A. 2023, MNRAS, 521, 6105 [NASA ADS] [CrossRef] [Google Scholar]
  118. Pensabene, A., Carniani, S., Perna, M., et al. 2020, A&A, 637, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Pensabene, A., Decarli, R., Bañados, E., et al. 2021, A&A, 652, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Petric, A. O., Ho, L. C., Flagey, N. J. M., & Scoville, N. Z. 2015, ApJS, 219, 22 [Google Scholar]
  121. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Pons, E., McMahon, R. G., Banerji, M., & Reed, S. L. 2020, MNRAS, 491, 3884 [Google Scholar]
  123. Pozzi, F., Calura, F., Zamorani, G., et al. 2020, MNRAS, 491, 5073 [NASA ADS] [CrossRef] [Google Scholar]
  124. Reed, S. L., McMahon, R. G., Martini, P., et al. 2017, MNRAS, 468, 4702 [NASA ADS] [CrossRef] [Google Scholar]
  125. Reed, S. L., Banerji, M., Becker, G. D., et al. 2019, MNRAS, 487, 1874 [Google Scholar]
  126. Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82 [NASA ADS] [CrossRef] [Google Scholar]
  127. Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, ApJ, 902, 78 [NASA ADS] [CrossRef] [Google Scholar]
  128. Riechers, D. A., Walter, F., Carilli, C. L., & Lewis, G. F. 2009, ApJ, 690, 463 [NASA ADS] [CrossRef] [Google Scholar]
  129. Rosario, D. J., Burtscher, L., Davies, R. I., et al. 2018, MNRAS, 473, 5658 [NASA ADS] [CrossRef] [Google Scholar]
  130. Saintonge, A., Kauffmann, G., Wang, J., et al. 2011, MNRAS, 415, 61 [NASA ADS] [CrossRef] [Google Scholar]
  131. Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2 [Google Scholar]
  132. Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22 [Google Scholar]
  133. Salak, D., Hashimoto, T., Inoue, A. K., et al. 2024, ApJ, 962, 1 [NASA ADS] [CrossRef] [Google Scholar]
  134. Salomé, Q., Krongold, Y., Longinotti, A. L., et al. 2023, MNRAS, 524, 3130 [CrossRef] [Google Scholar]
  135. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  136. Sanders, D. B., Mazzarella, J. M., Kim, D. C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607 [Google Scholar]
  137. Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Santini, P., Fontana, A., Castellano, M., et al. 2023, ApJ, 942, L27 [NASA ADS] [CrossRef] [Google Scholar]
  139. Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19 [NASA ADS] [CrossRef] [Google Scholar]
  140. Schneider, R., Bianchi, S., Valiante, R., Risaliti, G., & Salvadori, S. 2015, A&A, 579, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  142. Scoville, N., Faisst, A., Weaver, J., et al. 2023, ApJ, 943, 82 [NASA ADS] [CrossRef] [Google Scholar]
  143. Shangguan, J., Ho, L. C., Bauer, F. E., Wang, R., & Treister, E. 2020a, ApJS, 247, 15 [NASA ADS] [CrossRef] [Google Scholar]
  144. Shangguan, J., Ho, L. C., Bauer, F. E., Wang, R., & Treister, E. 2020b, ApJ, 899, 112 [NASA ADS] [CrossRef] [Google Scholar]
  145. Shao, Y., Wang, R., Jones, G. C., et al. 2017, ApJ, 845, 138 [NASA ADS] [CrossRef] [Google Scholar]
  146. Shao, Y., Wang, R., Carilli, C. L., et al. 2019, ApJ, 876, 99 [NASA ADS] [CrossRef] [Google Scholar]
  147. Shao, Y., Wang, R., Weiss, A., et al. 2022, A&A, 668, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  148. Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677 [NASA ADS] [CrossRef] [Google Scholar]
  149. Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 497, 956 [NASA ADS] [CrossRef] [Google Scholar]
  150. Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122 [NASA ADS] [CrossRef] [Google Scholar]
  151. Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
  152. Stefan, I. I., Carilli, C. L., Wagg, J., et al. 2015, MNRAS, 451, 1713 [NASA ADS] [CrossRef] [Google Scholar]
  153. Stone, M. A., Lyu, J., Rieke, G. H., Alberts, S., & Hainline, K. N. 2024, ApJ, 964, 90 [NASA ADS] [CrossRef] [Google Scholar]
  154. Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [Google Scholar]
  155. Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [NASA ADS] [CrossRef] [Google Scholar]
  156. Toyouchi, D., Inayoshi, K., Hosokawa, T., & Kuiper, R. 2021, ApJ, 907, 74 [NASA ADS] [CrossRef] [Google Scholar]
  157. Trakhtenbrot, B., Lira, P., Netzer, H., et al. 2017, ApJ, 836, 8 [NASA ADS] [CrossRef] [Google Scholar]
  158. Trinca, A., Schneider, R., Valiante, R., et al. 2022, MNRAS, 511, 616 [NASA ADS] [CrossRef] [Google Scholar]
  159. Tripodi, R., Feruglio, C., Fiore, F., et al. 2022, A&A, 665, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  160. Tripodi, R., Lelli, F., Feruglio, C., et al. 2023a, A&A, 671, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  161. Tripodi, R., Feruglio, C., Kemper, F., et al. 2023b, ApJ, 946, L45 [NASA ADS] [CrossRef] [Google Scholar]
  162. Tripodi, R., Scholtz, J., Maiolino, R., et al. 2024, A&A, 682, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  163. Tsukui, T., Wisnioski, E., Krumholz, M. R., & Battisti, A. 2023, MNRAS, 523, 4654 [NASA ADS] [CrossRef] [Google Scholar]
  164. Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, A&A, 677, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  165. Uzgil, B. D., Bradford, C. M., Hailey-Dunsheath, S., Maloney, P. R., & Aguirre, J. E. 2016, ApJ, 832, 209 [CrossRef] [Google Scholar]
  166. Valiante, R., Schneider, R., Salvadori, S., & Bianchi, S. 2011, MNRAS, 416, 1916 [Google Scholar]
  167. Vallini, L., Pallottini, A., Ferrara, A., et al. 2018, MNRAS, 473, 271 [NASA ADS] [CrossRef] [Google Scholar]
  168. Vallini, L., Witstok, J., Sommovigo, L., et al. 2024, MNRAS, 527, 10 [Google Scholar]
  169. Venemans, B. P., Bañados, E., Decarli, R., et al. 2015, ApJ, 801, L11 [Google Scholar]
  170. Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37 [Google Scholar]
  171. Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 837, 146 [NASA ADS] [CrossRef] [Google Scholar]
  172. Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130 [Google Scholar]
  173. Viero, M. P., Sun, G., Chung, D. T., Moncelsi, L., & Condon, S. S. 2022, MNRAS, 516, L30 [NASA ADS] [CrossRef] [Google Scholar]
  174. Vietri, G., Piconcelli, E., Bischetti, M., et al. 2018, A&A, 617, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  175. Volonteri, M. 2012, Science, 337, 544 [NASA ADS] [CrossRef] [Google Scholar]
  176. Walter, F., Bertoldi, F., Carilli, C., et al. 2003, Nature, 424, 406 [Google Scholar]
  177. Walter, F., Neeleman, M., Decarli, R., et al. 2022, ApJ, 927, 21 [NASA ADS] [CrossRef] [Google Scholar]
  178. Wang, R., Wagg, J., Carilli, C. L., et al. 2011, AJ, 142, 101 [NASA ADS] [CrossRef] [Google Scholar]
  179. Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44 [Google Scholar]
  180. Wang, F., Yang, J., Fan, X., et al. 2019a, ApJ, 884, 30 [Google Scholar]
  181. Wang, F., Wang, R., Fan, X., et al. 2019b, ApJ, 880, 2 [NASA ADS] [CrossRef] [Google Scholar]
  182. Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1 [Google Scholar]
  183. Wang, F., Yang, J., Hennawi, J. F., et al. 2024, ApJ, 962, L11 [NASA ADS] [CrossRef] [Google Scholar]
  184. Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  185. Willott, C. J., Omont, A., & Bergeron, J. 2013, ApJ, 770, 13 [NASA ADS] [CrossRef] [Google Scholar]
  186. Willott, C. J., Bergeron, J., & Omont, A. 2015, ApJ, 801, 123 [NASA ADS] [CrossRef] [Google Scholar]
  187. Willott, C. J., Bergeron, J., & Omont, A. 2017, ApJ, 850, 108 [Google Scholar]
  188. Witstok, J., Smit, R., Maiolino, R., et al. 2022, MNRAS, 515, 1751 [Google Scholar]
  189. Witstok, J., Jones, G. C., Maiolino, R., Smit, R., & Schneider, R. 2023, MNRAS, 523, 3119 [NASA ADS] [CrossRef] [Google Scholar]
  190. Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512 [Google Scholar]
  191. Yang, J., Venemans, B., Wang, F., et al. 2019, ApJ, 880, 153 [NASA ADS] [CrossRef] [Google Scholar]
  192. Yang, J., Wang, F., Fan, X., et al. 2023, ApJ, 951, L5 [NASA ADS] [CrossRef] [Google Scholar]
  193. Yue, M., Eilers, A.-C., Simcoe, R. A., et al. 2024, ApJ, 966, 176 [NASA ADS] [CrossRef] [Google Scholar]
  194. Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976 [Google Scholar]
  195. Zappacosta, L., Piconcelli, E., Fiore, F., et al. 2023, A&A, 678, A201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  196. Zhang, Y., Ouchi, M., Gebhardt, K., et al. 2023a, ApJ, 948, 103 [NASA ADS] [CrossRef] [Google Scholar]
  197. Zhang, H., Behroozi, P., Volonteri, M., et al. 2023b, MNRAS, 523, L69 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Individual objects

A.1. SDSS J010013.02+280225.8 (HYPERION)

For QSO SDSS J010013.02+280225.8 (hereafter J0100+2802) at z[CII] = 6.327 (Wang et al. 2019b), Wu et al. (2015) estimated a bolometric luminosity of Lbol = 4.29 × 1014 L and a BH mass of MBH = 1.24 × 1010 M, making it the most optically luminous QSO with the most massive SMBH known at z > 6. Both measurements have been recently confirmed by JWST (Eilers et al. 2023). Wang et al. (2019b) performed a multi-frequency analysis of the dust SED, but they could not obtain a precise determination of the dust properties, concluding that J0100+28 has either a high dust emissivity (β ≳ 2) or a high dust temperature (Tdust ≳ 60 K), or a combination of thereof.

A.2. PSO J036.5078+03.0498 (HYPERION)

PSO J036.5078+03.0498 (hereafter J036+03, Venemans et al. 2015) at z = 6.5405 was observed for the first time in the Pan-STARRS1 survey (Venemans et al. 2015), and it has a BH mass of MBH,MgII = (2.6 − 3.09) × 109 M from the analysis of MgII emission line and a bolometric luminosity of Lbol = (2.13 − 3.16) × 1047 erg s−1 (Zappacosta et al. 2023; Mazzucchelli et al. 2023). Decarli et al. (2022) studied the CO(6-5), CO(7-6) and [CI]2−1 emission lines of this QSO host with NOEMA observations, and they found LCO(6−5) = 12.7 × 109 K km s−1 pc2, LCO(7−6) = 10.7 × 109 K km s−1 pc2 and L[CI] = 5.7 × 109 K km s−1 pc2, implying a molecular gas mass of M H 2 , CO = 5 . 0 0.6 + 0.5 × 10 10 M $ M_{\rm H2,CO}=5.0^{+0.5}_{-0.6} \times 10^{10}\,\mathrm{M}_\odot $ and M H 2 , [ CI ] = 7 . 1 1.4 + 1.6 × 10 10 M $ M_{\rm H2,[CI]}=7.1^{+1.6}_{-1.4} \times 10^{10}\,\mathrm{M}_\odot $. The [CII] emission shows ordered motion, with a clear and regular velocity gradient in the moment-1 map and vrot/σ > 3, and it has a L[CII] = 3. × 109 L and a size of ∼ 2.4 × 1.6 kpc2 (Venemans et al. 2020; Neeleman et al. 2021). From the modeling of the velocity rotation they estimated a dynamical mass of M dyn = 2 . 9 0.7 + 1.1 × 10 10 M $ M_{\rm dyn}=2.9^{+1.1}_{-0.7}\times 10^{10}\,\mathrm{M}_\odot $, and from the [CII] underlying continuum emission they derived a gas mass of M H 2 , cont = 2 . 8 1.1 + 15 × 10 10 M $ M_{\rm H2,cont}=2.8^{+15}_{-1.1} \times 10^{10}\,\mathrm{M}_\odot $, assuming a gas-to-dust ratio of 100 and a molecular-to-total gas mass fraction of 0.75. Greiner et al. (2021) did not find any companion brighter than M1450(AB) < −26 mag within 0.1-3.3h−1 comoving Mpc search radius, using the simultaneous seven-channel Gamma-ray Burst Optical/Near-infrared Detector, confirming the results of Venemans et al. (2020), who did not detect any companion for this source using ALMA observation of [CII] emission.

A.3. VDES J022426.54-471129.4 (HYPERION)

VDESJ022426.54-471129.4 (hereafter J0224-4711, Reed et al. 2017) at z = 6.5222, firstly discovered by Reed et al. (2017), is one of the most X-ray luminous QSOs at z > 5.5 and the most X-ray luminous QSO at z > 6.5 (Pons et al. 2020; Zappacosta et al. 2023). It belongs to the HYPERION sample, the XQR-30 sample and the ASPIRE survey (Yang et al. 2023). It has a bolometric luminosity of Lbol = 3.47 × 1047 erg s−1, and BH mass of MBH,MgII = (1.30 − 2.29) × 109 M from the analysis of MgII emission line (Reed et al. 2019; Wang et al. 2021; Zappacosta et al. 2023; Mazzucchelli et al. 2023) and MBH,Hβ = 2.15 × 109 M, from the analysis of Hβ emission line (Yang et al. 2023). It has the most extreme broad and blueshifted [OIII] lines observed to date, even compared to observations of lower-redshift QSOs, with a velocity shift of −1690 km s−1 relative to the narrow [OIII], suggesting powerful ionized outflows (Yang et al. 2023).

A.4. PSO J231.6576-20.8335 (HYPERION)

PSO J231.6576-20.8335 (hereafter J231-20) at z = 6.587 has been discovered using the Pan-STARRS1 survey (Mazzucchelli et al. 2017), and it is one of the brightest objects at z > 6.5. Decarli et al. (2017) detected a [CII]-bright nearby companion at < 10 kpc separation. Pensabene et al. (2021) performed an extensive study of both the QSO and its companion, detecting [NII], CO(7-6), CO(10-9) emission lines, two OH transitions and their underlying continuum in both of them. Additionally, the CO(15-14), CO(16-15) and three transitions of H2O emission line have been detected for the central QSO. Analyzing the cold-dust SED of both the QSO and its companion, they derived Tdust = 54 K, Mdust = 5.1 × 108 M and Tdust = 35 − 46 K, Mdust = (2.3 − 3.4) × 108 M, respectively. The estimates on the dust temperature suffer of high uncertainties given that they lack of high-frequency observations to probe the peak of the SED. (Neeleman et al. 2021) found a dynamical mass Mdyn = 1.4 × 1010 M from the modeling of the velocity rotation curve, and reported a BH mass of MBH = 4.1 × 109 M. Finally, Bischetti et al. (2022) classified this as a broad absorption line (BAL) QSO, indicating that J231-20 may be caught in a phase of strong BH feedback.

A.5. ULAS J134208.10+092838.35 (HYPERION)

ULAS J134208.10+092838.35 (hereafter J1342+0928) at z = 7.54, the most distant QSO known to date, was discovered by Bañados et al. (2018) who reported an absolute AB magnitude M1450 = −26.8, bolometric luminosity of Lbold = 1013 L, and an SMBH mass of 8 × 108 M. It was followed-up with NOEMA resulting in the detection of bright [CII] emission and upper limits on several CO lines (Venemans et al. 2017). Novak et al. (2019) presented ALMA observations of the dust continuum and the ISM of the host galaxy J1342+0928. They well constrained the Rayleigh-Jeans tail of the dust SED, deriving a Mdust = 3.5 × 107 M and a SFR∼150 M yr−1, fixing the temperature at 47 K. They also detected many atomic fine structure line, such as [CII], [NII], [OIII], and limits on [CI], [OI] and multiple CO transitions (with a tentative stack detection).

A.6. SDSS J114816.64+525150.3 (HYPERION)

SDSS J114816.64+525150.3 (hereafter J1148+5251) at z = 6.42 was discovered by Fan et al. (2003) who reported an absolute magnitude of M1450 = −27.82. It is one of the most studied QSOs at high-z. It was observed by Subary Telescope (Iwamuro et al. 2004), Spitzer (Jiang et al. 2006; Hines et al. 2006) and Herschel (Leipski et al. 2013), and therefore this allowed a full modeling of the SED of this QSO (Li et al. 2008; Valiante et al. 2011; Schneider et al. 2015; Carniani et al. 2019). In particular, Schneider et al. (2015) derived that the AGN contribution to the dust heating in this QSO can be between 30% and 70%. Gallerani et al. (2014) detected an exceptionally strong CO(17−16) line in this QSO with the Plateau de Bure interferometer (PdBI) and performed a detailed analysis of the CO SLED using previously detected lower CO transitions (Bertoldi et al. 2003; Walter et al. 2003; Riechers et al. 2009).

A.7. SDSS J231038.88+185519.7

QSO SDSS J231038.88+185519.7 (hereafter J2310+1855), first discovered in SDSS (Jiang et al. 2006; Wang et al. 2013), is one of the most FIR-luminous QSOs and one of the brightest optical QSOs known at z ∼ 6, with Lbol = 9.3 × 1013 L. The redshift measured with the QSO rest-frame UV line emission is z = 6.00 ± 0.03 (Wang et al. 2013). Feruglio et al. (2018) detected and analyzed the CO(6-5) and [CII] emission lines and the submm continuum of J2310+1855, deriving a size of the dense molecular gas of 2.9 ± 0.5 kpc and of 1.4 ± 0.2 kpc for the 91.5 GHz dust continuum and a molecular gas mass of M(H2) = (3.2 ± 0.2)× 1010 M. They estimated a dynamical mass of M dyn = ( 4 . 1 0.5 + 9.5 ) × 10 10 M $ M_{\rm dyn} = (4.1^{+9.5}_{-0.5})\times 10^{10}\,\mathrm{M_\odot } $, measuring a disk inclination of i ∼ 50 deg. They also inferred the BH mass from the CIV emission line, measured in the X-shooter/VLT spectrum of the QSO, obtaining MBH = (1.8 ± 0.5) × 109 M. Shao et al. (2019) presented a detailed analysis of the FIR and submm SED and derived a dust temperature of T ∼ 40 K, a dust mass of Mdust = 1.6 × 109 M, a FIR luminosity L FIR 8 1000 μ m = 1.6 × 10 13 L $ L_{\rm FIR}^{8-1000\,\mu m}=1.6\times 10^{13}\,\mathrm{L_{\odot }} $, and an SFR= 2400 − 2700 M yr−1. D’Odorico et al. (2018) detected a very metal-poor, proximate damped Lyman α system (DLA) located at z= 5.938646 ± 0.000007 in the X-shooter/VLT spectrum of J2310, which was associated with a CO emitting source at z = 5.939. This source, called Serenity-18, was detected through its CO(6-5) emission line at [RA, DEC] = 23:10:38.44, 18:55:21.95.

A.8. SDSS J205406.49-000514.8

SDSS J205406.49-000514.8 (hereafter J2054-0005) at z = 6.39 was selected from SDSS stripe 82 with m1450 = 20.60, that is, about one magnitude fainter than the objects discovered from the SDSS main survey (Jiang et al. 2008). Leipski et al. (2014) reported observations in band z, y, J, H, K and with Herschel, however they were not able to fully study its SED due to lack of observations in the mm/submm regime. Wang et al. (2013) reported a BH mass of MBH = 8.6 × 108 M, later updated to MBH = 1.48 × 109 M using the MgII emission line (Neeleman et al. 2021). Neeleman et al. (2021) studied the rotation curve of J2054-0005 using a high-resolution ALMA observation of the [CII] emission of this object, and they determined a lower limit for the dynamical mass of Mdyn > 2.9 × 1010 M.

A.9. ULAS J131911.29+095051.4

ULAS J131911.29+095051.4 (hereafter J1319+0950) at z = 6.133 was discovered in the UKIRT Infrared Deep Sky Survey (UKIDSS) with m1450 = −19.65, which lies in the typical magnitude range of the optically bright z ∼ 6 quasars selected from the SDSS main survey (Mortlock et al. 2008). Wang et al. (2011) detected the CO(6-5) emission line and its underlying continuum, deriving a gas mass of Mgas = 1.5 × 1010 M. Later, Wang et al. (2013) analyzed the [CII] line emission and underlying continuum, deriving a dynamical mass of Mdyn = 12.5 × 1010 M (also confirmed by Shao et al. 2017). Shao et al. (2017) estimated a BH mass MBH = (2.7 ± 0.6) 109 M from the MgII line, which contributes 2% of the dynamical mass of the system. Carniani et al. (2019) performed a detailed study of the cold-dust SED of J1319+1950 deriving a dust mass of log(Mdust/M) = 8.8 ± 0.2 and a dust temperature of T dust = 66 10 + 15 $ T_{\rm dust}=66^{+15}_{-10} $ K. Herrera-Camus et al. (2020) tentatively detected the OH 119 μm doublet in absorption, which is blueshifted with a median velocity that suggests the presence of a molecular outflow, although characterized by a modest molecular mass loss rate of ∼ 200 M yr−1.

A.10. PSO J183.1124+05.0926

PSO J183.1124+05.0926 (hereafter J183+05) at z = 6.439 was discovered by Bañados et al. (2018) using color-color selections from the Pan-STARRS1 database (Chambers et al. 2016) and follow-up photometric and spectroscopic observations. The [CII] luminosity in this source is the highest among 27 quasars at z > 6 surveyed in (Decarli et al. 2018). J183+05 showed a clear velocity gradient however, given the resolution of the observation, only an upper limit to dynamical mass has been derived Mdyn > 1.3 × 1011 M from the study of its rotation curve (Neeleman et al. 2021). It has a BH mass of MBH = 3.0 × 109 M derived from MgII emission line. Recently, Decarli et al. (2023) analyzed the cold-dust SED of this QSO in detail, deriving Tdust = 47 ± 2 K, Mdust = (8.7 ± 1.1) × 108 M and SFR = 1330 M yr−1. They also presented a multi-line study of this object, reporting detections of [CII], [OIII], [NII], CO(7-6), OH and two H2O transitions.

Appendix B: Analysis of continuum emission from archival observations

B.1. QSO J036+03

We analyzed the three ALMA observations available in B6 and B7 for J036+03, and we derived continuum flux densities and sizes at 243.11 GHz, 260.53 GHz and 338.71 GHz, performing a 2D Gaussian fit with CASA for each observation. All continuum emissions at these frequencies are spatially resolved and the values for flux densities, peak fluxes and sizes are reported in Tab. C.1.

In order not to miss the fainter and more extended flux, we tapered the higher resolution observations at 243.1 GHz and 260.5 GHz, performing the imaging with uvtaper=[0.9 arcsec], reaching a resolution of ∼ 0.7 × 0.7 arcsec2 and of ∼ 0.88 × 0.83 arcsec2 respectively. The sources in the tapered maps were fitted with the CASA 2D Gaussian fit, and the new continuum fluxes are reported in bold within brackets in Tab. C.1. At 243.1 GHz, the higher resolution observation missed ∼ 20% of the flux. Similar flux losses were also seen by Wang et al. (2019b) when analyzing high resolution observations of QSO J0100+2802. We checked that there was no further flux gain if tapering at even lower resolution.

B.2. QSO J0224-4711

We analyzed the three ALMA observations available in B3 and B6 for J0224-4711, and we derived continuum flux densities and sizes at 95.33 GHz, 245.67 GHz, and 260.51 GHz, performing a 2D Gaussian fit with CASA for each observation. The two B6 continuum emissions are spatially resolved and the values for flux densities, peak fluxes and sizes are reported in Tab. C.1. The emission in B3 is not resolved (see left panel of Fig. 3), therefore we considered the peak flux as total flux of the source, which is 0.103 ± 0.009 mJy/beam, analogously to B8 and B9.

Also for this object, we tapered the higher resolution observations at 245.67 GHz and 260.51 GHz, performing the imaging with uvtaper=[0.7 arcsec], reaching a resolution of ∼ 0.7 × 0.7 arcsec2. The sources in the tapered maps were fitted with the CASA 2D Gaussian fit, and we found no gain in flux even when tapering both observations at even lower resolution.

B.3. QSO J2054-0005

We analyzed the continuum emission at different frequencies of J2054-0005 from the archival observations reported in Tab. C.1. The analyses were carried analogously to the ones of J036+03 and J0224-4711, and the info about the fluxes and sizes obtained are reported in Tab. C.1.

We tapered the higher resolution observations at 92.26 GHz, 262.6 GHz, and 488.31 GHz using uvtaper=[0.7 arcsec], in order to account for the more extended and fainter emission. We achieved a resolution of 1.0 × 0.89 arcsec2, 0.88 × 0.82 arcsec2, and 0.88 × 0.81 arcsec2 for the 92.26 GHz, 262.6 GHz and 488.31 GHz observation, respectively. We did not find any further emission for the lowest frequency observation, while we gained ∼ 5% of the flux in the other two observations.

B.4. QSO J231-20

Fig. B.1 shows the B8 continuum emission of QSO J231-20. Details on the analysis can be found in Sect. 4.1.4.

thumbnail Fig. B.1.

406 GHz dust continuum map of QSO J231-20 (levels −3, −2, 2, 3, 5, 10, and 15σ, σ = 0.5 mJy/beam). The clean beam (4.3 × 2.9 arcsec2, PA=-85.8°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak.

Appendix C: Tables

Table C.1.

Continuum emission at different frequencies for J036+03, J0224-4711, J231-20, and J2054-0005

Table C.2.

Details of the observations targeting the CO line emissions

Table C.3.

Results of the SED fitting

All Tables

Table 1.

General properties of the QSOs in our sample.

Table 2.

Line properties of QSOs J0224−4711, J1319+0950, and J2054−0005.

Table 3.

Properties of QSOs in our sample.

Table C.1.

Continuum emission at different frequencies for J036+03, J0224-4711, J231-20, and J2054-0005

Table C.2.

Details of the observations targeting the CO line emissions

Table C.3.

Results of the SED fitting

All Figures

thumbnail Fig. 1.

Dust continuum maps. Top left panel: 404.9 GHz dust continuum map of QSO J036+03 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.6 mJy/beam). The clean beam (2.10 × 1.43 arcsec2, PA = 37.40°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Top right panel: 670.9 GHz dust continuum map of QSO J036+03 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.5 mJy/beam). The clean beam (1.12 × 0.97 arcsec2, PA = 21.11°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Bottom left panel: 405.2 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.88 mJy/beam). The clean beam (3.87 × 2.79 arcsec2, PA = −76.45°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Bottom right panel: 670.9 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, 8, and 11σ, σ = 1.3 mJy/beam). The clean beam (1.08 × 0.95 arcsec2, PA = 14.05°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak.

In the text
thumbnail Fig. 2.

Dust continuum and CO(6−5) emission line maps of QSOs J1319+0950 and J2054−0005. Top left panel: 103 GHz dust continuum map of QSO J1319+0950 (levels −3, −2, 2, 3, 5, 10, and 20σ, σ = 0.005 mJy/beam). The clean beam (0.30 × 0.30 arcsec2, PA = −78.56°) is indicated in the lower left corner of the diagram. Top right panel: CO(6−5) map of QSO J1319+0950 (levels −3, −2, 2, 3, 5, 10, and 15σ, σ = 0.017 mJy/beam). The clean beam (0.32 × 0.31 arcsec2, PA = −78.56°) is indicated in the lower left corner of the diagram. Bottom left panel: 92 GHz dust continuum map of QSO J2054−0005 (levels −3, −2, 2, 3, 5, 7 and 10σ, σ = 0.006 mJy/beam). The clean beam (0.42 × 0.32 arcsec2, PA = −61.3°) is indicated in the lower left corner of the diagram. Bottom right panel: CO(6−5) map of QSO J2054−0005 (levels −3, −2, 2, 3, 5, 7 and 10σ, σ = 0.025 mJy/beam). The clean beam (0.39 × 0.30 arcsec2, PA = −61.4°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak for each source.

In the text
thumbnail Fig. 3.

Dust continuum and emission line maps of QSO J0224−4711. Left panel: 95.33 GHz dust continuum map of QSO J0224−4711 (levels −3, −2, 2, 3, and 5σ, σ = 0.016 mJy/beam). The clean beam (3.95 × 2.32 arcsec2, PA = −89.43°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak. Central panel: CO(7−6) emission line map of QSO J0224−4711 (levels −3, −2, 2, 3, 5, and 8σ, σ = 0.042 mJy/beam). The clean beam (3.49 × 2.05 arcsec2, PA = −89.06°) is indicated in the lower left corner of the diagram. The cross indicates the position of the B3 continuum peak. Right panel: [CI] emission line map of QSO J0224−4711 (levels −3, −2, 1, 2, and 3σ, σ = 0.054 mJy/beam). The clean beam (3.49 × 2.05 arcsec2, PA = −89.06°) is indicated in the lower left corner of the diagram. The cross indicates the position of the B3 continuum peak.

In the text
thumbnail Fig. 4.

Moment maps of the CO(7−6) emission line and spectrum of CO(7−6) and [CI] emission lines of J0224−4711. From left to right: Integrated flux, mean velocity map, velocity dispersion map, and continuum-subtracted spectra of CO(7−6) and [CI]. The clean beam is plotted in the lower left corner of the moment maps. The cross indicates the peak position of the Band 3 continuum emission. The spectrum was extracted from the region included within ≥2σ in the CO(7−6) map.

In the text
thumbnail Fig. 5.

Moment maps of the CO(6−5) emission lines and spectra of J1319+0950 (top row), J2054−0005 (bottom row). From left to right: integrated flux, mean velocity map, and velocity dispersion map, continuum-subtracted spectrum of CO(6−5). The clean beam is plotted in the lower left corner of the moment maps. The cross indicates the peak position of CO(6−5). The spectra were extracted from the region included within ≥2σ in the CO(6−5) map of each source.

In the text
thumbnail Fig. 6.

Observed SEDs and best-fitting results. Left panels: Observed SEDs of QSOs J036+03, J0224-4711, and J231−20 (top, central, and bottom row). Our new ALMA B8 and B9 data are shown as yellow stars, and other archival observations are plotted as cyan diamonds. The best-fitting curve is shown as the solid blue line. Right panels: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid orange lines indicate the best-fitting value for each parameter, and the dashed lines mark the 16th and 84th percentiles for each parameter. For J231−20, the empty diamond is the flux in B8, not corrected for the presence of the companion, and the dashed line is the best-fitting SED considering this flux (for details, see Sect. 5.1).

In the text
thumbnail Fig. 7.

Same as Fig. 6 for QSO J2054−0005.

In the text
thumbnail Fig. 8.

Observed SEDs and best-fitting results. Left panels: Observed SEDs of QSOs J183+05, J1342+0928 (top and bottom row) as cyan diamonds. The best-fitting curve is shown as a solid blue line. Right panels: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid orange lines indicate the best-fitting value for each parameter, and the dashed lines mark the 16th and 84th percentiles for each parameter.

In the text
thumbnail Fig. 9.

Compilation of cold-dust SEDs of the six QSOs analyzed in this work, normalized at 850 μm.

In the text
thumbnail Fig. 10.

Results for the observed mean cold-dust SED. Left panel: Observed mean fluxes and best-fitting function. The dashed gray lines mark the six bins we used to develop the mean SED. Right panel: Corner plot showing the posterior probability distributions of Tdust, Mdust, and β. The solid green lines indicate the best-fitting value for each parameter, while the dashed lines mark the 16th and 84th percentiles for each parameter.

In the text
thumbnail Fig. 11.

CO SLED of J036+03 and J2054−0005 compared with those of other QSOs at z > 6 and at lower redshift. The CO SLEDs for J036+03 and J2054−0005 are shown as pink circles and cyan markers, respectively. For J2054−0005, the star refers to the CO(6−5) line studied in this work and the circle to the CO(7−6) in Decarli et al. (2022). The results for J036+03 are taken from Decarli et al. (2022). All other sources are displayed as shadowed squares: J1007+2115 at z = 7.5419 in red (Feruglio et al. 2023); J0439+1634 at z = 6.511 in gray (Yang et al. 2019); J1148+5251 at z = 6.419 in brown (Riechers et al. 2009; Gallerani et al. 2014); J0100+2802 at z = 6.327 in purple (Wang et al. 2019b); J2310+1855 at z = 6.003 in blue (Li et al. 2020); APM08279+5255 at z = 3.911 in green (Papadopoulos et al. 2001; Weiß et al. 2007); and Cloverleaf at z = 2.560 in orange (Bradford et al. 2009; Uzgil et al. 2016).

In the text
thumbnail Fig. 12.

Dust temperature as a function of redshift. Our sample is shown as stars (red for Hyperion QSOs, green otherwise). Top panel: Comparison of our results with samples of local QSOs (PG, blue colormap with contours, Petric et al. 2015), local star-forming galaxies (JINGLE, mean value as pink dot, Lamperti et al. 2019), local ULIRGs (HERUS, mean value as orange dot, Clements et al. 2018), 2 < z < 7 star-forming galaxies (SPT, pink colormap with contours Reuter et al. 2020), 2 < z < 5 QSOs (WISSH, yellow squares, Duras et al. 2017), 4 < z < 7.5 SMGs and star-forming galaxies (light blue dots, Witstok et al. 2023), three z > 6.7 REBELS galaxies (light blue diamonds, Algera et al. 2023), an average-REBELS galaxy (light blue triangle, Sommovigo et al. 2022), and two galaxies with very low dust temperatures at z > 6 (black crosses, Witstok et al. 2022; Harikane et al. 2020). The observed trends inferred by Viero et al. (2022) and Schreiber et al. (2018) are shown as an orange and green lines and shaded regions, respectively. Our best-fitting curve is shown as a purple line with a shaded region. The theoretical relation found in Sommovigo et al. (2022) is shown as a dashed gray line. The solid black line marks the CMB temperature level. Bottom panel: Focus on the comparison to the samples of local PG QSOs (blue color map with contours) and 2 < z < 5 WISSH QSOs. Our best-fitting curve considering QSOs alone is shown as a purple line with a shaded region.

In the text
thumbnail Fig. 13.

Dust emissivity index, β, as a function of redshift. The symbols and colors are the same as in Fig. 12. For the PG, WISSH, and SPT samples, the symbols are empty to indicate that the value of β in these cases was fixed and not derived from SED analyses.

In the text
thumbnail Fig. 14.

Star formation rate as a function of dust mass. The symbols and colors are the same as in Fig. 12. The SFRs for the WISSH and our sample are corrected for a factor of 50% to take the AGN contribution to the dust heating into account.

In the text
thumbnail Fig. 15.

SFR as a function of the molecular gas mass, MH2. The symbols for our sample are the same as in Fig. 12. We compare our results with the WISSH sample (Bischetti et al. 2021), with 1 < z < 1.5 QSOs (pink dots), and with z ∼ 0 sources such as galaxies, ULIRGs, QSOs, AGNs, and narrow-line Seyfert 1 galaxies (symbols and colors in the legend; Saintonge et al. 2011, 2017; Salomé et al. 2023; Shangguan et al. 2020a,b; Rosario et al. 2018; Koss et al. 2021). The gray lines represent fixed values of the gas depletion time (i.e., the inverse of SFE), which is reported at the top of the line. The solid black line is the galaxy main sequence derived by Sargent et al. (2014) considering massive galaxies (M* > 1010 M) up to z ∼ 2.

In the text
thumbnail Fig. 16.

Redshift distribution of the GDR. The symbols and colors for the WISSH and our sample are the same as in Fig. 12. We compare our results with the WISSH sample (Bischetti et al. 2021), a sample of 2 < z < 5 star-forming galaxies hosting a heavily obscured AGN in the Chandra Deep Field-South (magenta dots, D’Amato et al. 2020), and a sample of z > 5.5 QSOs (blue squares Calura et al. 2014). The two vertical lines at the bottom right side are the systematic uncertainties induced by the choice of αCO and dust mass estimation (black line, ∼0.3 dex) and r65 or r76 in computing the gas mass (brown line, ∼0.1 dex).

In the text
thumbnail Fig. 17.

BH mass vs. dynamical mass for our sample (stars with green contours for HYPERION QSOs, and red contours otherwise), WISSH QSOs at z ∼ 2 − 4 (gray diamonds; Bischetti et al. 2021), and luminous z ∼ 4 − 7 QSOs (gray dots and gray squares; Venemans et al. 2016, 2017; Willott et al. 2013, 2015, 2017; Kimball et al. 2015; Trakhtenbrot et al. 2017; Feruglio et al. 2018; Mortlock et al. 2011; De Rosa et al. 2014; Kashikawa et al. 2015; Neeleman et al. 2021). The dashed black line (and shaded area) is the local relation from Kormendy & Ho (2013). The light green dots are the remaining HYPERION QSOs for which we were not yet able to perform a detailed study of the dust properties due to a lack of observations in the submm regime. The stars are color-coded based on the value of Ggal. The thin red arrows indicate upper limits on the dynamical mass.

In the text
thumbnail Fig. B.1.

406 GHz dust continuum map of QSO J231-20 (levels −3, −2, 2, 3, 5, 10, and 15σ, σ = 0.5 mJy/beam). The clean beam (4.3 × 2.9 arcsec2, PA=-85.8°) is indicated in the lower left corner of the diagram. The cross indicates the position of the continuum peak.

In the text

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

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

Initial download of the metrics may take a while.