Accurate Dust Temperature and Star Formation Rate in the Most Luminous z > 6 Quasar in the Hyperluminous Quasars at the Epoch of Reionization (HYPERION) Sample

In the past decade, the Atacama Large Millimeter/submillimeter Array (ALMA), along with the Northern Extended Millimeter Array, the Very Large Array, and Herschel, have probed the cold gas and dust of quasi-stellar object (QSO) host galaxies. The dust continuum was detected in many z ∼ 6 QSOs, with far-infrared (FIR) luminosities of L_FIR = 10^11–13 L_⊙ and dust masses of about M_dust = 10^7–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 relying only on single-frequency continuum detections. However, if multifrequency ALMA observations are available in the rest-frame FIR, probing 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), resulting high accuracy in the determination of the star formation rate (SFR).

In this Letter, we present ALMA Band 9 observations of the QSO SDSS J010013.02+280225.8 (hereafter J0100+28) at z_[CII] = 6.327 (Wang et al. 2019). Wu et al. (2015) estimated a bolometric luminosity of L_bol = 4.29 × 10¹⁴ L_⊙ and a black hole (BH) mass of M_BH = 1.24 × 10¹⁰ M_⊙ for J0100+28, 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. 2022). Wang et al. (2019) performed a multifrequency analysis of the dust SED, but they could not obtain a precise determination of the dust properties, concluding that J0100+28 has either high dust emissivity (β ≳ 2) or high dust temperature (T_dust ≳ 60 K) or a combination of thereof.

J0100+28 belongs to the HYPerluminous quasars at the Epoch of ReionizatION (HYPERION) sample of z > 6 QSOs, selected to include the QSOs with the most massive SMBH at their epoch, possibly resulting from an exceptionally fast mass growth during their accretion history. The sample consists of 17 QSOs with L_bol = 10^47.3 erg s⁻¹, SMBH masses in the range M_BH = 10⁹–10¹⁰ M_⊙, and Eddington ratios >0.3–0.4. The HYPERION sample and the details of its selection will be presented in L. Zappacosta et al. 2023 (in preparation).

The goal of this work is to derive the most accurate estimate of the dust mass, dust temperature, and therefore of the SFR for the host galaxy of J0100+28. We adopt a ΛCDM cosmology from Planck Collaboration et al. (2020): H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_m = 0.315, and Ω_Λ = 0.685. Thus, the angular scale is 5.66 kpc arcsec⁻¹ at z = 6.3.

We analyze the data set 2021.2.00151.S from the ALMA 7 m array, designed to detect the continuum emission at a frequency of 670.91 GHz in Band 9, and with a total integration time of 2.2 hr. The visibility calibration and imaging are performed through the Common Astronomy Software Applications (CASA; McMullin et al. 2007), version 5.1.1–5. We apply tclean using natural weighting and a 3σcleaning threshold. We image the continuum by collapsing all channels, selected by inspecting the visibilities in all spectral windows. We obtain a clean beam of (1.97'' × 1.17''), corresponding to a spatial resolution of ∼11 kpc, and an rms noise of 0.8 mJy beam⁻¹ in the continuum.

3.1. QSO Continuum Emission and Dust Properties

Figure 1 presents the 670.9 GHz continuum emission map of J0100+28. The peak flux density is 6.99 ± 0.71 mJy beam⁻¹. The source is not spatially resolved. We check that the other FIR and radio flux measurements were extracted from a region similar to our resolution. ¹⁶

Zoom In Zoom Out Reset image size

We perform an SED fitting using the continuum emission measured at ∼671 GHz, together with the emissions presented in Wang et al. (2019) from 32 to 353 GHz, and in Liu et al. (2022) at 1.5, 6, and 10 GHz. Although we are interested in the cold dust properties of the QSO host galaxy, such as T_dust, M_dust and β, we consider also the contribution of the lower frequency emission to the dust SED since in general it may not be negligible. Liu et al. (2022) noticed a time variability among their and previous measurements of the radio continuum in the range [6–10] GHz. For the sake of simplicity, we consider the most recent measurements of the radio continuum emission (i.e., those from Liu et al. 2022). However, we verified that our results do not change by considering all the measurements available in the radio band.

We model the dust continuum with a modified blackbody (MBB) function and the low-frequency radio emission using a power law with an exponential cutoff (PLCO). Details about the fitting functions and procedure can be found in the Appendix. The model has six fitting parameters: dust temperature (T_dust), dust mass (M_dust), dust emissivity index (β) entering in the MBB function, the normalization (n), the radio power-law spectral index (α), and the cutoff frequency (ν_cutoff) for the PLCO. It is worth to stress that the galaxy is likely characterized by a distribution of dust temperatures, which may reach hundreds of kelvin close to the AGN (see, e.g., Walter et al. 2022) and may decrease toward the outskirts; hence, the temperature derived from the SED fitting should be interpreted as an "effective" dust temperature. ¹⁷ The best-fit model has T_dust = 48.4 ± 2.3 K, M_dust = (2.29 ± 0.83) × 10⁷ M_⊙, β = 2.63 ± 0.23, n = 0.08 ± 0.01 mJy, α = 0.48 ± 0.09, and ν_cutoff = 235 ± 100 GHz. Figure 2 shows the observed SED fitted by our best-fit model (left panel) and the posterior distributions for the dust parameters (right panel). Posterior distributions of all parameters are reported in the Appendix. We find a gas-to-dust ratio (GDR) = 236 ± 155 based on our M_dust estimate and the molecular gas mass obtained by Wang et al. (2019; see Table 1). This is in agreement with GDRs found in ultraluminous QSOs at z ∼ 2–4 (Bischetti et al. 2021) and in local galaxies at solar metallicities (De Vis et al. 2019). This latter comparison would imply that J0100+28's host galaxy has already been highly enriched with metals.

Zoom In Zoom Out Reset image size

Table 1. Properties of SDSS J0100+2802

Tdust	[K]	48.4 ± 2.3
Mdust	[107 M⊙]	2.29 ± 0.83
β		2.63 ± 0.23
SFR a	[M⊙ yr−1]	265 ± 32
GDR		236 ± 155
Mgas b	[1010 M⊙]	0.54 ± 0.16
Mdyn(CV) c	[1010 M⊙]	∼890
Mdyn(VT) d	[1010 M⊙]	3.25 ± 0.46
MBH e	[1010 M⊙]	1.05

Notes. Quantities above the line are from this work, while the others are taken from Wang et al. (2019).

^a The SFR is corrected by a factor of 50%, accounting for the contribution of the QSO. ^b The gas mass is estimated within a diameter of ∼1.4 kpc (Wang et al. 2019). ^c Dynamical mass estimated using the circular velocity, assuming a disk inclination of 5° within 1.8 kpc radius (Wang et al. 2019). ^d Dynamical mass estimated using the virial theorem within 1.8 kpc radius (Wang et al. 2019). ^e BH mass from Mg ii emission line (C. Mazzucchelli et al. 2023, in preparation).

Download table as:  ASCIITypeset image

We estimate the total infrared (TIR) luminosity for the best-fit model by integrating from 8 to 1000 μm rest frame, obtaining L_TIR = 5.30 ± 0.64 × 10¹² L_⊙. This would imply an SFR of 530 ± 64 M_⊙ yr⁻¹, adopting a Chabrier initial mass function (Chabrier 2003). However, several observations and radiative transfer simulations suggested that the radiative output of luminous QSOs substantially contributes to dust heating on a kiloparsec scale (Schneider et al. 2015; Di Mascia et al. 2021; Walter et al. 2022). Duras et al. (2017) showed that on average, ∼50% of the total IR luminosity in QSOs with L_bol > 10⁴⁷ erg s⁻¹ is due to dust heated by the AGN radiation. Recently, Di Mascia et al. (2023) found a correction factor of one-thirtieth for the SFR of the brightest object in their sample. However, the L_bol of J0100 is above the range explored in their simulations, and moreover, J0100 has very peculiar properties in terms of UV magnitude and dust temperatures with respect to the simulated QSOs in Di Mascia et al. (2023). Therefore, we consider the average correction proposed by Duras et al. (2017) as the most appropriate choice, and we obtain SFR = 265 ± 32 M_⊙yr⁻¹.

By fitting the dust SED of the QSO J0100+28, we found that our observation in Band 9 favors a dust temperature lower than the effective dust temperatures found in simulations in bright (L_bol > 10¹³ L_⊙) quasar hosts (∼90 K, e.g., Di Mascia et al. 2023). This discrepancy can be due to different dust spatial distributions between the simulated objects and J0100+28 or to limits in the dust modeling and radiative transfer postprocessing (e.g., the absence of a dusty torus in Di Mascia et al. 2021). With the current unresolved observation we are able neither to constrain different temperature components nor to determine the temperature distribution across the galaxy. The value found for the emissivity index is higher than β = 1.6 found for the average SED of high-z QSOs (Beelen et al. 2006). However, as noted by Beelen et al. (2006), large variations of β are found when considering individual QSO SED. Only excluding the Band 9 data, we would obtain a good fit with β < 2, having then a factor of 2 higher temperature (see also Wang et al. 2019). This shows that the results can be misguided by relying only on the lowest frequency data points, while Band 9 is essential to reliably estimate dust parameters up to the highest redshifts. Indeed, Novak et al. (2019), using observations up to Band 8 at ∼404 GHz, could not constrain T_dust in the z = 7.54 quasar ULAS J1342+0928. This work provides the first measurement of the Band 9 continuum of a luminous QSO host galaxy. Similarly, Bakx et al. (2021) used Band 9 observations to constrain T_dust for a galaxy at z = 7.13.

It is not straightforward to assess the physical reasons of the 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 larger than two (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 or the different properties of the grain populations or both. The strong anticorrelation between T_dust and β can arise from the Markov Chain Monte Carlo (MCMC) fitting procedure itself, mainly if there is a massive amount of cold dust (i.e., T ≲ 10 K; Galliano et al. 2018). However, in our case, the accurate sampling of the SED from low to high frequency significantly relaxed the strength of T_dust − β anticorrelation (Figure 2). Therefore, we are confident that the estimates of β and T_dust are not strongly biased due to the effect of this anticorrelation. Finding a physical explanation for the high value of β would require studies of the properties of the dust grains and/or a detailed analysis of the temperature mixing and is beyond the scope of this Letter.

The left panel of Figure 3 shows the ratio of integrated [C ii] to TIR luminosity for J0100+28 and a compilation of high-redshift QSOs and galaxies (see caption). For J0100+28 we find L_[CII]/L_TIR ∼ 7 × 10⁻⁴, given L_[CII] = 3.7 × 10⁹ L_⊙ (Wang et al. 2019). This [C ii] deficit is also predicted for high-z galaxies by semi-analytical models of galaxy evolution (e.g., Lagache et al. 2018), where the [C ii] deficit arises from the high intensity of the interstellar radiation field. Our estimate of L_[CII]/L_TIR agrees well with their results at z ∼ 6 when we extrapolate their predictions at higher L_TIR. Carniani et al. (2018) found that the local L_[CII]−SFR relation for star-forming galaxies (see De Looze et al. 2014) is still valid at high-z but with a dispersion twice higher than observed locally. Our results agree within <1σ with the correlation of Carniani et al. (2018) and with the results for high-z galaxies of Lagache et al.(2018).

Zoom In Zoom Out Reset image size

We evaluate the evolutionary state of the SMBH—host galaxy system. Wang et al. (2019) provide two estimates for M_dyn, one based on the assumption that the [C ii] line arises from a rotating disk with very low inclination and the other using the virial theorem with the hypothesis that the gas is supported by random motion (see Table 1). Although both estimates carry large uncertainties, the latter is in agreement with the dynamical masses commonly measured for other QSOs at high-z (see right panel of Figure 3), and we adopt it as a lower limit. The molecular gas mass is M_gas = 5.4 × 10⁹ M_⊙, and this implies a molecular gas fraction μ = M_gas/M_* = M_gas/(M_dyn − M_BH − M_gas) = 0.4, which is remarkably lower than the typical gas fractions at z ∼ 3–4 (Tacconi et al. 2018). We define the growth efficiency of the galaxy as SFR/M_g+s, where M_g+s = M_gas + M_*, and the SFR is corrected for the QSO contribution. We use M_g+s, instead of M_dyn (as in Tripodi et al. 2022), since the BH mass is not negligible (M_BH ∼ 40%M_dyn; Wang et al. 2019) and can bias our results. ¹⁸ We obtain M_g+s = 2.01 × 10¹⁰ M_⊙; therefore, SFR/M_g+s = 1.3 × 10⁻⁸ yr⁻¹. On the other hand, we derive a BH growth efficiency, ¹⁹ $(1-\epsilon ){\dot{M}}_{\mathrm{BH}}/{M}_{\mathrm{BH}}\,=2.5\times {10}^{-8}\,{\mathrm{yr}}^{-1}$, where we use the BH mass derived from Mg ii (Table 1) and assume a radiative efficiency = 0.1 (e.g., Marconi et al. 2004). The right panel of Figure 3 shows M_BH versus M_dyn for J0100+28, J2310+18 (Tripodi et al. 2022) and a compilation of QSOs at different redshifts, comparing them with the local M_BH − M_dyn relation. The majority of QSOs, including J0100+28, lie above the local relation in the BH dominance regime (Volonteri 2012). However, this tension can be partially softened if accounting for the uncertainties on the dynamical mass estimates. These mainly depend on the determination of the disk inclination and can be significantly high for some QSOs (Valiante et al. 2014; Pensabene et al. 2020). This is also the case of J0100+28, whose M_dyn can be, in principle, as high as 10¹² M_⊙ (see Table 1). For J0100+28, we found $(1-\epsilon ){\dot{M}}_{\mathrm{BH}}/{M}_{\mathrm{BH}}\gt \mathrm{SFR}/{M}_{{\rm{g}}+{\rm{s}}}$, suggesting that the BH is still dominating the process of a BH galaxy growing in this QSO at z = 6.327. This result is still valid if considering lower SFR (adopting the correction by Di Mascia et al. 2023) and/or higher M_dyn (i.e., higher M_g+s at fixed BH mass) since the galaxy growth factor would be even smaller. Our results do not consider the gas inflow. However, we expect this term to contribute on average to both SFR and M_g+s, leaving their ratio mostly unaffected. On the other hand, in QSO J2310+18 at z ∼ 6, AGN feedback might be slowing down the accretion onto the SMBH, while the host galaxy grows fast (Bischetti et al. 2022; Tripodi et al. 2022). The different evolutionary states of J0100+28 and J2310+18, separated by only ∼60 Myr (i.e., Δz ∼ 0.3), arise mainly from the difference in their BH masses (i.e., M_BH,J0100+28 ∼ 2 × M_BH,J2310+18) and in the SFRs (i.e., SFR_J0100+28 ∼ 0.2 × SFR_J2310+18). In principle, SFR/M_* is a better probe of the galaxy growth; however, this is not available for most high-redshift QSOs. For J0100+28, using SFR/M_* would not affect our results since M_g+s ≈ M_*. For J2310, M_* ≪ M_g+s, therefore resulting in an even flatter slope. The BH-dominated growth of J0100+28 matches the qualitative expectations for the evolutionary state of the HYPERION QSOs that were selected to be the most luminous QSOs with the most massive SMBH at their epochs.

We compare our results with "zoom-in" simulations of QSOs using the moving-mesh code AREPO, following BH growth and feedback via energy-driven outflows (Costa et al. 2014, 2015). We found that simulations reproduce BHs only up to masses ∼10⁹ M_⊙ that have host galaxies with dynamical masses ∼10¹¹ M_⊙. These are considerably more massive than the J0100+28 host galaxy. The growth of the system is characterized by intermittent phases, where the BH and galaxy-dominated growth phases change on short timescales. Hence, the diagnostic power of the relation $(1-\epsilon ){\dot{M}}_{\mathrm{BH}}/{M}_{\mathrm{BH}}$ - SFR/M_g+s needs to be validated in larger samples that we plan to investigate in a forthcoming work.

In summary, this work allowed us to measure the SFR with the smallest statistical error, reaching an accuracy of ∼15% for QSO J0100+28. Such unprecedented constraints on the host-galaxy SFR (and T_dust) highlight the critical role of ALMA Band 9 to obtain a robust overview of the buildup of SMBHs and their massive host galaxies at the Epoch of Reionization. The systematic uncertainty on the SFR is still high due to the uncertainty in the estimate of the QSO contribution. Recently, Tsukui et al. (2023) estimated this contribution for a z = 4.4 QSO using high-resolution ALMA observations up to band 9. Given our current unresolved observation, we are not able to use the same approach. However, we plan to determine this correction factor in forthcoming works using both high-resolution observations and radiative transfer models to reproduce the observed SED.

We thank the anonymous referee for the useful suggestions that helped us improving this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2021.2.00151.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. R.T. acknowledges financial support from the University of Trieste. R.T., C.F., F.F., and M.B. 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. R.T. and M.B. acknowledge support from INAF grant 1.05.12.04.01 "Mini-feedback". F.K. acknowledges the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, financed by MCIN/AEI/10.13039/501100011033. S.C. is supported by the European Union (ERC, WINGS,101040227). R.M. acknowledges ERC Advanced Grant 695671 QUENCH, and support from the UK Science and Technology Facilities Council (STFC). R.M. also acknowledges funding from a research professorship from the Royal Society. L.Z. acknowledges 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" and from INAF "Progetti di Ricerca di Rilevante Interesse Nazionale" (PRIN), Bando 2019 "Piercing through the clouds: a multiwavelength study of obscured accretion in nearby supermassive black holes". F.C. is supported by funding of the NASA Physics of the Cosmos Program.

Facility: ALMA. -

Software: astropy, CASA (v5.1.1-5, McMullin et al. 2007).

We model the dust continuum with an MBB function given by

$$\begin{eqnarray}\begin{array}{rcl}{S}_{{\nu }_{\mathrm{obs}}}^{\mathrm{obs}} & = & {S}_{\nu /(1+z)}^{\mathrm{obs}}=\displaystyle \frac{{\rm{\Omega }}}{{\left(1+z\right)}^{3}}[{B}_{\nu }({T}_{\mathrm{dust}}(z))\\ & & -\,{B}_{\nu }({T}_{\mathrm{CMB}}(z))](1-{e}^{-{\tau }_{\nu }}),\end{array}\end{eqnarray} \tag{ A1 }$$

where ${\rm{\Omega }}={\left(1+z\right)}^{4}{A}_{\mathrm{gal}}{D}_{{\rm{L}}}^{-2}$ is the solid angle with A_gal and D_L is the surface area and luminosity distance of the galaxy, respectively (Carniani et al. 2019). The dust optical depth is

$$\begin{eqnarray}&&{\tau }_{\nu }=\displaystyle \frac{{M}_{\mathrm{dust}}}{{A}_{\mathrm{galaxy}}}{k}_{0}{\left(\displaystyle \frac{\nu }{250\ \mathrm{GHz}}\right)}^{\beta },\end{eqnarray} \tag{ A2 }$$

with β the emissivity index and k₀ = 0.45 cm² g⁻¹ the mass absorption coefficient (Beelen et al. 2006). The solid angle is estimated using the continuum mean size from resolved observations (Wang et al. 2019). The effect of the cosmic microwave background (CMB) on the dust temperature is given by

$$\begin{eqnarray}&&{T}_{\mathrm{dust}}(z)={\left({\left({T}_{\mathrm{dust}}\right)}^{4+\beta }+{T}_{0}^{4+\beta }[{\left(1+z\right)}^{4+\beta }-1]\right)}^{\tfrac{1}{4+\beta }},\end{eqnarray} \tag{ A3 }$$

with T₀ = 2.73 K. We also considered the contribution of the CMB emission given by B_ν (T_CMB(z) = T₀(1 + z)) (da Cunha et al. 2013).

We model the low-frequency radio emission using a PLCO that is

$$\begin{eqnarray}&&{F}_{{\nu }_{\mathrm{rest}}}=n\times {\left({\nu }_{\mathrm{rest}}/{\nu }_{0}\right)}^{-\alpha }\times \exp (-{\nu }_{\mathrm{rest}}/{\nu }_{\mathrm{cutoff}}),\end{eqnarray} \tag{ A4 }$$

where n is the normalization; α is the radio power-law spectral index; and ν₀ and ν_cutoff are the reference frequency and the cutoff frequency, respectively. We set ν₀ = 59 GHz to ease the fitting procedure and in order to minimize the covariance between n and α. This choice does not affect our results.

The total model has six fitting parameters: dust temperature (T_dust), dust mass (M_dust), β entering in the MBB function, and n, α, and ν_cutoff for the PLCO. We explore the six-dimensional parameter space using an MCMC algorithm implemented in the EMCEE package (Foreman-Mackey et al. 2013). We assume uniform priors for the fitting parameters: 10 K < T_dust < 300 K, 10⁵ M_⊙ < M_dust < 10⁹ M_⊙, 1.0 < β < 3.0, 0.001 mJy < n < 0.5 mJy, 0.01 < α < 1.0, and 75.0 GHz < ν_cutoff < 500 GHz. The best-fit model has T_dust = 48.4 ± 2.3 K, M_dust = (2.29 ± 0.83) × 10⁷ M_⊙, β = 2.63 ± 0.23, n = 0.08 ± 0.01 mJy, α = 0.48 ± 0.09, and ν_cutoff = 235 ± 100 GHz obtained from an MCMC with 60 chains, 6000 trials, and a burn-in phase of ∼100. Figure 4 shows the posterior distributions of the six fitting parameters T_dust, M_dust, β, n, α, and ν_cutoff.

Zoom In Zoom Out Reset image size