ALMA Reveals a Gas-rich, Maximum Starburst in the Hyperluminous, Dust-obscured Quasar W0533–3401 at z ∼ 2.9

We construct the multiwavelength SED of W0533−3401 by compiling the optical to millimeter broadband photometry from various catalogs available in the literature (see Table 2). Optical/near-infrared photometry in five broad bands, grizY, are retrieved from the first public data release of the Dark Energy Survey (DES DR1; Abbott et al. 2018).¹⁰ Near-infrared J band photometry has been obtained by SOAR/OSIRIS (Assef et al. 2015). The WISE W3 and W4 photometry of W0533−3401 are from the ALLWISE Data Release (Cutri 2013). W3 and W4 flux densities and uncertainties have been converted from catalog Vega magnitude by using zero-points of 29.04 and 8.284 Jy, respectively (Wright et al. 2010). While for WISE W1 and W2 photometry, we do the aperture photometry based on the unblurred coadded WISE images¹¹ (unWISE, Lang 2014; Meisner et al. 2017). The photometry errors have been estimated based on the inverse variance images. We also collect the FIR photometry of W0533−3401 obtained with Herschel (Pilbratt et al. 2010) PACS (Poglitsch et al. 2010) at 70 and 160 μm and SPIRE (Griffin et al. 2010) at 250, 350, and 500 μm in our previous work (Fan et al. 2016). Our ALMA observations show a marginal detection (∼4.2σ) of 3 mm dust continuum. The measured observed-frame continuum flux density is 0.140 ± 0.033 mJy (see Figure 1).

Zoom In Zoom Out Reset image size

For the multiwavelength SED analysis of W0533−3401, we use a forthcoming version (Y. Han et al. 2019, in preparation) of the Bayesian SED modeling and interpreting code bayesed¹² (Han & Han 2012, 2014, 2019). In the new version of bayesed, the stellar emission, dust attenuation, and dust emission can be consistently connected by assuming an energy balance, a technique similar to that employed in magphys (da Cunha et al. 2008) and cigale (Noll et al. 2009; Boquien et al. 2019). The stellar emission is modeled by using the Bruzual & Charlot (2003) SSP with a Chabrier (2003) initial mass function, an exponentially declining star formation history (SFH) and the Calzetti et al. (2000) dust attenuation law. The energy of stellar emission absorbed by dust is assumed to be totally re-emitted at the IR band, which is modeled by a graybody. The graybody model is defined as ${S}_{\lambda }\,\propto \left(1-{e}^{-{\left(\tfrac{{\lambda }_{0}}{\lambda }\right)}^{\beta }}\right){B}_{\lambda }({T}_{\mathrm{dust}})$, where B_λ is the Planck blackbody spectrum and λ₀ = 125 μm. Dust temperature T_dust and the emissivity index β are two free parameters. Finally, as in Fan et al. (2016), the active galactic nucleus (AGN) torus emission is modeled independently with the extensive database¹³ of 1,247,400 SEDs from the clumpy torus model (Nenkova et al. 2008a, 2008b). A k-dimensional tree-based nearest-neighbor searching technique has been employed to allow us to evaluate the clumpy torus model at any point in its 6D parameter space. Note that the clumpy model includes not only the torus dust emission, but also a part of the AGN accretion disk emission that is scattered into our line of sight or not absorbed by torus dust. Thus the clumpy model can provide a consistent description of AGN UV-to-millimeter SED. In total, the three-component model, including stars, AGN, and cold dust emissions, has 12 free parameters. The priors for them are summarized in Table 3.

Table 3.  The Model Parameters, Their Priors and Best-fitting Quantities for Multiwavelength SED Analysis

Name	Prior	Min	Max	Best-fitting Value
Stellar componenta
$\mathrm{log}({age}\,{\mathrm{yr}}^{-1})$	Uniform, ${age}\lt {{age}}_{{\rm{U}}}(z)$	5	10.3	${6.7}_{-0.1}^{+0.1}$
$\mathrm{log}(\tau \,{\mathrm{yr}}^{-1})$	Uniform	6	12	${9.1}_{-1.9}^{+1.8}$
$\mathrm{log}(Z/{{\rm{Z}}}_{\odot })$	Uniform	−2.3	0.7	${0.28}_{-0.38}^{+0.18}$
${A}_{{\rm{V}}}/\mathrm{mag}$	Uniform	0	4	${1.81}_{-0.08}^{+0.08}$
Graybody model
${T}_{\mathrm{dust}}/{\rm{K}}$	Uniform	10	100	${78.1}_{-5.2}^{+5.9}$
β	Uniform	1	3	${1.84}_{-0.13}^{+0.14}$
clumpy modelb
N0	Uniform	1	15	${5.9}_{-1.0}^{+1.3}$
Y	Uniform	5	100	${48.2}_{-28.2}^{+32.6}$
i	Uniform	0	90	${37.8}_{-23.8}^{+20.0}$
q	Uniform	0	3	${2.2}_{-0.4}^{+0.4}$
σ	Uniform	15	70	${62.3}_{-5.1}^{+8.5}$
${\tau }_{{\rm{V}}}$	Uniform	10	300	${22.2}_{-5.1}^{+5.4}$

Notes.

^aFour parameters of stellar components: age (age), the e-folding time for the exponentially declining SFH (τ), stellar metallicity (Z), and dust attenuation (${A}_{{\rm{V}}}$). ^bThe detailed description and definition of six free parameters for the clumpy model can be found in Nenkova et al. (2008a, 2008b).

Download table as:  ASCIITypeset image

In Figure 2, we show the best-fit three-component SED model (black solid line) with bayesed to the observed UV-to-millimeter SED of W0533−3401 (red points). With a newly developed function in the new version of bayesed, we can provide the confidence regions (CR) for our best-fit SED model. We plot the 68% and 95% CR with color-filled regions (cyan and purple, respectively). Absorbed stellar emission, AGN emission, and cold dust emission are shown by the green dotted line, the blue dashed line, and the gray dotted–dashed line, respectively. The derived properties have been shown in Table 4 and the results will be discussed in the next section.

Zoom In Zoom Out Reset image size

Table 4.  Derived Properties of W0533−3401

Parameter	Value	Unit	Note
${M}_{\star }\,$	(3.5 ± 0.9) × 1010	M⊙	(1)
SFR	6985 ± 3006	M⊙ yr−1	(2)
sSFR	200 ± 100	Gyr−1	(3)
LAGN	(7.0 ± 0.6) × 1013	L⊙	(4)
LGB	(3.5 ± 0.4) × 1013	L⊙	(5)
${L}_{\star }^{\mathrm{unabs}}$	(3.6 ± 0.4) × 1013	L⊙	(6)
${L}_{\star }^{\mathrm{abs}}$	(9.7 ± 1.0) × 1011	L⊙	(7)
${M}_{{{\rm{H}}}_{2}}$	(8.4 ± 0.5) × 1010	M⊙	(8)
${M}_{\mathrm{dyn}}$	1.6 × 1011	M⊙	(9)
${M}_{\mathrm{BH}}$	2.2 × 109	M⊙	(10)
Mdust	(8.9 ± 1.6) × 107	M⊙	(11)
Tdust	78.1 ± 5.6	K	(12)
β	1.84 ± 0.14		(13)
AV	1.81 ± 0.08	mag	(14)
${t}_{\mathrm{depl}}$	∼12–28	Myr	(15)

Note. (1): Stellar mass. (2): Star formation rate. (3): Specific star formation rate. (4): Bolometric luminosity of AGN component. (5): IR luminosity of cold dust. (6): Bolometric luminosity of attenuation-corrected stellar emission. (7): Bolometric luminosity of absorbed stellar emission. (8): Molecular gas mass. (9): Dynamical mass. (10): The central suppermassive black hole mass by assuming λ_Edd = 1. (11): Dust mass. (12): Dust temperature. (13): Dust emissivity index. (14): Dust attenuation. (15): The gas depletion timescales.

Download table as:  ASCIITypeset image

In Figures 3 and 4, we show the resulting line detection, both as spectrum and moment-maps. The CO $(3\mbox{--}2)$ emission line peaks at a frequency corresponding to a redshift of ${z}_{\mathrm{CO}(3\mbox{--}2)}\,=2.9026\pm 0.0003$. The detected CO $(3\mbox{--}2)$ line is spatially extended and shows a velocity gradient. We fit the velocity-integrated visibilities with an elliptical Gaussian function. The derived source size deconvolved from the beam is (0.73 ± 0.14) × (0.37 ± 0.14) arcsec² with P.A. = 126 ± 18. We measure the line flux ${I}_{\mathrm{CO}(3\mbox{--}2)}=2.01\pm 0.13$ Jy km s⁻¹. The rms is determined in the inner 10''of the moment map excluding the central part that contains the source itself. A single Gaussian fit to the continuum-subtracted CO $(3\mbox{--}2)$ spectrum gives CO $(3\mbox{--}2)$ FWHM = 566 ± 44 km s⁻¹. From the CO $(3\mbox{--}2)$ line flux, we can derive the CO line luminosity ${L}_{\mathrm{CO}(3-2)}^{{\prime} }\,=(8.4\pm 0.5)\times {10}^{10}$ K km s⁻¹ pc², by using Equation (3) in Solomon & Vanden Bout (2005). The continuum is marginally detected, as shown in Figure 1. The extend of the continuum emission is smaller than that of the CO $(3\mbox{--}2)$ line, which is likely due to a combination of the optical depth effects and the signal-to-noise ratio of continuum detection. The results are summarized in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can infer the molecular gas mass directly from the measured CO line luminosity. We adopt the excitation ratio between the CO $(3\mbox{--}2)$ and CO(1–0) line of r_32/10 = 0.8 as suggested by Banerji et al. (2017). The adopted ratio is intermediate between the typical values for submillimeter galaxies (SMGs, r_32/10 = 0.66) and optical quasars (r_32/10 = 0.97) from Carilli & Walter (2013), and consistent with the expectation that obscured quasars represent the transition phase from SMGs to unobscured quasars. The calculated CO(1–0) line luminosity is ${L}_{\mathrm{CO}(1\mbox{--}0)}^{{\prime} }=(10.5\pm 0.6)\times {10}^{10}$ K km s⁻¹ pc². We also adopt the CO-to-H₂ conversion factor, α_CO = 0.8 M_⊙(K km s⁻¹ pc²)⁻¹, which is suggested to be appropriate for starbursts and quasar hosts (Carilli & Walter 2013). Under these considerations, we obtain a molecular gas mass of ${M}_{{{\rm{H}}}_{2}}\,=(8.4\pm 0.5)\times {10}^{10}\,{M}_{\odot }$.

We show the CO $(3\mbox{--}2)$ velocity map in the right panel of Figure 3 and the position versus velocity (PV) diagram extracted along the major axis at P.A. = 330° in Figure 5. The CO $(3\mbox{--}2)$ emission line is marginally resolved. Both figures show the possible presence of a velocity gradient. Such a velocity gradient, which has also been observed in many other high-redshift quasars (e.g., Leung et al. 2017; Brusa et al. 2018; Feruglio et al. 2018; Talia et al. 2018), could be possibly the signature of a rotating disk of molecular gas, although the explanation is not unique. Especially, tentatively emission is seen at ∼3σ level extending to the west side of the source (see the left panel of Figure 3). This suggests that the gas-rich late-stage major merger and gas outflow are also possible for originating such a velocity gradient (e.g., Springel & Hernquist 2005; Hopkins et al. 2009; Ueda et al. 2014; Hung et al. 2015). However, due to limited depth and resolution, we cannot make a solid conclusion. Deeper observation would be required to distinguish them.

Zoom In Zoom Out Reset image size

Assuming that the gas is distributed in a rotating disk, we can investigate the kinematic properties of the molecular gas traced by the CO $(3\mbox{--}2)$ line, using ^3DBAROLO, a tool for fitting 3D tilted-ring models to emission-line data cubes (Di Teodoro & Fraternali 2015). We assume a disk model with three rings and a ring width of 02. We fix P.A. = 330° and adopt an inclination of 59° ± 14°, which is inferred from the observed ratio of minor to major axis. The derived rotation velocity (V_rot) and the intrinsic velocity dispersion σ are about ∼240 km s⁻¹ and 60 km s⁻¹, respectively. The ratio of rotation velocity to velocity dispersion V_rot/σ for W0533−3401 is about 4, which is close to the typical value ∼7 for the molecular gas in z > 1 star-forming galaxies (Tacconi et al. 2013). The large ratio V_rot/σ indicates that the molecular gas is turbulent, which is possibly due to a thick, dynamically hot disk in W0533−3401. A similar finding has been reported by Tadaki et al. (2018); they found a gravitationally unstable, rotating gas disk in an extreme starburst galaxy at z ∼ 4. The dynamical mass within the CO-emitting region can be estimated by applying the relation in Wang et al. (2013), ${M}_{\mathrm{dyn}}/{M}_{\odot }=1.16\times {10}^{5}\,\times {(0.75\times {\mathrm{FWHM}}_{\mathrm{CO}})}^{2}\times D/\sin i$, where D is the disk diameter in kiloparsecs from the CO $(3\mbox{--}2)$ measurement and i is the inclination angle. This gives M_dyn = 1.6 × 10¹¹ M_⊙.

Cold dust emission heated by stars has been modeled with a graybody model in our multiwavelength SED analysis (see Section 3.2 and Figure 2). The best-fit model gives the estimations of dust properties: IR luminosity L_GB, dust temperature T_dust, and the emissivity index β, which are listed in Table 4. The derived values are consistent with those in our previous work based on only IR SED decomposition (Fan et al. 2016). Dust mass can also be calculated by using the following equation:

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{D}_{{\rm{L}}}^{2}}{(1+z)}\times \displaystyle \frac{{S}_{{\nu }_{\mathrm{obs}}}}{{\kappa }_{{\nu }_{\mathrm{rest}}}B({\nu }_{\mathrm{rest}},{T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 1 }$$

where D_L is the luminosity distance, ${S}_{{\nu }_{\mathrm{obs}}}$ is the flux density at observed frequency ν_obs, ${\kappa }_{{\nu }_{\mathrm{rest}}}={\kappa }_{0}{(\nu /{\nu }_{0})}^{\beta }$ is the dust mass absorption coefficient at the rest frequency of the observed band, and $B({\nu }_{\mathrm{rest}},{T}_{\mathrm{dust}})$ is the Planck function at temperature T_dust. Adopting the best-fit values T_dust = 78.1 K, β = 1.84, and κ_{1 THz} = 20 cm² g⁻¹ which is the same as in Wu et al. (2014) and Fan et al. (2016), we derive the dust mass M_dust = (8.9 ± 0.5) × 10⁷ M_⊙. The small uncertainty of dust mass estimation only takes the uncertainties of the derived T_dust and β values into account. We note that the largest uncertainty can arise from the adopted ${\kappa }_{{\nu }_{\mathrm{rest}}}$ value, which can vary by over one order of magnitude at a certain frequency/wavelength. For instance, ${\kappa }_{850\mu {\rm{m}}}$ can vary from ∼0.4 to ∼11 cm² g⁻¹ in the literature (e.g., James et al. 2002; Draine 2003; Dunne et al. 2003; Siebenmorgen et al. 2014).

We derive the gas-to-dust mass ratio of W0533−3401, ${\delta }_{\mathrm{GDR}}=944\pm 77$, based on the estimations of ${M}_{{{\rm{H}}}_{2}}$ and M_dust. The derived δ_GDR value of W0533−3401 is thus about one order of magnitude higher than the typical value ∼50–150 derived for the Milky Way (Jenkins 2004), the local star-forming galaxies (SFGs) and ultraluminous IR galaxies (ULIRGs) at solar metallicity (Draine et al. 2007; Rémy-Ruyer et al. 2014) and high-redshift SMGs (Magnelli et al. 2012; Miettinen et al. 2017). The high δ_GDR value in W0533−3401 may be due to several possible reasons: the uncertainty of dust mass estimation, the low efficiency of dust formation, and/or the high efficiency of dust destruction. The dust mass derived by using a graybody model is about half of that derived by the Draine & Li (2007) model (Magdis et al. 2011). This will result in an overestimation of δ_GDR by a factor of two. The δ_GDR value is reported to increase with the decreasing metallicity and the increasing redshift (e.g., Rémy-Ruyer et al. 2014; Miettinen et al. 2017). It is possible that W0533−3401 has a low, subsolar metallicity. It is also possible that dust destruction in W0533−3401 may be efficient due to the strong radiation field from massive young stars and AGN, and the supernova shock waves which are expected to be frequent in this maximum-starburst galaxy (Jones 2004).

Based on our best-fit SED model presented in Section 3.2, we derive the stellar mass and star formation rate (SFR) of W0533−3401 and list them in Table 4. The stellar mass M_⋆  is derived by adopting an exponentially declining SFH, which is represented by a young stellar population. We check if there is an old stellar population which is omitted by the present SED fitting procedure. We add an SSP with an age of 1 Gyr to our multiwavelength SED model. The result is that the contribution of the old SSP to the derived stellar mass is not larger than 20%, though the constraint from the SED fitting is loose. Given M_⋆  = 3.5 × 10¹⁰ M_⊙ and ${M}_{{{\rm{H}}}_{2}}=8.4\times {10}^{10}\,{M}_{\odot }$, we can calculate the molecular gas fraction ${f}_{\mathrm{gas}}={M}_{{{\rm{H}}}_{2}}/({M}_{{{\rm{H}}}_{2}}\,+{M}_{\star }\,)=0.71$, which indicates that W0533−3401 is gas-rich. The derived SFR is ∼7000 M_⊙ yr⁻¹. The specific SFR or sSFR = 200 Gyr⁻¹ of W0533−3401 is over one order of magnitude higher than the star-forming main-sequence (MS) at z ∼ 3 (Speagle et al. 2014), suggesting that it is a maximum starburst. The uncertainty of SFR is large (∼3000 M_⊙ yr⁻¹) due to a wide range of possible SFH. By using CO line luminosity, we can estimate SFR in an individual way. First, we convert the measured ${L}_{\mathrm{CO}(3-2)}^{{\prime} }$ to ${L}_{\mathrm{CO}(5-4)}^{{\prime} }$ by taking ${L}_{\mathrm{CO}(5-4)}^{{\prime} }/{L}_{\mathrm{CO}(3-2)}^{{\prime} }=0.7$ which is intermediate between the typical values for SMGs and quasars (Carilli & Walter 2013). Then we estimate FIR luminosity using the high-J CO versus FIR luminosity relation presented in Liu et al. (2015). Using the relation SFR (M_⊙ yr⁻¹) = 4.5 × 10⁻⁴⁴ L_FIR (erg s⁻¹), we calculate SFR ∼ 3000 M_⊙ yr⁻¹ (Kennicutt 1998), which is generally consistent with the best-fit SED result. From the derived SFR and molecular gas mass, we can estimate the gas depletion timescale t_depl = M_gas/SFR. We infer t_depl ∼ 12–28 Myr using SFRs based on the best-fit SED result and the CO line luminosity, respectively. The gas depletion timescale of W0533−3401 is similar to other obscured quasars (e.g., Aravena et al. 2008; Brusa et al. 2018), but is much shorter than MS galaxies and SMGs (Bothwell et al. 2013; Sargent et al. 2014), indicating that W0533−3401 as an obscured quasar is consuming its residual gas more rapidly.

In Figure 6, we report the star formation efficiency (SFE), traced by the ratio of IR to CO(1–0) line luminosities, as a function of IR luminosity for W0533−3401 and the compiled samples of SMGs, unobscured and obscured quasars at z > 1 in Perna et al. (2018). For all three samples, SFE shows a positive correlation with IR luminosity. W0533−3401 shows a high SFE = ${L}_{\mathrm{IR}}/{L}_{\mathrm{CO}}^{{\prime} }=464\pm 56\ {L}_{\odot }/$(K km s⁻¹ pc²), which is well above the best-fit relation and 1σ scatter for massive MS galaxies (Sargent et al. 2014). Similar to W0533−3401, both unobscured and obscured quasars at z > 1 show a higher SFE at a given IR luminosity than MS galaxies and SMGs. Higher SFE in unobscured and obscured quasars relative to MS galaxies and SMGs is possibly due to the presence of starburst activity and/or the depletion of cold gas by AGN feedback in quasar host galaxies. The former enhances L_IR and the latter reduces ${L}_{\mathrm{CO}}^{{\prime} }$. It is possible that L_IR can be overestimated in quasars, as AGN-heated dust on kiloparsec scales can contribute significantly to IR luminosity (Duras et al. 2017; Symeonidis 2017).

Zoom In Zoom Out Reset image size

We derive the AGN bolometric luminosity L_AGN of W0533−3401 by integrating the AGN component of the best-fit UV-to-millimeter SED (see Section 3.2). Although we do not have a direct measurement of the black hole mass in W0533−3401, we can make a rough estimate from the derived L_AGN. By assuming Eddington ratios of 0.3, 1.0, and 3.0, the corresponding black hole masses are 7.3 × 10⁹, 2.2 × 10⁹, and 7.3 × 10⁸M_⊙, respectively. Wu et al. (2018) reported the black hole masses and Eddington ratios of five luminous obscured quasars at z ∼ 2, which are taken from the same parent sample of W0533−3401, based on broad H_α lines. They found that the average black hole mass is about 10⁹ M_⊙ and the derived Eddington ratios are close to unity. Tsai et al. (2018) also reported the measurement of M_BH and L_Edd for an extremely luminous, obscured quasar W2246−0526 at z = 4.6, which is selected with the same criteria as W0533−3401. They found that the central SMBH in it is growing rapidly by accreting at a super-Eddington ratio (λ_Edd = 2.8). Adopting the Eddington ratio λ_Edd = 1.0 for W0533−3401 will be a reasonable approximation.

From SED-based stellar mass M_⋆  and the estimated SMBH mass, we can infer the black hole to stellar mass ratio of W0533−3401. In Figure 7, we plot the observed M_BH/M_⋆  as a function of redshift for W0533−3401 and several other samples in the literature. If assuming λ_Edd = 1.0, the black hole to stellar mass ratio of W0533−3401 (M_BH/M_⋆  = 0.063) is close to that of an X-ray selected unobscured quasar, CID-947, which has the highest M_BH/M_⋆  (1/8) at z ∼ 3.3 known so far. A low Eddington ratio in W0533−3401, for instance λ_Edd < 0.3, seems not likely, otherwise M_BH/M_⋆  of W0533−3401 will be higher than 0.2. Even if taking λ_Edd = 1.0, the inferred M_BH/M_⋆  of W0533−3401 is not only over one order of magnitude higher than the typical values ∼0.0002–0.0005 in the local universe (Kormendy & Ho 2013), but also ∼5 times higher than the expected value by the evolutionary trend of M_BH/M_⋆  (McLure et al. 2006; Peng et al. 2006; Merloni et al. 2010; Targett et al. 2012; Bongiorno et al. 2014; Matsuoka et al. 2018). We derive the black hole mass growth rate ${\dot{M}}_{\mathrm{BH}}$ of W0533−3401 from the bolometric luminosity by using the relation ${L}_{\mathrm{AGN}}=(\eta {\dot{M}}_{\mathrm{BH}}{c}^{2})/(1-\eta )$ and adopting $\tfrac{\eta }{1-\eta }=0.1$. We infer ${\dot{M}}_{\mathrm{BH}}=49$ M_⊙ yr⁻¹, which suggests a rapid growth of the central SMBH. We then calculate the ratio between the black hole mass growth rate and SFR and obtain ${\dot{M}}_{\mathrm{BH}}$/SFR = 0.007, which is close to the local M_BH/M_⋆  values. The result indicates that W0533−3401 is evolving toward the evolutionary trend of M_BH/M_⋆  (the dashed line in Figure 7).

Zoom In Zoom Out Reset image size