© ESO 2021

1. Introduction

Evidence compiled over the past two decades suggests that there is a broad relationship between the growth of massive galaxies and their central super-massive black holes (SMBHs); namely, that there is a connection between the histories of star formation rate (SFR) in the former and the accretion rate of the latter. This is evidenced by the approximately synchronous co-evolution of these two galaxy properties over cosmic time (Madau & Dickinson 2014, and references therein), which suggests that some sort of self-regulating process must be at work. Indeed, the high-end of the galaxy mass function displays a steep slope (Muzzin et al. 2013; Davidzon et al. 2017) that requires a mechanism able to quench massive galaxies on short time scales. Feedback from luminous active galactic nuclei (AGN), also regularly labeled as “quasar feedback”, is currently the preferred physical process with which theoretical models of galaxy evolution and cosmological simulations attempt to reproduce observations (Croton et al. 2006; Anglés-Alcázar et al. 2017a; cf., Muzzin et al. 2012). Quasar-mode feedback, which manifests itself in the form of galactic gas outflows, has been identified in a large number of z ≳ 2 optically and infrared-selected quasars (OpQs and IRQs, respectively) via fast-moving (v_out ≥ 1000 km s⁻¹), ionized-gas ([O III]λλ4959,5007 Å) winds (Zakamska et al. 2016; LaMassa et al. 2017; Bischetti et al. 2017; Wu et al. 2018; Perrotta et al. 2019; Temple et al. 2019; Jun et al. 2020; Finnerty et al. 2020).

Access to the kinematics of the interstellar medium (ISM) is important not only to identify the potential evacuation of material due to powerful quasar feedback, but also to characterize the dynamical state of the host. Investigating how gas is transferred and (re)cycled in high-redshift galaxies, as it is accreted from the cosmic web and expelled by the central AGN, is key to establishing whether outflows are actually able to quench –or at least delay– further star formation. Moreover, it allows us to assess the role of turbulence as an additional heating mechanism of the ISM gas in dense environments, where galaxy collisions may be frequent.

The Atacama Large Millimeter/sub-millimeter Array (ALMA) enables these studies by providing access to a range of the electro-magnetic spectrum of galaxies that is not affected by dust extinction. Thus, the physical properties and kinematics of the ionized, neutral, and molecular phases of the ISM can be studied without the degeneracies introduced by self-absorption (with the exception of very extreme cases; e.g., Scoville et al. 2015). This has allowed the systematic detection and analysis of the ISM of high-redshift quasars and star-forming galaxies, mostly via the fine-structure [C II]λ157.7 μm (²P_3/2→²P_1/2) ([C II]₁₅₈) emission line and its underlying dust continuum, which are red-shifted into the ALMA highest frequencies bands at z ≳ 1 (e.g., Wang et al. 2013; Decarli et al. 2018; Venemans et al. 2018; Bischetti et al. 2018; Díaz-Santos et al. 2018; Hodge et al. 2019; Leung et al. 2019; Jones et al. 2020, among many others).

Within the general framework of galaxy evolution, nearby quasars are thought to be the product of gas-rich, massive galaxy mergers (Sanders et al. 1988; Hopkins et al. 2008) –a short-lived, ultra-luminous phase lasting a few tens of Myr (Hopkins et al. 2005; Trainor & Steidel 2013). This phase is subsequently divided into two stages, separated by the expulsion of the gas and dust that has been accreted onto the central engine: obscured quasars harbor SMBHs still embedded in a dusty cocoon that absorbs most of their UV and optical light and re-emits it in the IR, while optically visible quasars have already cleared out the surrounding medium through the action of SMBH feedback and are visible in the optical.

An important population of high-redshift IRQs are hot, dust-obscured galaxies (Hot DOGs). Hot DOGs were first discovered by the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), selected to be strongly detected at 12 and 22 μm, but weakly or not detected at 3.4 and 4.6 μm (Eisenhardt et al. 2012; Wu et al. 2012). Most of their extreme bolometric luminosity, L_bol ≳ 10¹³ L_⊙, is produced by accretion onto their central SMBH. The host galaxies are on average less massive than what would be expected from such hyper-luminous AGN, suggesting they are radiating at (or even above) the Eddington limit (Assef et al. 2015; Wu et al. 2018; Tsai et al. 2018).

Statistically, Hot DOGs seem to be located in galaxy over-densities (Jones et al. 2014; Assef et al. 2015) and host powerful ionized outflows (Wu et al. 2018; Jun et al. 2020, 2021; Finnerty et al. 2020), which is in agreement with the merger-driven model of quasars and their expected feedback, and also suggesting that they may be at the centers of protoclusters. Further evidence comes from ALMA [C II]₁₅₈ observations of the most luminous obscured quasar and Hot DOG known (L_bol = 3.5 × 10¹⁴ L_⊙; Tsai et al. 2015), namely WISE J224607.55–052634.9 (W2246–0526 hereafter; at z = 4.601), which revealed a highly turbulent ISM, with a full width at half maximum of FWHM_[CII] ≥ 500 km s⁻¹ spread across the entire host galaxy, over ∼2.5 kpc (Díaz-Santos et al. 2016). Furthermore, recent deep ALMA imaging of the dust continuum at λ_rest = 212 μm of W2246–0526 shows streamers and bridges of dust extending up to ∼35 kpc in projected size, which seem to connect at least three companion galaxies with the central Hot DOG. These results suggest the most luminous obscured quasars may be interacting systems –the result of ongoing merger-driven peaks of SMBH accretion and massive galaxy assembly in the early Universe (Díaz-Santos et al. 2018).

In this work, we study a sample of seven of the most luminous Hot DOGs, at redshifts of z ∼ 3.0−4.6, observed with ALMA in the [C II]₁₅₈ line. In Sect. 2, we introduce the ALMA observations, ancillary data, and comparison samples of other high-z OpQs and IRQs. In Sect. 3, we describe the analysis of the data, and in Sect. 4, we present results from the integrated, line and continuum measurements, including the assessment of the global ISM conditions. In Sect. 5, we introduce the results from the kinematic analysis and resolved morphology, and in Sect. 6 we place the Hot DOG hosts in the context of the star-formation galaxy’s main sequence. Finally, in Sect. 7, we speculate on the implications of these results and discuss the effect of the SMBH on the physics of the surrounding gas and dust, and we consider how Hot DOGs fit within the broader, merger-driven quasar paradigm. Our conclusions are summarized in Sect. 8. Throughout the paper, the following flat cosmology is applied: Ω_M = 0.28, Ω_Λ = 0.72, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Observations

2.1. The sample

Our sources are drawn from a catalog of over 2000 Hot DOGs, photometrically selected from the WISE All-Sky Data Release (Cutri et al. 2012) using the selection criteria presented by Eisenhardt et al. (2012) (see also Assef et al. 2015 for details). A subset of the sample have secure redshifts from rest-UV spectroscopic observations that will be presented by Eisenhardt et al. (in prep.). We obtained ALMA observations of seven targets selected from this catalog that have declinations accessible to ALMA and are at rest-UV redshifts that put the [C II]₁₅₈ emission line ([C II], hereafter) within the frequency range of ALMA’s band 8 and 7, in regions of high atmospheric transparency. This pre-selection implies the targets probe the high-redshift end of the Hot DOG population (i.e., concentrated in the 3 ≲ z ≲ 4.6 range; see Fig. 1 of Assef et al. 2015) and thus also the high end of the luminosity distribution (in the “extremely” luminous, EL, regime), as they all have L_bol ≥ 10¹⁴ L_⊙ (Tsai et al. 2015). Basic information about the ALMA Hot DOG sample can be found in Table 1. The observations for the most distant and most luminous Hot DOG currently known, W2246–0526, have already been presented in Díaz-Santos et al. (2016).

2.2. ALMA observations and data processing

ALMA targeted the [C II] emission line for the sample of seven Hot DOGs under programs 2013.1.00576.S (W0116–0505, W0134–2922, W0220+0137 and W2246–0526) and 2015.1.00612.S (W0236+0528, W0410–0913 and W0831+0140), both PI: R. Assef. The observations were carried out following the standard observatory calibration procedure under good weather conditions (i.e., a typical precipitable water vapor, PWV ≃ 0.4−0.6 mm). More details regarding these observations can be found in Table 2. For all sources except W0236+0528, the data were acquired in a single execution block (EB). For W0236+0528, two EBs were obtained, which were concatenated into a single measurement set.

Owing to their redshifts, W2246–0526 was observed in band 7 and the remaining galaxies in band 8. While the requested sensitivity of the observations was 1 mJy beam⁻¹ over 100 km s⁻¹ channels for the band 8 targets and 1 mJy beam⁻¹ over 20 km s⁻¹ channels for W2246–0526, the final data were deeper in most cases. Similarly, the requested angular resolution was 0.24″ for the band 8 targets and 0.32″ for W2246–0526, but higher resolutions (< 0.2″) were achieved for W0236+0528 and W0831+0140 (see Tables 3 and 4).

The center of the reference spectral window (SPW) was tuned to match the expected observed (redshifted) frequency of [C II] for each source. The redshifts were based on optical spectroscopic (rest-frame UV) observations, whose reduction and analysis will be presented and discussed in Eisenhardt et al. (in prep.). The spectral region around the Lyman αλ1216 Å line is shown in Fig. 1 for all seven galaxies except W0410–0913, for which Lyman αλ1216 Å landed in the gap between detector chips of the instrument (Keck/DEIMOS). The rest-UV redshfits, z_UV, were estimated based on the central wavelength(s) of the available line(s), which were fitted using a Gaussian profile. The lines used depend on the source spectrum, but for the Hot DOGs included here are generally based on Ly α, the Ly β/O IV blend, and N V. C IV was excluded except for W0410–0913, whose UV redshift is also based on the He II line (see discussion in Sect. 4.1). The reference SPW was normally the reddest of the two within the side-band, which in turn was the reddest of the two side-bands. However, adjustments to this configuration were made due to strong absorption lines and the cut-off of the atmospheric band. The SPWs not covering the [C II] line were used to sample the dust continuum.

Fig. 1. Optical spectra of the Hot DOG sample around the Ly-α emission line, whose rest-frame is set at zero velocity using zALMA (Table 1). The feature located at ∼5000 km s−1 is the N Vλ1240 Å emission line. The spectrum of W0410–0913 is not shown since Ly-α landed in a gap between detector chips.

The scripts provided by ALMA as part of the data delivery were executed using the appropriate version of the Common Astronomy Software Application (CASA; McMullin et al. 2007) to process each specific raw dataset and calibrate the visibilities. The resulting measurement sets were inspected to check for possible issues in the amplitude, band-pass and phase calibrations. All four SPWs were usable for all sources and no additional flagging was needed.

We also used CASA (v5.4.0-70) to process and clean the ALMA products. The cleaning algorithm was run using three different Briggs weighting schemes for the uv visibility plane: (1) robust parameter set to 0.5 (nearly half-way between uniform and natural weighting); (2) robust parameter set to 2 (similar to natural weighting); and (3) robust parameter set to 2 and a uv-taper with a resolution approximately twice that provided by scheme (2). The cleaning area was usually a circle of 1″ diameter centered on the object, but elliptical shapes of varying sizes were used when necessary. The angular size (FWHM) of the restored, synthesized beams at the frequency of the [C II] line range from ∼0.35″ for W2246–0526 to ∼0.15″ for W0236+0528 and W0831+0140 (considering the weighting scheme (2)). The spaxel size of the cubes was selected such that the beam is sampled across by 5−10 spaxels. The native channels were averaged to a width of ∼50 km s⁻¹ to maximize the signal-to-noise ratio (S/N) while keeping the final spectral elements sufficiently small to derive kinematic information. The line and continuum cubes were cleaned down to a level of 2σ_rms over the general background¹, where σ_rms is the average rms per channel of the respective SPW², calculated after masking out the object with the same aperture used later for the cleaning process. More technical details about the observations can be found in Table 2.

2.3. Comparison samples

2.3.1. High-redshift Quasars

For comparisons with high-redshift optically selected quasars (OpQs) we use the [C II] data compiled by Decarli et al. (2018), which includes originally published sources supplemented with OpQs from the literature, all at z ≳ 6, including Maiolino et al. (2005), Walter et al. (2009), Venemans et al. (2012, 2016, 2017), Wang et al. (2013, 2016), Willott et al. (2013, 2015, 2017), Bañados et al. (2015), and Mazzucchelli et al. (2017). In addition, we use the results from the dust continuum data analysis published in Venemans et al. (2020) for the same sample of OpQs.

We also use the sample of high-redshift IR quasars (IRQs) from Trakhtenbrot et al. (2017), which has been recently expanded and updated by Nguyen et al. (2020). We note that while these are also optically identified quasars found in the SDSS, before the sample was observed by ALMA, the sources were selected a priori to have significant (if not strong) far-IR continuum emission as detected by Herschel. Most of these IRQs are located at z ≃ 4.8.

2.3.2. High-redshift star-forming galaxies

To put the Hot DOG hosts in context with purely star-forming galaxies, we use the ALMA Large Program to INvestigate [C II] at Early times survey (ALPINE; Le Fèvre et al. 2020, and references therein), which observed more than a hundred galaxies at 4.4 ≲ z ≲ 5.9 in [C II] and dust continuum emission. ALPINE sources are representative of the star-forming galaxy main sequence (MS) at their respective redshifts.

2.3.3. GOALS: Nearby IR galaxies and AGN

We use the Great Observatories All-sky LIRG Survey (GOALS; Armus et al. 2009) as the reference sample for nearby galaxies. GOALS is a complete sample of luminous and ultra-luminous infrared galaxies ((U)LIRGs; L_IR = 10¹¹⁻¹² and ≥ 10¹² L_⊙, respectively) in the nearby Universe (z ≤ 0.1). The sample has extensive multi-frequency photometric and spectroscopic coverage from space- and ground-based facilities, from X-rays (Iwasawa et al. 2011; Torres-Albà et al. 2018), through ultraviolet (Howell et al. 2010), optical (Rich et al. 2015), near-IR (Inami et al. 2018), mid-IR (Stierwalt et al. 2013, 2014; Díaz-Santos et al. 2010, 2011), far-IR (Díaz-Santos et al. 2017), and sub-mm (Herrero-Illana et al. 2019), to radio wavelengths (Barcos-Muñoz et al. 2017; Linden et al. 2017). GOALS bridges the gap between MS galaxies and the most powerful starbursts and AGN in the local Universe, and therefore it is useful for comparisons with any IR galaxy population identified at high redshift.

2.3.4. Data homogenization

Throughout the paper, physical areas are calculated as , where R_{160 μm, eff} is the intrinsic (PSF-corrected) 160 μm dust continuum emission half-light radius (=FHWM/2, for a 2D Gaussian). These were obtained from Herschel photometry (Lutz et al. 2016; Chu et al. 2017) for the local (U)LIRG sample, and from ALMA observations for the IRQ sample of Nguyen et al. (2020) and the Hot DOG sample presented here.

Venemans et al. (2018) first reported dust continuum sizes for the OpQ sample of Decarli et al. (2018). However, most sources were unresolved or the data had limited S/N with which to calculate an accurate size of the emitting region. We thus use the updated values recently published in Venemans et al. (2020), which are based on higher angular resolution data. The continuum sizes reported in the latter work tend to be significantly smaller than those listed in the former (sometimes by a factor of ∼3), likely reflecting the improved resolution of the more recent observations and the fact that a fraction of the more extended continuum emission may have been resolved out. Due to a lack of data coverage at shorter IR wavelengths than the ALMA observations, only L_FIR values are available for some of the OpQs. In these cases, we use the average L_IR/L_FIR ≃ 1.7 ratio found for nearby (U)LIRGs to estimate their L_IR.

IR luminosities are available for only the 23 sources in the ALPINE survey that were detected in dust continuum (Béthermin et al. 2020). In order to plot the remainder of the sources on our figures, we estimate a “virtual” L_IR from their UV-based SFRs using the standard conversion given in Murphy et al. (2009). We emphasize that these are not IR luminosities per se, but rather represent what the L_IR of these galaxies would be if all their star formation was obscured. Fujimoto et al. (2020) do not provide sizes for the dust continuum detections, and therefore we use the [C II] sizes to calculate luminosity surface density quantities. Since these are mostly star-forming galaxies, no correction is performed to their [C II] sizes, such that the line and continuum sizes are assumed to be the same.

3. Analysis

3.1. Measurement of line and continuum fluxes

To calculate the line fluxes of our Hot DOGs, the cubes were continuum-subtracted in the uv-plane using representative frequency ranges from the continuum SPW adjacent to the reference SPW, and in some cases from the SPWs of the other side-band, when no slope due to the Rayleigh-Jeans continuum from dust was identified. Both the spectrum of the dust continuum and that of the continuum-subtracted line were extracted using an aperture of 1″ diameter centered at the peak of the galaxy emission (an aperture of 2.5″ was used for W0410–0913 due to its large extent). Figure 2 shows the spectra for all four SPWs.

Fig. 2. ALMA spectra of the luminous Hot DOG sample, including the four SPWs (black histogram). The average size of the synthesized beam of the cube is noted in the legend at the top left, as well as the radius of the aperture used for the extraction of the spectra, which is also shown as a dashed circumference in Fig. 6. The continuum spectrum (obtained from the line-subtracted cube) underlying the emission line in each source is shown in green, and the Gaussian fit to the [C II] line is shown in blue. The polynomial fit to the overall line-subtracted continuum as well as its ±1σ uncertainty are shown in red, with solid and dashed lines, respectively. The solid gray lines are the atmospheric transmission curve for a precipitable water vapor of 1 mm3, where the maximum transmission within the displayed spectral range has been normalized to the maximum y-axis value of the whole velocity range displayed. The vertical dashed line shows the expected redshift of the line based on the rest-UV spectrum, zUV.

The emission line profile³ was fitted by a Gaussian function with the central position, peak intensity and FWHM left as free parameters. The luminosity was obtained from the analytic expression of the Gaussian using the best-fit parameters, as well as by directly integrating the observed continuum-subtracted spectrum over a frequency range defined by ±2.5σ from the line’s peak frequency as obtained from the fit, where σ is the standard deviation of the Gaussian. Both line luminosities are reported in Table 3. In the case of W0410–0913, the line center is redshifted out of the frequency range covered by the observations and only the bluest wing of the line is detected (see Fig. 2). In this particular case, the fit to the line was performed using the available portion of the profile and fixing the central position of the Gaussian to the redshift obtained from the CO emission line (discussed in Sect. 4.1).

We also report two observed-frame continuum flux density measurements for each object. The first is calculated by obtaining the moment 0 of the entire cube after subtracting the emission line cube from the reference SPW cube. The second measurement is derived from a first-order polynomial fit to the extracted spectrum of the entire cube and interpolated to the central wavelength (see Fig. 2). These measurements can be found in Table 4 and are all consistent within the uncertainties.

The velocity and dispersion maps of the line were obtained by computing the moment 1 and 2 statistics of the continuum-subtracted cubes. For consistency and completeness, we re-analyze here the data of W2246–0526 that was first presented in Díaz-Santos et al. (2016). All the newly derived quantities are within the uncertainties of those given in the previous work. The maps of the continuum flux density and line flux, as well as those of the velocity field and dispersion of the line (when detected; see next section) are presented in Figs. 6, 7, 9, and 10, respectively.

3.2. Sizes and extended structure

The intrinsic sizes of the line and continuum emission were measured by fitting a 2D-Gaussian to the sources, from which the FWHM of the beam is subtracted in quadrature. The values are reported in Tables 3 and 4, respectively. We note that because of the limited number of degrees of freedom provided by a Gaussian function to describe the complex spatial variations seen in the flux distribution of some of our sources, these luminosity-weighted Gaussian fits are representative only of the “core emission” of the Hot DOGs. That is, even though the cores of some of the galaxies are compact or close to be unresolved, there may still be significant low surface brightness emission that could extend over a region several times larger than the core sizes derived from the Gaussian fitting (see Figs. 6 and 7).

In order to quantify the contribution of these potential extended structures surrounding or underlying the central, unresolved point source of each Hot DOG, we estimated the fraction of extended emission (f_ext) of both the line and the dust continuum. This procedure is performed in the plane of the sky only for those sources with robust line and continuum detections (see below) to take advantage of the high fidelity of the ALMA images. The f_ext is calculated by iteratively scaling the well-defined, clean PSF image of the observations to the [C II] and dust continuum images until finding the ratio that minimizes the dispersion in the central region of the residual map, within an aperture defined by an ellipse with semi-major and minor axes equal to those of the PSF. After the optimal scaling is found, the unresolved emission is subtracted from the source image and the f_ext is defined as the fraction of the flux remaining in the residual image, measured in a fixed circular aperture equal to that used in the original image. The aperture is set manually for each Hot DOG with the goal of encompassing the total flux (dashed circumferences in Fig. 6). The only exception is W2246–0526, for which we defined an aperture that purposefully avoids the companion galaxy located 1″ north-east of the central Hot DOG. The f_ext values are provided in Tables 3 and 4 for the line and continuum emission, respectively.

Finally, we also created radial profiles of the line and continuum flux surface densities for each source. The radial profiles are computed by measuring the emission within increasingly expanding ellipse-shaped rings with the same major-to-minor axis ratio and position angle (PA) as the best 2D-Gaussian fit to the source. Profiles are calculated for the source, PSF, scaled PSF, and residual images (see above). The profiles of the most compact and extended Hot DOGs in line and continuum emission are presented in Fig. 8.

4. Integrated emission properties and morphology

4.1. Line detection, velocity shifts and luminosities

The [C II] line is robustly detected (S/N ≥ 20) in four out of the seven Hot DOGs in the sample: W0134–2922, W0220+0137, W0831+0140, and W2246–0526, the latter of which was already reported in Díaz-Santos et al. (2016). Figure 2 shows the spectra of all sources in all four SPWs. We know that for two sources, W0116–0505 and W0410–0913, the line is redshifted completely and partially, respectively, out of the frequency range covered by the reference SPW (and side-band) used to target [C II]. This is based on additional ALMA observations later obtained for these two sources, which targeted the mid-J CO(4→3) and CO(6→5) transitions (González-López et al., in prep.), respectively, as well as on the shallower observations of CO(4→3) analyzed in Fan et al. (2018) for W0410–0913. In addition, JVLA observations of the CO(1→0) line of both sources have also been recently presented in Penney et al. (2020). The detection of the mid-J CO transitions with ALMA owes to the particular set up of these observations, in which the reference SPW dedicated to the line was not the reddest but instead the bluest of the side-band. Moreover, the lower frequency of these lines compared to [C II] enables a larger bandwidth in velocity space. For both galaxies, the mid-J CO lines were clearly detected in the reddest SPW, which was originally dedicated to probe the dust continuum, allowing for a new redshift determination based on these additional far-IR data, z_CO. The velocity offsets of z_CO with respect to z_UV for W0116–0505 and W0410–0913 are Δv ≃ 1250 ± 73 and 2487 ± 66 km s⁻¹, respectively. The redshifts derived from the [C II] line, z_[CII], for those galaxies in which the line is identified and the profile fully sampled have velocity offsets with respect to z_UV of Δv ≃ −37 ± 75, 989 ± 74, 49 ± 62 and −64 ± 55 km s⁻¹ for W0134–2922, W0220+0137, W0831+0140, and W2246–0526, respectively (see Table 3). Therefore, most of the newly derived redshifts for our Hot DOG sample indicate positive systemic velocities with respect to those inferred from the optical emission lines, with an average offset of Δv ≃ 779 ± 412 km s⁻¹.

We tentatively detect the [C II] line in the spectrum of W0236+0528 at roughly the frequency (−161 km s⁻¹) where it was expected to appear, based on the z_UV (see Fig. 2). However, there seems to be another spectral feature in emission detected at a S/N ≃ 4 in the adjacent, redder SPW, with a velocity offset of ≃950 km s⁻¹. Given the mostly positive velocity offsets found for other Hot DOGs in the sample, we cannot discard that this feature could actually be the [C II] line. Unfortunately there is no additional spectroscopic information for this source, and the noise in the entire side-band is large. We attempted to create an image of the line using two velocity ranges encompassing the frequencies of the two possible feature identifications, but no source could be found in either. However, using a uv-tapering of 0.5″ to clean the cube and generating the same line images delivers a S/N ≃ 3 detection at the SPW where the line was expected, while nothing can be found in the redder SPW. Given the large uncertainties, we thus consider this only as a tentative detection and provide the line luminosity in Table 3, but we refrain from including this source in the figures or discussion throughout the remainder of the paper.

Redshifted far-IR emission lines with respect to rest-frame UV emission features have also been observed in a significant fraction of high-z OpQs (see Vietri et al. 2018, for an optical line emission perspective). Decarli et al. (2018) showed that the average of the velocity offsets between the [C II] and optical lines including all samples is 620 ± 8 km s⁻¹ (i.e., [C II] is also frequently redshifted). We compare our velocity shifts with those of OpQs as well as with the IRQs from Nguyen et al. (2020) in Fig. 3. For reference, we also show the velocity offsets of normal star-forming galaxies at z ∼ 4−6. While the [C II] line is, on average, redshifted with respect to z_UV in all the quasar samples, Hot DOGs and the IRQ sample seem to lack sources that are significantly blueshifted, in contrast to what is found in a number of OpQs or in general in normal star-forming galaxies, which tend to be slightly blueshifted, with Δv between ≃0 and −500 km s⁻¹.

Fig. 3. Velocity shift of zALMA with respect to zUV, Δv, as a function of zALMA (see Table 3). Hot DOGs are displayed as red squares. The z ≃ 4.7 IRQs from Trakhtenbrot et al. (2017) and Nguyen et al. (2020) are represented as yellow pentagons. The compilation of z ≳ 6 OpQs from Decarli et al. (2018) are displayed as green stars when the optical redshift is based on the Mg II emission line, and as blue stars when other lines or methods are used. Star-forming galaxies at 4 < z < 6 from the ALPINE survey (Cassata et al. 2020) are shown as pink pentagons.

The integrated [C II] luminosities of the Hot DOG sample range between ∼3 × 10⁹ L_⊙ and ∼10¹⁰ L_⊙ (see Table 3)⁴. This is roughly similar to the range exhibited by OpQs at z ≳ 4.7 (e.g., Trakhtenbrot et al. 2017; Willott et al. 2017) and typical, UV-selected star-forming galaxies at z ∼ 5.1−5.7 (e.g., Capak et al. 2015; Le Fèvre et al. 2020), and it is several orders of magnitude larger than normal and sub-L^* star-forming galaxies at z ∼ 6−7 (Ouchi et al. 2013; Ota et al. 2014; Maiolino et al. 2015; Knudsen et al. 2016). The line FWHMs in the Hot DOG sample are always in excess of ∼500 km s⁻¹, which is usually larger than the FWHMs measured in other high-redshift sources (see discussion in Sect. 5), perhaps with the exception of a few OpQs (Bischetti et al. 2019; Izumi et al. 2019) and the massive dusty star-forming galaxy HXMM05 at z = 2.985 studied in Leung et al. (2019), which shows a comparable line width to the Hot DOG population with FWHM ∼ 600 km s⁻¹. The large, ∼730 km s⁻¹ FWHM measured in W0410–0913, even if estimated using a small portion of the blue wing of the line, seems to be in good agreement with that calculated by Fan et al. (2018), of ∼750 km s⁻¹, based on a clear detection of the CO(4→3) emission line.

4.2. [C II] luminosities and Ly-α equivalent widths

Harikane et al. (2018) reported an anti-correlation between the luminosity of [C II] and the equivalent width (EW) of the Ly-α line for a large sample of Lyman alpha emitters (LAEs) at z ∼ 5−7. While the authors do not provide a physical interpretation for this anti-correlation, Fig. 4 shows that luminous Hot DOGs lie well above the trend found for the LAEs (gray circles; the best fit is represented by a solid line), in a distinctively different region of the parameter space, displaying more than an order of magnitude higher L_[C II] for a given Ly-α EW (or vice versa). No Ly-α EW measurements are published for the OpQ and IRQ samples of Decarli et al. (2018) and Nguyen et al. (2020). W0220+0137 and W0831–0140 are the two most extreme Hot DOGs, where L_[C II] shows an excess of almost three orders of magnitude. Interestingly, Fig. 1 shows that these two Hot DOGs, in addition to a narrow Ly-α component with FWHM ∼ 500−1000 km s⁻¹, also present an extremely broad component with FWHM ∼ 7500−10 000 km s⁻¹ underlying the narrow component. W0116–0505 also displays a narrow+broad Ly-α profile, but unfortunately the [C II] line was likely red-shifted out of the SPW for this Hot DOG. W0236+0528, in which [C II] is tentatively detected, only displays a broad Ly-α component, while W2246–0526 and W0134–2922 only show narrow emission, with weak or no broad component.

Fig. 4. [C II] luminosity as a function of the rest-frame equivalent width (EW) of Ly-α. Hot DOGs are displayed as red squares. Normal star-forming galaxies at z ∼ 4−6 are shown as pink diamonds (Le Fèvre et al. 2020). The sample of Ly-α emitters (LAEs) at z ≃ 5−7 from Harikane et al. (2018) are displayed as gray circles, together with the fit to the correlation they found using the Ly-α EW corrected for inter-galactic medium absorption (their Eq. (25)). While we do not attempt such a correction in our sample of luminous Hot DOGs, the average increase in the EW of LAEs due to this effect is less than a factor of two. Hot DOGs deviate from the correlation by more than one order of magnitude, toward larger L[C II] for a given Ly-α (or vice versa). The arrows indicate a proportionally equal increase in the luminosity of both emission lines by a factor of 10 and 500 from the correlation toward W0134–2922 and W2246–052, and toward W0220+0137 and W0831–0140, respectively.

Figure 4 also shows that the star-forming galaxies of the ALPINE survey lie in the region of the parameter space between the LAEs and the Hot DOGs. This is interesting because of the redshift of each galaxy population. The LAE population is at z ≃ 5−7, the ALPINE survey is at z ≃ 4−6, and our luminous Hot DOG sample lies at z ≃ 3−4.6 (with only W2246–0526 being at z > 4). These differences could be related to the average escape fraction of Ly-α or ionizing photons (Harikane et al. 2018). Unfortunately there is no information regarding this quantity for the ALPINE galaxies or our Hot DOG sample and therefore no comparisons can be made among the samples.

4.3. Dynamical support

The |V_max|/σ_V ratio, where |V_max| is the maximum projected rotational velocity of the galaxy and σ_V its overall velocity dispersion, is a first order measurement that indicates whether a system can be gravitationally supported by rotation, when |V_max|/σ_V > 1, or by dispersion or turbulence, when |V_max|/σ_V < 1.

To calculate the quantities |V_max| and σ_V, we perform rigorous kinematic modeling of the two sources for which the angular resolution and quality of the data allow such a detailed characterization: W0220+0137 and W2246–0526. A description of the modeling, using the tool ^3DBarolo, together with the obtained results can be found in Appendix A. The |V_max| (=V_rot) and σ_V (=FWHM_mod/2.355) values obtained from this analysis are reported in Table 3. We also attempted to model the kinematics of W0134–2922 and W0831+0140, but the analysis resulted in an unconstrained morphological inclination and position angle for the former, and in position-velocity diagrams of the minor and major axes indicating clear non-rotation for the latter. Thus, for these sources we use |V_max| and σ_V values derived from visual inspection of the velocity field maps (Fig. 9) and from the FWHM measured directly in the spectra (Fig. 2), respectively. Specifically, |V_max| is calculated as (|V₊| + |V₋|)/2, where V₊ and V₋ are the maximum positive and negative projected mean velocities per resolution element identified in the moment 1 line maps, evaluated at the S/N = 3−4 contours; σ_V is calculated as FWHM_[CII]/2.355. These values are also tabulated in Table 3.

Considering this, we report |V_max|/σ_V ratios for each of the rings used in the kinematic modeling of W0220+0137 and W2246–0526. For W0220+0137, the set of {d_r, |V_max|/σ_V} values are d_r = [0.38, 1.14, 1.90] kpc and |V_max|/σ_V = [1.15 ± 0.24, 1.90 ± 0.35, 2.11± 0.71], where d_r is the distance of the ring to the kinematic center of the system. For W2246–0526, the set of values are d_r = [0.46, 1.37] kpc and |V_max|/σ_V = [0.88 ± 0.14, 0.71 ± 0.13]. The measured |V_max|/σ_V ratios for W0134–2922 and W0831+0140 are 0.30 ± 0.42 and 0.48 ± 0.27, respectively. We warn that since kinematic modeling was not possible for these two sources given the current data, beam smearing can be significant. As an example, while the FWHMs derived from the modeling of W0220+0137 and W2246–0526 are similar to those obtained directly from the spectra, the projected |V_max| measured in the moment 1 maps seem to be systematically underestimated with respect to the modeled rotation velocities by a factor of ∼2. Therefore, the |V_max|/σ_V ratios for W0134–2922 and W0831+0140 are likely to be at least ≳0.6.

Tan et al. (2019) observed 8 nearby (z < 0.19) IR QSOs with ALMA in the CO(1→0) emission line, a tracer of the cold molecular gas component in galaxies. They modeled the kinematics of the four isolated QSO host galaxies in the sample displaying velocity gradients and show that three of these sources are rotation dominated, with |V_max|/σ_V = 4−6, whereas the other one shows evidence for a disturbed disk with turbulent molecular gas. The remaining four QSOs display complex kinematics, three of them likely due to ongoing interactions with companion galaxies (see discussion in Sect. 7). The |V_max|/σ_V ratios found for the two Hot DOGs for which a kinematic modeling was possible are significantly lower, below ∼2. This range is also at the lower end of the distribution of values measured in star-forming galaxies at z ∼ 2 (Förster Schreiber et al. 2006, 2009), Lyman Break Galaxies (LBGs) at z ∼ 3−3.5 (Gnerucci et al. 2011), Hα-selected galaxies from z ∼ 0.4 to z ∼ 3.3 (Wisnioski et al. 2015; Gillman et al. 2019), and H-band selected galaxies in the CANDELS field at z ∼ 1.5−3.2 (Price et al. 2019). All these samples have average |V_max|/σ_V ratios in the range of ≃1−3, and, in particular, Gillman et al. (2019) and Price et al. (2019) showed that the ratio decreases as a function of redshift, which would agree with the relatively low values found for our Hot DOG sample.

4.4. Global ISM conditions

4.4.1. The [C II] deficit

By virtue of originating from multiple ISM phases (ionized and neutral gas), the [C II] line is the current workhorse used for the study of the gas content and kinematics of galaxies on or above the MS at high redshift. In combination with the IR luminosity, the so-called [C II] deficit is one of the most widely used far-IR diagnostics with which to infer the state of a galaxy’s ISM (e.g., Malhotra et al. 2001; Brauher et al. 2008; Stacey et al. 2010). Defined as the ratio of [C II]-to-IR or to-far-IR luminosity, low values of these ratios (large deficits) are a characteristic signature of nearby LIRGs and ULIRGs (Graciá-Carpio et al. 2011; Díaz-Santos et al. 2013; Farrah et al. 2013). While there is a clear trend for more IR luminous sources to display progressively lower L_[C II]/L_IR ratios (≲10⁻³) with respect to normal star-forming galaxies (≃5 × 10⁻³), the most clear and tight correlations of the [C II] deficit relate to the dust temperature and the IR luminosity surface density, Σ_IR (Díaz-Santos et al. 2013, 2017; Lutz et al. 2016; Smith et al. 2017), as well as the star formation efficiency (L_IR/M_H2; Graciá-Carpio et al. 2011).

The top left and right panels in Fig. 5 show the L_[C II]/L_IR ratio as a function of Σ_IR and L_IR, respectively, for our sample of Hot DOGs as well as for the comparison populations of high-z QSOs and the GOALS sample. In the top-left panel, high-redshift QSOs seem to follow the tight trend displayed by nearby (U)LIRGs, in which higher Σ_IR sources display lower L_[C II]/L_IR ratios regardless of their physical size (see color coding of the symbols) –a trend whose extrapolation to the highest Σ_IR intercepts the location where the Hot DOGs lie (see discussion in Sect. 7.1). Interestingly, while most OpQs and IRQs (star and diamonds symbols, respectively) seem to follow the local correlation within ∼2σ of the dispersion, they also tend to be systematically above it.

Fig. 5. Correlations between [C II] and IR luminosities, and far-IR sizes for the Hot DOG sample and the comparison high-z quasar and star-forming galaxy populations. Top-left: [C II] deficit (here defined as L[C II]/LIR) as a function of the IR luminosity surface density, ΣIR, color-coded as a function of the intrinsic 160 μm radius. Top-right: [C II] deficit as a function of LIR, color-coded as a function of ΣIR. Bottom-left: [C II] EW as a function of ΣIR, color-coded as a function of L[C II]. Bottom-right: [C II] surface density, Σ[CII], as a function of the intrinsic ∼160 μm radius (see Sect. 2.3), color-coded as a function of the IR luminosity. Hot DOGs are shown as squares, along with their names. The symbols of the comparison quasar samples are as in Fig. 3 (filled stars and diamonds). In addition, galaxy companions found in the neighborhood of some quasars in the sample of Decarli et al. (2018) are shown as open stars (from Neeleman et al. 2019). The local, z < 0.1 (U)LIRG sample from GOALS (Armus et al. 2009) is shown as small circles and limits (Díaz-Santos et al. 2017). Black crosses indicate local galaxies in which an AGN contributes ≥50% to the bolometric luminosity. In general, black symbols denote galaxies that do not have information regarding the property shown in the color bar. In all panels, the solid line is a fit to the local (U)LIRGs (the dotted lines show the 1σ standard deviation around the fit; see Díaz-Santos et al. 2017 for details), and the dotted-dashed line is the extrapolation toward lower and higher values of the x-axes. The dashed line in the top-left and bottom-right panels represents a constant Σ[CII] ≃ 2 × 109 L⊙ kpc−2, that seems typical of the Hot DOG sample, independently of their LIR. The two arrows in the top-left panel refer to the discussion in Sect. 7.1. We note that in this figure we are using LIR to show the [C II] deficits and luminosity surface densities, rather than the more traditional approach of using LFIR that is commonly found in the literature.

4.4.2. [C II] equivalent widths, Σ_[C II], and core sizes

The bottom-right panel of Fig. 5 shows that the dust continuum emission of Hot DOGs is as compact as in z ≃ 4.7 IRQs but neither population reach sizes as small as those displayed by some z ≳ 6 OpQs. The most extreme of these cases show physical sizes (FWHM = 2 × R_eff) as small as ∼0.5 kpc, more compact than the majority of local (U)LIRGs. It is interesting to note that, despite the wide range in dust continuum (and [C II]) physical sizes, our luminous Hot DOGs seem to exhibit consistently large concentrations of [C II] at any R_eff, with L_[C II] surface densities as large as Σ_[CII] = 2 × 10⁹ L_⊙ kpc⁻² (or ≃1 × 10⁹ L_⊙ kpc⁻² if [C II] sizes are used; see below). Moreover, the entire sample of Hot DOGs seem to be within ∼50% of that value. Some high-z OpQs and IRQs reach the same level of L_[C II] density, but the lower end of the distribution extends down to much smaller values, Σ_[CII] ≃ 10⁸ L_⊙ kpc⁻².

In Hot DOGs with secure [C II] detections, the size of the line emitting core (see Sect. 3.2) is clearly resolved in one case, W0831–0140, and marginally resolved in the rest (FWHM_int ≃ FHWM_beam; see Table 3). The circularized, intrinsic [C II] FWHMs are distributed around a narrow range of physical sizes, 1.8−2.4 kpc, with an average close to ≃2 kpc. In terms of the continuum emission, the sizes range from unresolved (W0116–0505), up through barely resolved (W0134–2922, W0220+0137 and W2246–0526) and up to securely resolved (W0236+0528, W0410–0913 and W0831+0140). The circularized dust sizes span a range between 1.1−2.6 kpc, with a mean of 1.6 kpc. In all cases, when both the line and continuum are detected, the size of the [C II] emission is larger than that of the continuum. Venemans et al. (2020) report [C II] and dust continuum sizes for the sample of 27 OpQs at z ≳ 6 introduced in Decarli et al. (2018). They find ranges between ≃0.7−6 kpc and 0.5−7 kpc, respectively. We note that the ratio of line-to-continuum sizes in Hot DOGs is very narrow, ≃1.35−1.85, with a mean of 1.61 ± 0.10, and remarkably similar to the averages found for the sample of z ≳ 6 OpQs, 1.68 ± 0.11, and for the z ≃ 4.7 IRQs, 1.56 ± 0.08. This suggests that the line and continuum sizes in high-redshift QSOs, regardless of their selection, could be linearly proportional, and may be pointing to a common physical origin for the size difference. However, to derive any statistically robust conclusion higher angular resolution observations are needed.

4.5. Spatial distribution and extended emission

As described in Sect. 3, the intrinsic sizes derived from the 2D Gaussian fits are only representative of the core emission of the galaxy and, as can be seen in Figs. 6 and 7, some sources that are barely resolved in line or continuum emission based on this measurement, still show faint structures with clumpy morphologies extended toward particular directions. The most luminous of these clumps are galaxy neighbors located at distances 1−5″ from the central Hot DOG host (e.g., W2246–0526; Díaz-Santos et al. 2016). A blind search for companion galaxies in the ALMA [C II] cubes will be presented in a separate work by González-López et al. (in prep.). The most extended Hot DOG in the sample is W0410–0913, which displays low surface brightness continuum emission extending up to 1−1.5″ north-east and south-west from the central quasar, equivalent to ≃10 kpc (projected). Unfortunately, the [C II] line was mostly missed in this object due to the redshift offset with respect to its z_UV, as described in Sect. 4.1. While fainter, W0236+0528 and W0831+0140 show a number of continuum emission knots surrounding the central galaxy. These clumps have flux density peaks with S/N ≃ 3 and are surrounded by large areas of overall positive emission, tentatively suggesting that they may be real.

Fig. 6. Dust continuum flux density maps of the Hot DOG sample. The beam is shown as a white ellipse on the bottom-left, together with its size and the depth (rms) of the collapsed, line-subtracted cube. The dashed white circumference denotes the aperture used to extract the spectra shown in Fig. 2. The dotted contour shows the zero-flux intensity, and the solid contours represent the [2.5, 3.0, 2(4 + n)/2]-σ levels, where σ is the rms of the map and n = [0, 1, 2…]. Symmetric, negative levels are shown as dashed contours. Field of views (FoVs) of 2″ and 4″ on a side are displayed depending on the galaxy. Roughly, 1″ ≃ 7 kpc.

Fig. 7. [C II] intensity (moment 0) maps of the Hot DOGs in which the line has been detected. The beam is shown as a white ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are displayed as in Fig. 6. The FoVs are also the same as in Fig. 6 to facilitate the comparison.

In addition to the sizes, Tables 3 and 4 present the fraction of extended emission (f_ext; line and continuum) underlying the central, unresolved point-source. In all Hot DOGs, the f_ext, [CII] is ≳50%, and up to ∼95% in the case of W0831+0140, indicating that most of the ionized gas is extended over several kpc and is likely associated with star formation in the host galaxy (see discussion in Sect. 7). On the other hand, the dust continuum shows f_ext, cont ranging from very compact (∼20%, W0116–0505) to very extended sources (∼80%, W0831+0140). These values are always smaller than f_ext, [CII], implying that the continuum emission is arising from a volume closer to the central AGN than that responsible for [C II]. Because we are measuring the rest-frame continuum at ∼160 μm, this implies that either: (a) the AGN could be contributing down to the coldest components of the distribution of dust temperatures in the galaxy (∼20 K), something rarely seen in the local Universe (e.g., not even in ultra-luminous quasars like Mrk 231), or (b) the temperature distribution of the dust heated by the AGN is typical, with the emission falling rapidly at λ ≳ 40 μm (e.g., Elvis et al. 1994; Nenkova et al. 2008; Mullaney et al. 2011), but its total luminosity completely overwhelms that of the underlying galaxy, and thus still contributes significantly to the dust continuum at longer wavelengths. We note that the physical resolution of the current ALMA data is ∼1.5−2 kpc, and therefore much higher resolution observations are needed to further isolate the gas and dust structures located under the sphere of influence of the SMBH; roughly within a radius of a few hundred pc for the estimated masses of these BHs (Wu et al. 2018; Tsai et al. 2018).

Figure 8 shows the radial density profiles of the continuum and line emissions, for the most compact and extended sources in the Hot DOG sample (see Sect. 3 for details). The continuum emission of W0116–0505 is very compact. Although the core of the galaxy is effectively unresolved (Table 4), about 20% of the total flux still arises from a very extended area of low surface brightness emission. W0831+0140 shows both the most extended continuum emission and the most extended line emissions, with f_ext ≳ 80% in both cases. However, we note that W0410–0913 has the largest continuum core, which suggests that W0831+0140 shows the largest f_ext, cont because the beam size of these observations is half that of W0410–0913. Indeed, if we use a uv-tapering of 0.5″ to place all the cubes on a common beam size, the f_ext, cont of W0410–0913 becomes the largest of all the Hot DOGs (≃50%) and that of W0831+0140 is reduced to ≃30%.

Fig. 8. Top two panels: Azimuthal profiles of the dust continuum flux density for the Hot DOGs with the smallest (top, W0116–0505) and largest (bottom, W0831+0140) fext, cont. Bottom two panels: Azimuthal profiles of the [C II] flux for the Hot DOGs with the smallest (top, W2246–0526) and largest (bottom, W0831+0140) fext, [CII]. The black line is the profile of the observed emission. The dashed line is the profile of the clean PSF. The dotted-dashed line is the scaled PSF that minimizes the scatter of the residuals within Bmaj/2 = 1 (see Sect. 3 for details and values in Tables 3 and 4). The dotted line is the profile of the residuals (observed – scaled-PSF). The PSF scaling value, fext and the radius used to perform the aperture photometry are noted at the top right of each panel (see dotted circumferences in Fig. 6).

5. Kinematics

5.1. Velocity fields and dispersion maps

Figure 9 presents the moment 1 maps of the Hot DOG sample for the sources in which the [C II] line has been clearly identified. As can be seen, the kinematics are very diverse: W0134–2922 shows a very small velocity gradient; W0220+0137 displays a smoothly varying structure more consistent with a rotating disk; W0831+0140 shows a complex velocity distribution; and W2246–0526 shows relatively slow rotation, with Δv ≃ 150 km s⁻¹ (cf., Díaz-Santos et al. 2016). While W0410–0913 is not shown, the detected, blue-shifted component of the line reaches velocities up to ∼600 km s⁻¹ (assuming the [C II] systemic velocity is the same as that derived from z_CO).

Fig. 9. [C II] projected radial velocity (moment 1) maps down to a 3σ level (see contours in Fig. 7). The beam is shown as a black ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are drawn at [0, ±50, ±100, ±200, ±300…] km s−1. The FoVs are the same as in Fig. 6 to facilitate the comparison.

The case of W0831+0140 is of particular interest. In this source, the complex velocity field is likely caused by the presence of multiple narrow-velocity components that can be seen in the spectral cube, suggesting that this source could actually be a system of clumps or galaxies in the process of merging. However, higher angular resolution and S/N observations are needed to accurately disentangle these components spatially and in velocity space. In fact, the sources are in general either not well resolved (W0134–2922) or do not have velocity fields smooth enough to perform rigorous kinematic modeling (W0831+0140). The exceptions are W0220+0137 and to a more limited extent, W2246–0526. We show the results from their kinematic modeling, using the tool ^3DBarolo, in Appendix A.

The moment 2 maps of the Hot DOG sample (see Fig. 10) show relatively uniform values of the velocity dispersion across the area where the [C II] line is detected. Moreover, the measured peak and integrated FWHM values are always in excess of 500 km s⁻¹, thus suggesting a highly turbulent gas component across the host galaxies. Histograms of the integrated FWHM_[CII] for our sample of EL Hot DOGs and for the comparison quasar samples are presented in Fig. 11. For reference, the distribution of MS galaxies at z ∼ 4−6 is shown as well. Quasars have on average higher [C II] velocity dispersions than star-forming galaxies, and in turn Hot DOGs have systematically larger dispersions than OpQs and IRQs. A two-sample Kolmogorov–Smirnov test performed over the Hot DOG and the combined quasar samples yields a value of the statistic D = 0.95 with a significance of p ≃ 1 × 10⁻⁴, implying that the two distributions are very likely to be drawn from different parent populations. Individual tests comparing the Hot DOGs with each OpQ and IRQ sample also provide D ≥ 0.9 and p ≤ 3 × 10⁻⁴. We explore the implications of this result in the discussion in Sect. 7.

Fig. 10. [C II] FWHM (moment 2) maps down to a 3σ level (see contours in Fig. 7). The beam is shown as a black ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are drawn at [0, 25, 50, 100, 200, 300…] km s−1. The FoVs are the same as in Fig. 6 to facilitate the comparison.

Fig. 11. Histogram of the [C II] FWHM for the sample of Hot DOGs (red line). The z ≃ 4.7 IRQs from Trakhtenbrot et al. (2017) and Nguyen et al. (2020) are shown in yellow. The compilation of z ≳ 6 OpQs from Decarli et al. (2018) is displayed as a green/blue dashed line. The 4 ≲ z ≲ 6 star-forming galaxies from Le Fèvre et al. (2020) are shown in pink.

5.2. Dynamical masses from modeling

The kinematic analysis described in Appendix A allows for a proper derivation of the dynamical masses of W0220+0137 and W2246–0526 host galaxies (see Eqs. (A.1) and (A.2)), which we can compare to the stellar masses estimated in Sect. 6 and to BH masses published in the literature (see Table 5).

The best-fit kinematic model for W0220+0137 yields a M_dyn = 2.8 × 10¹⁰ M_⊙ regardless of whether the galaxy is assumed to be rotation-supported (i.e., where most of the mass is a thin disk with ordered rotation) or dispersion-supported. In turn, W2246–0526 results in M_dyn = 5.6 × 10⁹ M_⊙ assuming the former, and M_dyn = 5.5 × 10¹⁰ M_⊙ considering the latter. Since the modeling of W2246–0526 strongly favors a dispersion-supported system, we refer from here on to the value derived from this assumption as the most likely M_dyn. This value is within the uncertainties of the M_dyn obtained by Díaz-Santos et al. (2018) based on simple considerations of the [C II] line width and the size of the emitting region, M_dyn = 8 ± 4 × 10¹⁰ M_⊙.

In Sect. 6 we derive the stellar masses of the Hot DOG hosts from a variety of assumed SFR histories. The values obtained for W0220+0137 and W2246–0526 are M_⋆= 1.2 ± 0.6 × 10¹¹ M_⊙, and M_⋆= 3.0 ± 1.6 × 10¹¹ M_⊙, respectively. Notwithstanding that the values of M_⋆ and M_dyn still agree at a ≃3σ level (given the large uncertainties in both measurements), the stellar mass of the Hot DOG hosts appears to exceed their dynamical mass by a factor of ∼5, suggesting M_⋆ might be significantly over-estimated. Indeed, it should be noted that the stellar masses are based on photometric data (for details, see the supplementary material in Díaz-Santos et al. 2018) that have much lower angular resolution than the ALMA observations used to derive the dynamical masses. Specifically, while the observed angular sizes measured from the [C II] line for both W0220+0137 and W2246–0526 are ≃0.3−0.4″, the Spitzer and WISE broad-band photometry have an angular resolutions of ≃2−6″; around an order of magnitude larger in projected area than the ALMA data. Thus, the near-IR photometry could include emission from companion galaxies and extended structures not captured by the ALMA-derived M_dyn. This caveat is not only important for Hot DOGs but also when comparing stellar and dynamical masses of any quasar or galaxy population at high redshift. Sub-arcsecond, rest-frame near-IR observations with JWST will be needed to estimate the stellar masses of Hot DOG’s hosts that allow for a direct comparison with the dynamical masses obtained from ALMA at similar physical scales.

Tsai et al. (2018) presented a rest-UV spectrum observed with Keck/OSIRIS of W2246–0526 and used the Mg II 2799 Å emission line to obtain an estimate of the BH mass. They use the calibration from Wang et al. (2009) to arrive to a M_• = 4 × 10⁹ M_⊙ and an Eddington ratio (λ_Edd = L_bol/L_Edd ∝ L_bol/M_•) of λ_Edd ≃ 2.8. Similarly, Finnerty et al. (2020) used rest-frame optical spectroscopy obtained with Keck/NIRES to derive an upper limit to the mass of the SMBH in W0220+0137. They assumed the M_• − σ_⋆ relation derived by Kormendy & Ho (2013), where σ_⋆ is determined from the narrowest detected emission feature in the spectrum, to obtain M_• ≤ 4 × 10¹⁰ M_⊙ and from it a λ_Edd ≥ 0.06.

Figure 12 shows the SMBH mass as a function of the system’s dynamical mass for W0220+0137 and W2246–0526 compared to other OpQs and IRQs in the literature, in the z ∼ 2−6 range (see Bischetti et al. 2021). The black hole mass of W2246–0526 is well within the dynamical mass of the system, with M_• representing a fractional ≃0.07 of M_dyn. Unfortunately, the black hole mass of W0220+0137 only allows us to place a rough upper limit to its contribution to the dynamical mass, which could be up to a 100%. However, if we were to assume the Eddington ratio of W0220+0137 to be similar to that of W2246–0526 (i.e., slightly super-Eddington), the contribution of M_• to M_dyn would decrease to a fraction of ≃0.03 (see long red arrow), more similar to the value obtained for the latter as well. Considering this, these two luminous Hot DOGs seem to host BHs that are at least an order of magnitude more massive than nearby elliptical and spheroid galaxies with comparable dynamical masses, in agreement with Assef et al. (2015). However, their location in the M_• − M_dyn parameter space is not significantly different from other quasar populations at similar, lower, or higher redshifts.

Fig. 12. SMBH mass, MBH, as a function of dynamical mass, Mdyn, for the Hot DOGs W0220+0137 and W2246–0526. Literature OpQ and IRQ samples, from z ∼ 2 to z ∼ 6, are displayed for comparison (Bischetti et al. 2021, and references therein). The black solid line is the MBH versus dynamical bulge mass relation from the nearby galaxy sample of Häring & Rix (2004), and the dashed represents the 1σ intrinsic scatter. The dotted line is the MBH/Mdyn = 1 ratio. The SMBH mass of W0220+0137 is an upper limit (Finnerty et al. 2020), and the large red arrow pointing downwards represents the change in MBH if a λEdd = 2.8 is assumed, equal to that derived for W2246–0526 (see text for details).

6. Hot DOGs on the main sequence

In order to estimate the stellar masses for the Hot DOG sample, we perform a spectral energy distribution (SED) fitting to the UV-through-mid-IR photometry available for each source following Assef et al. (2010, 2015). A complete description of the methodology can be found in the Supplementary Materials section of Díaz-Santos et al. (2018), which describes the example of W2246–0526. For additional clarity, we also re-describe the methodology in Appendix B, where we also present the photometry used for the fitting of each source, as well as the best combination of model templates that best reproduce their SED (see Fig. B.1). To obtain the SFR in the host, we derive SFR_IR estimates based on both the far-IR rest-frame dust continuum emission and the [C II] luminosity of each galaxy when available. Both observables are corrected by multiplying by the fractions of extended emission (f_ext, see Tables 3 and 4) to account only for the resolved component of the flux and thus avoid any potential AGN contamination. For the former, continuum-based method, we use two SED templates: one is the starburst galaxy M 82 and the other is a typical Sd spiral galaxy, which we scale to the rest-frame 158 μm dust continuum to obtain the SFR. For the latter, we use L_[C II] and two representative values of the [C II]-deficit seen in local, non AGN-dominated (U)LIRGs (L_[C II]/L_IR = [5,0.2] × 10⁻³, see Fig. 5) to estimate the star-forming L_IR, which in turn is converted to a SFR using the calibration from Murphy et al. (2011): SFR/L_IR = 1.48 × 10⁻¹⁰ [M_⊙ yr⁻¹]/[L_⊙], for a Kroupa (2001) IMF. Table 5 presents the average M_⋆ and SFRs for the Hot DOG sample.

Figure 13 shows the position of the Hot DOG hosts on the SFR–M_⋆ plane, together with the star-formation MS (from Schreiber et al. 2015) at the specific redshift of each galaxy (solid lines, color-coded as a function of redshift). The individual estimates of M_⋆ and SFR, based on the assumptions described above, are also shown (as small open squares). We note that the error bars around the average values (displayed as big solid squares) represent the full extent of the range of estimates and not the standard deviation. In addition, we also show the median values as open circles. As can be seen, EL Hot DOGs are spread over a large region of the parameter space. Three of the four z ∼ 3 Hot DOGs (purple/magenta squares) and W2246–0526 are consistent with being on the MS. W0831+0140 and especially W0410–0913, at z ∼ 3.9 and z ∼ 3.6, respectively, seem to be well above it. Interestingly, these two Hot DOGs have the largest dust continuum emission core sizes and show among the largest f_ext, which would agree with potential wide-spread star formation in the hosts. W0116–0505, on the other hand, is well below the MS.

Fig. 13. SFR as a function of M⋆ for the host galaxies of the Hot DOG sample, color-coded as a function of redshift. The galaxy star-formation MS at the specific redshift of each source is represented by solid lines (from Schreiber et al. 2015) also colored to match their redshifts. Mean values of SFR and M⋆ are shown as solid squares; median values are shown as open circles. Small, open squares represent the individual M⋆ and SFR estimates based on the SED modeling and the assumed SFR histories and scaling relations (see text and Appendix B). Note that the error bars associated to the mean values do not represent standard deviations but rather extend the full range of individual estimates.

Given the small number of sources and the large uncertainties inherited from the wide range of M_⋆ and SFR values obtained from the SED modeling and far-IR scaling relations, respectively, it is difficult to provide any detailed interpretation regarding the overall population of EL Hot DOGs or any particular object. Perhaps the only inference that can be made is that, when considered as a whole population, EL Hot DOGs do not all seem to be lying systematically below or above the star-formation MS. This is rather surprising considering the active accretion their SMBHs are experiencing, which should be accompanied by active star formation in the host, thus elevating the galaxies significantly above the MS if both the host and the SMBH were roughly co-evolving, as expected (Madau & Dickinson 2014). However, when considering sources individually, it is possible that they can deviate from the MS significantly, above or below, and therefore the average position of the population on the SFR–M_⋆ plane would not have much meaning. Riguccini et al. (2015, 2019) found a similar result for a sample of DOGs at z ∼ 1−3, which display a large scatter around the MS with no specific trend based on the galaxies being dominated by star formation, or having a significant contribution from a moderately obscured or Compton-thick AGN.

A larger number of Hot DOGs and a better sampling of their photometric SEDs are necessary to obtain more precise stellar masses and SFRs, and thus be able to investigate whether individual cases such as those of W0410–0913, W0831+0140 or W0116–0505 are significant.

7. Discussion

7.1. Disentangling the power-source(s) of IR emission in extremely luminous Hot DOGs

In Sect. 4.4 we described where our sample of Hot DOGs lie in a number of diagrams relating the luminosity of the [C II] line to that of the underlying dust continuum and its physical size (Fig. 5). In this section we discuss how the presence of an extremely luminous (EL) AGN (L_bol ≥ 10¹⁴ L_⊙) may be able to alter the regular ISM signatures of the host galaxy, and provide a possible interpretation for the position of the most luminous Hot DOGs in the parameter space delimited by the ALMA data.

The dust continuum emission signature of AGNs and, in particular, that of Hot DOGs, is most noticeable at rest-frame mid-IR wavelengths, where the contribution of the highest-temperature components of the dusty structure reprocessing most of the SMBH’s power dominates the SED (Wu et al. 2012; Tsai et al. 2015). However, as the AGN becomes extremely luminous and starts to bolometrically account for nearly the entire power output of a galaxy, it can also emerge as the dominant source even at λ ≳ 100 μm. This is all the more likely if the coldest component of the temperature distribution of the obscuring medium is still warmer than the overall dust temperature of the underlying galaxy (see, e.g., Wu et al. 2012, 2014), or if a significant amount of the AGN dust is optically thick in the far-IR. Were this to happen, the measured size of the dust continuum would progressively appear more compact, as it would be representative of the more dominant central point-source. These two effects (emission and size) would also likely come in hand with the fact that, assuming the powerful far-IR AGN is also a significant X-ray emitter (L_X ≳ 10¹² L_⊙), its hard ionizing spectrum could reduce the C⁺ abundance in the warm ionized medium (WIM) in the vicinity of the SMBH (≲few hundred pc), converting some fraction of C⁺ to C²⁺ and higher ionization states (Langer & Pineda 2015). This would have the effect of suppressing the surrounding [C II] luminosity up to almost two orders of magnitude in the most extreme cases, generating a cavity of higher ionization gas. All these effects combined would strongly drive galaxies that are bolometrically dominated by an AGN to have high Σ_IR and very low [C II] EWs, i.e., toward the lower-right region of the parameter space in the bottom-left panel of Fig. 5, significantly below the region where EL Hot DOGs actually lie, and below the extrapolation from local (U)LIRGs.

Indeed, there are a number of observables suggesting that some of these effects may actually not be as important as expected. The bottom-right panel in Fig. 5 shows the luminosity surface brightness of [C II], Σ_[CII], as a function of the intrinsic radius of the dust continuum emission. EL Hot DOGs do not seem to have particularly small far-IR sizes when compared to the most extreme, nearby systems. This implies that even if the unresolved AGN in Hot DOGs is two orders of magnitude more luminous than that in local (U)LIRGs, its relative contribution to the far-IR continuum, which would make the size smaller, is somehow compensated by a proportional increase of emission by the resolved host (or underlying extended gas structure in general), in order to keep the system physically resolved. This is in agreement with the extended component of the Hot DOGs accounting for a significant fraction of the dust continuum emission (from f_ext ∼ 20% up to ∼70−80% in the most extreme cases; see Table 4 and Fig. 8). The [C II] emission also has a substantial extended component fraction, larger than that of the continuum in all cases.

Instead of the size, the main difference with respect to the local, IR-luminous galaxy population is actually Σ_[CII], which in EL Hot DOGs is consistently around an order of magnitude higher than local (U)LIRGs for their given size range, and only matched by some of the most compact z ≳ 6 OpQs and a few z ≃ 4.7 IRQs (see bottom-right panel of Fig. 5). This excess in [C II] luminosity per area also seems to argue against the scenario in which an energetic AGN would deplete the surrounding abundance of C⁺ by pumping higher ionization states and thus suppressing the [C II] luminosity. Or alternatively, it could imply that Hot DOGs are not be powerful X-ray sources (Stern et al. 2014; Stern 2015).

The location of the EL Hot DOGs in the two bottom panels of Fig. 5 is then likely set by a combination of an IR and [C II] luminosity excess with respect to the extrapolation of local (U)LIRGs, and not (or to a much lesser degree) to a difference in galaxy size. While the Σ_IR excess can be easily explained by the contribution of the dominant AGN to the bolometric luminosity of the galaxy, the origin of the Σ_[CII] excess is less evident, and can be interpreted in at least two (non exclusive) ways:

First, at higher redshifts, the gas content of galaxies progressively accounts for a larger fraction of their baryonic mass (Scoville et al. 2017; Tacconi et al. 2018; Walter et al. 2020). If the total gas reservoir of galaxies, M_gas, roughly scales with L_[C II] (Zanella et al. 2018), then EL Hot DOGs should have an order of magnitude more gas mass than local (U)LIRGs of similar size (Fig. 5, bottom-right panel). In addition, in order to explain the narrow range of relatively low [C II] EWs (bottom-left panel), the luminosity of the continuum under the line needs to increase roughly proportionally, which would be consistent with the dust emission being also a tracer of the gas content in the host (or at least of the ionized gas probed by [C II]). In this case, the AGN has no impact on the properties of the Hot DOG other than boosting the total L_IR of the galaxy (at shorter than far-IR wavelengths). In other words, if we were to assume that the typical EL Hot DOG host is a galaxy similar to a maximal, star-forming local (U)LIRG (Σ_IR ≃ 5 × 10¹¹ L_⊙ kpc⁻² and R_eff ≃ 0.5−1 kpc; see Fig. 5, top-left panel) but containing an order of magnitude more gas mass (and thus a proportionally higher L_[C II]/L_IR ratio), when corrected by the two orders of magnitude increase in total luminosity due to the AGN (i.e., along the arrows in the same figure), the resulting galaxy would be put back on the extrapolation of the local IR population (dotted-dashed line). However, the actual ISM content and relative contribution of the energy sources (SFR vs. AGN) in the resulting Hot DOG may be very different from a galaxy that is based on the direct extrapolation of a purely star-forming ULIRG toward higher Σ_IR following the local correlation.

Secondly, an alternative possibility is the existence of an additional heating source, other than star formation, that enhances the [C II] emission. Evidence of excess C⁺ with respect to the maximum line cooling efficiency expected from pure photo-dissociation regions has been seen in purely star-forming galaxies at z ≃ 1−2 (Brisbin et al. 2015). In the local Universe, Seyfert galaxies hosting powerful AGN-driven radio jets and winds inject turbulence into the surrounding ISM, contributing to the excitation of a larger mass of the C⁺ reservoir (Guillard et al. 2015; Appleton et al. 2018; Smirnova-Pinchukova et al. 2019) (in contrast with the depletion scenario based on ionization described above, from Langer & Pineda 2015). In addition, low-velocity shocks driven by galaxy collisions can also turbulently heat the gas on large scales in dense environments (Appleton et al. 2013, 2017). Given that Hot DOG hosts do not seem to be forming stars at a much higher rate than normal MS galaxies at the same redshifts (see Fig. 13 and Sect. 6), significant turbulent heating could be provided by the energy input from the central SMBH (Díaz-Santos et al. 2016) or by the dynamical friction consequence of the accretion of neighbor galaxies (see next section). Both cases would be in concordance with the large velocity dispersion of the [C II] line seen in all EL Hot DOGs, FWHM ≳ 500 km s⁻¹ (Sect. 5), and with the fact that most seem to live in over-dense, merger-driven environments (e.g., Jones et al. 2014, 2017; Assef et al. 2015; Díaz-Santos et al. 2018).

Because of the size constraint, the most likely scenario is a combination of these two possibilities. That is, the typical EL Hot DOG galaxy host could be a scaled-up version of a local LIRG (∼10¹¹⁻¹² L_⊙) –that is, with up to an order of magnitude larger total L_IR (∼10¹²⁻¹³ L_⊙), continuum emission at ∼160 μm (νL_ν160 μm ∼ 10¹²⁻¹³ L_⊙), and L_[C II] (∼10⁹⁻¹⁰ L_⊙; see Table 3). This would be in agreement with the position of Hot DOGs in the SFR–M_⋆ plane, lying on or perhaps slightly above the star-forming galaxy MS at their respective redshifts (see Fig. 13). In this scenario, the host galaxy would contribute equally to the line and dust continuum emission under it, thus maintaining the EW (Fig. 5, bottom-left). In addition, the AGN could also account for a fraction of the continuum and line flux, the first originating from the coldest dust component heated by the SMBH, and the latter potentially by its injection of turbulence into the surrounding ISM, in addition to galaxy-neighbor collisions. And while contributing to both continuum and line, the AGN would dominate neither, nor would Hot DOGs appear as mostly unresolved, but rather a combination of a central source plus an extended host (see Σ_IR profiles presented in Sect. 4.5 and Fig. 8), which is kinematically and morphologically disturbed by the interaction with companion galaxies (as some of the moment 0 maps suggest in Figs. 6 and 7). Yet, because of the larger dynamic range and higher dust equilibrium temperatures, the AGN would boost at least another order of magnitude the total IR luminosity of the galaxy, dominating the emission at shorter wavelengths, and thus accounting for ≥90% of the L_IR. Hot DOGs are very strong emitters in the rest-frame mid-IR (by selection) and therefore if this picture is correct, future observations of z ≲ 1.5 Hot DOGs with the JWST should display mostly unresolved emission at observed-frame λ ≳ 20 μm.

7.2. An evolutionary perspective

Within the current paradigm of galaxy evolution, the most luminous quasars are a relatively short phase during the merger process of two massive, gas-rich, star-forming galaxies (e.g., Sanders et al. 1988; Hopkins et al. 2008). As the galaxies coalesce, gas and dust are piled up onto the nascent common nucleus, triggering both massive star formation and boosting SMBH accretion. This dusty, obscured phase is assumed to be traced by IRQs. Once the feedback from the quasar sets in, the radiation pressure exerted on the surrounding ISM by the central AGN clears out the dusty cocoon, making the quasar detectable in the optical (OpQ). However, while this scenario can account for ULIRGs in the nearby Universe, the physical processes and environment through which quasar activity proceeds at high redshift is likely not as straightforward. At high redshift, galaxies have significantly higher gas fractions (Scoville et al. 2017; Tacconi et al. 2018), and EL quasars, likely being hosted by galaxies at the knots of the cosmic web, may be subject to intermittent, yet sustained accretion from in-falling companion galaxies (Díaz-Santos et al. 2018; Bischetti et al. 2018, 2021), and not just to a single major-merger event.

In Sect. 5 we presented the velocity and dispersion maps of the Hot DOG sample. We found that EL Hot DOGs have very diverse dynamical states: from sources that are barely rotating, to resolved hosts with ordered, circular motions, to complex, disturbed systems that are likely the result of ongoing mergers. These findings are very similar to the diversity shown by ultra-luminous IR quasars in the local Universe. In a sample of eight IRQs, Tan et al. (2019) found rotating disks, disturbed velocities, and some sources with morphologies suggesting they are undergoing a merger process (albeit these results were drawn from observations of the CO(1→0) emission line rather than [C II], thus probing a different gas phase of the ISM). However, in contrast to nearby IRQs, which generally show moderate velocity dispersion (FWHM_CO(1−0) ≃ 250 km s⁻¹), EL Hot DOGs are characterized by FWHM_[CII] ≳ 500 km s⁻¹, which in a number of cases is uniform across most of the area where the underlying host is detected. This is also reflected in the fact that the |V_max|/σ_V ratios of rotation-dominated local IRQs are in the range of ∼4−6, while in Sect. 4.3 we found that our Hot DOG sample display a lower range, with |V_max|/σ_V ≲ 2.

In terms of the connection between the general galaxy evolution framework described above and the particular context of this study, we have reported on a population of EL, hot dust-obscured quasars that do not share a consistent or characteristic velocity field, but display a common kinematic signature in the form of a high velocity dispersion of the neutral and ionized gas as traced by [C II]. While a number of non-exclusive scenarios can be invoked to explain this diversity in morphology and dynamical states, they also need to account for the short-lived nature of the Hot DOG phenomenon, which lasts on average, and overall, for less than a few tens of Myr as implied by the limited number of such objects detected in the entire sky (≲0.01 deg⁻² for L_bol ≥ 10¹⁴ L_⊙; Assef et al. 2015). Here we put forward a few potentially important aspects to consider in order to characterize and explain the rise of Hot DOGs:

The nature of the activation mechanism(s). It is possible that there is more than one mechanism triggering the EL Hot DOG phenomenon. Large-volume cosmological simulations predict that smooth accretion of pristine gas from the cosmic web could fuel the SMBH of the central galaxy in a protocluster with enough gas to make the AGN become extremely luminous, as well as to trigger significant star formation in the host (Dekel et al. 2009; Fumagalli et al. 2011). On the other hand, mergers are also able to provide an intermittent supply of material in the form of gas that has been pre-processed in neighbor galaxy companions (Decarli et al. 2018; Díaz-Santos et al. 2018; Bischetti et al. 2021). Distinguishing between these two mechanisms is complicated, as the identification and measurement of low-metallicity gas reservoirs around galaxies at high redshift is very challenging (e.g., Rafelski et al. 2012). Companion galaxies have been found in around 70% of the Hot DOGs analyzed here (González-López et al., in prep.). In comparison, the fraction of OpQs at z ≳ 4.5 with detected nearby companions is significantly lower (≲30%; Decarli et al. 2018; Trakhtenbrot et al. 2017; Nguyen et al. 2020), although recent results by Bischetti et al. (2021) suggest that around ∼80% of those quasars at the brightest end of the QSO luminosity function are surrounded by one or more bright companion galaxies. Deeper and more uniform observations of both optical quasars and Hot DOGs are needed to investigate whether this is revealing a real evolution or it is simply a consequence of the selection bias (or probably both).

Time and spatial scales. Mergers may be the dominant mechanism triggering EL Hot DOGs, but there is a time sequence connecting the different steps in the accretion process that is correlated with the physical scales over which observables enable its tracking. That is, the identification of the merger nature of a system is usually done through the detection of companion galaxies (e.g., Díaz-Santos et al. 2018; Miller et al. 2018; Hodge et al. 2019); the dynamical state of the host galaxy is probed by the kinematics of the gas (e.g., Carniani et al. 2013; Kimball et al. 2015; Tadaki et al. 2018); and the presence of an accreting SMBH is revealed by the fact that an obscured AGN dominates the bolometric output of the source (Tsai et al. 2015). Each of these signatures has a characteristic spatial scale attached to it: neighbor galaxies in the vicinity of high-z quasars are identified at projected distances greater than at least a few kpc (Petitjean et al. 1996; McGreer et al. 2014; Trakhtenbrot et al. 2017; Decarli et al. 2018, 2019; Nguyen et al. 2020); the dynamical state of the host, which is driven by the interaction between the pre-existing bound gas and the in-falling companions, is characterized by physical scales of ∼1−2 kpc that, when tied to the gas velocity dispersion observed in the host galaxies (FWHM ≳ 500 km s⁻¹), provides a time-scale of ≲a few Myr; and the funneling of gas into the SMBH happens efficiently only inside of its sphere of influence (within a few hundred pc, considering BH mass estimates of M_• ∼ 10⁹ M_⊙; Wu et al. 2018; Tsai et al. 2018; Jun et al. 2020; Finnerty et al. 2020), where it can be transported further in by rotationally supported structures (Anglés-Alcázar et al. 2021). The lack of synchronicity and delay among all these phases (including with respect to the star-formation process; Hopkins 2012) can potentially generate the diversity of dynamical states displayed by Hot DOGs despite all sharing a common origin.

Recurrence. Independently of the mechanism dominating the accretion and the time/spatial-scales describing the accretion, this process could occur once, or it may be recurrent. The former scenario is usually simulated in the local Universe as a single-merger event of two massive, gas-rich galaxies leading to a quasar phase (Hopkins et al. 2008). This is supported by the fact that the majority of nearby, IR-selected ULIRGs and IRQs display signatures suggesting they are the product of major mergers in the latest stage of the interaction, immediately previous to (or right after) the coalescence of the two galactic nuclei (Veilleux et al. 2002, 2009; Haan et al. 2011; Kim et al. 2013; Stierwalt et al. 2013). However, at high-z, galaxies may have large enough gas reservoirs stored in clumpy, unstable disks (e.g., Decarli et al. 2016; Falgarone et al. 2017; Hodge et al. 2019) such that minor mergers could potentially be sufficient to trigger a number of short-lived extremely luminous, highly dissipative phases associated with rapid SMBH and galaxy mass growth. In the particular case of EL Hot DOGs, which seem to be living in over-dense environments (Jones et al. 2014; Assef et al. 2015; Fan et al. 2017; González-López et al., in prep.), a frequent in-fall of gas in the form of whole galaxies, large clumps, or tidally disrupted material would also help to explain the observed widespread turbulence of their ISM, in addition to the energy and momentum injected by the central AGN. Indeed, kpc-scale turbulence is a clear signature of the existing gravitationally bounded gas in a host galaxy attempting to relax while still being subjected to external accretion or internal AGN feedback (Appleton et al. 2017). If that is the case, EL Hot DOGs could be characterized as a recurrent phase (rather than a single event, like nearby ULIRGs/quasars) in which a galaxy at the center of a young protocluster (where the red galaxy sequence has not been formed yet; Penney et al. 2019) is experiencing chaotic episodes of accretion, each individually lasting for a time-scale shorter than a dynamical time. This would keep the gas turbulent in the host, promoting cloud-cloud collisions and the removal of angular momentum, which would make the gas more likely to fall under the gravitational influence of the SMBH. Some of the gas could be subsequently funneled to the SMBH accretion disk, a fraction of which could be expelled via feedback, further contributing to the turbulence in the host in a self-sustaining, intermittent process as long as external accretion continues.

This is a potential scenario that would overall agree with the evidence presented in this work. It is also important to note that within this hypothesis of recurrent accretion, and considering that EL Hot DOGs seem to have Eddington ratios close to unity or above (e.g., Tsai et al. 2018; Jun et al. 2020; Finnerty et al. 2020), each individual EL episode (L_bol ≳ 10¹⁴ L_⊙) would need to be significantly shorter than the characteristic timescale for Eddington-limited BH growth (also known as Salpeter time: t_Salp = M_•/Ṁ_Edd ≃ 40 (η/0.1) Myr; where η is the conversion efficiency of the accreted mass into electromagnetic radiation, and λ_Edd is assumed to be unity). That is, with BH masses on the order of ∼10⁹ M_⊙ or slightly above (see Table 5), a large number of EL events with typical efficiencies of η = 0.1 jointly lasting for a period significantly longer than t_Salp would lead to final BH masses in excess of those seen in SMBHs hosted by the brightest cluster galaxies (BCGs) today (M_• ≃ 10¹⁰ M_⊙; McConnell & Ma 2013). Therefore, if the EL Hot DOG phase is to be recurrent, each episode would need to be shorter than a few Myr (t_e ≲ t_Salp , where N_e is the number of episodes), which would be in agreement with the observed AGN flickering time, of about 10⁵ yr (Schawinski et al. 2015), and with recent ultra-high resolution simulations of the BH accretion rate down to the central sub-pc scale in high-redshift quasars at the center of protoclusters (Anglés-Alcázar et al. 2021).

We note that the picture put forward here may not be valid for lower luminosity Hot DOGs and it is possible that not all accretion events produce an EL, hot obscured phase. To investigate if this scenario holds in general, a larger number of Hot DOGs spanning a wider range in bolometric luminosity is needed –and in particular less powerful Hot DOGs (10¹³ ≲ L_bol/L_⊙ ≲ 10¹⁴). A larger sample is also required to perform statistical studies regarding the presence of galaxy companions, the time lags between the different accretion phases, and in general to accurately characterize the most luminous, dust-obscured population of galaxies so far detected in the Universe. Deeper observations are also critical to be able to perform a detailed modeling of their kinematics (see Appendix A), and investigate whether simulations of massive galaxy formation (e.g., Dekel et al. 2009; Faucher-Giguère & Quataert 2012; Gaspari et al. 2013; Anglés-Alcázar et al. 2017b, 2021; Lupi et al. 2019) are able to reproduce the different morphologies, velocity fields, and potential quasar-driven outflows of the most rapidly evolving galaxies at the nodes of the cosmic web. These simulations are critical, as they can be used to infer and predict the duration and recurrence of Hot DOG phase(s), thus allowing to quantitatively translate snap-shots of space into actual, distinctive physical stages of galaxy evolution.

8. Summary

We present ALMA observations of the [C II] emission line for a sample of seven of the most extremely luminous (EL; L_bol > 10¹⁴ L_⊙) IR quasars, at redshifts z ∼ 3−4.6, from the population of hot, dust-obscured galaxies (Hot DOGs) discovered by WISE. We have compared their physical properties to other high-redshift populations of IR and optically selected QSOs at z ≃ 4.7 (IRQs) and ≳6 (OpQs), respectively. The results of the study are summarized as follows:

• The [C II] line is robustly detected in four out of the seven EL Hot DOGs (W0134–2922, W0220+0137, W0831+0140 and W2246–0526). Based on other ALMA emission line datasets, we know that for two sources without an identification (W0116–0505 and W0410–0913) the line was likely red-shifted out of the frequency range covered by the reference spectral window and side-band. All Hot DOGs are detected in dust continuum emission at λ_rest ≃ 158 μm.

• In a significant number of sources, the [C II] or other far-IR emission lines are systematically red-shifted with respect to the systemic velocities inferred from rest-frame UV spectroscopy, with an average shift Δv ≃ 780 km s⁻¹.

• The EL Hot DOGs in the sample exhibit the highest IR luminosity surface densities, Σ_IR, and the lowest L_[C II]/L_IR ratios among galaxy populations at any redshift, with values commonly ≥10¹³ L_⊙ kpc⁻² and ≲10⁻⁴, respectively. They also tightly follow the extrapolation of the L_[C II]/L_IR vs. Σ_IR relation found in nearby (U)LIRGs.

• All Hot DOGs detected in [C II] show L_[C II] surface densities of Σ_[CII] = 1−2 × 10⁹ L_⊙ kpc⁻² with a very narrow scatter. Only some of the most compact high-z OpQs and IRQs reach such high [C II] surface densities, but these populations also extend to surface areas one order of magnitude smaller. The [C II] equivalent widths of the Hot DOG sample range between ≃0.4−0.8 μm.

• The core line emission is spatially resolved clearly in 1/4 of Hot DOGs and marginally resolved in 3/4. In dust continuum, 3/7 are securely resolved, 3/7 are marginally resolved, and 1/7 is unresolved. The average [C II] and continuum emission sizes are ≃2.1 kpc and ≃1.6 kpc, respectively. In all cases when both line and continuum are detected, the core size of the [C II] emission is larger than that of the continuum, and the range of line-to-continuum size ratios is very narrow, 1.35−1.85 with an average of 1.61 ± 0.10, remarkably similar to other high-z QSO populations, and therefore suggesting the sizes could be linearly proportional. Although the core emissions in most objects are very compact, nearly all Hot DOGs display large extended emission fractions, from f_ext ∼ 20−80% in continuum to ∼50−80% in [C II].

• The velocity fields of EL Hot DOGs are very diverse: from barely rotating sources, to resolved hosts with ordered, circular motions, to complex, disturbed systems that are likely the result of ongoing mergers. However, all Hot DOGs are characterized by large line-velocity dispersions, FWHM_[CII] ≳ 500 km s⁻¹, which on average are significantly larger than optically selected QSOs at similar or higher redshifts. This is a clear signature of a turbulent ISM that may be a characteristic feature common to the EL Hot DOG phase, despite their diversity in rotation properties.

• As a population, the Hot DOG hosts seem to be centered around the star-formation MS. Individually, three of the most nearby Hot DOGs, at z ≃ 3, and W2246–0526 are consistent with being roughly on the MS, W0410–0913 and W0831+0140 seem to be lying above, and W0116–0505 below. However, uncertainties are large, and additional observations are needed to pin down the position of individual sources on the SFR–M_⋆ plane and better locate the population as a whole.

In summary, there is no unique morphological or dynamical state that describes the extremely luminous, AGN-powered Hot DOG phase: not all sources have galaxy companions identified in their fields; the morphology of their line and continuum emissions show relatively compact cores on top of very extended structures; and a wide variety of velocity fields are observed. The only common kinematic signature shared among the EL Hot DOG population is the large turbulence observed in the gas phase traced by [C II], which is likely related to their systematically large Σ_[CII] with respect to high-z OpQs and local (U)LIRGs with similar physical sizes. In the discussion section we explored what may be the source(s) powering their turbulent ISM, and investigated how different scenarios may agree with the particular location of the Hot DOG sample in a number of diagnostic diagrams enabled by the ALMA data. We also tried to place EL Hot DOGs within the merger-driven galaxy evolution framework and speculate on their origin and nature, proposing a possible hypothesis by which the most luminous Hot DOGs could be characterized as recurrent, yet short-lived periods (≲a few Myr; significantly shorter than the Eddington time) in which the galaxy at the center of a protocluster is undergoing rapid accretion in the form of companion galaxies, clumps or tidal structures, which frequently replenish the gas reservoir in the host and maintains its turbulence. This likely chaotic in-fall and lack of organized motions would promote the removal of angular momentum, allowing the SMBH to trap gas under its sphere of influence. There, matter can be subsequently funneled to the accretion disk, simultaneously keeping the quasar obscured. Some of the accreted gas can be expelled away in the form of feedback, further contributing to the turbulence in the host in a self-sustained process. We note, however, that other interpretations are also conceivable. In addition, the proposed scenario may not be valid for lower luminosity Hot DOGs and it is possible that not all the accretion events produce an EL hot, obscured phase.

The confirmation of the proposed picture thus requires higher signal-to-noise ALMA observations of a larger, statistically significant number of Hot DOGs for which a proper kinematic modeling of the galaxies can be performed and compared with detailed, physics-driven simulations. In addition, future JWST data will help to elucidate what is the main heating mechanism of the dust at different spatial scales in this outstanding galaxy population.

Acknowledgments

We thank the referee for their constructive comments and useful suggestions, which significantly improved the paper. We would like to thank Roberto Decarli and Bram P. Venemans for insightful discussions and early access to ALMA datasets from high redshift quasars. We would like to thank Jordan I. Penney for comments on early drafts of the paper. T.D-S. acknowledges support from the CASSACA and CONICYT fund CAS-CONICYT Call 2018. R.J.A. was supported by FONDECYT grant No. 1191124. H.D.J. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2017R1A6A3A04005158). G.C.J. acknowledges ERC Advanced Grant 695671 “QUENCH” and support by the Science and Technology Facilities Council (STFC). M.A. and R.J.A. have been supported by the grant CONICYT + PCI + Instituto Max Planck de Astronomia MPG190030. C-W.T. was supported by a grant from the NSFC (No. 11973051). S.E.L. was funded by the Chinese Academy of Sciences President’s International Fellowship Initiative. Grant No.2019PM0017/2021PM0076. J.W. acknowledges support from MSTC through grant 2016YFA0400702. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00576.S and ADS/JAO.ALMA#2015.1.00612.S. ALMA is a partnership of ESO (representing its member states), NSF (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. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Portions of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration, and of NASA’s Astrophysics Data System (ADS) abstract service.

References

Appendix A: Kinematic modeling

The continuum-subtracted [C II] cubes of each source (considering natural weighting –robust parameter equal to 2– and no uv-tappering), containing only the line SPW, were spatially trimmed to the central 200 × 200 spaxels, centered on the galaxy. After trimming, these cubes were fit with tilted ring models using the code ^3DBarolo (hereafter BBarolo; Di Teodoro & Fraternali 2015) in order to determine the dynamical masses of each source. First, the line emission in each cube was identified using the SEARCH algorithm in BBarolo, which is based on the code DUCHAMP (Whiting 2012). After automatically determining the rms noise level of the cube, this algorithm searches for pixels with intensities above a user-provided S/N ratio (S/N_upper). A three-dimensional (i.e., RA, DEC, velocity) search is then conducted around these peaks for emission above a second user-provided S/N ratio (i.e., S/N_lower), thus creating a 3-D signal mask.

Each cube is then collapsed over all channels containing line emission using the CASA toolkit task im.moments. The resulting integrated intensity (i.e., moment 0) map is fit with a 2D Gaussian using the CASA toolkit task im.fitcomponents, yielding a peak intensity, integrated flux density, FWHM of the major and minor axes, position angle, and central position. In addition, a velocity dispersion (i.e., moment 2) map is created using the same task, including only pixels identified through the SEARCH algorithm. With these maps in hand, the main function of BBarolo (i.e., 3DFIT) is used to fit a tilted ring model to the line emission identified through the SEARCH routine. While this task is able to automatically predict many parameters, and is well-tested for low-S/N and low-resolution observations, its performance is dependent on the initial parameter estimates and the overall model geometry. For our fits, the width of each ring is set to the FWHM of the minor axis of the restoring beam divided by 2.5, while the maximum radius of the model is set to the morphological major axis FWHM. The central position of all model rings is fixed to the morphological central position, while the initial guess for the velocity dispersion is the maximum value in the moment 2 map, and a thin disk is assumed (z_° = 0.01″). The inclination is allowed to vary between 10−80°, with an initial guess of 45°. BBarolo provides an initial guess for the rotational velocity and position angle.

Using these initial morphological and kinematic estimates, BBarolo creates a model of the first ring by populating a physical volume with discrete clouds, such that the rotational velocity, position angle, etc. of the estimates are recreated. This model is then converted into an observational (i.e., RA, Dec, velocity) cube, convolved with the synthesized beam of the data cube, normalized by setting the azimuthally-averaged flux in the model and data cubes equal, and compared to the model cube. Each parameter is varied, until the absolute residual (i.e., |model–observation|) is minimized. This process is then repeated with the next ring, and so on.

In the first run, a number of parameters were fit (i.e., rotational velocity, velocity dispersion, systemic velocity, inclination, and position angle). The three latter variables are then set to the average values between all rings, and the fit algorithm is run again, with only the rotational velocity and velocity dispersion as variables.

The final outputs of this process are morphological parameters (e.g., central position, integrated flux) and kinematic parameters (e.g, inclination, position angle, velocity dispersion and rotation curves, systemic velocity/redshift). We use these outputs to generate two estimates of the dynamical mass of the system. First, we consider a rotation-dominated system:

(A.1)

where r and v are the radius and velocity of the outermost ring, respectively, and G is the gravitational constant. Second, we consider a dispersion-dominated system:

(A.2)

where σ_v is the velocity dispersion of the outermost ring and α is a constant which typically takes a range of values from 2.5 to 8.5 depending on geometry. Here we follow Stott et al. (2016) and adopt a value of α = 3.4. From Figures A.1 and A.2, it is clear that the two systems analyzed here seem to be dispersion-dominated and therefore we use Equation A.2 to estimate their dynamical masses.

As shown in Figure A.1, the model is able to reproduce the kinematic properties of W0220+0137 to an excellent level, with very small moment 0 residuals. The only significant residuals appear in the moment 1 and 2 maps at the outskirts of the galaxy, where the S/N is low. The spectral residuals (bottom-right panel) are small as well. The galaxy has a well-defined rotation (moment 1 model panel), with circular velocities that range between [−100, 300] km s⁻¹. Despite the strong coherent rotation, the velocity dispersion is still large, with a peak of around FWHM ≃ 500 km s⁻¹. The morphological and kinematic position angles agree (41 ± 16° and 43 ± 15°, respectively), as well as the inclination on the plane of the sky (51 ± 14° and 48 ± 13°, also respectively).

Figure A.2 shows the modeling results for W2246–0526. The intrinsic kinematics of the galaxy are not well recovered by the model, as evidenced by the large residuals in all moment maps, including the moment 0, where the model cannot reach the peak intensity of the spectrum (bottom-right panel). The primary evidence for rotation is only visible in the position-velocity diagram along the visual major axis of the galaxy, and the velocity gradient, clear in the moment 1 map, is not captured in the modeled rotation. The model, instead, seems to favor the existence of a large intrinsic velocity dispersion over the entire galaxy of FWHM ≃ 500 km s⁻¹, in agreement with the results from Díaz-Santos et al. (2016). The inclination of the galaxy is constrained between 31 ± 17° and 47 ± 17°, and the position angle between 150 ± 50° and 163 ± 8°, corresponding to the morphological and kinematic estimates, respectively.

Fig. A.1. Dynamical modeling of W0220+0137 using BBarolo. The first three rows show the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field) for the data and model cubes (created using the CASA toolkit task im.moments), as well as the difference between them. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. In order to avoid the low-signal portion of the model, the model moment maps are only shown over the spaxels where the moment zero intensity is ≥1% of the maximum value. The bottom row shows (from left to right) the major axis postion-velocity diagram (PVD), minor axis PVD, and extracted spectra. For each PVD, the data is shown by the background color, while the overlaid contours represent the model. The data, model, and residual spectra in the bottom-right panle are depicted by the red, blue, and green lines, respectively.

Fig. A.2. Dynamical modeling of W2246–0526 using BBarolo. Same details as Figure A.1.

Appendix B: SED modeling

Fig. B.1. Results from the SED fitting. The photometric data of the sources are shown as black open circles. The AGN component is shown as a pink, solid line, while the E, Sbc and Im galaxy templates are displayed as red, green and blue dashed lines, respectively. The SED of the host galaxy (E + Sbc + Im) is shown as a yellow solid line, while the total SED of the source (including the AGN component) is displayed as a black solid line, with the synthetic photometry obtained from it (filter-convolved, color-corrected) shown as black open diamonds. The H-band flux of the host, used to obtain its stellar mass via scaling with a [M/L]H ratio (see text), is shown as a yellow open square.

This section describes the methodology used to infer the star formation rate (SFR) and stellar masses (M_⋆) of the quasar-host galaxies used in Section 6 to place the EL Hot DOGs in context with the MS of star-forming galaxies at the same redshifts. The analysis starts with the construction of the spectral energy distributions (SEDs) of galaxies from the optical to the mid-IR, based on HST, ground-based optical, Spitzer and WISE broad-band imaging. The SED coverage depends on the availability of photometric data for each source. Details about the telescopes and filters are provided in Table B.1, together with measurements and associated uncertainties. In all cases we collected the photometry with the smallest aperture available to avoid including emission from the neighbor galaxies. While the angular resolution of HST allows us to easily separate all companions, the Spitzer and WISE angular resolutions are coarser, and thus the fluxes could be contaminated. We use the AllWISE catalog to investigate the mid-IR WISE photometry, which is extracted via a point-spread-function (PSF) fitting procedure that also provides the goodness of fit of the model. In all WISE bands the galaxies in the sample are unresolved and no deblending was needed. Therefore, any contamination of extended structures beyond 6″ to the mid-IR fluxes is negligible. Regarding the near-IR bands, the Spitzer photometry was obtained with a 6″-diameter circular aperture using Sextractor (Griffith et al. 2012). While Sextractor performs a deblending process in crowded fields, extended emission could potentially contaminate the photometry. We note, however, that no other source of emission is evident in the IRAC imaging within the apertures. Minding these caveats, errors associated to uncertainties in the photometry are propagated via the Monte Carlo analysis performed to estimate the uncertainty in the H-band photometry of the host (see below).

We follow the approach of Assef et al. (2015) and use the low spectral resolution templates from Assef et al. (2010), spanning the wavelength range from 0.03 to 30 μm, to fit the photometric SEDs of the Hot DOGs and obtain an estimate of the host galaxy luminosity after removing the AGN contribution. There are four empirical templates in the suite: one used to fit the AGN component and three representative of galaxy hosts dominated by old (E template), intermediate (Sbc template), and young (Im template) stellar populations. The three galaxy templates allow for enough spectral versatility such that it is not necessary to account for additional dust extinction at the time of fitting, other than the intrinsic reddening the templates already have. On the other hand, the AGN template is allowed to have additional obscuration, which is modeled using a dust attenuation law that is a combination of a Small Magellanic Cloud (SMC) reddening law at shorter wavelengths and a Milky Way reddening law at longer wavelengths. Intergalactic medium absorption is also considered. Figure B.1 shows the best-fit SED model for each Hot DOG in the sample. The galaxy hosts dominate the rest-frame optical wavelengths in most cases, while the AGN becomes relevant above ≳1 μm. We note that, although broad and narrow optical emission lines are detected in the spectra of many Hot DOGs, the photometry is unlikely to be severely biased by these features, as the equivalent widths of rest-frame optical emission lines in Hot DOGs are small (Jun et al. 2020; Finnerty et al. 2020) compared to the width of the broad-band filters. Using HST imaging, Assef et al. (2020) has recently identified a few Hot DOGs displaying a rather blue UV/optical slope, most noticeable at λ_rest ≲ 3000 Å. The origin of this excess blue light is likely emission from the accretion disk of the central SMBH, which has been scattered into our line of sight by gas or dust. Two of the sources in the study of Assef et al. (2020) are in the Hot DOG sample analyzed in this work (W0116–0505 and W0220+0137). The contribution of this component to the rest-frame near-IR emission is, however, small (≲10% at H-band; see their Fig. 1) and therefore has a limited impact on the estimation of the stellar masses of Hot DOGs given the significantly larger uncertainties associated to the SED modeling assumptions.

After the best SED fit is found, the three galaxy templates are added together to obtain the rest-frame H-band flux density of the host, free from AGN contamination. In order to estimate the uncertainty associated with this value, we run 1000 Monte Carlo simulations for each Hot DOG, sampling from the errors of the photometric data-points assuming they follow a normal distribution, and calculate the H-band flux density for each instance. The uncertainty is obtained as the standard deviation of the simulated results. Because the galaxy (and AGN) templates are purely empirical, the SED fitting procedure does not provide information about the SFR history (SFRH) or any other intrinsic physical property of the stellar populations of the source. Therefore, in order to obtain the stellar masses of the host galaxies, an additional step is required; that is, a [M/L]_H ratio. To derive it, we use the code ezgal (Mancone & Gonzalez 2012) and the Galaxev stellar population SED models (Bruzual & Charlot 2003), which include a contribution from thermally pulsing-asymptotic giant branch (TP-AGB) stars, with a metallicity Z = 0.008 (≡0.4 Z_⊙, where Z_⊙ is the metallicity of the Sun) and a Chabrier intial mass function (IMF) (Chabrier 2003). We consider exponentially declining SFRHs that depend only on two parameters: the time of formation, t_f, defined with respect to the current age of the galaxy, t, where Δt = t − t_f, and the characteristic decaying time-scale of the SFR, τ. The SFRH is then parametrized as SFR(t′) ∝ e^{−(t′−t_f)/τ} with t′ ∈ [t_f, t]. The integral of the SFRH provides the M_⋆ of the model together with the resulting SED, from which the [M/L]_H is derived. In turn, the scaling with the data yields the normalization factor with which to obtain an absolute value of SFR(t).

Because it is extremely difficult to infer any information regarding the SFRH of galaxies based on sparse photometric data, for this work we decide to consider four simple limiting cases that are selected to represent reasonable lower and upper boundaries for Δt and τ, given the redshifts of the Hot DOGs, which correspond to between t ≃ 1.3 and 2.2 Gyr after the Big Bang. Specifically, we assume that star formation began Δt = {500, 100} Myr before the age of each galaxy, t, and two time-scales, of τ = {0.1, 1} Gyr. Note that the combination of Δt = 500 Myr and τ = 0.1 Gyr closely mimics a past starburst event, while the combination of Δt = 100 Myr and τ = 1 Gyr is equivalent to having an almost constant SFR. The M_⋆ values derived based on these four combination of parameters are used in Section 6 to place the Hot DOGs along the x-axis in the star-formation main sequence (Figure 13). In turn, the SFRs obtained from these SFRHs and scalings are not used in the calculation of the average SFR of each source in order to avoid the intrinsic correlation existing between SFR and M_⋆ for each assumed SFRH.

All Tables

All Figures

	Fig. 1. Optical spectra of the Hot DOG sample around the Ly-α emission line, whose rest-frame is set at zero velocity using zALMA (Table 1). The feature located at ∼5000 km s−1 is the N Vλ1240 Å emission line. The spectrum of W0410–0913 is not shown since Ly-α landed in a gap between detector chips.
In the text

Fig. 2. ALMA spectra of the luminous Hot DOG sample, including the four SPWs (black histogram). The average size of the synthesized beam of the cube is noted in the legend at the top left, as well as the radius of the aperture used for the extraction of the spectra, which is also shown as a dashed circumference in Fig. 6. The continuum spectrum (obtained from the line-subtracted cube) underlying the emission line in each source is shown in green, and the Gaussian fit to the [C II] line is shown in blue. The polynomial fit to the overall line-subtracted continuum as well as its ±1σ uncertainty are shown in red, with solid and dashed lines, respectively. The solid gray lines are the atmospheric transmission curve for a precipitable water vapor of 1 mm3, where the maximum transmission within the displayed spectral range has been normalized to the maximum y-axis value of the whole velocity range displayed. The vertical dashed line shows the expected redshift of the line based on the rest-UV spectrum, zUV.

In the text

Fig. 3. Velocity shift of zALMA with respect to zUV, Δv, as a function of zALMA (see Table 3). Hot DOGs are displayed as red squares. The z ≃ 4.7 IRQs from Trakhtenbrot et al. (2017) and Nguyen et al. (2020) are represented as yellow pentagons. The compilation of z ≳ 6 OpQs from Decarli et al. (2018) are displayed as green stars when the optical redshift is based on the Mg II emission line, and as blue stars when other lines or methods are used. Star-forming galaxies at 4 < z < 6 from the ALPINE survey (Cassata et al. 2020) are shown as pink pentagons.

In the text

Fig. 4. [C II] luminosity as a function of the rest-frame equivalent width (EW) of Ly-α. Hot DOGs are displayed as red squares. Normal star-forming galaxies at z ∼ 4−6 are shown as pink diamonds (Le Fèvre et al. 2020). The sample of Ly-α emitters (LAEs) at z ≃ 5−7 from Harikane et al. (2018) are displayed as gray circles, together with the fit to the correlation they found using the Ly-α EW corrected for inter-galactic medium absorption (their Eq. (25)). While we do not attempt such a correction in our sample of luminous Hot DOGs, the average increase in the EW of LAEs due to this effect is less than a factor of two. Hot DOGs deviate from the correlation by more than one order of magnitude, toward larger L[C II] for a given Ly-α (or vice versa). The arrows indicate a proportionally equal increase in the luminosity of both emission lines by a factor of 10 and 500 from the correlation toward W0134–2922 and W2246–052, and toward W0220+0137 and W0831–0140, respectively.

In the text

Fig. 5. Correlations between [C II] and IR luminosities, and far-IR sizes for the Hot DOG sample and the comparison high-z quasar and star-forming galaxy populations. Top-left: [C II] deficit (here defined as L[C II]/LIR) as a function of the IR luminosity surface density, ΣIR, color-coded as a function of the intrinsic 160 μm radius. Top-right: [C II] deficit as a function of LIR, color-coded as a function of ΣIR. Bottom-left: [C II] EW as a function of ΣIR, color-coded as a function of L[C II]. Bottom-right: [C II] surface density, Σ[CII], as a function of the intrinsic ∼160 μm radius (see Sect. 2.3), color-coded as a function of the IR luminosity. Hot DOGs are shown as squares, along with their names. The symbols of the comparison quasar samples are as in Fig. 3 (filled stars and diamonds). In addition, galaxy companions found in the neighborhood of some quasars in the sample of Decarli et al. (2018) are shown as open stars (from Neeleman et al. 2019). The local, z < 0.1 (U)LIRG sample from GOALS (Armus et al. 2009) is shown as small circles and limits (Díaz-Santos et al. 2017). Black crosses indicate local galaxies in which an AGN contributes ≥50% to the bolometric luminosity. In general, black symbols denote galaxies that do not have information regarding the property shown in the color bar. In all panels, the solid line is a fit to the local (U)LIRGs (the dotted lines show the 1σ standard deviation around the fit; see Díaz-Santos et al. 2017 for details), and the dotted-dashed line is the extrapolation toward lower and higher values of the x-axes. The dashed line in the top-left and bottom-right panels represents a constant Σ[CII] ≃ 2 × 109 L⊙ kpc−2, that seems typical of the Hot DOG sample, independently of their LIR. The two arrows in the top-left panel refer to the discussion in Sect. 7.1. We note that in this figure we are using LIR to show the [C II] deficits and luminosity surface densities, rather than the more traditional approach of using LFIR that is commonly found in the literature.

In the text

Fig. 6. Dust continuum flux density maps of the Hot DOG sample. The beam is shown as a white ellipse on the bottom-left, together with its size and the depth (rms) of the collapsed, line-subtracted cube. The dashed white circumference denotes the aperture used to extract the spectra shown in Fig. 2. The dotted contour shows the zero-flux intensity, and the solid contours represent the [2.5, 3.0, 2(4 + n)/2]-σ levels, where σ is the rms of the map and n = [0, 1, 2…]. Symmetric, negative levels are shown as dashed contours. Field of views (FoVs) of 2″ and 4″ on a side are displayed depending on the galaxy. Roughly, 1″ ≃ 7 kpc.

In the text

	Fig. 7. [C II] intensity (moment 0) maps of the Hot DOGs in which the line has been detected. The beam is shown as a white ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are displayed as in Fig. 6. The FoVs are also the same as in Fig. 6 to facilitate the comparison.
In the text

Fig. 8. Top two panels: Azimuthal profiles of the dust continuum flux density for the Hot DOGs with the smallest (top, W0116–0505) and largest (bottom, W0831+0140) fext, cont. Bottom two panels: Azimuthal profiles of the [C II] flux for the Hot DOGs with the smallest (top, W2246–0526) and largest (bottom, W0831+0140) fext, [CII]. The black line is the profile of the observed emission. The dashed line is the profile of the clean PSF. The dotted-dashed line is the scaled PSF that minimizes the scatter of the residuals within Bmaj/2 = 1 (see Sect. 3 for details and values in Tables 3 and 4). The dotted line is the profile of the residuals (observed – scaled-PSF). The PSF scaling value, fext and the radius used to perform the aperture photometry are noted at the top right of each panel (see dotted circumferences in Fig. 6).

In the text

	Fig. 9. [C II] projected radial velocity (moment 1) maps down to a 3σ level (see contours in Fig. 7). The beam is shown as a black ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are drawn at [0, ±50, ±100, ±200, ±300…] km s−1. The FoVs are the same as in Fig. 6 to facilitate the comparison.
In the text

	Fig. 10. [C II] FWHM (moment 2) maps down to a 3σ level (see contours in Fig. 7). The beam is shown as a black ellipse on the bottom-left, together with its size and the depth (rms) of the reference SPW used to create the map. Contour levels are drawn at [0, 25, 50, 100, 200, 300…] km s−1. The FoVs are the same as in Fig. 6 to facilitate the comparison.
In the text

	Fig. 11. Histogram of the [C II] FWHM for the sample of Hot DOGs (red line). The z ≃ 4.7 IRQs from Trakhtenbrot et al. (2017) and Nguyen et al. (2020) are shown in yellow. The compilation of z ≳ 6 OpQs from Decarli et al. (2018) is displayed as a green/blue dashed line. The 4 ≲ z ≲ 6 star-forming galaxies from Le Fèvre et al. (2020) are shown in pink.
In the text

Fig. 12. SMBH mass, MBH, as a function of dynamical mass, Mdyn, for the Hot DOGs W0220+0137 and W2246–0526. Literature OpQ and IRQ samples, from z ∼ 2 to z ∼ 6, are displayed for comparison (Bischetti et al. 2021, and references therein). The black solid line is the MBH versus dynamical bulge mass relation from the nearby galaxy sample of Häring & Rix (2004), and the dashed represents the 1σ intrinsic scatter. The dotted line is the MBH/Mdyn = 1 ratio. The SMBH mass of W0220+0137 is an upper limit (Finnerty et al. 2020), and the large red arrow pointing downwards represents the change in MBH if a λEdd = 2.8 is assumed, equal to that derived for W2246–0526 (see text for details).

In the text

Fig. 13. SFR as a function of M⋆ for the host galaxies of the Hot DOG sample, color-coded as a function of redshift. The galaxy star-formation MS at the specific redshift of each source is represented by solid lines (from Schreiber et al. 2015) also colored to match their redshifts. Mean values of SFR and M⋆ are shown as solid squares; median values are shown as open circles. Small, open squares represent the individual M⋆ and SFR estimates based on the SED modeling and the assumed SFR histories and scaling relations (see text and Appendix B). Note that the error bars associated to the mean values do not represent standard deviations but rather extend the full range of individual estimates.

In the text

Fig. A.1. Dynamical modeling of W0220+0137 using BBarolo. The first three rows show the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field) for the data and model cubes (created using the CASA toolkit task im.moments), as well as the difference between them. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. In order to avoid the low-signal portion of the model, the model moment maps are only shown over the spaxels where the moment zero intensity is ≥1% of the maximum value. The bottom row shows (from left to right) the major axis postion-velocity diagram (PVD), minor axis PVD, and extracted spectra. For each PVD, the data is shown by the background color, while the overlaid contours represent the model. The data, model, and residual spectra in the bottom-right panle are depicted by the red, blue, and green lines, respectively.

In the text

	Fig. A.2. Dynamical modeling of W2246–0526 using BBarolo. Same details as Figure A.1.
In the text

Fig. B.1. Results from the SED fitting. The photometric data of the sources are shown as black open circles. The AGN component is shown as a pink, solid line, while the E, Sbc and Im galaxy templates are displayed as red, green and blue dashed lines, respectively. The SED of the host galaxy (E + Sbc + Im) is shown as a yellow solid line, while the total SED of the source (including the AGN component) is displayed as a black solid line, with the synthetic photometry obtained from it (filter-convolved, color-corrected) shown as black open diamonds. The H-band flux of the host, used to obtain its stellar mass via scaling with a [M/L]H ratio (see text), is shown as a yellow open square.

In the text