Physical Properties of Hyperluminous, Dust-obscured Quasars at z ∼ 3: Multiwavelength Spectral Energy Distribution Analysis and Cold Gas Content Revealed by ALMA

We present a UV to millimeter spectral energy distribution (SED) analysis of 16 hyperluminous, dust-obscured quasars at z ∼ 3, selected by the Wide-field Infrared Survey Explorer. We aim to investigate the physical properties of these quasars, with a focus on their molecular gas content. We decompose the SEDs into three components: stellar, cold dust, and active galactic nucleus (AGN). By doing so, we are able to derive and analyze the relevant properties of each component. We determine the molecular gas mass from CO line emission based on Atacama Large Millimeter/submillimeter Array (ALMA) observations. By including ALMA observations in the multiwavelength SED analysis, we derive the molecular gas fractions, gas depletion timescales, and star formation efficiencies (SFEs). Their sample median and 16th–84th quartile ranges are fgas∼0.33−0.17+0.33 , t depl ∼ 39−28+85 Myr, and SFE ∼ 297−195+659 K km s−1 pc−2, respectively. Compared to main-sequence galaxies, they have lower molecular gas content and higher SFEs, similar to quasars in the literature. This suggests that the gas in these quasars is rapidly depleted, likely as the result of intense starburst activity and AGN feedback. The observed correlations between these properties and the AGN luminosities further support this scenario. Additionally, we infer the black hole to stellar mass ratio and black hole mass growth rate, which indicate significant central black hole mass assembly over short timescales. Our results are consistent with the scenario that our sample represents a short transition phase toward unobscured quasars.


INTRODUCTION
The coevolution between the central supermassive black hole (SMBH) and host galaxy is now widely ac-Corresponding author: Lulu Fan llfan@ustc.edu.cnknowledged (Kormendy & Ho 2013).This is evidenced by the tight correlation between the mass of central SMBHs and the stellar bulge masses in galaxies (e.g., Magorrian et al. 1998;Ferrarese & Ford 2005).In one of the most popular coevolution scenarios, galaxy gas-rich major galaxy mergers trigger intense starbursts, provide the fuel for central SMBH accretion, and trigger active galactic nucleus (AGN) activity delayed after the triggering of starbursts (Hopkins et al. 2008).This stage of AGN-starburst composite systems often leads to the observation of galaxies as dust-obscured quasars.The galaxies will evolve to unobscured quasars after the accreting SMBH experienced a "feedback"phase, which clears the dust and gas in the galaxy in the form of powerful outflowing winds (see Fabian 2012;Somerville & Davé 2015, for recent reviews).During this evolutionary sequence first presented by Sanders et al. (1988), dust-obscured quasars have been identified as essential links between starbursts and unobscured quasars.They serve a critical role in the rapid assembly of both the SMBH and galaxy mass, as well as in AGN feedback (Hickox & Alexander 2018).Evidence from both observations and theoretical models have suggested that the efficiency of galaxy-scale outflows increases with quasar bolometric luminosity (see, e.g., Heckman & Best 2014;Hopkins et al. 2016;Fiore et al. 2017;García-Burillo et al. 2021).Therefore, luminous, obscured quasars are good candidates for investigating the interplay between host galaxies and their central SMBHs.
The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) has revealed an important population of luminous, dust-obscured galaxies at z ∼ 3 by selecting sources that are strongly detected at 12 and 22 µm, but weakly or not detected at 3.4 and 4.6 µm (Eisenhardt et al. 2012).A variety of follow-up studies utilizing different techniques have been carried out.Multiband spectral energy distribution (SED) analyses have played an essential role in the study of these highredshift objects.Through the construction of median IR SEDs, it has been revealed that these galaxies exhibit a high mid-infrared (MIR) to submillimeter luminosity ratio, elevated dust temperatures, and extraordinary bolometric luminosities over 10 13 L ⊙ (Wu et al. 2012).The number density of these luminous, DOGs is comparable to that of equally luminous type 1 quasars (Assef et al. 2015).These galaxies represent an exceedingly rare class of DOGs, now commonly referred to as hot, dust-obscured Galaxies (Hot DOGs).Further investigations have shown that the IR SEDs of the most luminous Hot DOGs are dominated by hot dust at a temperature exceeding 450 K (Tsai et al. 2015), and their IR SEDs can be decomposed through two component, AGN-starburst SED fitting (Fan et al. 2016b(Fan et al. , 2017(Fan et al. , 2018a)).UV-optical spectral analyses of Hot DOGs found black hole masses ∼ 10 9 M ⊙ which are accreting near or above the Eddington limit and also host powerful ionized outflows (Wu et al. 2012;Tsai et al. 2018;Wu et al. 2018;Finnerty et al. 2020;Jun et al. 2020).Millimeter interferometric observations like Atacama Large Millimeter/submillimeter Array (ALMA) of several Hot DOGs revealed a highly turbulent ISM and also provide evidence of possible molecular outflows (Wu et al. 2014;Díaz-Santos et al. 2016, 2018, 2021;Fan et al. 2018bFan et al. , 2019;;Penney et al. 2020).The X-ray observations of Hot DOGs consistently find high column densities close to Compton thick (Stern et al. 2014;Piconcelli et al. 2015;Assef et al. 2016Assef et al. , 2020;;Ricci et al. 2017;Vito et al. 2018;Zappacosta et al. 2018).The environments where Hot DOGs reside were found to be significantly overdense (Jones et al. 2015(Jones et al. , 2017;;Silva et al. 2015;Penney et al. 2019;Ginolfi et al. 2022;Luo et al. 2022;Zewdie et al. 2023).All results are generally consistent with the merger-driven coevolution scenario.
SED fitting is an effective method for decomposing emissions between star formation and AGN (Sokol et al. 2023).Some widely used SED fitting codes, including SED3FIT (Berta et al. 2013), CIGALE (Boquien et al. 2019;Yang et al. 2022), and BayeSED (Han & Han 2012, 2014, 2019), have been proven to be efficient tools for exploring AGN-starburst systems, and AGN models are constantly evolving and becoming more and more accurate and realistic (Fritz et al. 2006;Nenkova et al. 2008a,b;Hönig & Kishimoto 2010, 2017;Siebenmorgen et al. 2015;Stalevski et al. 2016).At high redshift, it is hard to resolve the central AGN emission from the host galaxy.SED observations, modeling, and fitting are indispensable to investigate the physical properties of these high-redshift AGNs and their host galaxies (e.g., Merloni et al. 2010;Bongiorno et al. 2014;Suh et al. 2019;López et al. 2023).
Molecular gas, predominantly traced by carbon monoxide (CO) emission lines, acts as the fuel for both star formation and black hole accretion.Additionally, it plays an important role in energy feedback from AGN (e.g., Feruglio et al. 2017;Bischetti et al. 2019;Fluetsch et al. 2019).Investigations of the molecular gas content (M H2 ) in quasars, coupled with other observables such as SED-derived dust masses (M dust ), stellar masses (M ⋆ ), and star formation rates (SFRs), provide valuable insights into the physical processes driving the coevolution of galaxies and SMBHs (e.g., Brusa et al. 2015Brusa et al. , 2018;;Banerji et al. 2017;Kakkad et al. 2017;Perna et al. 2018;Bischetti et al. 2021).However, for Cosmic noon (z ∼ 2-3) when both star formation and black hole accretion activity in the Universe peaked (Shapley 2011), there is a lack of systematic investigation into the molecular gas properties of luminous quasars, primarily limited to analyses of individual sources or relatively small sample sizes.
To investigate systematically the physical properties of Hot DOGs at z ∼ 3, particularly their cold gas con-  (2012) and Tsai et al. (2015).( 5) and ( 6): The IR luminosities of AGN torus and cold dust emission derived by IR decomposition as reported in Fan et al. (2016b).
tent, we conduct a comprehensive UV to millimeter SED analysis of a subsample of 16 Hot DOGs selected from Fan et al. (2016b).Table 1 lists the properties of these Hot DOGs, as reported in Fan et al. (2016b).This subsample has ALMA observations of CO emission lines, which were used to constrain their molecular gas content.This study represents the largest sample of high-redshift, luminous obscured quasars so far in this kind of study.We derived the properties of stellar, cold dust, AGN, and gas components and calculated their relative values (e.g., molecular gas fraction f gas = M H2 /(M H2 + M ⋆ ) and SFE), further testing the role of Hot DOGs in massive galaxy and SMBH coevolution.In Section 2, we present details of ALMA observations and the subsequent data analysis.Section 3 covers the construction of our multiwavelength SED and the SED modeling approach we used.Results and discussions are described in Section 4 and Section 5, respectively.Finally, in Section 6, we summarize our findings and draw a conclusion.Throughout this work, we assume a Lambda cold dark matter (ΛCDM) cosmology (see Komatsu et al. 2011) with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

ALMA OBSERVATIONS AND DATA ANALYSIS
Observations were carried out with ALMA using the Band-3 and Band-4 receivers during Cycle 5. Most of our sample sources were observed in our project 2017.1.00441.S (PI: L. Fan).A few sources were allocated to a different project and available in the ALMA archive (2017.1.00358.S), and we also include these observations. 1Table 2 summarizes the details of the observations, including a list of the calibrators.We note that while the observations in 2017.1.00441.S used a spectral setup for the spectral window (spw) of the sideband that was predicted to include the emission line (and continuum mode for the spectral windows of the other sideband), project 2017.1.00358.S used a continuum setup for all spectral windows.
Reduction, calibration, and imaging were done using casa (Common Astronomy Software Application 2 ; Mc-Mullin et al. 2007).The pipeline-reduced data delivered from the observatory was of sufficient quality such that no additional flagging and further calibration were necessary.The pipeline includes the steps required for a standard reduction and calibration, such as flagging, bandpass calibration, as well as flux and gain calibration.A conservative estimate of the uncertainty of the absolute flux calibration is 10%.
The data was imaged both as continuum and spectral cube using natural weighting.The casa task uvcontsub was used to subtract the continuum from uv data for sources for which the continuum was detected.A continuum image was produced combining all spectral windows, while a spectral cube was constructed for the spectral windows tuned to the redshifted CO line.The rms sensitivity and the synthesized beam size achieved by imaging with a natural weighting scheme are given in Table 2.

MULTIWAVELENGTH DATA AND SPECTRAL
ENERGY DISTRIBUTION FITTING

UV to Millimeter Spectral Energy Distribution Data
To decompose the host galaxy emission from the central AGN and estimate their physical properties, such as stellar mass (M ⋆ ) and SFR, we constructed UV to millimeter SEDs for all objects in our sample.Various  2014).The IR broadband photometry at wavelengths ranging from 12 µm to the millimeter band were directly collected from their parent samples as reported in Fan et al. (2016b).The ALMA continuum observations at rest-frame 3 mm have been included to constrain the cold dust component more accurately.

SED Analysis
The UV to millimeter SED fitting in our study was conducted with the latest version of BayeSED (Han & Han 2012, 2014, 2019) 10 , namely BayeSED V3.0.This new version has improved the accuracy and speed of the stellar population synthesis algorithm and has been tested with mock galaxies to show good quality and speed for parameter estimation of galaxies (Han et al. 2023).For SED models given as a SED library, principal component analysis is employed to reduce memory usage.Then, interpolation between the SED models is conducted with artificial neural networks or K-nearest neighbors to allow a fast and continuous sampling of high-dimensional parameter space.Finally, the Multi-Nest algorithm is employed to sample the parameter space and calculate the posterior probability distribution of the parameters.
The SED of each galaxy in our sample was decomposed into three components: stellar, cold dust, and central AGN.The stellar emission was modeled by using the Bruzual & Charlot (2003) simple stellar population library with a Chabrier ( 2003) initial mass function (IMF), and an exponential declining star formation history (SFH).The Calzetti et al. (2000) dust attenuation law was used, and the energy of stellar emission absorbed by cold dust was assumed to be totally reemitted in the IR band.This assumption can break the degeneracy in the UV and optical bands by the complement of IR photometry (da Cunha et al. 2008;Buat et al. 2014), and stellar population properties can be more robustly constrained.The cold dust emission was modeled conventionally as a graybody, which was defined as S λ ∝ (1 − e −( λ 0 λ ) β )B λ (T dust ), where λ 0 =125µm, and the dust temperature T dust and the emissivity index β are two free parameters in the SED fitting.The CLUMPY torus model (Nenkova et al. 2008a,b) 11 has significantly advanced the modeling of IR emission in various AGN samples (Ramos Almeida & Ricci 2017), and has been utilized to model the UV-IR emission of the central AGN of our samples.Six parameters are employed in the CLUMPY model to describe the geometry and physical properties of the CLUMPY torus.The SED analysis method employed in this study is consistent with that used in a previous work by Fan et al. (2019), who presented a comprehensive analysis of W0533-3401.However, we have restricted the T dust parameter to a range of 30-50K to address AGN contamination in the FIR, as discussed in Section 5.4.The 12 free parameters used in the fitting process are listed in Table 3 of Fan et al. (2019).For the three Hot DOGs without optical-NIR photometry, we excluded the stellar population component to prevent overfitting.Thus, the stellar population parameters, such as M ⋆ , have not been estimated for these three galaxies.The upper flux density limit was taken into account by setting the corresponding flux density error to a negative value according to the convention in BayeSED.

CO Emission Lines and Continuum
We detect significant CO(3-2) or CO(4-3) emission lines for all but four galaxies in our sample.Among the four galaxies, W0248+2705 and W1814+3412 have a tentative detection, while no emission lines were found for the other two sources.The emission lines were fitted with a single Gaussian line profile to estimate the redshift, peak flux density, line width, and integrated intensity.The results are given in Table 4.For the four nondetections, 3σ upper limits of the peak flux density and integrated intensity (assuming a line width of 300 km s −1 ) are reported in Table 4. Galaxy centers were estimated from the moment-0 maps.Figure 1 shows the CO line spectra and moment 0 maps of our galaxies.
The continuum was detected in all but four galaxies.The continuum flux density and size were measured by fitting 2D Gaussian profile to the maps.We note that W0126-0529 has no line detection but has a clear detection of its continuum.The relevant measurements of continuum are reported in Table 4, where nondetections are listed as 5σ upper limits.
The line luminosities L ′ CO in K km s −1 pc 2 were calculated by using an equation from Solomon & Vanden Bout (2005): where S CO ∆v is the CO line flux in Jy km s −1 , D L is the luminosity distance in Mpc, and ν obs is the observed frequency in GHz.The line emission data for W0149+2305, W0220+0137, and W0410-0913 were collected from Fan et al. (2018b), and the line emission data for W2246-0526 were obtained from Díaz-Santos et al. ( 2018).W0126-0529 has been re-observed recently and a redshift of z = 0.832 (Vito et al. 2018) has been reported, suggesting that it may be classified as a low-redshift Hot DOG, similar to W1904+4853 in Li et al. (2023).An uncertain redshift might be responsible for its nondetection in our observations.Thus, we chose to discard this source from our analysis.Utilizing line ratios of L 87, as recommended for QSOs (Carilli & Walter 2013), we derived the CO(1-0) line luminosities.We note that a intermediate value of L ′ CO(3−2) / L ′ CO(1−0) = 0.8 between SMGs and QSOs was adopted in Banerji et al. (2017) and Fan et al. (2019).However, our results in Section 4.2 indicate that the bolometric luminosities of our sample are dominated by central AGN emission, favoring the line ratios typical of QSOs.The CO-to-H 2 conversion factor, α CO , relates the CO(1-0) luminosity to total molecular gas mass.For the Milk Way, we have α CO ∼ 4.6 M ⊙ (K km s −1 pc 2 ) −1 (Bolatto et al. 2013).However, for starburst galaxies, α CO is significantly lower, with a value of 0.8 M ⊙ (K km s −1 pc 2 ) −1 (Carilli & Walter 2013).We adopted α CO = 0.8 M ⊙ (K km s −1 pc 2 ) −1 for our galaxies, given the starburst nature (Fan et al. 2019) and the high merger fraction (Fan et al. 2016a) of Hot DOGs.The logarithmic molecular gas mass log M H2 [M ⊙ ] of our sample ranges from 9.55 to 11.59, with a median value of 10.56.Our estimations of molecular gas mass based on CO line observations is 0.56 dex higher than the values predicted by dust mass assuming a Milky Way dust-to-gas ratio of 0.01 in Fan et al. (2016b).The line luminosities and molecular gas mass of our sample are listed in Table 5

Results of the SED Fitting and Dust Properties
The best-fit SEDs are shown in Figure 2. Thanks to high AGN obscuration, the host galaxies of our Hot DOGs sample are easily observable, and the stellar emission can be separated out so that we can estimate the physical properties of the host galaxies.We adopted the median of the posterior probability distribution of each parameter as the fiducial value, and the uncertainties are reported as the 68% confidence intervals around the fiducial values.The derived properties are listed in Table 6.Based on our UV-millimeter three-component SED modeling described in Section 3, the stellar population parameters, including stellar masses and SFRs, as well as the cold dust properties from the cold dust component and AGN luminosities have been obtained.
The cold dust infrared luminosity L IR was calculated by integrating the cold dust component from 8 to 1000 µm.Our Hot DOGs exhibit L IR ∼ 10 13 L ⊙ .The estimated dust temperatures range from 42K to 48K, with a median of 45 K, consistent with T dust -L IR relation of SMGs (e.g., Magnelli et al. 2012;Roseboom et al. 2012).The parameter β denotes the power-law index of the optical depth, with τ λ = ( λ0 λ ) β .It ranges from 2.0 to 2.8,  with a median of 2.5.The measured T dust and β of Hot DOGs is similar to the most luminous quasar at z = 6.327 (Tripodi et al. 2023).A relatively high β value indicates optically thick dust in IR bands, which has also been reported in other compact starburst galaxies (e.g., Scoville et al. 2017).With the cold dust temperature T dust and the emissivity index β, we derived the dust mass with the formula: where D L is the luminosity distance, S ν obs is the flux density at observed frequency ν obs , κ νrest = κ 0 (ν/ν 0 ) β is the absorption coefficient at the corresponding restframe frequency, and B(ν rest , T dust ) is the Planck function per unit frequency at temperature T dust .We adopted κ 850µm = 3.8 cm 2 g −1 by following Wu et al. (2014) and Fan et al. (2016bFan et al. ( , 2019)).The estimated dust masses of our sample are consistent with previous IR SED decomposition results in Fan et al. (2016b), with a sample median of 8.0 × 10 7 M ⊙ .We note that β was fixed to 1.6 in Fan et al. (2016b) to avoid degeneracy.The uncertainties in the dust masses shown in Table 6 will be larger if we consider the adopted κ 850µm value which can vary from ∼0.4 to ∼11 cm 2 g −1 (e.g., James et al. 2002;Draine 2003;Dunne et al. 2003;Siebenmor-gen et al. 2014).All these dust properties are listed in Table 6.
Together with the CO measurements, these parameters allow us to derive the molecular gas fractions and star formation efficiencies (Table 7).

Gas-to-dust ratio
The estimated molecular gas mass M H2 is plotted as a function of M dust in Figure 3.The two dashed lines represent gas-to-dust ratios δ GDR = 50 and 150, which cover the typical values derived for the Milky Way (Jenkins 2004), local star-forming galaxies (Draine et al. 2007;Rémy-Ruyer et al. 2014), and high-redshift SMGs (Magnelli et al. 2012;Miettinen et al. 2017).Most of our Hot DOGs exhibit a high δ GDR value ∼ 474 +711 −324 , with a median uncertainty of 43% (uncertainties propagated from M H2 and M dust ).Bischetti et al. (2021) found a median δ GDR value of 180 for their nine hyperluminous, Type I quasars at z ∼ 2 − 4, slightly higher than the typical values, and they attributed this to an increasing δ GDR with redshift (e.g., Miettinen et al. 2017).We note that our sample is selected to have Herschel PACS and SPIRE observations, and have either SPIRE 500 µm or SCUBA-2 850 µm detection, and therefore may be biased toward the most intense starbursting systems, where the supernova-shock-driven dust heating and destruction may be more significant (Jones 2004).The dust mass decreases with increasing dust temperature (Fan et al. 2016b, and references there in).The estimated T dust ∼ 45 K may trace a warmer dust component associated with photodissociation regions by starbursts, instead of the diffuse ISM with temperature below 30 K which represents the bulk of the dust mass (Draine & Li 2007;Liang et al. 2019;Sommovigo et al. 2020;Pozzi et al. 2021).The dust mass estimated from a simplistic graybody model may be underestimated by up to a factor of 3 compared to that derived by the dust models (for example, Draine & Li 2007) which include more parameters and adequately describe the multicomponent dust properties (Conroy 2013, and references therein).When we take a dust temperature of 30 K we found the dust mass increases by ∼ 0.5 dex and the corresponding δ GDR approximately decreases by a factor of 3.

Stellar mass and molecular gas fraction
The logarithmic stellar masses range from 10.3 to 12.1, with a median of 10.8, indicating the majority of our Hot DOGs are massive galaxies.However, Díaz-Santos et al. ( 2021) discovered higher SED-based stellar masses than dynamical masses based on ALMA [C II] observations.They attributed this overestimation to the lower angular resolution of optical/NIR data than that of interferometric [C II] data.We defer a comparative analysis of the stellar and dynamical masses of our sample to a future work.With M ⋆ derived from SED fitting and molecular gas mass M H2 inferred from CO line observations, we calculated the molecular gas fraction, which is defined as f gas = M H2 /(M H2 + M ⋆ ), and listed the results in Table 7. Due to the lack of optical-NIR photometry data for W0248+2705, W0615-5716 and W1248-2154, we cannot obtain SED-based stellar masses M ⋆ for these three Hot DOGs.To put an upper limit on the molecular gas fraction for these three Hot DOGs, we adopted M ⋆ = 10 10.3 M ⊙ which is the minimum M ⋆ estimated for the other galaxies in our sample.
Our Hot DOGs have comparable redshift and stellar mass ranges to the SMGs studied in Miettinen et al. (2017) (z ∼ 2.3 and log M ⋆ [M ⊙ ] = 11.09+0.41  −0.53 , respectively).However, the molecular gas fraction of our Hot DOGs (0.33 +0.33  −0.17 ) is much lower than the SMGs in Miettinen et al. (2017) (0.62 +0.27  −0.23 ).In Figure 4. we show the molecular gas fraction f gas as a function of redshift for our Hot DOGs, as well as literature samples of SMGs, obscured quasars at z > 1 from Perna et al. (2018) and Palomar-Green quasars in the local Universe from Shangguan et al. (2020).We note that the median uncertainty of f gas of our Hot DOG sample is approximately 60%, which results from the uncertainties in M H2 and M ⋆ and is shown as gray in the bottom-right corner in Figure 4.The sample of high-redshift SMGs from Miettinen et al. (2017) is also represented in this figure.These literature samples all have M ⋆ ∼ 10 11 M ⊙ , which is comparable to our Hot DOGs.To demonstrate how the molecular gas fractions of these galaxies compare with MS galaxies, we present the gas fraction evolutionary trend of MS galaxies of M ⋆ = 10 11 M ⊙ , M ⋆ = 2 × 10 10 M ⊙ , and M ⋆ = 5 × 10 11 M ⊙ predicted from the 2-SFM model (Sargent et al. 2014)  forming galaxies.The sample of SMGs lacking AGN exhibits a higher molecular gas fraction compared to MS galaxies.The discrepancy in f gas between active and normal galaxies is likely due to the depletion of cold gas by the AGN feedback.We color coded our Hot DOGs in the figure according to their AGN bolometric luminosity(L AGN )derived by integrating the AGN component of the best-fit UV to millimeter SED.Most Hot DOGs with an f gas lower than MS galaxies exhibit a relatively higher L AGN , which is consistent with recent findings of a positive correlation between the efficiency of AGN feedback traced by the mass outflow rate and L AGN (see, e.g.Hopkins et al. 2016;Fiore et al. 2017;García-Burillo et al. 2021).AGN-driven outflows deposit energy and momentum into the surrounding gas and affect the evolution of the host galaxy by heating and ejecting the ISM (e.g., Weymann et al. 1991;Faucher-Giguère & Quataert 2012;Marasco et al. 2020).Our Hot DOGs stand out from the other samples by their extremely high L AGN , which may explain their relatively large deviation from the MS relation.

Star formation rate and star formation efficiency
The SFRs of Hot DOGs were derived by averaging the modeled exponential SFH in the last 100 Myr.The cold dust component is modeled by adding a grey body component, whose energy budget is identical to the attenuated luminosity in the UV-optical band, i.e., the so called energy balance assumption.Figure 5 6.By invoking the MS   7) demonstrate that most of our Hot DOGs are extreme starburst systems, which could be triggered by gas-rich galaxy mergers (e.g., Noguchi & Ishibashi 1986;Mihos & Hernquist 1996).With the molecular gas masses and SFRs, we derived the gas depletion timescale t depl = M H2 / SFR.Similar to obscured quasars (e.g., Aravena et al. 2008;Brusa et al. 2018), our Hot DOGs exhibit a short gas depletion timescale.In contrast, typical starburst galaxies have a t depl value of several 100 hundred megayears (e.g., Genzel et al. 2010;Bothwell et al. 2013;Miettinen et al. 2017).Hot DOGs have been discovered to exhibit outflow mass loss rates of several thousand solar masses per year (Finnerty et al. 2020).Considering the powerful AGN outflows, the gas depletion timescale could be shorter.The molecular gas in Hot DOGs may be depleted within several tens of megayears, resulting in the lower molecular gas fraction in Hot DOGs.The few Hot DOGs with a high gas fraction may represent a relatively earlier stage after the AGN activity has been triggered.
In Figure 6, we show the SFEs (= SFR / M H2 ) as a function of IR luminosity L IR for our Hot DOGs, as well as for the compiled samples of SMGs and unobscured and obscured quasars from Perna et al. (2018).We use the corresponding observables, L IR obtained from the cold dust component of our SED decomposition and L ′ CO(1−0) acquired through ALMA observations, to calculate the SFE as L IR / L ′ CO(1−0) .The SFEs exhibit a tight correlation with L IR for both MS galaxies (Sargent et al. 2014) and SMGs.Our Hot DOGs, along with obscured and unobscured quasars at z > 1, display high SFEs that are well above the relation for MS galaxies.This is likely due to the rapid depletion of cold gas by AGN feedback and star formation.We plot the relation between SFE and L AGN for Hot DOGs and obscured quasars in Figure 7.We also included WISE-SDSS selected hyper-luminous (WISSH) quasars from Bischetti et al. (2021) and PG QSOs from Shangguan et al. (2020).The AGN bolometric lumi-  nosities of the obscured quasars from Perna et al. ( 2018) is calculated from their X-ray luminosities, assuming a luminosity-dependent bolometric correction from Duras et al. (2020).Across all four quasar samples, there is a positive correlation between SFE and L AGN .The Spearman's rank correlation coefficient ρ = 0.532 (p = 3.55e-8), where p is the probability of the null hypothesis that a correlation does not exist.When excluding the PG QSOs given that they are in local Universe, the ρ deceases to 0.404 with p = 2.21e-3.This correlation can be explained by the enhanced outflow rates in galaxies with high AGN luminosity, and is consistent with their low gas fractions (Figure 4).
With optical-NIR detections Without optical-NIR detections  Cyan square symbols represent a sample of hyperluminous, type I quasars at z ∼ 1-4 in Bischetti et al. (2021).The downward arrows mark upper limits.

Influence of AGN contamination to FIR
Initially, we fitted our sample with BayeSED by setting the T dust parameter range to default values (10-100K).However, we obtained dust temperatures ∼ 60 − 90K, which are higher than DOGs and submillimeter galaxies (SMGs), but are consistent with previous studies employing a similar fitting methodology (Wu et al. 2012;Fan et al. 2016bFan et al. , 2019)).We refer to this SED fitting as T dust -unconstrained fitting and denote the resulting dust temperature and cold dust luminosity as T dust,0 and L IR,0 , listed in Table 6.The L IR,0 is consistent with those given in Table 1.For the T dust -unconstrained fitting, the FIR emission is entirely attributed to cold dust heated by massive stars.
However, many studies suggested that AGNs could heat dust up to kiloparsec scales, significantly contributing to their FIR emission (see section 7 in Netzer 2015, for a recent review).For example, Schneider et al. (2015) conducted radiative transfer modeling and reported that AGN-heated dust contributes from 30% up to 70% of the FIR luminosity for a high-redshift quasar host galaxy.Duras et al. (2017) employed the same technique and found 40% to 60% quasar contribution to the FIR emission for their sample.Tsukui et al. (2023) determined an AGN contribution of approximately 53% based on image decomposition of spatially resolved ALMA continuum observations.The high cold dust temperature estimated for Hot DOGs may result from contribution by AGN-heated warmer dust, similar to the findings of Tsukui et al. (2023).The hot dust MIR emission from the central torus may be reprocessed to the FIR by optical thick dust in nuclear region (e.g., Scoville et al. 2017;Sokol et al. 2023).Additionally, the FIR contribution of the central AGN may also originate from emission by narrow-line region (NLR) polar dust, which has been proven to be heterogeneous in nearby AGNs (Netzer 2015, and references therein).Given the high AGN luminosities of our Hot DOGs, the NLR dust is expected to extend to kiloparsec scales.González-Martín et al. (2019) improved SED fitting for higher-luminosity AGN by employing two-component models proposed by Hönig & Kishimoto (2017), which incorporate a clumpy disk and a polar clumpy wind, in contrast to the CLUMPY model.There is evidence suggesting that the extended polar dust emission is likely associated with AGN-driven dusty outflows (e.g., Alonso-Herrero et al. 2021).Lyu & Rieke (2018) utilized their AGN SED libraries, which include an extended polar dust component, to model the SEDs of Hot DOGs.Their results attribute a significant portion of the FIR emission to the polar dust component.
It is crucial to consider accurately the AGN contribution to the FIR emission and to estimate the cold dust luminosity appropriately.Based on the flux density ratios of 350 µm and 1.1 mm continuum of several Hot DOGs, Wu et al. (2012Wu et al. ( , 2014) ) suggested that the cold dust component heated by star formation is not very different from those in starburst galaxies, characterized by temperatures ranging between 30 and 50 K (e.g., Magnelli et al. 2012).We constrained the T dust parameter of the graybody component to be within the range 30-50 K and refitted the SEDs in our sample.This fitting approach is denoted as T dust -constrained fitting.We show an example in Figure 8 which compare the best-fit SEDs obtained via T dust -constrained fitting (solid line) and T dust -unconstrained fitting (dash line).For the T dustconstrained fitting, the AGN emission dominates up to rest-frame 50 µm, while for the T dust -unconstrained fitting, the AGN emission only dominates up to rest-frame 25 µm.The reduced χ 2 only improve by a value of 0.10 after we apply our prior knowledge of the T dust for W0533-3401.We calculate the Bayes factor, defined as the difference in the Bayesian evidence between the two fitting methods, to be in the range of 0.1-3.6,suggesting no strong evidence in favor of the T dust -unconstrained fitting (Han & Han 2014).In the earliest torus model, Pier & Krolik (1993) enlarged the torus to account for the FIR emission by a ∼100 pc scale torus.Similarly, for the CLUMPY model of Nenkova et al. (2008a,b), the AGN FIR emission can come from an extended torus with a large torus outer radius R out = Y R in (Drouart et al. 2014), where Y is the radial extent, one of the six free parameters used to define the CLUMPY model.R in is the inner radius set by the location of the dust at the sublimation temperature, and is computed using the AGN bolometric luminosity L AGN by the equation: R in = 0.4 ( L AGN 10 45 erg s −1 ) 0.5 pc. (4) Based on the T dust -constrained fitting, we obtain L AGN ∼ 10 47.6 and Y ∼ 65 for our sample, which give R out ∼ 0.5 kpc.The L IR is smaller by 0.16-0.47dex compared to the L IR,0 shown in Hot DOGs are usually under heavy dust obscuration, and optically seen as type 2 AGN (Wu et al. 2012).However, Assef et al. (2016) discovered a subpopulation of eight Hot DOGs which exhibit a blue UV-optical SED similar to blue bump by the AGN accretion disk.They named them BHDs.Their rest-frame UV spectra are of typical type 1 quasars, showing broad emission lines (Assef et al. 2020).They confirmed that the UV-optical SEDs of this subpopulation Hot DOGs are due to 1% scattered light from the highly obscured, hyperluminous AGN into our line of sight, using X-ray and imaging polarization observations (Assef et al. 2016(Assef et al. , 2020(Assef et al. , 2022)).One Hot DOG in our sample, namely W0220+0137, has been identified as a BHD in Assef et al. (2016Assef et al. ( , 2020)).
The CLUMPY model of Nenkova et al. (2008a,b) is a geometrical torus model that assumes a certain geometry and dust composition and conducting radiative transfer modeling.The absorption and scattering coefficients given in Ossenkopf et al. (1992) are used, where the UV-optical emission is dominated by the AGNscattered radiation (Nenkova et al. 2008a).Ichikawa et al. (2015) employed the CLUMPY model to fit the SEDs of type 2 Seyferts with scattered light or a hidden broad-line region (HBLR), which is similar to BHDs, and they identified differences in the modeled torus geometry of HBLRs compared to type 2 Seyferts lacking HBLR signatures.As shown in Figure 2, our SED modeling of W0220+0137 is consistent with the results of Assef et al. (2016Assef et al. ( , 2020)), where the UV-optical SED is dominated by an AGN component.For almost all other Hot DOGs with an optical detection, we found that the UV light is dominated by AGN component.There are three Hot DOGs, namely W0134-2922, W1603+2745 and W2054+0207, whose optical-NIR band also exhibits more AGN emission than stellar emission.These objects are also likely to be compatible with BHDs, albeit to a relatively less extent.For these four BHDs, their stellar mass estimations are more uncertain.For example, Merloni et al. (2010) decomposed the UV to MIR SEDs of 89 type 1 AGNs with two-component SED fitting, and they assigned an upper limit to the stellar mass for galaxies whose contribution of stellar light in the K band is less than 5%.None of our Hot DOGs exhibit such low stellar contribution, and the cold dust FIR emission from the host galaxy serves as an additional constraint to the stellar component.We have highlighted these four galaxies in Table 6 and Figure 4. We also note that W2246-0526 and W2305-0039 have nearly equal contributions from AGN and stellar components in the optical-NIR band.

Estimation of the Central Supermassive Black Hole Mass and Growth Rate
Recent studies of Hot DOGs have consistently revealed that their Eddington ratios are near or above the Eddington limit (Wu et al. 2018;Finnerty et al. 2020;Jun et al. 2020).Tsai et al. (2018) reported the measurement of M BH of W2246-0526.By using the L AGN of W2246-0526 from our SED fitting, we estimated a super-Eddington ratio λ Edd = L AGN /L Edd = 1.7,where L Edd /L ⊙ = 3.28 × 10 4 (M BH /M ⊙ ).Given these studies, we assume an Eddington ratio λ Edd = 1.0 for our Hot DOG sample and infer the mass of the SMBH M BH and the black hole to stellar mass ratio M BH /M ⋆ .The results are listed in Table 7.
In Figure 9, we plot M BH /M ⋆ as a function of redshift for Hot DOGs and other samples of AGN for different redshift and AGN bolometric luminosity ranges.Similar to other luminous quasars at redshift 2-4 (e.g., Targett et al. 2012;Trakhtenbrot et al. 2015;Matsuoka et al. 2018), the inferred M BH /M ⋆ of the Hot DOG sample is about 2 times higher than the evolutionary trend of M BH /M ⋆ seen by McLure et al. (2006).Bischetti et al. (2021) also revealed an extremely high ratio between SMBH mass and dynamical mass M BH /M dyn in WISSH quasars.These WISSH quasars also exhibit a close or super unit Eddington ratio.In contrast, obscured and unobscured AGN with moderate AGN bolometric luminosities at redshift 1-2 (e.g., Merloni et al. 2010;Bongiorno et al. 2014) exhibit a relatively low M BH /M ⋆ value.Using the equation L AGN = (η ṀBH c 2 )/ (1 − η) and adopting η 1−η = 0.1, we derived high black hole growth rates of ṀBH ∼ 65 +39 −29 M ⊙ yr −1 for our sample.These results suggest rapid black hole growth in Hot DOGs, consistent with Wu et al. (2018).During this high-accretion phase, the majority of the black hole mass could be assembled within the Salpeter timescale, which is consistent with the gas depletion timescale and the high luminosity state timescale measured in Tsai et al. (2015).

SUMMARY AND CONCLUSION
We present a UV to millimeter SED analysis and molecular gas content measurements in a sample of 16 WISE-selected, hyperluminous dust-obscured quasars at z ∼ 3. We inferred the physical properties of this sample, such as gas-to-dust ratio, molecular gas fraction, and SFE.This study represents the largest sample to date in which a systematic investigation of the cold gas content in hyperluminous quasars at Cosmic Noon has been conducted.The main results can be summarized as follows: 1. Based on ALMA observations of the  or CO(4-3) lines, we have calculated the molecular gas mass M H2 = 3.69 +9.92 −2.61 × 10 10 M ⊙ of our sample by adopting L  2. We modeled the observed UV to millimeter SEDs of our sample using an updated version of BayeSED.The median values and 16th-84th quartiles of all sources in our sample are given as follows.For the cold dust emission represented by the graybody function, we estimated a typical T dust ∼ 45 +2 −1 K, β ∼ 2.5 +0.2 −0.2 , L IR ∼ 1.1 +1.3 −0.3 × 10 13 L ⊙ , M dust = 8.0 +4.0 −4.8 × 10 7 M ⊙ .For the stellar component, we infered typical M ⋆ ∼ 5.8 +12.3  −2.3 × 10 10 M ⊙ and SFR ∼ 1008 +597 −409 M ⊙ yr −1 .For the AGN component, we obtained typical L AGN ∼ 9.7 +5.8 −4.3 × 10 13 L ⊙ .
3. We estimated the gas-to-dust ratio, finding δ GDR ∼ 474 +711 −324 .Most Hot DOGs in our sample exhibited a higher δ GDR than the typical values for the Milky Way, local star-forming galaxies, and high-redshift SMGs.This discrepancy can potentially be attributed to the intense radiation field and high dust temperatures caused by starbursts and potentially the central AGNs, which may result in an underestimation of M dust .4. We inferred the molecular gas fraction, finding f gas ∼ 0.33 +0.33 −0.17 .By comparing our findings with SMGs, obscured quasars from the literature, and the f gas -z relation for MS galaxies, we found that Hot DOGs exhibit a relatively low gas content.This lower gas content is likely attributed to the depletion of gas caused by AGN-driven outflows.
5. The remarkable offset of our Hot DOGs from the MS, ∆ MS ∼ 6.12 +5.1 −2.9 suggests that the majority of our Hot DOGs are extreme starburst systems.The gas depletion timescales, 39 +85 −28 Myr are remarkably short.When comparing the average SFE of ∼ 297 +659 −195 K km s −1 pc −2 with those of SMGs, obscured and unobscured quasars, as well as MS galaxies, we found that Hot DOGs exhibit higher SFEs similar to optically luminous quasars and obscured quasars, rather than SMGs and MS galaxies.Moreover, we discovered a positive correlation between SFE and AGN bolometric luminosity.
6. Based on AGN bolometric luminosity, we inferred a typical black hole growth rate ∼ 65 +39 −29 M ⊙ yr −1 and a typical black hole mass ∼ 3.0 +1.8 −1.3 × 10 9 M ⊙ by adopting η 1−η = 0.1 and λ Edd = 1.0, respectively.These results suggest that the majority of the black hole mass can be assembled within a Salpeter timescale, which is consistent with the gas depletion timescale and the high luminosity state timescale of Hot DOGs suggested by Tsai et al. (2015).The observed black hole to stellar mass ratio ∼ 0.042 +0.029 −0.026 is similar to other highredshift luminous quasars.
We conclude that our results are consistent with the scenario that our sample represents a phase when both star formation and AGN activity are at their peak, leading to rapid depletion of gas and dust, ultimately transiting the galaxies toward unobscured quasars.

Figure 1 .
Figure 1.The CO emission line moment 0 maps and continuum-subtracted spectra.

Figure 2 .
Figure 2. Best-fit model SEDs by T dust -constrained fitting for 16 Hot DOGs in our sample.The source ID, redshift, χ 2 /N d (N d is the photometric data points), and the systematical errors (see Han & Han 2014) are shown in each panel.The red points with error bars represent the observed photometric data, and those with downward arrows mark flux density upper limits.The solid blue, orange, and green lines represent, respectively, the stellar, cold dust, and AGN components.The solid black line represents the total SED.

Figure 3 .
Figure3.Molecular gas mass MH 2 vs. cold dust mass M dust .The red points represent our sample values, while the downward arrows mark the upper limit molecular gas mass for W2210-3507.The typical uncertainties are shown as a gray point with an error bar.The two dashed lines represent a gas-to-dust ratio δGDR value of 50-150, which is the typical value derived for the Milky Way(Jenkins 2004), local star forming galaxies(Draine et al. 2007;Rémy-Ruyer et al. 2014) and high-redshift SMGs(Magnelli et al. 2012;Miettinen et al. 2017).The solid line represents a gas-todust ratio δGDR value of 521, which is typical for our sample based on SED analysis and ALMA observations.

Figure 4 .
Figure 4. Molecular gas fraction fgas as a function of redshift.The 16 Hot DOGs represented by circle symbols in our sample are color coded according to their AGN bolometric luminosity, and their typical uncertainties are plotted in the bottom-right corner as a gray point with error bars.The three sources that did not have an SED-based stellar mass and were estimated to have a stellar mass of 10 10.3 M⊙ are additionally highlighted with red circles, and four sources whose optical emission are dominated by an AGN component (Section 5.4) are additionally highlighted with orange circles.The purple square and olive diamond symbols represent the samples of SMGs and obscured quasars, respectively, compiled in Perna et al. (2018).The black squares represent the sample of SMGs in Miettinen et al. (2017).Note that Miettinen et al. (2017) have excluded from their sample the SMGs with evidence of AGN, whereas a large fraction of SMGs studied by Perna et al. (2018) have AGN activities.The blue pentagon represents the sample of local PG quasars in Shangguan et al. (2020).The median and 16th-84th ranges of our Hot DOGs are shown as a larger filled red circle with an error bar, and the sample median and 16th-84th quantile ranges of the literature samples are shown correspondingly with larger symbols with error bars.The redshift evolutionary trend predicted from the 2-SFM model (Sargent et al. 2014) of MS galaxies of M⋆ = 10 11 M⊙ (typical for our sample) is shown as a solid curve, while MS galaxies of M⋆ = 2 × 10 10 M⊙ and M⋆ = 5 × 10 11 M⊙ are labeled with dashed lines.

Figure 5 .Figure 6 .
Figure 5. SFR traced by the IR luminosity.The sources with optical-NIR detections are shown as blue circles, and they can be represented by the relation as labeled with a solid curve in the figure, which is about 0.10 dex lower than the Kennicutt (1998) relation calibrated by Chabrier IMF.The SFRs of three Hot DOGs without optical-NIR photometry, namely W0248+2705, W0615-5716, and W1248-2154, are calculated based on this relation.Their typical uncertainties are shown as a gray point with an error bar.

Figure 7 .
Figure 7. Positive correlation between SFE and AGN bolometric luminosity.Red circles with error bars represent our Hot DOGs.Green diamond symbols represent the sample of obscured quasars in Perna et al. (2018).The blue pentagons represent the sample of local PG quasars in Shangguan et al.Cyan square symbols represent a sample of hyperluminous, type I quasars at z ∼ 1-4 inBischetti et al. (2021).The downward arrows mark upper limits.

Figure 8 .
Figure 8. Best-fit model SED comparison between T dustconstrained fitting (solid line) and T dust -unconstrained fitting (dashed line) for W0533-3401.The source ID, redshift, χ 2 /N d offset, Bayes factor offset, and the systematical error offset (see Han & Han 2014) are shown on the panel.The color of the legend is identical to that of Figure 2.

Figure 9 .
Figure 9. Redshift evolution of black hole to stellar mass ratio (MBH/M⋆).The MBH/M⋆ ratios of Hot DOGs are shown as red circles, and the sample median and 16th-84th quantile ranges are shown as a larger symbol with error bars.Blue triangles represent a sample of 89 moderately luminous broad-line AGN in the redshift range 1 < z < 2.2 (Merloni et al. 2010).A sample of 21 X-ray obscured, red AGNs with moderate luminosity in Bongiorno et al. (2014) is labeled with a cyan diamond symbol.For luminous quasars with AGN bolometric luminosity > 10 46 erg s −1 , two luminous SDSS quasars at z ∼ 4(Targett et al. 2012) and an extremely red dust-obscured quasar(Toba, & Nagao 2016;Matsuoka et al. 2018) are denoted with purple filled triangles and gold filled square respectively.The green-filled diamond denotes an X-ray selected luminous unobscured quasar, CID-947 at z ∼ 3.3, which has an extremely high black hole to stellar mass ratio MBH/M⋆ = 1/8(Trakhtenbrot et al. 2015).The solid line and two dashed lines show the best-fit evolutionary trend of MBH/M⋆ and 1σ errors at z < 2(McLure et al. 2006).The gray area shows the typical range of MBH/M⋆ ∼ 0.002-0.005 in the local Universe(Kormendy & Ho 2013).

Table 1 .
Sample properties a W0126−0529 has been excluded from our sample for its ambiguous redshift identification.(2)and(3): The WISE coordinates from the AllWISE database.(4):The spectroscopic redshift fromWu et al.

Table 2 .
Summary of the ALMA Observations

Table 4 .
CO and Millimeter Continuum Measurements Based on ALMA Observations.

Table 5 .
CO Line Luminosities and Molecular Gas Masses.
Díaz-Santos et al. (2018), andPenney et al. (2020) of the SMGs compiled inPerna et al. (2018)exhibit AGN activity, while the SMG sample inMiettinen et al. (2017)has excluded SMGs that demonstrate evidence of hosting an AGN.Similar to the PG QSOs fromShangguan et al. (2020)and the SMGs and obscured QSOs fromPerna et al. (2018), the f gas of most of our Hot DOGs is below the relation for MS galaxies with comparable M ⋆ .This is consistent withDíaz-Santos et al. (2018), andPenney et al. (2020), which concluded that Hot DOGs may have lower cold molecular gas content than ordinary star- sys=0.33

Table 6 .
Summary of the results of SED fitting.

Table 7 .
Summary of the physical properties.