© ESO 2020

1. Introduction

Active galactic nuclei (AGN) are characterized by accretion onto super massive black holes (SMBHs) which powers their emission over twenty decades in frequency. Almost all massive galaxies are thought to host a central SMBH which has passed (and possibly will pass) through active accretion phases during their life. Indeed, there is increasing evidence for a close connection between the star formation rate (SFR) and BH accretion history, suggesting a coevolution of the BH and its host galaxy; in particular, a crucial phase of such evolution occurs at z = 2 − 3, when both the SFR and BH accretion reach their maximum over cosmic time (e.g., Gruppioni et al. 2013; Madau & Dickinson 2014). The unified models postulate that nuclear absorption is due to a parsec-scale medium distributed in a toroidal geometry around the central engine (Urry & Padovani 1995). Nevertheless, the large amount of gas required to sustain the vigorous star formation (SF) of these objects during the epoch of their activity peak suggests that the UV and X-ray emission from the very central region of the galaxy can be obscured by the interstellar medium (ISM) located up to kiloparsec-scale in the host galaxy (Hickox & Alexander 2018). In addition, the high-z (3 <  z ≲ 5) heavily obscured (column density log N_H >  23) AGN fraction is estimated to be ∼70% (Vito et al. 2018), which is significantly higher than in the local Universe (intrinsic distribution ∼43% for the same column density, Burlon et al. 2011). This suggests an evolution of the luminous obscured AGN fraction, which may be driven by the larger gas content observed in high-z galaxies (e.g., Carilli & Walter 2013).

Several concurring physical processes are responsible for the formation of galaxies and the build-up of their mass: for example, major mergers between gas-rich galaxies (e.g., Hopkins & Beacom 2006), mass inflows from the cosmic web (Dekel et al. 2009; Bournaud et al. 2011), or fast collapse in early halos during which most of the mass of the system is accreted (e.g., Lapi et al. 2014, 2018); in all these scenarios the central regions are fueled with gas, triggering SF and BH accretion. Recent data from Atacama Large Millimeter/submillimeter Array (ALMA), Institut de radioastronomie millimétrique (IRAM) and Herschel indicate short gas depletion timescales (some 10⁸ yr, see Genzel et al. 2015; Tacconi et al. 2018). The submillimeter galaxies (SMGs, defined as submilllimeter objects having flux densities ≳1 mJy at 850 μm, Carilli & Walter 2013; see also Circosta et al. 2019 and references therein) are often indicated as those objects that trace the brief phase at the peak of SFR in the galaxy evolution. They are usually detected at median redshift z ∼ 2 − 3 (Simpson et al. 2014), and are characterized by extremely dense ISM, showing molecular masses up to 10¹⁰ M_⊙ and dust masses up to few ×10⁸ M_⊙ distributed on kiloparsec scales (Swinbank et al. 2014; Barro et al. 2014; Ikarashi et al. 2015; Massardi et al. 2018; Lapi et al. 2018). The estimated SFR of these objects is ∼ 100 − 1000 M_⊙ yr⁻¹ (Blain et al. 2002; Swinbank et al. 2014). These properties make SMGs the best targets to investigate the contribution of the host galaxy in obscuring the AGN. Moreover, many observations revealed that these sources have morphologies showing complex gas motions and that the star-forming gas can lie in compact, rotating disk structures (Tacconi et al. 2008; Hodge et al. 2012; Ivison et al. 2013; De Breuck et al. 2014). Given their high gas and dust density and small size, SMGs are also considered among the best candidates as progenitors of compact quiescent galaxies (cQGs; Cimatti et al. 2008), a class of quiescent galaxies that are considerably smaller and denser than local galaxy counterparts of similar mass, showing stellar half-light radii of ∼0.5 − 2 kpc and typical stellar masses of ∼10¹⁰ M_⊙ (Cassata et al. 2011).

Estimating the column density N_H associated with the ISM of the host galaxy requires a multiwavelength approach exploiting data from the X-rays to the submillimeter: on the one hand, the X-rays directly probe the innermost region of the sources and can be used to evaluate the line of sight absorption to the central engine emission. On the other hand, the mm band is crucial to estimate the ISM properties such as the mass and geometrical distribution of the gas and dust components. In particular, measurement of the gas mass, dominated by the molecular component at z ≳ 2 (Carilli & Walter 2013), can be achieved with different methods: the most common exploits the low-J transition of CO as a tracer of the total molecular gas mass. An alternative method exploits the monochromatic continuum dust luminosity at 850 μm (L_850 μm): Scoville et al. (2016) showed an empirical relation between L_850 μm and the CO(1–0) emission line luminosity – valid for local ultraluminous infrared galaxies (ULIRGs, L_{8 − 1000 μm} ≳ 10¹² L_⊙), z ∼ 2 SMGs and local star-forming galaxies – from which they derived a direct relation between L_850 μm and the molecular gas mass, assuming a Galactic CO-to-H₂ conversion factor. Once the gas mass and size are measured through high resolution observation, the gas column density can be derived assuming a geometrical configuration.

Gilli et al. (2014) found that the column densities derived from the X-ray spectral analysis and the ISM mass are comparable for a SMG at z = 4.755 hosting a heavily obscured AGN (XID 539, included in our sample, Table 1). This source is also included in the sample studied by Circosta et al. (2019), who performed an in-depth analysis on seven X-ray selected, far infrared (FIR) detected AGN with spectroscopic redshift z >  2.5 in the 7 Ms Chandra Deep Field-South (CDF-S), in order to assess whether and to what extent the host ISM contributes to the AGN obscuration. They performed X-ray spectral fitting to derive the column density responsible for the central obscuration, assuming different spectral models. Then they exploited the available UV-to-FIR photometric data to derive the 850 μm flux from spectral energy distribution (SED) fitting, and they converted it into gas mass assuming the $L_{1-0}^{\prime }\propto L_{850\mu \mathrm{m}}$ Scoville et al. (2016) empirical relation and the CO-to-H₂ conversion factor α_CO = 0.8 M_⊙ (K km s⁻¹ pc²)⁻¹. On the basis of a simple assumption of spherical gas distribution with uniform density, they derived column densities of the host ISM that are comparable with those extracted from the X-ray fitting, and concluded that the ISM likely plays a significant role in the obscuration.

In this work we analyze for the first time the millimetric emission of the high-J CO emission and the dust continuum of six heavily obscured AGN (X-ray column densities N_H, X >  10²³ cm⁻²) in the CDF-S at spectroscopic z >  2.5, included in the sample of Circosta et al. (2019) and observed with ALMA. In Sect. 2 we illustrate the criteria adopted to select the targets and summarize the sample properties. In Sect. 3 we describe the ALMA observations and the data reduction process. In Sect. 4 we perform the imaging and obtain the flux density, the angular size and line shape of the sources from both the spectrum and visibility fitting. We exploit these quantities in Sect. 5 to derive the mass and physical size of both the dust and molecular gas content. As for the gas mass, we explore several cases with different assumptions on the involved conversion factors, to address the large uncertainties in deriving such mass from both high-J CO transitions and dust long-wavelength emission. Then we derive the column densities for a simple spherical geometry and compare them with those measured from the X-ray spectra by Circosta et al. (2019). We also calculate the dust mass in the optically-thin regime and the dynamical mass of the two sources for which the rotating disk model has been tested. Finally, in Sect. 6 we summarize and discuss the results and future perspectives.

Throughout the paper we adopt a standard ΛCDM cosmology: H₀ = 69.6, Ω_M = 0.286, Ω_Λ = 0.714 (Bennett et al. 2014). As for the reported velocities we adopt the radio definition v = (ν₀ − ν)/ν₀, where ν₀ and ν are the rest and observed frequency, respectively, in a Barycentric (i.e., referred to the Sun and Earth center of mass) rest-frame.

2. Sample selection

The Chandra Deep Field South (CDF-S, ∼480 arcmin², Luo et al. 2017) provides the best X-ray spectral information currently available for distant AGN, thanks to its 7 Ms exposure time (the average flux limit within 1 arcmin from the center of the field is F_{0.5 − 7 keV} ≈ 2 × 10⁻¹⁷ erg s⁻¹ cm⁻², Luo et al. 2017). In addition, this field features an exceptional multiwavelength coverage that allows the characterization of the overall SED of most X-ray sources with extraordinary accuracy. In particular, deep IR observations have been carried out with the Hubble Space Telescope (HST), Herschel and Spitzer observatories. Moreover, many observations in the submillimeter band were taken with instruments such as the Submillimeter Common-User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope and the Large Apex BOlometer CAmera (LABOCA) on the Apex Telescope. The central area of the CDF-S (about 1/3 of the whole field) is covered by the Great Observatories Origins Deep Survey (GOODS-S, Giavalisco et al. 2004), in which the Cosmic Near-IR Deep Extragalactic Legacy Survey (CANDELS, Grogin et al. 2011), a powerful imaging survey of the distant universe, has been carried out with HST. CANDELS consists of two different areas and depth portions, the deeper and smaller CANDELS/Deep and the shallower and more extended CANDELS/Wide. The CDF-S was also observed by the Galaxy Evolution from Morphology and SEDs (GEMS, Rix et al. 2004) survey, a large-area (800 arcmin²) two-color (F606W and F850LP) imaging survey with the Advanced Camera for Surveys (ACS) on the HST.

In order to select our sources, we combined the sample from Vito et al. (2013), consisting in 34 X-ray selected AGN at z >  3 in the 4 Ms CDF-S with eight heavily obscured quasars (QSOs) in the redshift range 1.1–3.7, X-ray selected in the 1 Ms CDF-S exposure and observed at 850 μm with SCUBA (Rigopoulou et al. 2009). To select our targets from this combined sample we applied the following criteria. Firstly we required a secure (errors < 300 km s⁻¹) spectroscopic redshift > 2.5 as given by Luo et al. (2017), in order to avoid the large uncertainties related to the photometric redshift that would propagated in the derived fundamental quantities for our study (such as the intrinsic luminosity) and properly select the CO emission lines to target in the ALMA observations. The redshift threshold at z = 2.5 was chosen as a compromise between a sufficiently large gas content in the galaxies (expected to be larger at high-z, such that the host galaxy obscuration might be important) and the size of the sample. Secondly we required at least one > 3σ detection at λ_obs >  100 μm, in order to characterize the FIR-SED and maximize the detection probability with ALMA. Finally we included among the sample only the sources with gas column density log N_H >  23 (derived from the X-ray spectral fitting), in order to consider only heavily obscured objects and trace the evolutionary stage in which the AGN is still deeply buried and the SF is high.

This selection returns six targets spanning the redshift range ∼2.5–4.8. This sample includes six of the seven sources studied by Circosta et al. (2019), who adopted our same selection criteria except for the requirement of heavy obscuration. As pointed out by those authors, the fact that the FIR detection criterion is sufficient to select only obscured AGN indicates a connection between the SF-heated dust in the ISM that is responsible for the FIR to (sub)millimeter emission and the nuclear absorption. Despite the small number of sources (owing to the stringent selection criteria), this sample represents a complete and unique collection of objects: they are characterized by a relatively high number of photons (from ∼260 to ∼2000 in the 7 Ms exposure) in the observed 0.5–7 keV band which allows an accurate estimate of the X-ray absorption and by good IR coverage. The IDs, coordinates and redshift are reported in Table 1. Throughout the paper we adopt for our targets the Luo et al. (2017) X-ray Identification numbers (XIDs).

The selected objects are all included in CANDELS/Deep area, except for XID 539 that lies in CANDELS/Wide and XID 42 that is outside CANDELS, but it is included in the GEMS catalog, in the multiwavelength survey by Yale-Chile (MUSYC) catalog (consisting of deep optical medium-band and broadband optical and near-IR imaging of the CFDF-S, Gawiser & MUSYC Collaboration 2003), and it has been also observed in the H₁₆₀ band by HST wide field camera 3 (WFC3; Chen et al. 2015).

Circosta et al. (2019, see Table 4) derived their basic quantities (such as SFR, IR luminosity and stellar mass) by the SED decomposition technique, finding that the targets are massive (stellar mass ∼10¹¹ M_⊙) and very luminous (infrared luminosity in the 8–1000 μm band ≳10¹² L_⊙ yr⁻¹) galaxies. They derived SFRs of ∼200 − 900 M_⊙ yr⁻¹ (except for XID 42 which has ∼1680 M_⊙ yr⁻¹), which places all sources within 1σ from the best fit main sequence relation found by Schreiber et al. (2015) at the corresponding redshift (2σ for XID 337 which is below the main sequence). Three of them are well known heavily obscured AGN: XID 337 was firstly detected in the early 1 Ms CDF-S exposure and classified as a type II QSO object in the spectroscopic follow-up (Szokoly et al. 2004). This result was later confirmed by Mainieri et al. (2005) who also reported a submillimeter detection at 850 μm using SCUBA. Exploiting high resolution (∼0.2″) ALMA observations, Circosta et al. (2019) found that the single object associated with XID 337 in the Spitzer/MIPS, Herschel and SCUBA maps is actually a blend of two sources, whose emission is dominated by a bright QSO at 3.5″ separation; thus, the MIR/FIR photometric datapoints of XID 337 were converted into upper limits. XID 551 was the first type II AGN ever detected in the CDF-S 1 Ms exposure at z >  3 (Norman et al. 2002), for which Comastri et al. (2011) found a large column density (log N_H >  23) performing an X-ray spectrum analysis using 3.3 Ms XMM-Newton survey data. XID 539 (ALESS73.1, Hodge et al. 2013) is the most distant Compton-thick AGN known (Gilli et al. 2011), and one of the widely studied sources in the CDF-S (Coppin et al. 2010; Nagao et al. 2012; De Breuck et al. 2014). In particular, Gilli et al. (2014) have performed an in-depth analysis through an observation of the continuum at ∼1.3 mm, combining it with continuum observations from ALESS and from Nagao et al. (2012). They built a full SED including Herschel data and obtained dust temperature and size and the SFR. XID 666 (BzK 4892) is a Compton-thick QSO studied by Feruglio et al. (2011), Del Moro et al. (2016), Perna et al. (2018). XID 42 (ALESS57.1) has been recently studied by Calistro Rivera et al. (2018) an by Rybak et al. (2019), who resolved the CO(3–2) and [C II] emission, respectively. XID 746 is the less studied objects, despite being known as Type II QSO since the first 1 Ms exposure in the CDF-S (Szokoly et al. 2004).

3. Observations

ALMA dramatically improved our knowledge of SMGs, providing us with extremely deep and high-resolution images and spectra, which can be used to probe the dust mass and temperature of these objects, their star formation activity, and their gas mass.

We performed ALMA observations (Cycle 4, Band 4, ID: 2015.1.01205.S, PI: R. Gilli) of the continuum at ∼2.1 mm (142.78 GHz) and of one high-J CO transition for each of the targets, reported in Table 2 along with their rest- and observed-frame frequencies. The requested sensitivity was ∼180 μJy beam⁻¹ for the line emission (considering a line width of 200 km s⁻¹) that corresponds to 20 μJy beam⁻¹ for the continuum.

The observations were split in two blocks, each containing three targets: the 2 × 4 GHz spectral windows (spws) cover the range 126.56–130.33 GHz and 138.56–142.33 GHz to observe XID 666, XID 746 and XID 539, and the range 133.07–136.93 GHz and 145.17–148.93 GHz for XID 42, XID 337 and XID 551. Each spw is sampled in 240 channels of 7.812 MHz width corresponding, at the median frequency of the spws, to 17.5 km s⁻¹ in velocity for the first block and to 16.6 km s⁻¹ for the second one. The project was accepted in C priority and was only partially observed at the end of the Cycle 4. In particular, the first block was observed on the 23rd September 2016 for about 88 min, and the second block for about 46 min on the 29th September 2016. QSO J0334−4008 was used to calibrate the flux and bandpass, while the QSO J0348−2749 served as phase calibrator¹.

Three out of six sources (namely XID 42, XID 539 and XID 666) were detected in both the line and continuum emission. Concerning the undetected sources, two of them (XID 337 and XID 551) were observed for half of the time needed to achieve the requested sensitivity, while XID 746 is included in the execution block that was observed twice and reached the requested sensitivity.

We performed data-reduction using the ALMA calibration pipeline in CASA version 4.7.0-1. The pipeline procedure automatically flags bad data and performs the calibration; it also produces several diagnostic plots that are useful to identify possible calibration failures and discern additional needed flags. After the inspection of both the diagnostic plots and calibrated data, the calibration resulted to be satisfactory.

4. Data analysis

We performed a standard imaging of the sources by using the CLEAN task of CASA in “velocity” mode for the line emission and “mfs” mode for the continuum. The imaging parameters were chosen as best trade-off between sensitivity and angular resolution (see Sects. 4.1 and 4.3 for details, Table 2 for the resolution and sensitivity of each observation).

4.1. CO-line imaging

Since the line is not clearly visible from the visibility amplitude-frequency plane, we cannot unambiguously assess the channels containing it, thus we cannot a priori exclude them when we image the continuum emission. Thus we decided to image separately the spw containing the line and the other spws for the continuum, in order to avoid any possible contamination. To image the continuum-subtracted line we preliminarly ran the UVCONTSUB task, that performs a polynomial fitting of the continuum level and then subtracts it from the data. We fit the continuum using only the spws that do not contain the line, and subtract the resulting fit from the spw in which the line is expected to be. For XID 337, XID 551, and XID 746 it was not possible to determine unambiguously a channel in which the brightness at the source position was significantly above the noise level. For these sources no detection was obtained even applying natural weighting, tapering (0.45″) and a large channel width (120 km s⁻¹) to improve the chance of achieving a significant detection. For XID 42, XID 539 and XID 666 imaging we used a robust weighting equal to 0.5, 1, 1, and a channel width of 40, 120 and 40 km s⁻¹, respectively. We clearly identified the channels corresponding to the peak of the emission, whose images are reported in the top row of Fig. 1. We detect XID 42, XID 539 and XID 666 with a signal-to-noise ratio (S/N) equal to 6.2, 4.7 and 4.0/4.9, respectively. The S/N is defined as the ratio between the value of the brightest pixel and the root mean square (rms) derived averaging the rms measured in four regions placed at ∼10″ from any source. Since for XID 666 we detected a double peaked line, the two values refer to the blue- and red-shifted peaks. As a confirmation of the detection, we note the central position of the peak in the images (offset with respect to expected source positions < 0.5″, beam sizes in Table 2) and the line feature in the spectra visible in Fig. 2.

Fig. 1. Top row: CO-line peak channel of the three detected sources. The observed CO transition and the velocity corresponding to the center of the channel (defined as the velocity shift with respect to the systemic redshift) are reported in the up-left and bottom-right corner of the images, respectively. The blue ellipses are the regions from which the spectra in Fig. 2 were extracted. In the XID 666 CO-line peak image (right panel) we show the contours of the blue-shifted peak channel centered at −206.5 km s−1 (2, 3, 4σ, black solid lines) over the red-shifted peak channel (see Sect. 4 for details). Bottom row: dust-continuum image of the three detected sources. In all the images the black ellipse in the left-bottom corner represents the clean beam.

Fig. 2. Spectrum of XID 42, XID 539 and XID 666 (green shaded histogram), expressed as the mean brightness within the regions contained in the blue ellipses of Fig. 1. The black solid line represents the Gaussian line fitting. Due to its double-peaked feature, we fit the line of XID 666 with two Gaussian components, indicated in the plot with the blue and red dashed lines. The channel width Δv is reported for each spectrum.

We note that only for XID 539 the velocity of the peak channel (reported on the image) is consistent with zero, being the channel central velocity −111.6 km s⁻¹ and its width is 120 km s⁻¹. XID 666 has two peaks shifted respectively by ∼​​​​−200 km s⁻¹ and ∼150 km s⁻¹ centered on the rest-frame velocity. This double-peaked line, also clearly visible from the spectrum (Fig. 2), is interpreted as a signature of a rotating structure. This interpretation is also supported by the velocity map (Fig. 3) that shows a continuous gradient along one direction. However, in Sect. 5.5 we briefly explore the possibility that this source is actually a blend of two nonresolved, merging system. The spectroscopic redshift of XID 42 (z = 2.940) was derived by Szokoly et al. (2004) from the FOcal Reducer and low dispersion Spectrograph (FORS) of the Very Large Telescope (VLT). Re-analyzing the same spectrum using the high excitation lines CIV and NV, and the recombination line HeII (that is a robust redshift estimator) we found a better estimate z = 2.9369 ± 0.0003, assumed as the rest-frame hereafter (against the value z = 2.939 obtained using the Ly-α line, which is known to be challenging as exact redshift estimator due to the resonant scattering). From the CO(5–4) line we derived a redshift z = 2.94355 ± 0.00001, which is shifted by ∼500 km s⁻¹ with respect to the rest-frame velocity. We note that Calistro Rivera et al. (2018) derived a redshift z = 2.943 ± 0.002 from the CO(3–2) line, which is compatible within 3σ with that derived both from the optical spectrum and from the CO(5–4) transition. This redshift is also consistent with the [C II] line detected by Rybak et al. (2019). The rest frequency displacement between the CO and [C II] lines and the UV emission lines could be explained as a difference in radial velocity. Such as spectral displacements have been often reported in the literature (Pentericci et al. 2016; Matthee et al. 2017; Carniani et al. 2018), and interpreted as a consequence of the clumpy structure of high-z galaxies, in which the different components move with respect to each other. High excitation UV-rest-frame lines and the CO (and [C II]) lines indeed trace different gas phases and may arise from different regions.

Fig. 3. Total flux (moment 0), velocity map (moment 1) and velocity dispersion map (moment 2) (from top to bottom) extracted for XID 42, XID 539 and XID 666 (from left to right). Only the pixels above the 3σ level are shown. In the top right corner of the moment 0 images the angular scale is reported. In all the images the black ellipse in the left-bottom corner represents the clean beam.

4.2. Moments

By using the IMMOMENTS task we computed the moments 0, 1, 2 for the detected sources, defined as the frequency-integrated flux, the velocity field and the velocity dispersion, respectively. We performed a Gaussian fitting of the line spectrum and derived their central velocity v₀ (with respect to the rest-frame velocity at the given redshift) and Full Width at Half Maximum (FWHM); to calculate the moments we considered a velocity range given by v₀ ± FWHM, in order to include few channels aside the line and avoid the loss of signal from its wings. For XID 666, to calculate the moments of the entire source, we consider the channels spanning the velocity range from v₀ − FWHM of the blue component to v₀ + FWHM of the red one. In Fig. 2 we show the spectra of the sources and their fits.

We fit the lines of XID 42 and XID 539 with a single Gaussian component, since there are not evident double-peaked features in their spectra. As for XID 666, we fit two Gaussian components on the whole source region (blue ellipse in Fig. 1). Since the two spectral peaks of XID 666 are quite close and likely contaminated by each other, we additionally performed fitting with a single Gaussian component on the spectra of the two spatially separated regions corresponding to the peaks, in order to estimate such contamination and derive the most accurate line width. Both the line width and central velocity of the single blue and red shifted components are fully consistent (within 1σ) with those derived from the spectrum over the entire source region. In the following, we use the line width and central velocity derived from the two Gaussian components fit on the whole region. The fit results are reported in Table 3. We note that the brightness of the peak is not used in any of the following considerations and derived physical quantities.

We extracted the moments using as a threshold a value of 3σ calculated by the measured rms on the image. Moments are reported in Fig. 3.

4.3. Dust continuum imaging

For each source we imaged all the spws but those containing the CO line. As for the line imaging, for XID 337, XID 551, and XID 746 we apply natural weighting and a 0.45″ tapering, without achieving a significant detection. We detected the continuum emission for XID 42, XID 539 and XID 666 with a S/N equal to 6.8, 12.8 and 7.0, respectively, using a robust weighting equal to 1 for all of them. We show the continuum images of the detected sources in the bottom row of Fig. 1.

4.4. Visibility fitting

Since the reconstructed images depend on the deconvolution algorithm and its input values, they do not provide a unique and direct representation of the data but rather a nonlinear mapping from the Fourier domain into the sky plane. Hence, we firstly perform an image fitting with the sole purpose of using the obtained values as starting input parameters to perform the fits directly on the visibilities, which are unaffected by the reconstruction algorithm and pixel correlation as images are. In addition, if the sources are not or barely spatially resolved in the images, the visibility fitting provides an estimate of their angular size even on scales smaller than the synthesized beam.

As for the image fitting, to derive the information about the dust component, we fit the continuum image, while for the gas we fit the moment 0 of the CO-line since it represents the total integrated emission. We performed the fitting using the task IMFIT, which relies on the algorithm developed by Condon (1997) using as a model a spatial elliptical Gaussian 2D function. To select the region in which the fit is performed, for the continuum we chose a region relatively large (∼1.5″ diameter) around the sources to include all the signal, setting a 3σ threshold on the pixels to include. Given the high S/N this selection guarantees that only the source will be fit. For moment 0 this is not an issue because all the included pixels are above the 3σ significance level and are attributed to the source. For the line fitting we included the same channels (v₀ ± FWHM) used to calculate the moments.

Since the task UVMODELFIT, which performs the visibility fitting, takes the flux density as input, while IMFIT returns the velocity-integrated flux, in order to obtain the starting value for the visibility fitting of the CO line emission we divided the integrated flux obtained from the image fitting by the FWHM of the line reported in Table 3. We included the same native channel range corresponding to v₀ ± FWHM (on the continuum subtracted dataset) for the line fitting, and the rest of the band for the continuum. We initially left all the parameters free in order to check whether the central position of the Gaussian matches that derived from the image fitting. Since the offset was within 1σ, we rerun a second fit keeping the position fixed in order to reduce the uncertainties related to the flux densities and the sizes. Final fitting results are reported in Table 4.

The large error associated with the position angle is due to the fact that the sources are barely resolved and it is difficult to unambiguously derive an orientation. However, this does not affect our analysis since we are interested in the size of the sources and not in their orientation with respect to the celestial axes. For what concerns the dust, the continuum axial ratio (R = b/a, where b is the minor axis and a is the major axis) of XID 666 is consistent with 0. This is attributed to the very small size of the component that cannot be resolved in the minor axis direction.

5. Results and discussion

5.1. Source size

Given the source redshifts, we derived the physical source size from the resulting fit. For XID 666, the continuum axial ratio could not be estimated either from the visibilities (Table 4) or from the deconvoluted image fitting. However, from the images we obtained the size of a and b convoluted with the clean beam. Hence, for XID 666, a fair assumption on the value of R is that derived from the convoluted a and b, that is R = 0.8. This ratio can be affected by the convolution given that, for Gaussian beams and sources, sizes add in quadrature. However, the derived ratio is consistent within 1σ with those derived for the other sources, and given the large uncertainties this assumption does not significantly affect our results. We do not propagate the errors on the axes as obtained from the image fitting to estimate the uncertainty in the continuum axial ratio of XID 666, since the error on R derived from the visibility fitting is generally larger, as visible from the fits of the other sources in Table 4. These (relative) errors span the range 30%–50%. Thus, using a conservative approach, we assigned a relative error of 50% to R for XID 666, corresponding to the maximum error obtained from the other fits. In Table 5 the physical sizes of each source are reported for both the molecular gas and dust components. In order to derive an upper limit of the gas density, for the undetected sources (XID 337, XID 551 and XID 746) we assume as “typical” major and minor axes the mean of those of the detected ones. Several analyses have shown that this kind of objects have all comparable sizes (e.g., Swinbank et al. 2010; Tacconi et al. 2008). In addition, the mean redshift z = 3.42 of the detected sources is similar to that of the undetected ones (z = 3.47), thus we may expect that they have a similar size. We still use a conservative approach assigning a 30% relative error to a and 50% to b, corresponding to the maximum uncertainties obtained from the visibility fits. The assumed a and b of XID 337, XID 551, XID 746 are 2.7 ± 0.8 kpc and 1.6 ± 0.8 kpc for the CO component, and 2.0 ± 0.6 kpc and 1.2 ± 0.6 kpc for the dust component, respectively.

We found that the sources are extremely compact objects, all having a major axis of the Gaussian fit ≲3 kpc. These results are in agreement with the hypothesis that these sources may be identified as the likely progenitors of the compact Quiescent Galaxies (cQGs), a class of extremely compact objects (≲1 kpc, Cimatti et al. 2008; Talia et al. 2018) that have experienced a brief phase of intense SF, eventually quenched by the SF and AGN feedback. We note that the molecular gas size obtained from the CO(5–4) line fitting for XID 42 is nearly half of that derived by Calistro Rivera et al. (2018) on the basis of the CO(3–2) observation, and more similar to both the optical HST H₁₆₀ and [C II] size, whose emission is dominated by a compact (radius < 1 kpc) component (Rybak et al. 2019). The difference between the CO(3–2) and CO(5–4) region sizes could be attributed to the fact that the CO(5–4) emission traces the innermost, denser and more excited phase of the gas. Indeed, Engel et al. (2010) showed how the extent of the CO emission may decrease for increasing J-transitions. However, we cannot exclude that the size of the CO(5–4) component is biased by the sensitivity limit of the observation. As for XID 539, the dust size is in agreement with De Breuck et al. (2014), who constrained the diameter of the 1 mm continuum emission to a region < 2 kpc. They also resolved the [C II] line emission, which has been spatially fit with a Gaussian, obtaining a size of 2.4 ± 0.2 kpc half-light radius, corresponding to about twice the CO(7–6) radius detected in this work (defined as the half of the mean between the major and minor axes). However, we note that – due to the large error on the CO(7–6) size – the two measurements are still consistent within ∼2σ.

We also found that the dust component is confined in a smaller region, (∼70% radius for XID 539 and XID 666, ∼90% for XID 42) than the molecular gas component. Our results are in agreement with the typical size of the dust in SMGs that has been shown to be ≲2 kpc (e.g., Ikarashi et al. 2015; Hodge et al. 2016). Calistro Rivera et al. (2018) showed how the (apparent) compactness of the dust can be explained by a radial decrease of the gas column density and dust temperature toward the outskirts of the sources, on the basis of a radiative transfer model presented by Weiß et al. (2007). It is interesting to note that Casasola et al. (2017), analysing the gas and dust distribution of a sample of 18 nearby face-on galaxies, found that – on average – the exponential scale-length of the dust-mass surface-density is ∼2.3 times higher than that of the molecular gas, that is the dust is distributed at larger radii than the molecular gas, contrary to what we found at high-z. However, we point out that given the sensitivity limit of our observation it is possible that we are tracing only the innermost, warmer phase of dust, whose emission flux increases linearly with the density and temperature (in the Rayleigh–Jeans regime, assuming a constant dust mass).

5.2. Dust mass

We derive the mass of the dust M_d from its emission in the optically thin regime:

$$\begin{array}{c}\hfill M_{\mathrm{d}}=\frac{D_{\mathrm{L}}^2{S}_{\mathrm{obs}}}{k_{\nu }{B}_{\nu }(T_{\mathrm{d}})(1+z)}\end{array}$$

where z is the redshift, D_L is the luminosity distance, S_obs is the flux density at the observed frequency ν_obs, B_ν(T_d) is the Planck function of the dust emission at the temperature T_d and k_ν is the dust opacity per mass unit (g⁻¹ cm²), assumed to scale with the frequency as k_ν = 4[ν/(1.2 THz)]^β (Draine & Li 2007). The index β is set equal to 2.0 (e.g., Magnelli et al. 2012; Gilli et al. 2014). For the undetected sources, an upper limit on S_obs is assumed equal to three times the rms measured on the continuum image. The rest-frame frequency ν = ν_obs(1 + z) is computed considering as ν_obs the median observed frequency of the spectral windows in which the imaging and fitting of the continuum have been performed. The temperatures T_d of the dust have been derived from the SED fit with a single modified black body (MBB) component S_ν ∝ B_ν(T_d)(1 − e^−τ_ν) by means of the code developed by Feltre et al. (2012), see also Circosta et al. (2019). The fiducial error associated with such temperatures is on the order of ≈ ± 5 K. The photometric datapoints used in this work to perform the SED fitting and derive the dust temperatures are the same of the SEDs reported in Circosta et al. (2019), taken from the UV- to mid-IR (MIR) provided by Hsu et al. (2014), complemented with the FIR data from Herchel/PACS and SPIRE and reported in the Magnelli et al. (2013) and Oliver et al. (2012) catalogs, respectively. Furthermore, we also added the measurements reported in this work at 2.1 mm. In addition to this, the other main difference with respect to Circosta et al. (2019) is that they fit the FIR SED using several templates accounting for the stellar, star formation and AGN contributions, rather than a single MBB component. We discuss the consequences of these differences in Sect. 5.3. The temperatures and dust masses are reported in Table 6. In fact, the high temperatures obtained by fitting only a single component to the SED could indicate that we are tracing only the warmer (brighter) phase of dust, and a second component would be required to characterize the colder (fainter) phase. In this sense, we point out that these temperatures (and masses) are luminosity-weighted. Unfortunately, the data coverage of the SED does not allow us to properly constrain a double-component MBB. According to Santini et al. (2014), the attempt of reproducing both the Wien and the Rayleigh–Jeans side of the MBB spectrum produces an overestimate of the dust temperature, and in turn an underestimate of the dust mass. However, our results are in agreement with Rybak et al. (2019) who measured a $M_{\mathrm{d}}=4.1_{-0.6}^{+0.5}\times {10}^8{M}_{\odot }$ for XID 42 by means of SED fitting, and with Gilli et al. (2014) who derived for XID 539 a M_d = 4.9 ± 0.7 × 10⁸ M_⊙ on the basis of the continuum emission at 1.3 mm. We stress that the dust temperatures do not affect the gas masses in our treatment, since we derive them in the RJ regime (see Sect. 5.3 for the discussion).

5.3. Gas mass

The mass of the molecular hydrogen (expressed in M_⊙) in distant galaxies is usually calculated as $M_{{\mathrm{H}}_2}={\alpha }_{\mathrm{CO}}{L}_{1-0}^{\prime }$, where $L_{1-0}^{\prime }$ is the luminosity of the CO(1–0) transition and α_CO is the assumed conversion factor in units (omitted in the following) of M_⊙ (K km s⁻¹ pc²)⁻¹. Carilli & Walter (2013) suggest a standard value 0.8 (Solomon & Vanden Bout 2005) usually assumed for the star forming nuclei of classical SMGs, and ∼4 for normal star forming galaxies (SFGs), accordingly to what is found for the giant molecular clouds (GMCs) in the Milky Way by Bolatto et al. (2013). The low value assumed for SMGs would imply that there is more CO emission per unit molecular gas mass. Papadopoulos et al. (2012) suggest that this is likely due to the high SFR that induces heating processes, resulting in an increase of the temperature and excitation level of the molecular gas with respect to normal galaxies. However, the average excitation level in SMGs may be less extreme than in QSOs (Weiß et al. 2007), likely due to the fact that in some SMGs the star formation occurs on a larger scale than in compact QSO hosts. Recent works based on the dust and CO comparison showed how a Galactic α_CO may hold also for sources up to almost an order of magnitude above the main sequence (Genzel et al. 2015; Tacconi et al. 2018).

As for the galaxies in our sample, they are characterized by extreme compactness (i.e., density), a SFR of several hundreds (up to 1600 for XID 42) M_⊙ yr⁻¹, and powerful nuclear activity (L_{2 − 10 KeV} up to tens 10⁴⁴ erg s⁻¹) that contributes up ∼15% of the total FIR (8–1000 μm) luminosity (Circosta et al. 2019), that likely boosts the excitation level of the gas, implying a low value of α_CO. Magdis et al. (2012) derived an estimate of α_CO for a few samples of both (local and distant) main sequence and starburst galaxies: they derived the molecular mass exploiting the observed M_gas/M_dust ∝ Z relation (where Z is the metallicity), and then obtained α_CO measuring $L_{1-0}^{\prime }$. For local ULIRGs and high-z SMGs they found α_CO ∼ 1. Using the Magdis et al. (2012) relation, for one of our sources (XID 539, which features nearly solar metallicity as shown by Nagao et al. 2012) Gilli et al. (2014) constrained the α_CO to a value of 0.8–1 from the measured M_dust (which is consistent within 1σ with M_dust derived in this work), in agreement with the standard α_CO = 0.8 suggested by Carilli & Walter (2013). We note that α_CO ≈ 1 is also adopted in all past works on our targets: Coppin et al. (2010) and Gilli et al. (2014) for XID 403, Calistro Rivera et al. (2018) for XID 42, and Circosta et al. (2019) for the entire sample. In particular, we aim to compare the gas masses and column densities with those derived by Circosta et al. (2019, assumed α_CO = 0.8). In addition, Talia et al. (2018) derived the gas mass of the GMASS 0953 from the CO(3–2) and CO(6–5) transitions, adopting α_CO = 0.8 (motivated by the compactness and the high SFR of the object), finding a good agreement with the gas mass derived from the [C I] emission; we note that GMASS 0953 is a z ∼ 2.2 SFG very similar to those in our sample, hosting a heavily obscured AGN (log N_H >  24) and with gas size ∼1 kpc, M_* ∼ 10¹¹ M_⊙ and SFR∼210 M_⊙ yr⁻¹. Given the aforementioned premises, we adopted α_CO = 0.8 to obtain the molecular masses used in Sect. 5.5 to derive the column densities. Then, because of the existing uncertainties in the CO luminosity to gas mass conversion factors, we also calculate the molecular gas masses assuming a Galactic α_CO = 6.5 (as in Scoville et al. 2016). This means that both the masses and column densities simply increase by a factor of 6.5/0.8 ≈ 8.1, since all the relations are linear. Assuming α_CO> 0.8 then strengthens the hypothesis that the ISM in the host galaxy contributes significantly to the AGN obscuration.

We calculate the gas masses from both the line and continuum emission, to test whether they provide consistent results and address the large uncertainty on the involved conversion factors. We firstly adopt α_CO = 0.8 to obtain M_H2 from $L_{1-0}^{\prime }$; we derive $L_{1-0}^{\prime }$ from the detected high-J CO lines using different CO Spectral Line Energy Distribution (SLED). As for the continuum emission, we derive L_850 μm (rest-frame luminosity at 850 μm, that is ≈353 GHz) from the observed flux density S_obs (at λ ≈ 2.1 mm) adopting a MBB model, and then we exploit the empirical $L_{1-0}^{\prime }\propto L_{850\mu \mathrm{m}}$ relation from Scoville et al. (2016) to obtain the CO(1–0) luminosity. We note that this is a pure luminosity relation that does not imply any assumption on the α_CO factor. Similar trends between (both monochromatic and integrated) FIR and CO luminosity are often found in the literature, showing that they are either nonlinear (e.g., Saintonge et al. 2018) or even point to different relations for SFGs and ULIRGs/SMGs (Daddi et al. 2010; Genzel et al. 2010; see also Carilli & Walter 2013). These results suggest that assuming a single constant $L_{\mathrm{CO}}^{\prime }/L_{\mathrm{FIR}}$ ratio might be too simplistic. However, all these relations show a large scatter (at least ∼1 dex for the FIR luminosity); for simplicity, we then assume the single functional form $L_{1-0}^{\prime }\propto L_{850\mu \mathrm{m}}$ from Scoville et al. (2016), also to compare our results with those of Circosta et al. (2019), who have assumed the same relation. Scoville et al. (2016) then adopt α_CO  =  6.5 to calibrate directly the M_H2 ∝ L_850 μm relation exploiting both the $L_{1-0}^{\prime }$ measurements and L_850 μm derived from SPIRE and SCUBA fluxes, finding a ratio α_850 μm  =  L_850 μm/M_H2² = $6.7\times {10}^{19}{M}_{\odot }^{-1}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{Hz}}^{-1}$. They chose a single value for all galaxies in their samples (i.e., SFGs, local ULIRGs and SMGs) arguing that a lower factor would be inappropriate for a globally distributed molecular gas. However, they point out that a α_CO value several times smaller is often assumed for ULIRGs and SMGs, and that such a low α_CO would be more appropriate to describe the very nuclear regions. In addition, in Kaasinen et al. (2019), where the Scoville et al. (2016) relation is exploited to test the use of long-wavelength dust emission as a gas mass tracer, the authors conclude that the main contribution to the uncertainty in the derived molecular masses (from either CO(1–0) or dust continuum) remains the α_CO factor. Given this large uncertainty, we additionally derive the gas masses assuming the same Galactic α_CO of Scoville et al. (2016), in order to provide a range of masses that fairly takes this uncertainty into account. We calculate the gas masses from the detected lines assuming α_CO  =  6.5 and then from the continuum using α_850 μm, and finally compare them.

We firstly estimate the molecular mass from the luminosity of the detected lines. We calculate the line flux F_CO as the flux density (derived from the uv-plane fitting, Table 4) times the line FWHM reported in Table 3. For XID 666 we found a double-peaked line, thus we calculate the flux multiplying the flux density by the sum of the two line widths derived from the two Gaussian components of the fit, which is 636 ± 72 km s⁻¹. For the undetected sources we consider an upper limit to the flux density equal to three times the rms measured in the peak channel and as line width we assumed a typical value of 400 km s⁻¹ (e.g., Carilli & Walter 2013), consistent within 1σ with the two single-peaked detected lines. We report F_CO in Col. 2 of Table 7. We estimate the CO luminosity $L_{\mathrm{CO}}^{\prime }$ in K km s⁻¹ pc² from F_CO, given the redshift, by means of the well known relation reported in Solomon et al. (1992): $L_{\mathrm{CO}}^{\prime }=3.25·{10}^7{\nu }_{\mathrm{obs}}^{-2}{(1+z)}^{-3}{D}_{\mathrm{L}}^2{F}_{\mathrm{CO}}$ where ν_obs is the observed frequency (in GHz) and D_L is the luminosity distance (in Mpc).

The luminosity of the observed transition has been converted into $L_{1-0}^{\prime }$ through the line ratio derived from a CO-SLED. This ratio depends on the excitation level of the molecular gas, that can significantly differ for the various astronomical objects (e.g., Carilli & Walter 2013). Historically, many sources have been classified as QSOs while others as SMGs, mainly because of the different selection techniques used; however, there is a significant overlap in their properties (Carilli & Walter 2013). Despite the main sequence nature of the galaxies in our sample, recent works (Pozzi et al. 2017; Mingozzi et al. 2018) have shown how the AGN contribution to the high-J CO transitions is not negligible even in low-luminosity QSOs (L_{2 − 10 keV} ∼ 10^{42 − 43} erg s⁻¹, for comparison our sources have L_{2 − 10 keV} ∼ 10⁴⁴ erg s⁻¹). In addition, in compact high-z galaxies the AGN boost to the high-J transitions can affect the whole source (Vallini et al. 2019), making its effect more important than in lower redshift objects. We also stress that the MIR (5–15 μm) emission of our sources is by far (> 50%) dominated by the AGN component (as shown in the SED decomposition by Circosta et al. 2019). Thus, we argue that a main-sequence CO-SLED could not correctly fit the physical properties of the gas in these objects, due to their compactness and high AGN luminosity; however, Daddi et al. (2015) showed that, in a sample of z = 1.5 normal SFGs, the high-J emission (at least from the CO(5–4) transition) is likely more similar to that of SMGs, rather than local SFGs (while for low-J transitions – up to CO(3–2) – they show a MW-like CO-SLED). At any rate, we point out that assuming a lower excitation CO ladder that SMG would result in an increase of the gas mass and, as a consequence, in a higher contribution of the host galaxy ISM to the AGN obscuration. Given these premises, we use two different AGN CO-SLEDs suitable for QSOs and SMGs, in order to encompass a large range of possible values. This, coupled with the additional continuum-based estimate, partly mediates the large uncertainties in measuring the ISM masses from high-J CO transitions, primarily due to the assumed excitation ladders. We use ratios between the observed and CO(1–0) lines reported in Table 2 of Carilli & Walter (2013), who derived averaged values for different selected object types (including QSO and SMG) based on all the available literature at the time. The CO(7–6)/CO(1–0) ratio (needed for XID 539) and CO(6–5)/CO(1–0) ratio (needed for XID 337 and XID 551) are not reported, so we have extrapolated them from the median of the CO-SLEDs presented in the same work. They are equal to 0.16 and 0.48 (for SMG and QSO, respectively) for the CO(7–6)/CO(1–0) ratio, and equal to 0.34 and 0.64 for the CO(6–5)/CO(1–0) ratio. These values for both the transitions are in good agreement (within 1σ) with those reported by Kirkpatrick et al. (2019, see Table 4) for AGN-dominated galaxies. We assume in the following the same designation of Carilli & Walter (2013) for the two adopted conversion factors, namely SMG CO-SLED and QSO CO-SLED.

We also estimated the molecular gas mass from the continuum luminosity. We firstly exploit the empirical relation between the luminosity of the CO(1–0) transition $L_{1-0}^{\prime }$ and L_850 μm (in units of erg s⁻¹ Hz⁻¹) shown in the left panel of Fig. 1 in Scoville et al. (2016), and valid for normal star forming galaxies, local ULIRGs, and SMGs: $L_{1-0}^{\prime }=3.02\times {10}^{-21}{L}_{850\mu \mathrm{m}}\mathrm{K}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^2$. Given the luminosity distance D_L of a source, L_850 μm can be derived from the red-shifted flux density at 850 μm $S_{\frac{{\nu }_{850\mu \mathrm{m}}}{1+z}}$ as $4\pi D_{\mathrm{L}}^2{S}_{\frac{{\nu }_{850\mu \mathrm{m}}}{1+z}}{(1+\mathrm{z})}^{-1}$. Thus, we firstly need to rescale the measured flux densities S_obs to $S_{\frac{{\nu }_{850\mu \mathrm{m}}}{1+z}}$. We assume a MBB law in the optically thin regime S_ν ∝ B_ν(T_d)(1 − e^−τ)≈B_ν(T_d)τ_ν, where τ_ν = (ν/ν₀)^β is the dust optical depth, with β = 2 and ν₀ = 1.5 THz, defined as the frequency at which the dust becomes optically thick (e.g., Conley et al. 2011; Rangwala et al. 2011; Gilli et al. 2014); at this wavelength the Planck function can be approximated by the Rayleigh–Jeans law B_ν, RJ(T_d) = 2k_BT_dν²/c², thus S_ν ∝ ν^β + 2. Hence, from the ratio between $S_{\frac{{\nu }_{850\mu \mathrm{m}}}{1+z}}\propto {[353\mathrm{GHz}/(1+z)]}^{\beta +2}$ and $S_{\mathrm{obs}}\propto {\nu }_{\mathrm{obs}}^{\beta +2}$ we derive $S_{\frac{{\nu }_{850\mu \mathrm{m}}}{1+z}}=S_{\mathrm{obs}}{(353\mathrm{GHz}/{\nu }_{\mathrm{obs}})}^{\beta +2}{(1+\mathrm{z})}^{-(\beta +2)}$. These flux densities are reported in the fourth column of Table 7. We compared these values with those obtained from both the template (performed by Circosta et al. 2019) and MBB (β = 2, Sect. 5.2 of this work) SED fitting. In general, the fluxes from the template fitting are a factor ∼1.7 (XID 539)–4.9 (XID 746) higher than those derived from the MBB SED fitting. This could be attributed to the intrinsic differences between the models, to the dependence of the flux density on the RJ slope (i.e., the assumed β) and to the uncertainties due to the lack of data in the RJ tail. In this respect, the MBB SED fitting in Sect. 5.2 does include the additional points (or upper limits) at 2.1 mm of this work. We found that the MBB SED fitting derived flux densities are in excellent agreement (1σ) we those derived in this section and reported in Table 7.

Finally, we can write the 850 μm luminosity as:

$$\begin{array}{c}\hfill L_{850\mu \mathrm{m}}=4\pi D_{\mathrm{L}}^2{\left(\frac{353\mathrm{GHz}}{{\nu }_{\mathrm{obs}}}\right)}^{\beta +2}{S}_{\mathrm{obs}}{(1+z)}^{-(\beta +3)}.\end{array}$$

We note that this luminosity, as well as the derived gas masses and column densities, does not depend on the dust temperature, due to the RJ regime assumption. In order to obtain the total gas mass we also need to account for the atomic hydrogen mass M_HI. We assume a typical ratio M_H2/M_HI ≈ 5 as found by Calura et al. (2014) on the basis of IR observations of a large sample of high-redshift QSOs (z >  5.7).

In Table 7 we report the luminosities of the observed transitions (Col. 3), the CO(1–0) luminosities (Col. 6) and gas masses (Col. 7) derived from both the CO-SLEDs and continuum (as indicated in Col. 5). In Fig. 4 we show the gas masses derived from the detected CO lines, assuming the two different CO-SLEDs, vs. the gas masses derived from the dust continuum emission. We also report the masses obtained by Circosta et al. (2019) from the SED fitting performed in the infrared band.

Fig. 4. Gas mass derived from the dust emission (x-axis) vs. gas mass derived from the CO line (y-axis). We used αCO = 0.8 to calculate the molecular gas mass from both the line and continuum emission, which have been previously converted into $L_{1-0}^{\prime }$ assuming different CO-SLEDs (SMG, red filled circles and QSO, blue filled circles) for the line, and the $L_{1-0}^{\prime }\propto L_{850\mu \mathrm{m}}$ empirical relation from Scoville et al. (2016, see left panel of their Fig. 1) for the continuum. We accounted for the atomic mass assuming it is 1/5 of the molecular mass. We also report on the y-axis (black open squares) the gas mass derived from the SED fitting by Circosta et al. (2019). The source IDs are reported in orange. For clarity, we slightly shifted the datapoints of XID 551 toward the left, as they would overlap with those of XID 337 (black shaded symbols). In addition, we report the horizontal error bars of the gas mass derived from the dust emission only on the black open squares. The black dashed line represents the 1:1 relation.

For the detected sources (XID 42, XID 539 and XID 666) the gas masses obtained from the dust emission are fully compatible (within 1σ) with those derived from both the SMG CO-SLED conversion factor and IR-SED fitting (Circosta et al. 2019). Moreover, for XID 42, Calistro Rivera et al. (2018) obtained a gas mass of 5 ± 2 × 10¹⁰ M_⊙ from the CO(3–2) line, which is in agreement with our SMG CO-SLED and continuum measurements (3σ) and SED fitting derived mass (1σ). XID 539 shows all masses compatible within 3σ with the molecular mass of ≃1.6 × 10¹⁰ M_⊙ found by Coppin et al. (2010) on the basis of CO(2–1) transition analysis. The lowest masses are the those derived using the QSO CO-SLED, as a consequence of the higher line ratios (Carilli & Walter 2013), which in turn lead to lower $L_{1-0}^{\prime }$ values. Finally we note that for the undetected sources the gas masses derived by Circosta et al. (2019) from the SED fitting are in agreement (within 1σ) with the upper limits derived in this work.

In order take in account also a Galactic α_CO, we calculate the gas mass (still accounting for the atomic gas mass as 1/5 of the molecular gas mass) from the CO emission assuming α_CO = 6.5, as used by Scoville et al. (2016) to calibrate their relation (this α_CO also includes a factor of 1.36 with respect to the molecular hydrogen gas mass, accounting for the mass of heavier elements). Then, we calculate the molecular gas mass directly from L_850 μm via the α_850 μm factor (we note that β = 1.8 assumed by Scoville et al. (2016) to derive the L_850 μm – and used to calibrate α_850 μm – slightly differs from that assumed in this work, that is β = 2; this difference is negligible considering the large uncertainty on the other assumed conversion factors). The results are shown in Fig. 5. Because of the linear scaling of the gas mass with α_CO at a given $L_{1-0}^{\prime }$, the CO-derived gas masses in Fig. 5 are a factor 6.5/0.8 = 8.1 times higher than those shown in Fig. 4. For the detected sources, we found that the gas mass derived from the continuum is comprised in the range of values encompassed by the gas masses derived assuming the two CO-SLEDs, and in general that they are all compatible within 3σ. The consistency between these measurements is further increased if one considers that an additional ∼25% uncertainty is present in the α_850 μm conversion factor.

Fig. 5. Gas mass derived from the dust emission (x-axis) vs. gas mass derived from the CO line (y-axis). We used αCO = 6.5 to calculate the molecular gas mass from the line, which has been previously converted into $L_{1-0}^{\prime }$ assuming different CO-SLEDs (SMG, red filled circles and QSO, blue filled circles), and α850 μm from Scoville et al. (2016, see right panel of their Fig. 1) to derive the molecular gas mass from the continuum. We accounted for the atomic mass assuming it is 1/5 of the molecular mass. The source IDs are reported in orange. For clarity, we slightly shifted the datapoints of XID 551 toward the left, as they would overlap with those of XID 337 (black shaded symbols). In addition, we report the horizontal error bars of the gas mass derived from the dust emission only on the black open squares. The black dashed line represents the 1:1 relation.

5.4. Notes on the morphology of the sources

XID 42 shows a gradient along one direction in the velocity map, at least with a beam separation, suggesting the presence of a differential kinematic structure. In addition, its velocity dispersion is mostly concentrated in the central core of the galaxy. Calistro Rivera et al. (2018) have inferred its inclination angle (40 ° ±50°) from the axial ratio of the CO(3–2) emission assuming a disk-like geometry, while Rybak et al. (2019) have fit the [C II] emission with a rotating disk model. We investigated the observation of the source carried out by the HST in the V band (GEMS tile ID: 27), corresponding to the rest-frame UV emission. In Fig. 6 we report the r-g-b image in which the red, green and blue channels correspond to the continuum dust emission, CO moment 0 and optical image, respectively. Originally the optical image was shifted by ∼0.2 arcsec from the center of the dust (and CO) emission in the northwest direction, in agreement with the astrometric error between ALMA and HST images shown by Dunlop et al. (2017). Thus, we corrected the HST image astrometry using the Gaia position measurements of the sources in the same field (Gaia Collaboration 2016). After the correction the bulks of the stellar, dust and molecular gas emission coincide, within the beam resolution. The optical image shows a very irregular shape that, coupled with the significant shift of the CO(5–4) line with respect to the optical rest-frame (Fig. 2), confirms the peculiar properties of this source. Therefore, we point out that XID 42 may not be securely associated to a regular rotating disk morphology, and further high-resolution observations are required to draw definitive conclusions on its nature.

Fig. 6. Composed r-g-b image of XID 42: the red, green and blue channels are the dust continuum, CO integrated flux and HST V600 image, respectively. In the bottom-right corner the physical scale is reported.

Indications of a rotating system for XID 666 come from the observed velocity gradient that is clearly distributed along a defined direction (moments 1, Fig. 3), and from the line that shows a double-peaked feature (Fig. 2), with the blue and red-shifted peaks that actually correspond to two different emitting regions. However, the very asymmetric distribution of the two peaks suggests the possibility that XID 666 could actually be a blend of two merging systems. As a support to this interpretation, the moment 2 map (Fig. 3) shows that the velocity dispersion is mostly concentrated in two clumps, a bigger one located in the innermost region of the source, and a smaller one located southeast from the main clump. Even improving the resolution during the imaging (i.e., lowering the robust weighting), we could not unambiguously distinguish the presence of two different sources, and future observations at higher resolution will help to reject or confirm this interesting prospect. For XID 539 it would be critical to assess the presence of a rotating system, due to both the spectral shape (Fig. 2) and low S/N per channel (< 5 in the peak). However, the velocity maps suggest the presence of a differential kinematic structure in one direction, at least with a beam separation between the red and blue-shifted region, and the velocity dispersion is confined in a relatively small region concentrated at the center of the source (moments 1 and 2, Fig. 3). In addition, De Breuck et al. (2014), using ALMA data with a resolution of ∼0.5″, demonstrated that the [C II] emission is well described by a rotating disk model with a similar position angle.

Despite these indications of rotating systems, we could not properly constrain a realistic disk-like model for the sources, due to the large uncertainties on the sizes. Instead, in order to derive the column densities in Sect. 5.5, we consider only a simple uniform spherical geometry. However, we point out that assuming a more compact morphology would simply result in compressing the same amount of gas reservoir in a smaller volume, implying higher densities. Future observations at higher resolution will dramatically reduce the uncertainties on the sizes and in turn on all the size-dependent quantities such as the ISM density.

5.5. Column density associated with the ISM

The obtained masses and sizes reveal that a large reservoir of gas is confined in a few kiloparsec region, in agreement with the literature for SMG in the same redshift range (e.g., Carilli & Walter 2013; Bothwell et al. 2013). We derived the ISM density as M_gas/V, where M_gas is the gas mass derived either from the CO line (using CO-SLEDs) or the continuum emission (using the Scoville et al. (2016) $L_{850\mu \mathrm{m}}\propto L_{1-0}^{\prime }$ relation), assuming α_CO = 0.8 in both the cases, and where V is the volume of a given geometrical model of the source. In order to obtain the number density of the hydrogen atoms, we further divide by the proton mass. Then, we calculated the column density toward the nucleus multiplying the number density by the length of the optical path, from the center of the source toward the galaxy outskirts.

We present a homogeneous spherical model for every source, which is the simplest assumption because it is symmetrical in any direction and independent on the inclination angle. Given the 3σ upper limit to the gas mass for the undetected sources, we derived a 3σ upper limit on the column densities as well. As the source diameter we considered the mean value between the major and minor axis of the Gaussian ellipse fit reported in the top panel of Table 5; as path length we assumed half of the diameter. In Table 7 we report the gas column densities N_HISM (Col. 8) for the spherical model using the gas mass derived with the various methods indicated in Col. 5 of the same table. In Col. 11 we report the X-ray column density N_HX derived by Circosta et al. (2019). The ISM vs. X-ray column densities plot is shown in Fig. 7.

Fig. 7. X-ray derived column densities (x-axis) vs. ISM column densities obtained with a uniform spherical model (y-axis). The values related to the SMG CO-SLED, the QSO CO-SLED, and the continuum emission method are indicated by the red circles, blue circles, and green diamonds, respectively. We also report on the y-axis (black open squares) the column densities derived from the SED fitting by Circosta et al. (2019). The source IDs are reported in orange. The black dashed line represents the NHX = NHISM curve. The error bars correspond to the 1σ significance level, while the upper limits are given at the 3σ level. For clarity the horizontal error bars of the X-ray derived column densities are reported only on the black open squares.

All the column densities derived through the different methods show values similar to N_HX (up to 10²⁴ cm⁻²), suggesting that the host galaxy can provide a significant amount of the column density toward the central AGN. These results are also in agreement with the column densities found by Circosta et al. (2019) from the masses obtained by the IR-SED fitting. We stress that if we consider the gas masses derived from the high-J CO lines and from the continuum emission using α_CO = 6.5 and α_850 μm, respectively, the ISM column density would increase by a factor of 8.1.

The spherical model likely underestimates the source density, since it distributes the mass on a larger volume than models based on more compact and realistic geometries (e.g., disk-like). Given the derived column densities (log N_H∼23–24) for the spherical model, we infer that a more compact morphology can produce even larger column densities, depending on the inclination angle and gas mass profile. This means that the host galaxy ISM may contribute to the nuclear obscuration up to the Compton-thick regime, especially adopting a high value for α_CO.

5.6. Dynamical mass

Given the indications of a rotating structure for the detected sources (Sect. 5.4), we derive the dynamical masses in the idealized case of a simple uniform thin disk geometry, where the diameter corresponds to the major axis a_G and the inclination angle i (defined with respect to the perpendicular to the line of sight, north to south) is derived from the axial ratio. We assume i = 32°, 37°, 27° for XID 42, XID 539 and XID 666, respectively. For XID 42 and XID 539 we derive the dynamical masses M_dyn from the FWHM (Δv) of the CO lines obtained from the spectral fitting reported in Table 3, exploiting the following relation adapted from Wang et al. (2013), which assumes a ratio of 0.75 between the maximum projected circular velocity v_c, proj and Δv:

$$\begin{array}{c}\hfill M_{\mathrm{dyn}}{sin}^{2}i=6.5·{10}^4{\left(\frac{\mathrm{\Delta }v}{{\mathrm{km}\mathrm{s}}^{-1}}\right)}^2\left(\frac{a_G}{\mathrm{kpc}}\right)M_{\odot }.\end{array}$$

For XID 666 we resolve the double peak of the line, thus we directly derive the dynamical mass from v_c, proj = 195 ± 14 km s⁻¹, calculated as the half of the peak separation, assuming that a_G is the distance at which the circular velocity flattens:

$$\begin{array}{c}\hfill M_{\mathrm{dyn}}{sin}^{2}i=1.16·{10}^5{\left(\frac{v_{\mathrm{c},\mathrm{proj}}}{{\mathrm{km}\mathrm{s}}^{-1}}\right)}^2\left(\frac{a_G}{\mathrm{kpc}}\right)M_{\odot }.\end{array}$$

We found $M_{\mathrm{dyn}}=3.7_{-0.7}^{+0.8}\times {10}^{10}{M}_{\odot }$, $2.9_{-1.5}^{+1.8}\times {10}^{10}{M}_{\odot }$, $0.8_{-0.3}^{+0.8}\times {10}^{10}{M}_{\odot }$ for XID 42, XID 539 and XID 666, respectively. These values are about one order of magnitude lower than the values previously found for SMG at z ∼ 2 based on the study of resolved and unresolved CO kinematics (Tacconi et al. 2008; Swinbank et al. 2010; Bothwell et al. 2013) and than the baryonic mass calculated as M_bar = M_* + M_H2 + M_HI, where M_* are the stellar masses obtained by the SED fitting (Circosta et al. 2019) which are equal to M_* = 2.16 ± 0.65 × 10¹¹ M_⊙, 2.15 ± 0.64 × 10¹¹ M_⊙, 4.41 ± 1.32 × 10¹¹ M_⊙ for XID 42, XID 539 and XID 666, respectively. We stress that, for XID 539, De Breuck et al. (2014) and Coppin et al. (2010) found M_dyn in agreement with the value derived in this work. De Breuck et al. (2014) have assumed different inclination angles (29°, 50° and 53°), while in Coppin et al. (2010) an angle i = 30° was assumed, as appropriate for randomly inclined disks in a sample of galaxies. Such discrepancy could be attributed either to an overestimate of the stellar mass by the SED fitting or to an underestimate of Δv and a_G, given the sensitivity limitation of the observations. As for XID 666, we assumed that v_c, proj is equal to the half of the distance of the peaks and that the disk diameter is the distance at which the velocity profile flattens, but due to the large uncertainty on the source geometry and size it is not possible to securely determine this distance. It is also likely that the stellar mass has been derived considering a different emitting area with respect to the molecular gas. Finally we note that the 30% error on the stellar masses reported in Circosta et al. (2019) is indicative only of the statistical error given the assumed Salpeter Initial Mass Function (IMF, Salpeter 1955), and the systematic error due to different choices of IMF are likely much larger. By assuming a Chabrier IMF (Chabrier 2003) the resulting stellar masses would be a factor ∼1.7 smaller, leading to a lower discrepancy in the M_bar/M_dyn ratio.

6. Summary and conclusions

Exploiting ALMA data we analyzed the high-J CO and the dust emission of a sample of six star forming galaxies hosting a heavily obscured (log N_HX >  23) AGN in the redshift range ∼2.5 <  z <  4.7 in the CDF-S. We detect both the line and continuum emissions for three sources, for which we extracted the flux densities and sizes in order to characterize their masses, gas densities and column densities, and evaluate to what extent the host galaxy ISM can contribute to the central AGN obscuration. Our results are summarized as follows.

– We fit the CO line and continuum emission in both the brightness and visibility spaces with a 2D Gaussian model, finding that all sources are extremely compact objects. The CO emission, tracing the molecular gas component, is contained in regions with major axis in the range 2.1–3 kpc. The continuum emission, tracing the dust component, is produced in somewhat smaller regions, with major axes that are ∼10–30% smaller than those derived from the CO emission. This difference in size can be explained by a radial decrease of the gas column density and dust temperature (Calistro Rivera et al. 2018), but we cannot exclude that it could be attributed to the sensitivity limit of the observations.

– We analyzed the spectra of the detected sources and performed a Gaussian fit to the lines. For one of these sources (XID 666) we found a double-peaked feature that we ascribe to a rotating disk-like morphology. However, given the strong asymmetry of the line shape (Fig. 2, bottom panel) and the dispersion velocity mainly distributed in two separated clumps, we cannot exclude that XID 666 is actually an unresolved merger of two different systems. XID 42 also shows indication of a rotating structure, even if it cannot be securely associated to a regular disk shape: the line is red-shifted by ∼500 km s⁻¹ with respect to the optical rest-frame velocity and the HST image (V-band, Fig. 6), that traces the UV rest-frame stellar emission, shows a very irregular morphology. These features make XID 42 the most peculiar target among the sample, and observations at higher resolution (∼0.1″) would be necessary to cover the source with a sufficient (∼5–6) number of beams, and unveil whether its kinematics points to relaxed system or to a chaotic, merging, structure.

– We calculated the gas mass as M_H2 + M_HI, assuming M_HI = M_H2/5 (Calura et al. 2014). The molecular gas mass was firstly derived from the high-J CO detected lines assuming two different CO-SLEDs typical for SMGs and QSOs to obtain $L_{1-0}^{\prime }$, and then assuming a $L_{1-0}^{\prime }$–M_H2 conversion factor α_CO = 0.8. We obtained masses in the range (0.5 − 2.7)×10¹⁰ M_⊙, in agreement with those derived by Circosta et al. (2019) from the SED fitting. Secondly, we derived the molecular gas masses from the continuum emission on the basis of the $L_{1-0}^{\prime }$–L_850 μm empirical relation found by Scoville et al. (2016), using a modified black-body model to calculate L_850 μm and then again α_CO = 0.8 to obtain the gas masses. We found that these masses are fully in agreement (within 1σ) with those derived from the CO line emission. In order to encompass the uncertainty on α_CO, we also calculated the gas masses from the CO lines using α_CO = 6.5, and then compared these values with those obtained from the direct Scoville et al. (2016) relation M_H2 and L_850 μm, which is based on the assumption of the same α_CO, finding that they are in good agreement. This corresponds to an increase of the gas masses (a factor of ∼8.1) and, as a consequence, in a more pronounced contribution of the ISM to the nuclear obscuration.

– We investigated the contribution of the host galaxy to the obscuration of the AGN by comparing the column density N_H, ISM derived from the ISM masses with that (N_H, X) derived from the X-ray spectral fitting (Circosta et al. 2019). As for the source geometry, we assumed a uniform spherical model for all the sources. We found that N_H, ISM is similar to N_H, X (∼10^{23 − 24} cm⁻²), and conclude that the host galaxy can significantly affect the obscuration of the central engine. All the N_H, ISM are also comparable to the column densities derived by Circosta et al. (2019) from the IR-SED fitting. We point out that despite the large number of assumptions, the N_H, ISM values derived from three independent measurements (i.e., CO line, dust continuum and SED fitting), are in agreement with each other and with N_H, X. This result strongly supports a scenario in which a substantial fraction of the obscuration may be produced on large (kpc) scales, at least for these FIR-bright objects characterized by a large gas reservoir confined in a very compact volume. Since the spherical model likely underestimates the ISM density, we point out that a more realistic (disk-like) geometry, compressing the gas mass in a smaller volume, would result in larger density; depending on the inclination angle of the disk and the gas distribution, this may produce column densities that can significantly affect the obscuration even in heavily Compton-thick AGN (≳10²⁴ cm⁻²).

– From the continuum emission we derived the dust mass, assuming a modified black body law in the optically thin regime. We found that the dust masses are in the range ∼(0.9 − 4.9)×10⁸ M_⊙.

– We derived the dynamical masses for the three detected sources, whose values are in the range (3.7 − 0.8)×10¹⁰ M_⊙. These masses are one order of magnitude lower than the sum of the stellar (Circosta et al. 2019) and gas masses, and consistent with those derived in previous work for XID 539, assuming different inclination angles (Coppin et al. 2010; De Breuck et al. 2014). Such discrepancy could be attributed to an overestimate of the stellar mass by the SED fitting or to an underestimate of the FWHM of the line and the gas component size (given the sensitivity limitation of the observations). Moreover, the stellar mass may be derived considering a different emitting area with respect to the molecular gas. Finally, we note that the scatter of the stellar masses derived by Circosta et al. (2019) does not take into account different assumptions for the IMF. Indeed, by assuming a Chabrier IMF the stellar masses reduce of a factor ∼1.7, leading to a smaller discrepancy between M_bar and M_dyn.

From the analysis of our sample, we conclude that SMGs environment is characterized by an extremely dense ISM confined in few kiloparsec. In such conditions, the ISM can produce column densities that contribute to the obscuration of the central AGN up to the Compton-thick regime. In addition, the sources’ velocity gradient indicates the presence of a rotating system, even though some of their morphological features suggest a more complex structure, possibly associated with mergers.

Acknowledgments

QD acknowledges A. Peca and G. Sabatini for all the support during the long journey, and R. Decarli, M. Mignoli and G. Iorio for the very helpful discussions. CC acknowledges support from the IMPRS on Astrophysics at the LMU (Munich). FV acknowledges financial support from CONICYT and CASSACA through the Fourth call for tenders of the CAS-CONICYT Fund and from CONICYT grants Basal-CATA AFB-170002. This paper makes use of the following ALMA data: 2015.1.01205.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.

References

All Tables

All Figures

Fig. 1. Top row: CO-line peak channel of the three detected sources. The observed CO transition and the velocity corresponding to the center of the channel (defined as the velocity shift with respect to the systemic redshift) are reported in the up-left and bottom-right corner of the images, respectively. The blue ellipses are the regions from which the spectra in Fig. 2 were extracted. In the XID 666 CO-line peak image (right panel) we show the contours of the blue-shifted peak channel centered at −206.5 km s−1 (2, 3, 4σ, black solid lines) over the red-shifted peak channel (see Sect. 4 for details). Bottom row: dust-continuum image of the three detected sources. In all the images the black ellipse in the left-bottom corner represents the clean beam.

In the text

Fig. 2. Spectrum of XID 42, XID 539 and XID 666 (green shaded histogram), expressed as the mean brightness within the regions contained in the blue ellipses of Fig. 1. The black solid line represents the Gaussian line fitting. Due to its double-peaked feature, we fit the line of XID 666 with two Gaussian components, indicated in the plot with the blue and red dashed lines. The channel width Δv is reported for each spectrum.

In the text

	Fig. 3. Total flux (moment 0), velocity map (moment 1) and velocity dispersion map (moment 2) (from top to bottom) extracted for XID 42, XID 539 and XID 666 (from left to right). Only the pixels above the 3σ level are shown. In the top right corner of the moment 0 images the angular scale is reported. In all the images the black ellipse in the left-bottom corner represents the clean beam.
In the text

Fig. 4. Gas mass derived from the dust emission (x-axis) vs. gas mass derived from the CO line (y-axis). We used αCO = 0.8 to calculate the molecular gas mass from both the line and continuum emission, which have been previously converted into $L_{1-0}^{\prime }$ assuming different CO-SLEDs (SMG, red filled circles and QSO, blue filled circles) for the line, and the $L_{1-0}^{\prime }\propto L_{850\mu \mathrm{m}}$ empirical relation from Scoville et al. (2016, see left panel of their Fig. 1) for the continuum. We accounted for the atomic mass assuming it is 1/5 of the molecular mass. We also report on the y-axis (black open squares) the gas mass derived from the SED fitting by Circosta et al. (2019). The source IDs are reported in orange. For clarity, we slightly shifted the datapoints of XID 551 toward the left, as they would overlap with those of XID 337 (black shaded symbols). In addition, we report the horizontal error bars of the gas mass derived from the dust emission only on the black open squares. The black dashed line represents the 1:1 relation.

In the text

Fig. 5. Gas mass derived from the dust emission (x-axis) vs. gas mass derived from the CO line (y-axis). We used αCO = 6.5 to calculate the molecular gas mass from the line, which has been previously converted into $L_{1-0}^{\prime }$ assuming different CO-SLEDs (SMG, red filled circles and QSO, blue filled circles), and α850 μm from Scoville et al. (2016, see right panel of their Fig. 1) to derive the molecular gas mass from the continuum. We accounted for the atomic mass assuming it is 1/5 of the molecular mass. The source IDs are reported in orange. For clarity, we slightly shifted the datapoints of XID 551 toward the left, as they would overlap with those of XID 337 (black shaded symbols). In addition, we report the horizontal error bars of the gas mass derived from the dust emission only on the black open squares. The black dashed line represents the 1:1 relation.

In the text

	Fig. 6. Composed r-g-b image of XID 42: the red, green and blue channels are the dust continuum, CO integrated flux and HST V600 image, respectively. In the bottom-right corner the physical scale is reported.
In the text

Fig. 7. X-ray derived column densities (x-axis) vs. ISM column densities obtained with a uniform spherical model (y-axis). The values related to the SMG CO-SLED, the QSO CO-SLED, and the continuum emission method are indicated by the red circles, blue circles, and green diamonds, respectively. We also report on the y-axis (black open squares) the column densities derived from the SED fitting by Circosta et al. (2019). The source IDs are reported in orange. The black dashed line represents the NHX = NHISM curve. The error bars correspond to the 1σ significance level, while the upper limits are given at the 3σ level. For clarity the horizontal error bars of the X-ray derived column densities are reported only on the black open squares.

In the text