EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A23 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A23
Number of page(s) 15
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202453226
Published online 27 February 2025

© The Authors 2025

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

The significant role of quasars in reionizing the Universe at z > 6 has been known for many years (Gunn & Peterson 1965). However, our understanding has advanced considerably over the past couple of decades, paralleling the growth in the number of known quasars, which now includes approximately 300 identified at the epoch of reionization (EoR; Fan et al. 2023). This has been achieved thanks to the combination of near-infrared photometric surveys and spectroscopic follow-up observations, which have proven to be reliable tools for discovering the quasar population at high redshifts (z > 5; e.g., Bañados et al. 2016; Shen et al. 2019).

Luminous quasars at the EoR are also crucial probes for testing the coevolution between supermassive black holes (SMBHs) and their host galaxies. In this regard, modern interferometric facilities have been essential for confirming the redshift of the objects and constraining the star formation activity and properties of the cold phase of the interstellar medium (ISM) of their host galaxies (e.g., Riechers et al. 2006, 2007, 2009; Wang et al. 2007, 2010). Far-infrared (FIR)/submillimeter observations have revealed the presence of highly star-forming host galaxies, with star formation rates (SFRs) of up to 1000–3000 M⊙ yr−1, and copious amount of dust (Mdust > 108 M⊙; Maiolino et al. 2005; Wang et al. 2013; Feruglio et al. 2018; Pensabene et al. 2022; Venemans et al. 2017a, 2020). These massive quasar hosts are crucial for benchmarking gas and dust content and star formation efficiencies (SFEs) during the mid-point and terminal stages of the EoR, but observations are still limited to a few tens of quasar hosts at these epochs. Furthermore, the scaling relation between gas and dust is not clear yet, and the methods that are currently adopted to derive molecular gas masses can be improved with more data. Given its brightness, the 158 μm emission line of the singly ionized carbon atom ([C II] hereafter) is the most used tracer of the ISM at z > 4. While the [C II] emission is undoubtedly useful for dynamical studies (e.g., Wang et al. 2013, 2024; Decarli et al. 2017; Venemans et al. 2019, 2020; Neeleman et al. 2021), it may not effectively trace the dense and cold components of the molecular gas, nor the star-forming regions (e.g., Pineda et al. 2013; Herrera-Camus et al. 2018; Neeleman et al. 2019; Bischetti et al. 2024; Izumi et al. 2024). On the other hand, the rotational excited transitions of carbon monoxide (CO) are much more reliable probes of the dense cold H2 and are pivotal when investigating the star formation process (Kaasinen et al. 2024). Despite being fainter than [C II], CO lines have been successfully detected in massive quasar hosts at z ∼ 6 (e.g., Wang et al. 2013; Decarli et al. 2022); however, they remain poorly sampled at the highest redshifts (z > 7) of quasars, with only one confirmed individual detection (Feruglio et al. 2023) and a possible detection suggested by the stacking of multiple lines (Venemans et al. 2017a; Novak et al. 2019) at z ∼ 7.1 − 7.5.

In this work we present new observations targeting CO(6–5) and (7–6) emission lines and the underlying continuum in five quasars at z > 7, obtained with the Northern Extended Millimeter Array (NOEMA), thus providing a complete census of the molecular gas in the population of z > 7 quasars known to date (Fan et al. 2023). We also investigated the star formation activity and dust properties for the five quasars by modeling the FIR spectral energy distribution (SED) using NOEMA and archival observations of the Atacama Large Millimeter/submillimeter Array (ALMA). Based on these measurements, we studied the growth of host galaxies and the implications for their coevolution with SMBHs in the reionization era.

This paper is organized as follows. Section 2 describes the observations and data reduction. Section 3 describes the procedures adopted to model the FIR SED and, in turn, measure the dust and cold molecular gas content, including a new calibration using the [CII] luminosity. In Sect. 4 we discuss the results and their implication for the coevolution of the host galaxy of quasars at the EoR, and compare them with objects at lower redshifts and the results from semi-analytic models (SAMs). We summarize the paper in Sect. 5.

Throughout the paper, we adopt a standard flat Λ cold dark matter cosmology with the matter density parameter ΩM = 0.30, the dark energy density parameter ΩΛ = 0.70, and the Hubble constant H0 = 70 km s−1.

2. Sample and observations

We focused on the quasars with z ≥ 7, which are eight in total according to Fan et al. (2023), aiming to obtain a complete census of the molecular gas content in these objects. The five quasars for which we acquired new NOEMA observations under project S23CX are listed in Table 1, while the entire list of z > 7 quasars with their names and identifications are reported in Table 2.

Table 1.

NOEMA observations.

Table 2.

Quasar properties.

Regarding project S23CS, receivers were tuned to detect the CO(6–5) and CO(7–6) and the underlying continuum in the lower and upper side band (LSB and USB), respectively. We calibrated the visibilities using the CLIC pipeline of the GILDAS software1. Imaging was performed with MAPPING in GILDAS, using a natural weighting scheme for visibilities, with a detection threshold equal to the noise of the pre-imaging visibilities. The observation targeting quasar J0313 was executed in bad weather conditions, so it did not reach the requested sensitivity to detect either continuum emission or CO lines in both LSB and USB. Given the low declination of J0038 (coordinates are listed in Table 1), the observation of this target was limited by shadowing, which reduced the total observing time on-source, hence the final sensitivity. In these two cases, we merged LSB and USB to increase the signal-to-noise ratio, and we used the merged cube to extract an upper limit on the continuum emission. Continuum-subtracted cubes were binned to spectral resolutions of 20 MHz, corresponding to ∼60 km s−1, to optimize the sensitivity.

To extend the sampling of the FIR emission of the quasar host galaxies, we collected archival ALMA observations for four of the five quasars2, which cover sky frequencies ∼220 − 240 GHz, where the [C II] emission line falls. ALMA-calibrated visibilities of observations covering [C II] and the underlying continuum were retrieved from the science archive (see the acknowledgements for the list of ALMA project IDs). The imaging was performed through the Common Astronomy Software Applications (CASA; McMullin et al. 2007), version 6.5.5–21. We ran tclean procedures using natural weighting, a 3σ cleaning threshold, and a channel width of 20 MHz. Continuum-subtracted cubes and continuum maps were created using the imfit procedure in CASA, fitting a constant function to all line-free channels.

The synthesized beams, root-mean-square (rms) noise levels, and representative frequencies of the continuum maps are reported in Table 1. Continuum maps of ALMA and NOEMA observations are shown in Fig. B.1.

3. Results

3.1. Dust properties and star formation rates

The continuum maps of the five targets are presented in Fig. B.1. We measured continuum flux densities of NOEMA maps with a fit of the visibilities with a point source model. The associated error accounts for both the statistical uncertainty of the fit and the calibration error (≲10%3). For J0252, J1243, and J2356, we also measured the continuum flux by integrating the 2D Gaussian function that fits the cleaned continuum emission map. The two approaches provide flux estimates that are consistent within the uncertainties. In the case of the NOEMA observations, the observing setup did not allow us to spatially resolve the target extension, all best-fit functions are consistent with a point-like source having a size comparable with the beam (see Table 1). Flux upper limits are estimated at a 3σ significance level using the rms noise and assuming a size equal to the beam size.

As visible in Fig. B.1, ALMA observations of J0038, J0313, J0252, and J1243 spatially resolve the extent of the continuum emission. This allowed us to measure the size of the continuum emitting region in quasar host galaxies by modeling the continuum maps with a 2D Gaussian function. In Table 1, we report the integrated flux and host-galaxy size, which is assumed to be equal to the full width half maximum (FWHM) of the best-fit 2D Gaussian function.

Using EOS-DUSTFIT4, we fitted the cold dust SED of the target galaxies. In EOS-DUSTFIT, following the prescription by (Draine & Li 2007, see also Carniani et al. 2019; Tripodi et al. 2024a), the cold dust continuum emission is modeled with a modified blackbody (MBB) function, assuming an optically thin regime. The MBB model can have up to three free parameters, namely dust temperature (Td), dust mass (Md), and the emissivity index (β). However, for the five sources of our sample, the limited amount of data prevents us from modeling the FIR SED with three free parameters. In particular, the peak of the SED is not sampled by current data, making it challenging to constrain the dust temperature. Consequently, we adopted a fixed temperature, Td = 55 K, for all five quasars based on the mean dust temperature in a sample of 10 z > 6 quasars (Tripodi et al. 2024a). This assumption is further supported by recent literature results on z > 6 quasars (e.g., Decarli et al. 2022, Sommovigo et al. 2022, Shao et al. 2022, Witstok et al. 2023, and Tripodi et al. 2023) and lower-z quasars with measured Td (e.g., Bischetti et al. 2021), as it would be expected in the case of compact, and dense star-forming regions (dust sizes of few kiloparsecs; e.g., Venemans et al. 2017a, Decarli et al. 2018, 2022, and Tripodi et al. 2024a). Their intense star formation still stands even after accounting for the active galactic nucleus (AGN) contribution to the dust heating, which can be substantial, especially in the nuclear region (Duras et al. 2017, Di Mascia et al. 2021, Walter et al. 2022, and Tsukui et al. 2023).

Combining NOEMA and ALMA detections in the observed-frame ∼80–230 GHz range, we estimated the dust emissivity parameter, β, for two objects in our sample (J0252 and J1243). For the others, we fixed β = 1.6, the local value (Beelen et al. 2006 and Witstok et al. 2023) and similar to that measured in J0252 and J1243. Since NOEMA observations do not allow us to resolve the source size, we assumed that the size does not vary significantly with frequency in the observed range ∼80 GHz to ∼230 GHz. Quasar J2356 only has an unresolved observation at 3 mm (Fig. B.1). Then, we assumed a size of 0.75 × 0.64 arcsec2, based on the median value of the FWHM dimensions measured for the four quasars of our sample with ALMA observations.

EOS-Dustfit explores the parameter space for each SED using a Markov chain Monte Carlo algorithm implemented in the emcee package (Foreman-Mackey et al. 2013). A uniform distribution for priors is assumed for fitting parameters in the range: 5 < log(Md/M⊙) < 9 for all the quasars, and 0.5 < β < 3.0 for both J0252 and J1243. We ran 40 chains, with 3000 trials and a burn-in phase of 150 steps for each dataset; we added in quadrature a 10% calibration uncertainty to the continuum flux errors. The upper limits of the continuum emission were included in the fit as a 1σ detection with a 2σ error. There is not yet a standard approach for dealing with upper limits in the fitting procedure, one could either treat them as detections with large errors (as we do; see also, e.g., Witstok et al. 2022 and Ronconi et al. 2024) or change the fitting code properly to ensure that above the upper limit level, the likelihood is zero. We adopt the first method given that it is commonly employed in many fitting routines (e.g., GalaPy and MERCURIUS) and gives us reasonable results. We adopted the 50th percentile of the posterior distribution as the best-fit value, while the errors are calculated considering the 16th and 84th percentiles (corner plots of the posterior distribution of free parameters are shown in Fig. A.2). In the case of the modeling of FIR SED assuming a fixed β, the uncertainties on Md are relatively small (see below) when compared to the fit where β is left free to vary. This is because of the marginalization over the distribution of β. To determine with accuracy to which extent the uncertainty on Md and LFIR may be underestimated due to marginalization on β, we performed additional fits for J0252 and J1243, this time setting β = 1.6. We then obtained uncertainties on Md and LFIR lower by ∼0.05 dex. A similar conclusion can be reached for the assumption of Td, leaving it free to vary over a limited range of temperature (Td = 40 − 70 K). In Sect. 3 we take care of this systematic by including an additional contribution to the error on Md and LFIR. The results from the SED fitting are reported in Table 2, while the best-fit models are shown in Fig. A.1. We measured the FIR luminosity (LFIR) by integrating the 40–1000 μm emission of the best-fit model produced by EOS-DUSTFIT in the five quasars. SFR is computed with the relation by Kennicutt (1998), SFR/(M⊙ yr−1) = 10−10 LFIR/L⊙, assuming a Chabrier initial mass function (IMF; Chabrier 2003). The same prescription is commonly assumed in high-z quasars (e.g., Duras et al. 2017, Bischetti et al. 2021, and Bertola et al. 2024), and is roughly consistent with the SFR derived assuming a Kroupa (2001) IMF (within 7%). A Chabrier IMF was also adopted for J1007 by Feruglio et al. (2023) and J1342 by Novak et al. (2019), while Venemans et al. (2017a) assumed a Kroupa IMF for J1120.

In two cases we modeled the dust emissivity index, for J1243 the best-fit value derived is fully consistent with the typical value observed in high-redshift quasars (Beelen et al. 2006, Venemans et al. 2017a, 2020, and Tripodi et al. 2024a). For J0252, we derived β = 0 . 93 − 0.21 + 0.20 $ \beta=0.93^{+0.20}_{-0.21} $, which is significantly lower (4σ) than β values measured in high-z quasars. A potential reason for the flat FIR SED in J0252 could be the presence of contaminant sources within the beam of the NOEMA observation (see Fig. B.1). Indeed, the ALMA continuum map of J0252 (see Fig. B.1) shows a marginally resolved source surrounded by elongated emission features, each detected with > 2σ significance level. These emission features may either hide the presence of close companion galaxies or interloper sources. However, the relatively low β value obtained from the fit does not significantly affect the Md estimate. Assuming β = 1.6 for J0252, we obtained a Md fully consistent with the one presented in Table 2.

We tested the consistency of our assumptions (namely, Td, β) on the results of the SED fitting. Using a higher (lower) dust temperature such as Td = 70 K (40 K) would have led to an average decrease (increase) in Md by a factor of ∼1.5 (∼0.18 dex), while the value of LFIR would have been higher (lower) by a factor of ∼2 (∼0.3 dex). Regarding β, assuming a steeper emissivity index (β = 2) for the targets with β fixed (namely, J0038, J0313, and J2356) would have resulted in a larger Md (lower LFIR) by a factor of ∼2. The results of the modeling assuming different β and Td values are shown in Fig. A.3. Conversely, choosing a flatter β (β = 1.2) results in a mean increase (decrease) in the value of Md (LFIR) by a factor of ≤1.2. Eventually, we tested the impact of assuming a different source size: increasing the source area up to a factor of 10 provided results consistent within the uncertainties with those listed in Table 2; if the source area was overestimated by a factor of 10 or more, this would have led to Md and LFIR being lower by a factor of ≥2.

The five quasars in our sample show Md in the range 0.6 − 4 × 108 M⊙ (see Fig. 1 and Table 2), consistent with the dust content derived in the population of luminous quasars at the EoR (e.g., Izumi et al. 2021; Witstok et al. 2023; Tripodi et al. 2024a) and at different cosmic epochs (e.g., Duras et al. 2017, Bischetti et al. 2021, and Bertola et al. 2024), and star-forming galaxies at lower redshifts (e.g., Mancini et al. 2015, Leśniewska & Michałowski 2019, Pozzi et al. 2021, Hygate et al. 2023, Algera et al. 2024). As shown in Fig. 1, comparing quasars at z > 7 and sources at later epochs, there is no clear evolution of Md as a function of the redshift. The quasars in our sample have already built large dust masses in a relatively short time (the Hubble time at z ∼ 7 is ∼800 Myr), suggesting that in high-z quasar hosts the physical processes that drive the formation of dust grains are very efficient and overtake those processes destroying dust particles (e.g., Popping et al. 2017). We continue the discussion on the dust formation processes in Sect. 4.1.

thumbnail Fig. 1.

Dust mass and molecular gas mass vs. redshift. Upper panel:Md vs. redshift. Red circles are the z > 7 quasars presented in this work. Green circles are for J1120 from Venemans et al. (2017a), J1007 from Feruglio et al. (2023), and J1342 from Novak et al. (2019). Gray circles are the 5 < z < 7 quasars from Venemans et al. (2017b), Izumi et al. (2021), Decarli et al. (2022), and Tripodi et al. (2024a), while blue dots are for quasars at intermediate redshifts (1 < z < 5) from Bischetti et al. (2021) and Salvestrini et al. (in prep.). Lower panel:MH2 vs. redshift. Data are color-coded as in the upper panel. We also include the quasars at intermediate redshifts (1 < z < 5) from Bertola et al. (2024).

3.2. Molecular gas mass

For four out of five quasars observed with NOEMA we derived 3σ upper limits on the cold H2 mass using the rms at the expected frequency of CO(6–5) and (7–6) emission lines. We assumed a line width of FWHM = 300 km/s, as expected for z > 6 quasars (e.g., Decarli et al. 2022 and Feruglio et al. 2023), and a galaxy size equal to the beam of the observations. We adopted the CO spectral line energy distribution (SLED) correction from Kaasinen et al. (2024), that is, r6, 1 = 0.92 for CO(6–5) and r7, 1 = 0.65. Assuming a luminosity to mass conversion factor αco = 0.8 M⊙ (K km s−1 pc2)−1 (see Bolatto et al. 2013 and Carilli & Walter 2013), we derive an upper limit MH2 in the range ∼0.9 − 1.5 × 1010 M⊙ (see Table 2). The same αco and CO SLED correction factor are adopted for the upper limits and estimate of the H2 masses of J1007, J1120, J1342, for which we collected CO luminosity from the literature and included in Table 2. Of the known quasars at z ≳ 7, five have only upper limits on their H2 masses as traced by CO emission lines, with J1007 being the sole source with a statistically significant detection (> 6σ) of molecular gas mass (Feruglio et al. 2023). The molecular gas mass and the upper limits for the z > 7 quasars (see Table 2 and Fig. 1) are consistent with the mean value measured in the population of 5 < z < 7 quasars (log(MH2/M⊙)∼10 ± 0.3; e.g., Venemans et al. 2017b; Decarli et al. 2022; Kaasinen et al. 2024).

To test the consistency of our MH2 upper limits, we compared the CO-derived measurement with that obtained from the [CII]158 μm ([CII], hereafter) luminosity. L[CII] is a viable way to derive molecular gas masses for z ≳ 4 non-AGN galaxies given its brightness and scaling relation between L[CII] and MH2 calibrated on lower redshift samples. [CII] is expected to trace multiple gas phases (e.g., Maio et al. 2022, Maio & Viel 2023, and Casavecchia et al. 2025), but the corresponding emission is mostly due to the star formation activity in the host. As shown in the left panel of Fig. 2, L[CII] is tightly correlated (Pearson correlation coefficient ∼0.71) with the SFR in quasar host galaxies at z > 6. The data points shown in Fig. 2 include all quasar at z > 7 and a collection of objects at z ≳ 6 from the literature with both CO and [CII] measurements (see Table C.1 for the full list of objects and references). For this reason, several proposed relations use the L[CII] to trace the molecular gas content based on the Kennicutt-Schmidt relation (Kennicutt 1998). As an example, the relation by Zanella et al. (2018) calibrated on MS and starburst galaxies at redshift up to z ∼ 6, with the bulk of the sample at z ∼ 2. Alternatively, the relation by Madden et al. (2020) allows us to measure the total MH2, including the potential contribution due to CO-dark clouds, and is calibrated on local dwarf galaxies. As shown in the right panel of Fig. 2, both relations overestimate significantly the MH2 measurement and upper limit derived from L′CO in quasars at z > 5. Indeed, several works in the literature studying the ISM properties of z > 6 quasars (e.g., Neeleman et al. 2021, Decarli et al. 2022, Kaasinen et al. 2024, and Tripodi et al. 2024a) observed a large discrepancy between MH2 derived with the prescription by Zanella et al. (2018) and Madden et al. (2020) and CO-derived ones, with [CII]-based MH2 that overestimates CO-based MH2 from a factor of a few to one order of magnitude (e.g., Kaasinen et al. 2024). The discrepancies in the determination of MH2 from L[CII] are likely due to the different physical properties of the ISM between the sample collected by Zanella et al. (2018) and Madden et al. (2020) and the host galaxies of high-z quasars. Indeed, when considering galaxies at high redshifts (z ≳ 4) with an intense star formation activity (Rizzo et al. 2020, 2021), the measured MH2-to-L[CII] ratio is considerably lower (∼4 − 8 M⊙/L⊙) than the value ∼30 M⊙/L⊙ measured by Zanella et al. (2018). This is also supported by the results of post-processed FIR emission lines from the SIMBA cosmological hydrodynamical simulations (Vizgan et al. 2022), suggesting that the MH2-to-L[CII] ratio should be lowered by up to a factor of ∼2 for synthetic star-forming galaxies at z ≃ 6.

thumbnail Fig. 2.

SFR and molecular gas mass vs. [CII] luminosity. Left panel:z > 7 quasars with L[CII] measurements: J1234, J0038, and J0252 (red circles) and J1120, J1342, and J1007 (green circle). Gray dots are the z ≳ 6 quasars from the literature (see Table C.1). The best-fit relation and the relative 68% confidence interval are represented as the solid black line and green-shaded region, respectively. The best-fit parameters of the relation are: slope 1 . 15 − 0.27 + 0.28 $ 1.15^{+0.28}_{-0.27} $, normalization − 8 . 16 − 0.06 + 0.06 $ -8.16^{+0.06}_{-0.06} $, and intrinsic dispersion 0 . 07 − 0.02 + 0.03 $ 0.07^{+0.03}_{-0.02} $ dex. The dotted line is the relation from De Looze et al. (2014) for low-metallicity dwarf galaxies in the local Universe. The dashed line is the best-fit relation from Romano et al. (2022) for galaxies from the ALPINE survey (4 < z < 6), including both detections and stacked non-detections. Right panel:z > 7 quasars with L[CII] measurements and MH2 upper limits: J1234, J0038, and J0252 (red circles) and J1342 (green circle); J1007 has both L[CII] and MH2 measurements. MH2 measurements and upper limits for z ≳ 6 quasars from the literature have been homogenized to those presented in this work by assuming a common αCO of 0.8 M⊙ pc−2 (K km s−1)−1 and a CO SLED. The best-fit relation and the relative 68% confidence interval are represented as the solid black line and gray-shaded region, respectively. The best-fit intrinsic dispersion of the relation is δ intr = 0 . 10 − 0.03 + 0.06 $ \delta_{\mathrm{intr}}=0.10^{+0.06}_{-0.03} $ dex. The dashed line is the scaling relation derived by Zanella et al. (2018) for main-sequence (MS) and luminous infrared galaxies, up to z ∼ 6. We also report the scaling relation provided by Madden et al. (2020) for local dwarf galaxies (dotted line).

Given that, we combined L[CII] and MH2 measurements and upper limits for the z > 7 quasars presented in Table 2 with measurements for 21 z ≳ 6 quasars from the literature (see Table C.1) to provide a new calibration for the MH2–L[CII] relation. We fitted data points of the resulting quasar sample of 25 objects with a linear regression based on a Bayesian approach, using the Python package linmix (Kelly 2007). This package allows us to consider errors on both variables, L[CII] and MH2, and upper limits on MH2. The best-fit relation is

log ( M H 2 / M ⊙ ) = ( 0 . 75 − 0.31 + 0.31 log ( L [ CII ] ) + ( 2 . 87 − 0.07 + 0.07 ) $$ \begin{aligned} \log (M_{\rm H2}/M_{\odot }) = (0.75^{+0.31}_{-0.31} \log (L_{\rm [CII]}) + (2.87^{+0.07}_{-0.07}) \end{aligned} $$(1)

and is shown in Fig. 2. The same relation corresponds roughly to a MH2-to-L[CII] ratio of ∼ 3 . 9 − 2.3 + 1.7 $ {\sim}3.9^{+1.7}_{-2.3} $ in the L[CII] regime 109 − 10 L⊙, almost one order of magnitude lower than that predicted by Zanella et al. (2018). We note that even if we considered a Milky Way like CO–H2 conversion factor (αCO = 4.3 M⊙ pc−2 (K km s−1)−1) to derive MH2 for high-z quasars, the corresponding MH2 would be lower by a factor of 1.5–2 than those derived assuming the calibration of Zanella et al. (2018). This relation provides MH2 estimates for luminous quasars (Lbol > 1046 erg/s) at z > 5 that are way more accurate than those provided by scaling relation calibrated from star-forming or dwarf galaxies at lower redshifts.

While applying this relation to high-redshift quasars, it is crucial to consider the following factors: (i) large-scale ionized winds driven by the accreting SMBH could significantly boost the observed L[CII] (e.g., Bischetti et al. 2024; ii) tidally stripped gas due to merging companions close to the quasar (e.g., Tripodi et al. 2024b; Decarli et al. 2024) could enhance L[CII] measurements, especially in observations with limited angular resolution. Both (i) and (ii) cases could significantly overestimate the molecular gas masses in the quasar’s host galaxy.

4. Discussion

4.1. Dust enrichment at high redshifts

The origin of large dust reservoirs in the early Universe is widely debated (see Schneider & Maiolino 2024 for a comprehensive review on this topic), since there is no consensus on the dominant formation mechanisms. However, among the different channels for forming dust in the first few hundred million years of the Universe, asymptotic giant branch (AGB) stars and supernova (SN) ejecta are expected to contribute significantly even at z > 7 (e.g., Sommovigo et al. 2022). Even considering the maximally efficient yield from SN ejecta, this process can produce dust masses up to a few times 107 M⊙ at the redshift of our targets (Mancini et al. 2015; Valiante et al. 2011). Nonetheless, we must also consider that a significant fraction of the dust created by SN ejecta (from 20% up to the total amount, depending on model assumptions; e.g., Micelotta et al. 2018 and Kirchschlager et al. 2019) could be destroyed by the reverse shocks that may follow the initial explosion (Schneider & Maiolino 2024). Regarding AGB stars, their contribution to dust production depends mostly on the gas metallicity and the IMF (e.g., Ventura et al. 2018; Dell’Agli et al. 2019), nevertheless it can reach up to 40% of the total dust at z ≳ 6 (e.g., Schneider & Maiolino 2024).

Moreover, two more processes are expected to significantly contribute to dust production in a pre-enriched ISM: aggregation and evolution of grains within the ISM, and quasar-driven winds. The aggregation of dust particles within the ISM usually requires timescales on the order of 1 Gyr and is observed to contribute up to 20–50% of the total dust production in the local Universe (Saintonge et al. 2018; Galliano et al. 2021). However, the extremely dense and turbulent ISM in the host galaxy of quasars (Gallerani et al. 2010, Valiante et al. 2011, 2014, Mancini et al. 2015, Decarli et al. 2023) could favor the formation of massive dust particles in a fraction of that time, making it the dominant channel in massive systems (e.g., Mancini et al. 2015). Indeed, local studies suggest that the aggregation of dust grains in the ISM primarily depends on gas density and temperature (e.g., Draine 2003, 2009, and Galliano et al. 2021). Given the luminosity of the quasar at the EoR included in this study, we cannot exclude that a non-negligible amount of dust grains that are produced and continue to grow within the outflowing clouds that are ejected due to the SMBH feedback (Elvis et al. 2002). Even if the role of this production channel is still debated (e.g., Pipino et al. 2011), recent results from hydrodynamical simulations (Sarangi et al. 2019) suggest that this mechanism can be quite efficient (up to few M⊙ yr−1 for a MBH = 108 M⊙) in luminous quasars such as our targets. This means this channel can contribute to at least a few percent of the total dust budget, considering an AGN cycle of a million years. Although this makes dust production in quasar winds less significant than the other mechanisms just described, it could still play a crucial role in dust enrichment in nuclear regions. To conclude, building dust masses that exceed 108 M⊙ in the host galaxies of luminous quasars at redshift z > 7 requires a combination of the physical processes described above.

As discussed above, the amount of dust traces the metal enrichment of the ISM. For this reason, the gas-to-dust ratio (GDR = MH2/Md) is a crucial parameter to understand the rapid evolution and growth of galaxies at EoR. Here, we adopted the MH2 measurement and upper limits listed in Table 2 to derive the GDR for 7 out of 8 quasars at z > 7, among which J1007 is the only object with a GDR estimate. We similarly derived GDRs for the quasars listed in Table C.1 and quasars with similar luminosity at Cosmic Noon from Bischetti et al. (2021). Since the high-z quasar population with dust mass estimates are biased toward luminous objects at FIR wavelength, we also considered a sample of local Seyfert galaxies (LFIR ∼ 1010 − 12 L⊙) to check for any potential selection bias. In Fig. 3, the resulting GDRs are shown as a function of Md. Globally, the population of quasars at the EoR shows GDRs that are consistent (GDRmean = 80 ± 40) with the local value GDR = 100 (Draine & Li 2007), similar to that measured in local star-forming galaxies (e.g., Casasola et al. 2020). This local value is also commonly adopted at high redshifts to derive the molecular gas mass given the dust mass (e.g., Neeleman et al. 2021). Quasars at z > 7 show upper limits on GDR that are consistent or just below the local value. J0252 shows the tighter constraint of GDR among z > 7 targets, but it is still consistent with GDR values observed in similarly luminous quasars. We recall that assuming a higher (lower) Td value would result in a higher (lower) GDR up to a factor of ∼2 with Td = 65 K (45 K). Looking at Fig. 3, GDR decreases at increasing Md, with a similar trend in quasars at different epochs and Seyfert galaxies. Indeed, by fitting the GDR versus Md with LINMIX, we find a negative slope (=−0.41 ± 0.16) when considering only quasars at the EoR. Consistent results can be obtained by including either lower-redshift quasars (Bischetti et al. 2021) or local Seyfert galaxies (Salvestrini et al. 2022; in the former case, there is a shift in the normalization). The decreasing trend is likely driven by the level of metal enrichment of the ISM, which is higher at higher dust reservoirs. This confirms that dust-rich systems are already enriched in metals, which are key ingredients of dust particles (Rémy-Ruyer et al. 2014). However, we recall that our MH2 estimates are derived with a constant αCO factor, while an assumption of a αCO dependence on metallicity (e.g., Amorín et al. 2016) would have balanced out the GDR anticorrelation with Md. We conclude that quasars at the EoR exhibit GDR values (both detection and upper limits) consistent with those measured in quasars at lower redshifts and in local objects. This confirms that assuming a local GDR value (Draine & Li 2007) is a reasonable choice for quasar host galaxies up to z ∼ 7.5.

thumbnail Fig. 3.

GDR vs. dust mass. The horizontal black line represents GDR = 100. Red circles are the z > 7 quasars presented in this work with MH2 upper limits. Green circles are for J1007 from Feruglio et al. (2023), J1120 from Venemans et al. (2017a), and J1342 from Novak et al. (2019). Gray circles are the 6 ≲ z < 7 quasars from Table C.1, while blue dots are for quasars at lower redshifts from Bischetti et al. (2021), Decarli et al. (2022), and Salvestrini et al. (in prep.). Local Seyfert galaxies are shown as red diamonds.

4.2. Gas depletion and star formation efficiency

In the top-left panel of Fig. 4, we show the depletion time (tdep = MH2/SFR) as a function of redshift. As a comparison, we include a collection of quasars spanning a wide range of redshift, namely: quasars at 6 ≲ z < 7 (see Table C.1), luminous objects at Cosmic Noon from Bischetti et al. (2018, 2021), Decarli et al. (2022), Bertola et al. (2024), and Salvestrini et al. in prep., and in the nearby Universe (Bischetti et al. 2019a; Shangguan et al. 2020, PDS456 and the PG survey, respectively). We also extend the comparison to AGN at intermediate luminosities (Lbol ∼ 1044 − 45 erg/s), with samples of local AGN and Seyfert galaxies from Koss et al. (2021) and Salvestrini et al. (2022) and highly accreting AGN hosts from Salomé et al. (2023). Eventually, we added star-forming and luminous infrared galaxies with CO detections from Tacconi et al. (2018), which cover a wide range of redshift, from the nearby Universe up to redshift z ∼ 4. As is clearly visible, the quasars at the EoR show relatively low tdep values (≲0.1 Gyr), with few objects reaching depletion times of a few tens of megayears. If we limit the analysis to the EoR, there is no evidence of the evolution of tdep in quasar host galaxies.

thumbnail Fig. 4.

Molecular gas consumption timescales and SFEs. Upper-left panel: Depletion time as a function of (1 + z), shown in logarithmic scale. Sources are color-coded as follows: the z > 7 quasars are red (this work) and green (Venemans et al. 2017a, Novak et al. 2019, and Feruglio et al. 2023). Literature samples include: z < 4 star-forming galaxies (pink-shaded region; Tacconi et al. 2018); local Seyfert galaxies and quasars (red diamonds; Bischetti et al. 2019a, Shangguan et al. 2020, Salvestrini et al. 2020, 2022, and Salomé et al. 2023); quasars at Cosmic Noon (blue dots; Bischetti et al. 2021 and Bertola et al. 2024); and quasars at 6 ≲ z < 7 (gray circles; see Table C.1). The parameter space occupied by the simulated quasars from GAEA is shown as an orange-shaded region, with 10% isodensity lines. The dashed line represents the evolution of the depletion time with redshift (Sommovigo et al. 2022). Upper-right panel: SFE displayed as a function of LIR. The best-fit lines are shown for illustration purposes only and were obtained from a linear fit of the data shown in the legend. Data points are color-coded as in the upper-left panel. Lower-left panel: SFE vs. Lbol. The black diamonds represent the median SFE value in each of the equally spaced Lbol bins. The vertical lines represent the distance between the 84th and 16th percentile of the SFE distribution in each bin. Data points are color-coded as in the upper-left panel. Lower-right panel: The histograms of the SFE distribution of each sample. The dashed line represents each distribution’s median value.

To further investigate this scenario, we contrasted our findings against theoretical predictions from the GAlaxy Evolution and Assembly (GAEA) model. In particular, the latest rendition of the model (De Lucia et al. 2024) combines the explicit partitioning of the cold gas into its molecular and neutral phases (the molecular gas ratio depends on the mid-plane pressure Blitz & Rosolowsky 2006, more details in Xie et al. 2017), with improved modeling of AGN activity (both in terms of cold gas accretion onto SMBHs and AGN-driven feedback Fontanot et al. 2020). The GAEA SAM has recently been coupled (Fontanot et al. in prep.) to the P-MILLENNIUM simulation (Baugh et al. 2019), which spans a volume of 8003 Mpc3, with a mass resolution of mp = 1.06 × 108 M⊙, and assuming cosmological parameters consistent with the first year results from the Planck satellite (Planck Collaboration XVI 2014), to predict the evolution of galaxy properties from z ∼ 20 to the local Universe (Cantarella et al. in prep.). We extracted all predicted AGN galaxies from the simulated box with a Lbol > 3 × 1046 erg s−1 at z > 6, resulting in a sample of 42 individual sources. For each of these sources, we studied their mean molecular gas content and SFR by considering the evolution of these quantities over a timescale compatible with the Eddington time (∼4 × 107 yr), before the peak of the predicted AGN activity. The synthetic objects show a distribution of tdep that covers a larger range of values, with the bulk of the population of 42 GAEA quasars that have depletion times of a few hundred million years. This can be due to the prescription used in the GAEA SAM that instantaneously triggers the quasar phase when a sufficient amount of gas is concentrated in a certain radius in the host galaxy. This favors the emergence of bright quasars, as well as the case of less massive objects, with respect to our current picture of bright quasars. Indeed, GAEA quasars show MBH that is roughly two orders of magnitude smaller (log(MBH/M⊙)∼7.40 ± 0.23) than that observed in the quasar at redshift z > 6 (log(MBH/M⊙)∼9.33 ± 0.24). This is likely due to the flat prescription adopted for black hole seeding, which does not consider the hypothesis that massive (MBH > 106 M⊙) can form via the direct collapse of giant gas clouds in the early Universe (for an alternative approach, see Trinca et al. 2022). Furthermore, the different distribution between synthetic and observed quasars can be partially due to an observational bias because we only detect the brightest sources at the EoR. Indeed, considering the quasars at Cosmic Noon that span a larger regime Lbol and LFIR (see the upper-right and lower-left panels of Fig. 4), quasars from the SUPER (SINFONI Survey for Unveiling the Physics and Effect of Radiative feedback) and KASHz (KMOS AGN Survey at High redshift) samples from Bertola et al. (2024) show depletion times that are consistent with the bulk of GAEA objects. Eventually, AGN and star-forming galaxies at lower redshifts show a wider range of tdep value, but the bulk of the distributions peak close to tdep ∼ 1 Gyr.

An interpretation of this result requires that the feedback from the SMBH in high-redshift and luminous quasars can remove or heat the cold gas reservoir from their host (Brusa et al. 2018; Fiore et al. 2017; Fluetsch et al. 2019). However, several studies of SMBH-driven winds in quasars at different redshifts have shown that they are not powerful enough to impact the star formation (Bischetti et al. 2019b, Tripodi et al. 2023, and Novak et al. 2019), at least in the cold phase of the gas. Simulations suggest that the quenching effect of the accreting SMBH is not instantaneously effective and, hence, that accretion onto the SMBH and massive star formation coexist (Costa et al. 2020, Costa et al. 2022, and Valentini et al. 2021).

A different interpretation proposes that the relatively low depletion times are likely due to the highly efficient star formation of the quasar host that is favored by the concentration of cold gas in a relatively compact size. To examine this scenario, we then derived the gas SFE, defined as SFEgas = 1/tdep = SFR/MH2, which is represented as a function of LIR in the upper-right panel of Fig. 4. The trend between SFE and LFIR shown in the upper-right panel of Fig. 4 suggests that the galaxies that are brighter at FIR wavelengths (hence have higher SFR), are more efficient at forming stars compared to galaxies at lower LFIR. This is true, irrespective of the presence of an AGN at their center. Nonetheless, quasar hosts exceed up to one order of magnitude the SFE measured in star-forming and luminous galaxies (Tacconi et al. 2018; dash-doted and dotted lines in the upper panel of Fig. 4) at high luminosities (LFIR > 1012 L⊙). In particular, quasars at the EoR show the highest efficiencies among the objects considered in our analysis, except for the hyper-luminous quasars at Cosmic Noon from the WISE/SDSS-selected hyper-luminous quasar (WISSH) sample (Bischetti et al. 2021 and Salvestrini et al. in prep.). Regarding z > 7 massive quasar hosts, we find lower SFE limits of > 10−8 yr−1, slightly higher than the SFE measured in J1007: SFE = 9.7 ± 2.5 × 10−9 yr−1.

As discussed in Sect. 3.2, a highly star-forming host galaxy is a common characteristic among bright (Lbol > 1012.5 L⊙) quasars at high redshifts (e.g., Venemans et al. 2017c, Decarli et al. 2022, Izumi et al. 2021, and Tripodi et al. 2024a). However, this effect could result from a combination of observational bias and the evolution of the MS (Renzini & Peng 2015) at high redshifts. Regarding the first point, we still lack a complete sampling, especially at (sub)millimeter wavelengths, of galaxy populations at intermediate luminosities (LFIR < 1012 L⊙). Concerning the second point, the SFR is expected to increase for MS galaxies with redshift as a result of the higher cosmological accretion rate at early times (0 < z < 6; Tacconi et al. 2020, Walter et al. 2020, and Popesso et al. 2023). This is confirmed by observational studies (e.g., Tacconi et al. 2020 and Sommovigo et al. 2022) that found that the depletion time decreases with redshift, with values of ∼0.01 Gyr at z ∼ 7, consistent with those observed in quasars at the EoR (see the upper-left panel of Fig. 4). These two points should be taken into account when interpreting the lower-right panel of Fig. 4. In each histogram, AGN and quasar sources are divided by redshift, while the luminous and star-forming galaxies from Tacconi et al. (2018) are shown for comparison in the lower panel. Quasars at the EoR represent only the population of the brightest sources at this epoch, and we are likely missing the rich population of intermediate luminosity AGN recently discovered with James Webb Space Telescope (JWST; e.g., Harikane et al. (2023), Maiolino et al. 2024a, and Maiolino et al. 2024b). In this regard, the results from the GAEA simulations predict a broader SFE distribution, with a peak at lower SFEs with respect to observations of quasars at similar redshifts.

Thanks to JWST, it is now possible to detect the stellar light from the quasar host (Ding et al. 2023), and to model the corresponding emission to derive stellar masses (e.g., Ding et al. 2023 and Yue et al. 2024). We took the estimate and upper limit of the stellar mass from Yue et al. (2024) for four quasars at the EoR (J1120 and three other quasars), which are listed in Table C.1. We then calculated the specific SFR (sSFR =SFR/M⋆) and molecular gas fraction (fgas = MH2/M⋆. The lower limit and estimate of the sSFR (> 10−8 yr−1) are way above the expected value for MS galaxies at the EoR. Their molecular gas fraction (fgas ∼ 0.05 − 0.5) is lower than what observed in local quasars and AGN (Shangguan et al. 2020 and Salvestrini et al. 2022), but consistent with that derived from dynamical measurements in luminous quasars at Cosmic Noon (Bischetti et al. 2021). This confirms that these objects are experiencing an intense starburst phase that is coeval with the bright phase of the quasars, justifying the shift between the median value of the SFE distribution between AGN and non-AGN host galaxies in the histograms of Fig. 4.

In the lower-left panel of Fig. 4, we show the SFE as a function of the bolometric luminosity of the accreting SMBH. The solid black line connects the median value of the SFE of quasars and AGN divided into four bins of Lbol. The increasing trend between the SFE and Lbol5 shown in the lower-left panel of Fig. 4 suggests that the processes that convey a significant amount of gas toward the nuclear region (e.g., disk instabilities) are more efficient than the effect of the feedback (e.g., outflows) produced by the accreting SMBH, which is expected to be proportional to the quasar luminosity (Fiore et al. 2017; Fluetsch et al. 2019). We then conclude that among the known population of bright quasars at the EoR there is no clear evidence of an efficient quenching of the host-galaxy star formation due to quasar feedback. On the contrary, at such early epochs, a bright quasar phase is likely coeval with an intense buildup phase in the host galaxy. The coexistence between SMBH growth and intense star formation in the host is visible at all redshifts, when considering luminous quasars, suggesting that quasar activity benefits from the gas concentration in the nuclear region of the galaxy. An alternative scenario can be that the quasar activity triggers the star formation in the host, by compressing the gas in the surrounding medium (Cresci & Maiolino 2018).

5. Conclusions

In this work we present new NOEMA observations targeting CO emission lines for five z > 7 quasars, completing the survey of molecular gas properties for all eight known quasars at this epoch. These observations represent the highest-redshift investigation of cold dust and gas to date in a sample of quasar host galaxies at the EoR. By modeling the FIR emission with a MBB using EOS-DUSTFIT, we derived dust properties and SFR estimates for our sample. Although no statistically significant CO emission lines were detected in the five targets, we derived upper limits on the molecular gas mass. We compared the molecular gas and dust properties of all known z > 7 quasars with a compilation of quasar hosts and star-forming galaxies at different redshifts. Our results are:

  • Among the eight known quasars at z > 7, only one object (J1007; Feruglio et al. 2023) has a statistically significant molecular gas mass estimate. We find no massive gas reservoirs (MH2< a few times 1010 M⊙) at z > 7.

  • Combining the new observations presented in this work with measurements from the literature, we provide a new calibration to derive MH2 from L[CII] for z ≳ 6 quasars.

  • Quasar host galaxies at z > 7 had already accumulated large dust reservoirs (Md ∼ 108 M⊙) in a relatively short time (a few hundred megayears) after the Big Bang. This suggests that the physical processes responsible for dust enrichment are very efficient and may include contributions from evolved stellar populations, SN, and reprocessing within the ISM. The GDR estimate and upper limits for z > 7 quasars presented here align with the mean GDR for luminous quasars and AGN at later epochs, which is consistent with the local value (GDR = 100; Draine & Li 2007).

  • Quasars at the EoR are hosted in galaxies undergoing intense starburst phases, with SFEs up to an order of magnitude higher than those expected for non-AGN host galaxies. This suggests that the emergence of a luminous quasar phase is coeval with the rapid buildup of the host galaxy.

  • Semi-analytic models of quasar host galaxies from GAEA at the EoR also support the idea that the quasar phase is triggered during periods of efficient star formation.

2

Project IDs of the dataset analyzed are listed in the Acknowledgements. J2356 observation with project ID 2023.1.00443.S is still in a proprietary period at the moment of this analysis.

4

EOS-DUSTFIT is a publicly available tool for fitting the cold dust SED of galaxies (https://github.com/roberta96/EOS-Dustfit), which has been used in this work and in Tripodi et al. (2024a). Details about the modeling can be found on the GitHub page and in Tripodi et al. (2024a).

5

The Pearson correlation coefficient for the SFE-Lbol relation is 0.53 and 0.46 considering all quasars at z > 1 or the all AGN and quasars shown in the plot, respectively.

Acknowledgments

FS and CF thank J.M. Winters for his help in data reduction and calibration. FS thanks V. D’Odorico for the useful discussion and suggestions. 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]. This work is based on observations carried out under project number S23CX with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). LZ, EP, AT, and FS acknowledge financial support from the Bando Ricerca Fondamentale INAF 2022 Large Grant “Toward an holistic view of the Titans: multi-band observations of z > 6 QSOs powered by greedy supermassive black holes”. CF, and FS acknowledge financial support from PRIN MUR 2022 2022TKPB2P – BIG-z and Ricerca Fondamentale INAF 2023 Data Analysis grant “ARCHIE ARchive Cosmic HI & ISM Evolution”. M.B. acknowledges support from INAF under project 1.05.12.04.01 – 431 MINI-GRANTS di RSN1 “Mini-feedback” and support from UniTs under project DF-microgrants23 “Hyper-Gal”. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01188.S, #2019.A.00017.S, #2019.1.00074.S, #2019.1.01025.S, #2021.1.00934.S, #2021.1.01262.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. Software: ASTROPY (Astropy Collaboration 2013, 2018, 2022); CASA (CASA Team et al. 2022); EMCEE (Foreman-Mackey et al. 2013, 2019); GILDAS https://www.iram.fr/IRAMFR/GILDAS/; LINMIX (Kelly 2007); MATPLOTLIB (Hunter 2007); NUMPY (Harris et al. 2020); SEABORN (Waskom 2021); SCIPY (Virtanen et al. 2020).

References

  1. Algera, H. S. B., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 6867 [Google Scholar]
  2. Amorín, R., Muñoz-Tuñón, C., Aguerri, J. A. L., & Planesas, P. 2016, A&A, 588, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bañados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11 [Google Scholar]
  7. Baugh, C. M., Gonzalez-Perez, V., Lagos, C. D. P., et al. 2019, MNRAS, 483, 4922 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bertola, E., Circosta, C., Ginolfi, M., et al. 2024, A&A, 691, A178 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bischetti, M., Piconcelli, E., Feruglio, C., et al. 2018, A&A, 617, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bischetti, M., Piconcelli, E., Feruglio, C., et al. 2019a, A&A, 628, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bischetti, M., Maiolino, R., Carniani, S., et al. 2019b, A&A, 630, A59 [NASA ADS] [CrossRef] [EDP Sciences] [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., Choi, H., Fiore, F., et al. 2024, ApJ, 970, 9 [NASA ADS] [CrossRef] [Google Scholar]
  15. Blitz, L., & Rosolowsky, E. 2006, ApJ, 650, 933 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [CrossRef] [Google Scholar]
  17. Brusa, M., Cresci, G., Daddi, E., et al. 2018, A&A, 612, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105 [NASA ADS] [CrossRef] [Google Scholar]
  19. Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939 [NASA ADS] [Google Scholar]
  20. CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501 [NASA ADS] [CrossRef] [Google Scholar]
  21. Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Casavecchia, B., Maio, U., Péroux, C., & Ciardi, B. 2025, A&A, 693, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  24. Costa, T., Pakmor, R., & Springel, V. 2020, MNRAS, 497, 5229 [Google Scholar]
  25. Costa, T., Arrigoni Battaia, F., Farina, E. P., et al. 2022, MNRAS, 517, 1767 [NASA ADS] [CrossRef] [Google Scholar]
  26. Cresci, G., & Maiolino, R. 2018, Nat. Astron., 2, 179 [Google Scholar]
  27. Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Nature, 545, 457 [Google Scholar]
  28. Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97 [Google Scholar]
  29. Decarli, R., Pensabene, A., Venemans, B., et al. 2022, A&A, 662, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Decarli, R., Pensabene, A., Diaz-Santos, T., et al. 2023, A&A, 673, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Decarli, R., Loiacono, F., Farina, E. P., et al. 2024, A&A, 689, A219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. De Lucia, G., Fontanot, F., Xie, L., & Hirschmann, M. 2024, A&A, 687, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Dell’Agli, F., Valiante, R., Kamath, D., Ventura, P., & García-Hernández, D. A. 2019, MNRAS, 486, 4738 [Google Scholar]
  35. Dietrich, M., & Hamann, F. 2004, ApJ, 611, 761 [NASA ADS] [CrossRef] [Google Scholar]
  36. Di Mascia, F., Gallerani, S., Behrens, C., et al. 2021, MNRAS, 503, 2349 [NASA ADS] [CrossRef] [Google Scholar]
  37. Ding, X., Onoue, M., Silverman, J. D., et al. 2023, Nature, 621, 51 [NASA ADS] [CrossRef] [Google Scholar]
  38. Draine, B. T. 2003, ARA&A, 41, 241 [NASA ADS] [CrossRef] [Google Scholar]
  39. Draine, B. T. 2009, in Cosmic Dust – Near and Far, eds. T. Henning, E. Grün, & J. Steinacker, ASP Conf. Ser., 414, 453 [Google Scholar]
  40. Draine, B. T., & Li, A. 2007, ApJ, 657, 810 [CrossRef] [Google Scholar]
  41. Duras, F., Bongiorno, A., Piconcelli, E., et al. 2017, A&A, 604, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Elvis, M., Marengo, M., & Karovska, M. 2002, ApJ, 567, L107 [NASA ADS] [CrossRef] [Google Scholar]
  43. Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373 [NASA ADS] [CrossRef] [Google Scholar]
  44. Feruglio, C., Fiore, F., Carniani, S., et al. 2018, A&A, 619, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Feruglio, C., Maio, U., Tripodi, R., et al. 2023, ApJ, 954, L10 [NASA ADS] [CrossRef] [Google Scholar]
  46. Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586 [NASA ADS] [Google Scholar]
  48. Fontanot, F., De Lucia, G., Hirschmann, M., et al. 2020, MNRAS, 496, 3943 [CrossRef] [Google Scholar]
  49. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  50. Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, J. Open Source Software, 4, 1864 [CrossRef] [Google Scholar]
  51. Gallerani, S., Maiolino, R., Juarez, Y., et al. 2010, A&A, 523, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18 [EDP Sciences] [Google Scholar]
  53. Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633 [Google Scholar]
  54. Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39 [NASA ADS] [CrossRef] [Google Scholar]
  55. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  56. Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2018, ApJ, 861, 95 [Google Scholar]
  57. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hygate, A. P. S., Hodge, J. A., da Cunha, E., et al. 2023, MNRAS, 524, 1775 [NASA ADS] [CrossRef] [Google Scholar]
  59. Izumi, T., Matsuoka, Y., Fujimoto, S., et al. 2021, ApJ, 914, 36 [CrossRef] [Google Scholar]
  60. Izumi, T., Matsuoka, Y., Onoue, M., et al. 2024, ApJ, 972, 116 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kaasinen, M., Venemans, B., Harrington, K. C., et al. 2024, A&A, 684, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Kashikawa, N., Ishizaki, Y., Willott, C. J., et al. 2015, ApJ, 798, 28 [Google Scholar]
  63. Kelly, B. C. 2007, ApJ, 665, 1489 [Google Scholar]
  64. Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189 [Google Scholar]
  65. Kirchschlager, F., Schmidt, F. D., Barlow, M. J., et al. 2019, MNRAS, 489, 4465 [CrossRef] [Google Scholar]
  66. Koss, M. J., Strittmatter, B., Lamperti, I., et al. 2021, ApJS, 252, 29 [CrossRef] [Google Scholar]
  67. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  68. Leśniewska, A., & Michałowski, M. J. 2019, A&A, 624, L13 [Google Scholar]
  69. Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Maio, U., & Viel, M. 2023, A&A, 672, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Maio, U., Péroux, C., & Ciardi, B. 2022, A&A, 657, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Maiolino, R., Cox, P., Caselli, P., et al. 2005, A&A, 440, L51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024a, A&A, 691, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Maiolino, R., Risaliti, G., Signorini, M., et al. 2024b, MNRAS, resubmitted [arXiv:2405.00504] [Google Scholar]
  75. Mancini, M., Schneider, R., Graziani, L., et al. 2015, MNRAS, 451, L70 [CrossRef] [Google Scholar]
  76. Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, ApJ, 849, 91 [Google Scholar]
  77. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127 [Google Scholar]
  78. Micelotta, E. R., Matsuura, M., & Sarangi, A. 2018, Space Sci. Rev., 214, 53 [NASA ADS] [CrossRef] [Google Scholar]
  79. Neeleman, M., Bañados, E., Walter, F., et al. 2019, ApJ, 882, 10 [Google Scholar]
  80. Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141 [NASA ADS] [CrossRef] [Google Scholar]
  81. Novak, M., Bañados, E., Decarli, R., et al. 2019, ApJ, 881, 63 [NASA ADS] [CrossRef] [Google Scholar]
  82. Pensabene, A., van der Werf, P., Decarli, R., et al. 2022, A&A, 667, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103 [CrossRef] [EDP Sciences] [Google Scholar]
  84. Pipino, A., Fan, X. L., Matteucci, F., et al. 2011, A&A, 525, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526 [Google Scholar]
  87. Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152 [NASA ADS] [CrossRef] [Google Scholar]
  88. Pozzi, F., Calura, F., Fudamoto, Y., et al. 2021, A&A, 653, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31 [Google Scholar]
  90. Renzini, A., & Peng, Y.-J. 2015, ApJ, 801, L29 [Google Scholar]
  91. Riechers, D. A., Walter, F., Carilli, C. L., et al. 2006, ApJ, 650, 604 [Google Scholar]
  92. Riechers, D. A., Walter, F., Carilli, C. L., & Bertoldi, F. 2007, ApJ, 671, L13 [NASA ADS] [CrossRef] [Google Scholar]
  93. Riechers, D. A., Walter, F., Bertoldi, F., et al. 2009, ApJ, 703, 1338 [NASA ADS] [CrossRef] [Google Scholar]
  94. Rizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201 [Google Scholar]
  95. Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952 [NASA ADS] [CrossRef] [Google Scholar]
  96. Romano, M., Morselli, L., Cassata, P., et al. 2022, A&A, 660, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Ronconi, T., Lapi, A., Torsello, M., et al. 2024, A&A, 685, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Saintonge, A., Wilson, C. D., Xiao, T., et al. 2018, MNRAS, 481, 3497 [NASA ADS] [CrossRef] [Google Scholar]
  99. Salomé, Q., Krongold, Y., Longinotti, A. L., et al. 2023, MNRAS, 524, 3130 [CrossRef] [Google Scholar]
  100. Salvestrini, F., Gruppioni, C., Pozzi, F., et al. 2020, A&A, 641, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Salvestrini, F., Gruppioni, C., Hatziminaoglou, E., et al. 2022, A&A, 663, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Sarangi, A., Dwek, E., & Kazanas, D. 2019, ApJ, 885, 126 [Google Scholar]
  103. Schneider, R., & Maiolino, R. 2024, A&ARv, 32, 2 [NASA ADS] [CrossRef] [Google Scholar]
  104. Shangguan, J., Ho, L. C., Bauer, F. E., Wang, R., & Treister, E. 2020, ApJ, 899, 112 [NASA ADS] [CrossRef] [Google Scholar]
  105. Shao, Y., Wang, R., Weiss, A., et al. 2022, A&A, 668, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45 [Google Scholar]
  107. Shen, Y., Wu, J., Jiang, L., et al. 2019, ApJ, 873, 35 [Google Scholar]
  108. Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122 [NASA ADS] [CrossRef] [Google Scholar]
  109. Stefan, I. I., Carilli, C. L., Wagg, J., et al. 2015, MNRAS, 451, 1713 [NASA ADS] [CrossRef] [Google Scholar]
  110. Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [Google Scholar]
  111. Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [NASA ADS] [CrossRef] [Google Scholar]
  112. Trinca, A., Schneider, R., Valiante, R., et al. 2022, MNRAS, 511, 616 [NASA ADS] [CrossRef] [Google Scholar]
  113. Tripodi, R., Feruglio, C., Kemper, F., et al. 2023, ApJ, 946, L45 [NASA ADS] [CrossRef] [Google Scholar]
  114. Tripodi, R., Feruglio, C., Fiore, F., et al. 2024a, A&A, 689, A220 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Tripodi, R., Scholtz, J., Maiolino, R., et al. 2024b, A&A, 682, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Tsukui, T., Wisnioski, E., Krumholz, M. R., & Battisti, A. 2023, MNRAS, 523, 4654 [NASA ADS] [CrossRef] [Google Scholar]
  117. Valentini, M., Gallerani, S., & Ferrara, A. 2021, MNRAS, 507, 1 [NASA ADS] [CrossRef] [Google Scholar]
  118. Valiante, R., Schneider, R., Salvadori, S., & Bianchi, S. 2011, MNRAS, 416, 1916 [Google Scholar]
  119. Valiante, R., Schneider, R., Salvadori, S., & Gallerani, S. 2014, MNRAS, 444, 2442 [NASA ADS] [CrossRef] [Google Scholar]
  120. Venemans, B. P., Walter, F., Decarli, R., et al. 2017a, ApJ, 837, 146 [NASA ADS] [CrossRef] [Google Scholar]
  121. Venemans, B. P., Walter, F., Decarli, R., et al. 2017b, ApJ, 845, 154 [Google Scholar]
  122. Venemans, B. P., Walter, F., Decarli, R., et al. 2017c, ApJ, 851, L8 [NASA ADS] [CrossRef] [Google Scholar]
  123. Venemans, B. P., Neeleman, M., Walter, F., et al. 2019, ApJ, 874, L30 [NASA ADS] [CrossRef] [Google Scholar]
  124. Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130 [Google Scholar]
  125. Ventura, P., Karakas, A., Dell’Agli, F., García-Hernández, D. A., & Guzman-Ramirez, L. 2018, MNRAS, 475, 2282 [NASA ADS] [Google Scholar]
  126. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  127. Vizgan, D., Greve, T. R., Olsen, K. P., et al. 2022, ApJ, 929, 92 [NASA ADS] [CrossRef] [Google Scholar]
  128. Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111 [Google Scholar]
  129. Walter, F., Neeleman, M., Decarli, R., et al. 2022, ApJ, 927, 21 [NASA ADS] [CrossRef] [Google Scholar]
  130. Wang, R., Carilli, C. L., Beelen, A., et al. 2007, AJ, 134, 617 [NASA ADS] [CrossRef] [Google Scholar]
  131. Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699 [Google Scholar]
  132. Wang, R., Wagg, J., Carilli, C. L., et al. 2011, AJ, 142, 101 [NASA ADS] [CrossRef] [Google Scholar]
  133. Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44 [Google Scholar]
  134. Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1 [Google Scholar]
  135. Wang, F., Yang, J., Fan, X., et al. 2024, ApJ, 968, 9 [NASA ADS] [CrossRef] [Google Scholar]
  136. Waskom, M. L. 2021, J. Open Source Software, 6, 3021 [NASA ADS] [CrossRef] [Google Scholar]
  137. Witstok, J., Smit, R., Maiolino, R., et al. 2022, MNRAS, 515, 1751 [Google Scholar]
  138. Witstok, J., Jones, G. C., Maiolino, R., Smit, R., & Schneider, R. 2023, MNRAS, 523, 3119 [NASA ADS] [CrossRef] [Google Scholar]
  139. Xie, L., De Lucia, G., Hirschmann, M., Fontanot, F., & Zoldan, A. 2017, MNRAS, 469, 968 [Google Scholar]
  140. Yang, J., Wang, F., Fan, X., et al. 2020, ApJ, 897, L14 [Google Scholar]
  141. Yang, J., Wang, F., Fan, X., et al. 2021, ApJ, 923, 262 [NASA ADS] [CrossRef] [Google Scholar]
  142. Yue, M., Eilers, A.-C., Simcoe, R. A., et al. 2024, ApJ, 966, 176 [NASA ADS] [CrossRef] [Google Scholar]
  143. Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976 [Google Scholar]

Appendix A: Far-infrared SED models

In Fig. A.1 we present the best-fit models of the FIR SED of the five quasars at z > 7 in our sample. The corresponding corner plots of the free parameters obtained with the Bayesian analysis of the MBB with the EOS-DUSTFIT code are shown in Fig. A.2. Details of the analysis are reported in Sect. 3.1.

In Fig. A.3 we show the results of the modeling of the continuum emission when assuming different values for β and Td for each of the five quasars. The baseline model (i.e., the one shown in Fig. A.1) is represented with a solid black line in each panel, while observed fluxes and upper limits are black diamonds and triangles, respectively. The impact of the assumption of a different dust temperature is shown as the best-fit model obtained with a Td = 70 (40) K, which is displayed with a red (orange) curve. In the case of J1243 and J0252, β was left free to vary, while it is fixed to 1.6 for the remaining targets.

For J0313, J0038, and J2356, we also evaluate the impact of the assumption of β = 2 and β = 1.2 with a blue and solid green line, respectively. For J1243 and J0252, the solid blue line represents the case of β = 1.6 fixed. When we varied β, we left Td = 55 K fixed. Relative residuals are shown in the lower box of each panel, color-coded accordingly to the best-fit lines of the upper box.

thumbnail Fig. A.1.

FIR SED. Observed fluxes and upper limits are shown as green diamonds and triangles, respectively. The best-fit model is the black line, and the corresponding parameters (namely β, Td, and Md) are reported in each panel. The gray-shaded region represents the 68% confidence interval of the best-fit parameters. Relative residuals are shown in the lower box of each panel.

thumbnail Fig. A.2.

Corner plots showing the posterior probability distributions of the best fit parameters: Md for each source, and β only for J0252 and J1243. Solid cyan lines indicate the best-fitting value for each parameter, while the dashed lines mark each parameter’s 16th and 84th percentiles.

thumbnail Fig. A.3.

FIR SED obtained with different parameter values. Observed fluxes and upper limits are shown as diamonds and triangles, respectively. The best-fit parameters and corresponding uncertainties are reported in the legend. Best-fit models and relative residuals in the lower box are color-coded accordingly.

Appendix B: Continuum maps

In Fig. B.1 we present the continuum maps produced from the interferometric data collected and analyzed as described in Sect. 2. Flux densities were extracted from the continuum maps as described in Sect. 2.

thumbnail Fig. B.1.

Continuum maps. The source name and the observed frequency are in the upper-left corner. The optical position of the quasar is reported with a black cross. Contour levels are -3, -2, 2, 3, 4, and 5 σ significance; the RMS is listed in Table 1. The clean beam for each observation is shown in the lower-left corner of the diagram, and the corresponding size is reported in Table 1.

Appendix C: Comparison sample of high-z quasars

In Table C.1 we report the information collected from the literature for 21 quasars at redshift z ≳ 6 and a quasar at z ∼ 5. Quasars have been selected to have observations of CO emission lines, and properties derived from the FIR SED fitting. In addition, 19 out of the 21 quasars at z ≳ 6 have L[CII] measurement. We report L ′ C O ( 1 − − 0 ) $ L\prime_{CO(1--0)} $ derived from L ′ C O ( J + 1 − − J ) $ L\prime_{CO(J+1--J)} $ available in the literature (see column Jup), assuming for all objects the same fixed CO SLED that we used in Sect. 3. In particular, we assumed r J + 1 , J = L ′ C O ( J + 1 − − J ) / L ′ C O ( 1 − − 0 ) $ r_{J+1,J}=L\prime_{CO(J+1--J)}/L\prime_{CO(1--0)} $ as follows: r7, 6 = 0.65, r6, 5 = 0.9, r2, 1 = 1. The host-galaxy SFR is then derived from LFIR by adopting a Chabrier IMF. Moreover, we report only the measurement uncertainties for the sources for which these quantities were reported in the literature (see references). During calculations and fitting procedures described in Sect. 3, for each measurement without uncertainty, we assigned an error that is equal to the median of the relative uncertainties available.

Table C.1.

Quasar properties. From left to right, columns include source name, coordinates, redshift, logarithm of bolometric luminosity, logarithm of the BH mass, [CII] luminosity, logarithm of the total infrared luminosity, CO(1-0) luminosity, upper J of the CO transition considered, logarithm of the dust mass. CO(1-0) luminosities are derived by assuming a CO SLED from CO(7-6) emission line, apart from a and b, where we used CO(2-1) and (6-5) transitions. References for optical and FIR measurements are: [1] Shen et al. (2011); [2] Mazzucchelli et al. (2017); [3] Wang et al. (2021); [4] Dietrich & Hamann (2004); [5] Shen et al. (2019); [6] Yang et al. (2021); [7] Kashikawa et al. (2015); [8] Decarli et al. (2018); [9] Venemans et al. (2017b); [10] Tripodi et al. (2024a); [11] Decarli et al. (2022); [12] Wang et al. (2011); [13] Decarli et al. (2023); [14] Pensabene et al. (2022); [15] Venemans et al. (2017a); [16] Kaasinen et al. (2024); [17] Stefan et al. (2015); [18] Feruglio et al. (2018).

All Tables

Table 1.

NOEMA observations.

In the text
Table 2.

Quasar properties.

In the text
Table C.1.

Quasar properties. From left to right, columns include source name, coordinates, redshift, logarithm of bolometric luminosity, logarithm of the BH mass, [CII] luminosity, logarithm of the total infrared luminosity, CO(1-0) luminosity, upper J of the CO transition considered, logarithm of the dust mass. CO(1-0) luminosities are derived by assuming a CO SLED from CO(7-6) emission line, apart from a and b, where we used CO(2-1) and (6-5) transitions. References for optical and FIR measurements are: [1] Shen et al. (2011); [2] Mazzucchelli et al. (2017); [3] Wang et al. (2021); [4] Dietrich & Hamann (2004); [5] Shen et al. (2019); [6] Yang et al. (2021); [7] Kashikawa et al. (2015); [8] Decarli et al. (2018); [9] Venemans et al. (2017b); [10] Tripodi et al. (2024a); [11] Decarli et al. (2022); [12] Wang et al. (2011); [13] Decarli et al. (2023); [14] Pensabene et al. (2022); [15] Venemans et al. (2017a); [16] Kaasinen et al. (2024); [17] Stefan et al. (2015); [18] Feruglio et al. (2018).

In the text

All Figures

thumbnail Fig. 1.

Dust mass and molecular gas mass vs. redshift. Upper panel:Md vs. redshift. Red circles are the z > 7 quasars presented in this work. Green circles are for J1120 from Venemans et al. (2017a), J1007 from Feruglio et al. (2023), and J1342 from Novak et al. (2019). Gray circles are the 5 < z < 7 quasars from Venemans et al. (2017b), Izumi et al. (2021), Decarli et al. (2022), and Tripodi et al. (2024a), while blue dots are for quasars at intermediate redshifts (1 < z < 5) from Bischetti et al. (2021) and Salvestrini et al. (in prep.). Lower panel:MH2 vs. redshift. Data are color-coded as in the upper panel. We also include the quasars at intermediate redshifts (1 < z < 5) from Bertola et al. (2024).
In the text
thumbnail Fig. 2.

SFR and molecular gas mass vs. [CII] luminosity. Left panel:z > 7 quasars with L[CII] measurements: J1234, J0038, and J0252 (red circles) and J1120, J1342, and J1007 (green circle). Gray dots are the z ≳ 6 quasars from the literature (see Table C.1). The best-fit relation and the relative 68% confidence interval are represented as the solid black line and green-shaded region, respectively. The best-fit parameters of the relation are: slope 1 . 15 − 0.27 + 0.28 $ 1.15^{+0.28}_{-0.27} $, normalization − 8 . 16 − 0.06 + 0.06 $ -8.16^{+0.06}_{-0.06} $, and intrinsic dispersion 0 . 07 − 0.02 + 0.03 $ 0.07^{+0.03}_{-0.02} $ dex. The dotted line is the relation from De Looze et al. (2014) for low-metallicity dwarf galaxies in the local Universe. The dashed line is the best-fit relation from Romano et al. (2022) for galaxies from the ALPINE survey (4 < z < 6), including both detections and stacked non-detections. Right panel:z > 7 quasars with L[CII] measurements and MH2 upper limits: J1234, J0038, and J0252 (red circles) and J1342 (green circle); J1007 has both L[CII] and MH2 measurements. MH2 measurements and upper limits for z ≳ 6 quasars from the literature have been homogenized to those presented in this work by assuming a common αCO of 0.8 M⊙ pc−2 (K km s−1)−1 and a CO SLED. The best-fit relation and the relative 68% confidence interval are represented as the solid black line and gray-shaded region, respectively. The best-fit intrinsic dispersion of the relation is δ intr = 0 . 10 − 0.03 + 0.06 $ \delta_{\mathrm{intr}}=0.10^{+0.06}_{-0.03} $ dex. The dashed line is the scaling relation derived by Zanella et al. (2018) for main-sequence (MS) and luminous infrared galaxies, up to z ∼ 6. We also report the scaling relation provided by Madden et al. (2020) for local dwarf galaxies (dotted line).
In the text
thumbnail Fig. 3.

GDR vs. dust mass. The horizontal black line represents GDR = 100. Red circles are the z > 7 quasars presented in this work with MH2 upper limits. Green circles are for J1007 from Feruglio et al. (2023), J1120 from Venemans et al. (2017a), and J1342 from Novak et al. (2019). Gray circles are the 6 ≲ z < 7 quasars from Table C.1, while blue dots are for quasars at lower redshifts from Bischetti et al. (2021), Decarli et al. (2022), and Salvestrini et al. (in prep.). Local Seyfert galaxies are shown as red diamonds.
In the text
thumbnail Fig. 4.

Molecular gas consumption timescales and SFEs. Upper-left panel: Depletion time as a function of (1 + z), shown in logarithmic scale. Sources are color-coded as follows: the z > 7 quasars are red (this work) and green (Venemans et al. 2017a, Novak et al. 2019, and Feruglio et al. 2023). Literature samples include: z < 4 star-forming galaxies (pink-shaded region; Tacconi et al. 2018); local Seyfert galaxies and quasars (red diamonds; Bischetti et al. 2019a, Shangguan et al. 2020, Salvestrini et al. 2020, 2022, and Salomé et al. 2023); quasars at Cosmic Noon (blue dots; Bischetti et al. 2021 and Bertola et al. 2024); and quasars at 6 ≲ z < 7 (gray circles; see Table C.1). The parameter space occupied by the simulated quasars from GAEA is shown as an orange-shaded region, with 10% isodensity lines. The dashed line represents the evolution of the depletion time with redshift (Sommovigo et al. 2022). Upper-right panel: SFE displayed as a function of LIR. The best-fit lines are shown for illustration purposes only and were obtained from a linear fit of the data shown in the legend. Data points are color-coded as in the upper-left panel. Lower-left panel: SFE vs. Lbol. The black diamonds represent the median SFE value in each of the equally spaced Lbol bins. The vertical lines represent the distance between the 84th and 16th percentile of the SFE distribution in each bin. Data points are color-coded as in the upper-left panel. Lower-right panel: The histograms of the SFE distribution of each sample. The dashed line represents each distribution’s median value.
In the text
thumbnail Fig. A.1.

FIR SED. Observed fluxes and upper limits are shown as green diamonds and triangles, respectively. The best-fit model is the black line, and the corresponding parameters (namely β, Td, and Md) are reported in each panel. The gray-shaded region represents the 68% confidence interval of the best-fit parameters. Relative residuals are shown in the lower box of each panel.
In the text
thumbnail Fig. A.2.

Corner plots showing the posterior probability distributions of the best fit parameters: Md for each source, and β only for J0252 and J1243. Solid cyan lines indicate the best-fitting value for each parameter, while the dashed lines mark each parameter’s 16th and 84th percentiles.
In the text
thumbnail Fig. A.3.

FIR SED obtained with different parameter values. Observed fluxes and upper limits are shown as diamonds and triangles, respectively. The best-fit parameters and corresponding uncertainties are reported in the legend. Best-fit models and relative residuals in the lower box are color-coded accordingly.
In the text
thumbnail Fig. B.1.

Continuum maps. The source name and the observed frequency are in the upper-left corner. The optical position of the quasar is reported with a black cross. Contour levels are -3, -2, 2, 3, 4, and 5 σ significance; the RMS is listed in Table 1. The clean beam for each observation is shown in the lower-left corner of the diagram, and the corresponding size is reported in Table 1.
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.