Behind the dust veil: A panchromatic view of an optically dark galaxy at z=4.82

XS55 was originally selected in the COSMOS-XS (Algera et al., 2020; van der Vlugt et al., 2021) catalog with an ID=55, hence we dub it XS55. It is detected in the ultra-deep COSMOS-XS S-band image with $S_{\rm 3GHz}$ = $6.35\pm 0.96\,{\rm\mu Jy}$, but undetected in the less deep COSMOS 3 GHz map (${\rm rms=2.5\,\mu Jy}$, Smolčić et al., 2017), and detected in the MeerKAT image with $S_{{\rm 1.3GHz}}=10.9\pm 2.1\,{\rm\mu Jy}$ (Jarvis et al., 2016; Heywood et al., 2022; Hale et al., 2024). XS55 has no optical counterpart (i.e., optically dark) and is not included in the COSMOS2020 catalogue (Weaver et al., 2022). It drops out in JWST F115W and F150W images, but is detected in IRAC 4.5$\mu$m (Fig. 1) and tentatively detected ($\sim 3\sigma$) in the ALMA 2 mm MORA map (Casey et al. 2021). By performing the super-deblending technique (Jin et al., 2018; Liu et al., 2018) with the radio prior, we measure the deblended Herschel and SCUBA-2 photometry of XS55. Interestingly, it is not detected in Herschel images ($3\sigma$ limiting depths: 250 ${\rm\mu}$m=$5.3\,{\rm mJy}$,350 ${\rm\mu}$m=$8.0\,{\rm mJy}$, 500 ${\rm\mu}$m=$8.7\,{\rm mJy}$, Jin et al. 2018), but is well detected in two SCUBA-2 bands; 450$\mu$m (Gao et al., 2024) and 850$\mu$m (Simpson et al., 2019) with a red 450$\mu$m/850$\mu$m colour ($S_{\rm 450\,\mu m}=5.6\pm 1.4\,{\rm mJy},S_{\rm 850\,\mu m}=5.7\pm 0.8\,{\rm mJy}$). Assuming typical dust templates from Magdis et al. (2012), the red SCUBA-2 colour suggests a FIR photometric redshift of $z>6$. Consequently, XS55 was followed up by two ALMA 3 mm line scan projects in Cycle 9 (ID: 2022.1.00884, PI: R. Gobat; ID: 2022.1.00863.S, PI: J. Hodge).

The two ALMA programmes were observed for a total of 1.9 hours on source. The frequency setups are identical in the two programs, adopting the same setups as in Jin et al. (2019) and covering 84–108 GHz with three tunings. We produce measurement sets with the Common Astronomy Software Applications (CASA, McMullin et al. 2007) pipeline for each observation programme. Following the methods presented in Jin et al. (2022), Zhou et al. (2024), and Sillassen et al. (2024), the calibrated data are converted to $uv$ table format and analysed with the GILDAS software package in $uv$ space. To enhance the signals, the $uv$ tables from the two programs are combined using the task $uv\_merge$. The final products reach a continuum sensitivity of 9.6$\,{\rm\mu Jy/beam}$ with a spatial resolution of $1\aas@@fstack{\prime\prime}6$, and a line sensitivity of 12$\,{\rm mJy/beam}$ over 500 km/s at $\sim 99\,{\rm GHz}$. As shown in see Figs. 2 and 3.1, We robustly detect continuum ($\sim 14\sigma$) and two lines ($\sim 10\sigma$ and $\sim 4\sigma$). The dust continuum is well-fitted by a point-source model using GILDAS uvfit, while fitting with an elliptical Gaussian does not yield useful constraints. Therefore, the dust continuum is unresolved, and we place an upper limit on the continuum size using Eq. (2) from Gómez-Guijarro et al. (2022) (Table 2).

XS55 was observed with JWST NIRCam in F115W, F150W, F277W and F444W bands, as a part of the COSMOS-Web survey (Casey et al., 2023). We use the image product versions from the Dawn JWST Archive (DJA¹¹1https://dawn-cph.github.io/dja/index.html, Valentino et al., 2023), and further, have verified these are consistent with the COSMOS-Web team’s map (Shuntov et al. in prep.). XS55 is well detected in both F277W ($\sim 10\sigma$) and F444W ($\sim 36\sigma$), but not detected in both F115W and F150W ($<3\sigma$), consistent with the H-dropout selection from Wang et al. (2019).

The COSMOS field has been fully observed in soft (0.5-2.0 keV) and hard (2.0-10 keV) X-rays with both XMM-Newton (50 ks per pointing, PI: G. Hasinger; Hasinger et al., 2007) and Chandra ($\sim$180 ks exposure) as part of the Chandra COSMOS (C-COSMOS, PI: M. Elvis; Elvis et al., 2009) and Chandra COSMOS Legacy (PI: F. Civano; Civano et al., 2016) surveys. As shown in Figs. 3 and 1, XS55 is tentatively detected in soft X-ray band of Chandra with $2.1\sigma$ significance, and detected in the stacked soft, medium and hard X-ray XMM-Newton images with $3.1\sigma$.

After extracting the ALMA spectra in all spectral windows (SPWs), we combine them in a single 1D spectrum (Fig. 2), and run a line-searching algorithm as in Jin et al. (2019) to search for emission line features with the highest significance. The continuum is fit with a power-law with fixed slope of 3.7 in frequency (assuming $\beta\approx 1.7$, Magdis et al., 2012; Sillassen et al., 2024), masking out the channels of significant emission lines. The detected emission lines are fitted with a Gaussian profile in the continuum subtracted spectrum. XS55 is detected in 3mm continuum at $\sim 14\sigma$ with a flux of $134.2\pm 9.3\,{\rm\mu Jy}$ at $\sim 96\,{\rm GHz}$. One line is detected at $98.99\,{\rm GHz}$ at $9.8\sigma$ (see Fig. 2), and we search for other lines in the spectrum consistent with this detection. We find a $3.9\sigma$ detection at $84.58\,{\rm GHz}$. The two lines are consistent with CO(5-4) and [CI](1-0) at $z=4.8214\pm 0.0004$. We note that there is a slight velocity offset between the CO and CI line peaks of $\sim 141\,{\rm km/s}$ (Fig. 2). By defining the [CI] line range using the velocity range of the CO line, the [CI] S/N is $3.3\sigma$, yielding a low $P_{\rm chance}=0.4\%$, where $P_{\rm chance}$ is the chance probability of finding a spurious second line (Jin et al., 2019). Further, we compare the redshift solution with the NIR-SED fitted redshift probability distribution PDF(z) of LePhare in the COSMOS-Web catalog ($z_{\rm phot}=4.73^{+0.52}_{-0.64}$, Shuntov et al. in prep.), and find they are in excellent agreement (Fig. 4). For a sanity check, we tested the redshift solution of $z=3.66$ in the case of the bright line being CO(4-3), as this redshift is seemingly consistent with the second peak of the NIR PDF(z). However, we found $z=3.66$ is very unlikely, because (1) the $z=3.66$ [CI](1-0) is not detected at the expected frequency 105.6 GHz ($<1\sigma$); (2) the SED fitting at $z=3.66$ yields an abnormally high dust mass to stellar mass fraction $0.09<{M_{\rm dust}/M_{\ast}}<0.13$ that is $>10\times$ above typical values (e.g., Donevski et al., 2020), again disfavouring the $z=3.66$ solution. Therefore, the multitude of evidence confirms the redshift of XS55 to be $z=4.8214\pm 0.0004$.

We model the F277W and F444W morphology of XS55 using Galfit (Peng et al., 2010). We adopt three separate components in both F277W and F444W: a compact central emission, a diffuse central emission, and a companion component at the south-east. The NIRCam PSFs were obtained using Webbpsf (Perrin et al., 2014) with a pixel scale of $0.05{\rm\aas@@fstack{\prime\prime}/pix}$. We consider two cases; one case using a point source (compact component) together with two Sérsic profiles (diffuse component and companion), and another case with three separate Sérsic profiles. In the case using a point source and two Sérsic models in F444W, there is a ring of emission left in the residual (Fig. 6-middle-right), suggesting the brightest part of the galaxy is marginally resolved, or the PSF modelling is imperfect. In the case of three Sérsic profiles, there are no clear structures, as would be expected with random noise (Fig. 6-top). In the three Sérsic profile model of F444W, the compact component has a size of $0\aas@@fstack{\prime\prime}061\pm 0\aas@@fstack{\prime\prime}007$ corresponding to $R_{\rm e}=0.40\pm 0.05\,{\rm kpc}$ with Sérsic index $n<1.2$, while the diffuse emission and the companion have sizes and Sérsic indices of $R_{\rm e}=1.7\pm 0.1\,{\rm kpc}$, $n=0.37\pm 0.08$ and $R_{\rm e}=1.1\pm 0.1\,{\rm kpc}$, $n=0.61\pm 0.33$ respectively (Table 3). In F444W, $33\pm 3$% of the total flux is coming from the compact component, while $53\pm 3$% is coming from the extended diffuse component, and $14\pm 3$% is coming from the companion. In JWST/F277W, the fit of the compact component yields a bulge-like $n=2.6\pm 0.8$ model with $R_{\rm e}=0.58\pm 0.07\,{\rm kpc}$ (Table 3), providing $87\pm 2\%$ of the total flux. The diffuse component cannot be fit, and the fit to the companion yields a $n=2.0\pm 1.4$ model with $R_{\rm e}=1.7\pm 0.5\,{\rm kpc}$, providing $13\pm 2\%$ of the total flux (Sections 3.2 and 3). At this redshift, [OIII]${\lambda\lambda 4959,5007}$ and ${\rm H\beta}$ fall within the JWST/F277W filter, consequently the flux is possibly boosted by line emission (Fig. 4).

To estimate the dust mass ($M_{\rm dust}$) and temperature ($T_{\rm dust}$), we fit the FIR and (sub)mm photometry of XS55 with a modified black-body (MBB) model (Magdis et al., 2012) using the code mercurius (Witstok et al., 2022). We explore two cases; one assuming optically thin dust, and the other assuming ‘self-consistent’ optically thick dust. For the optically thick case, we place an upper limit on the emitting area based on the half light radius of the galaxy in JWST/F444W ($9.08\,{\rm kpc^{2}}$; Section 3.2). The resulting fit and corresponding parameters are shown in Figs. 7 and 2. In both cases, the IR spectral index $\beta_{\rm IR}$ is consistent with a weighted average of $\beta_{\rm IR}=2.02\,{\pm\,0.12}$. For the optically thin case, we recover a low dust temperature of $T_{\rm dust}=27.6^{+3.1}_{-2.7}\,{\rm K}$, with an accompanying high dust mass of $M_{\rm dust}=2.4^{+1.4}_{-0.9}\times 10^{9}\,{\rm M_{\odot}}$. On the other hand, the optically thick dust yields a higher $T_{\rm dust}=32.5^{+2.2}_{-1.9}\,{\rm K}$ and a lower $M_{\rm dust}=1.7^{+0.5}_{-0.4}\times 10^{9}\,{\rm M_{\odot}}$.

To construct a panchromatic SED of XS55, we fit the stellar, AGN, and dust components using STARDUST (Kokorev et al., 2021) and available photometry from optical to radio wavelength, yielding an IR luminosity of $L_{\rm IR}=(5.43\pm 1.77)\times 10^{12}\,{\rm L_{\odot}}$ and a ${\rm SFR}_{\rm IR}=543\pm 177\,{\rm M_{\odot}\,yr^{-1}}$. For the stellar component, we recover the stellar mass $\log(M_{\ast}/{\rm M_{\odot}})=10.7\pm 0.1$, attenuated by $A_{v}=2.2\pm 0.3\,{\rm mag}$. As an additional check, we also fit the stellar part of the spectrum (up to F444W) using Bagpipes²²2https://bagpipes.readthedocs.io/en/latest/ (Carnall et al., 2018), using the same parameters as in Jin et al. (2024), yielding physical properties similar to those from STARDUST, $\log(M_{\rm\ast,bagpipes}/{\rm M_{\odot}})=10.7\pm 0.1$, $A_{\rm v,bagpipes}=2.1\pm 0.3$. The physical parameters calculated with LePhare in the COSMOS-Web catalog (Shuntov et al. in prep.) also agree with our results ($\log(M_{\rm\ast,LePhare}/{\rm M_{\odot}})=10.8\pm 0.3$, $E({\rm B-V})_{\rm LePhare}=1.1$). The resulting panchromatic SED and corresponding properties from STARDUST, and MBB fit from mercurius, are shown in Figs. 4 and 2.

With the derived dust mass, we infer a molecular gas mass ($M_{\rm mol}$) using the standard gas-to-dust mass ratio of star-forming galaxies ($\delta_{\rm gdr}=100$, assuming solar metallicity, Magdis et al., 2012), yielding $M_{\rm mol,thin}=2.6^{+0.5}_{-0.4}\times 10^{11}\,{\rm M_{\odot}}$ and $M_{\rm mol,thick}=1.7^{+0.5}_{-0.4}\times 10^{11}\,{\rm M_{\odot}}$. Furthermore, based on the detection of the [CI](1-0) line (Table 1), we infer the molecular gas mass $M_{\rm mol}/{\alpha_{\rm CI}}=(1.8\pm 0.5)\times 10^{11}/{\alpha_{\rm CI}}\,{\rm M_{\odot}}$, where $\alpha_{\rm C{\tiny I}}=17.0\pm 0.3\,{\rm M_{\odot}\,K^{-1}\,km^{-1}\,s\,pc^{-2}}$ (Dunne et al., 2022) that is consistent with other calibrations (Valentino et al., 2018; Heintz & Watson, 2020). Adopting instead $\alpha_{\rm CI}=4.1\pm 1.4\,{\rm M_{\odot}\,K^{-1}\,km^{-1}\,s\,pc^{-2}}$ from Frias Castillo et al. (2024) for high-$z$ SMGs, it yields a $M_{\rm mol}=(4.3\pm 1.9)\times 10^{10}\,{\rm M_{\odot}}$. Interestingly, the CO(5-4) emission is marginally resolved along the NE-SW direction (see Figs. 2 and 3), as the minor axis of the CO emission as well as the continuum source are unresolved, we obtain upper limits on their sizes using Eq. (2) of Gómez-Guijarro et al. (2022). We measure the size of the integrated CO(5-4) emission by fitting an elliptical Gaussian in $uv$-space with GILDAS/uvfit, and obtain $R_{e}=1.79\pm 1.22\,{\rm kpc}$ with ${\rm PA}=28\pm 20\,{\rm deg}$. We constructed the CO(5-4) moment-1 map using Cube Analysis and Rendering Tool for Astronomy (CARTA; Comrie et al., 2021, Fig. 3-right). In the moment-1 map there is a clear velocity gradient, with the highest redshifted velocity of $70\,{\rm km/s}$ and the highest blueshifted velocity of $-130\,{\rm km/s}$. This velocity gradient is comparable to that of the radio selected NIR-dark galaxies from Gentile et al. (2024b). We discuss the possible interpretations of the extended CO(5-4) emission in Sect. 4.2.

XS55 is tentatively detected in X-ray with both Chandra and XMM-Newton (Fig. 3, 1).Using a 3” radius aperture that corresponds to the mean PSF of the Chandra COSMOS Legacy Survey (Civano et al., 2016), we detect 3 counts in the 0.5-2 keV band of Chandra. This corresponds to a flux of $f_{\rm[0.5-2keV]}=(1.15\pm 0.54)\times 10^{-16}\,{\rm erg\,s^{-1}\,cm^{-2}}$ adopting the conversion rate from Civano et al. (2016). Since XS55 is extremely dusty, this flux should be considered a lower limit, yielding a soft X-ray luminosity of $L_{\rm[0.5-2keV]}>2.81\times 10^{43}\,{\rm erg\,s^{-1}}$. Using the soft X-ray to bolometric luminosity correction from Lusso et al. (2012) we obtained a bolometric luminosity of $L_{\rm bol}>5.35\times 10^{44}\,{\rm erg\,s^{-1}}$. XS55 would be an X-ray selected AGN using the criterion from Riccio et al. (2023) (i.e., $L_{0.2-2.3{\rm keV}}\geq 3\times 10^{42}\,{\rm erg\,s^{-1}}$). Due to the low significance of the X-ray detection, we cannot exclude that the detection could be spurious, or is originating from star-formation. XS55 does not show excess radio emission compared to the infrared radio correlation (IRRC) of Delvecchio et al. (2021), however, this could be due to the large scatter of IR-radio correlation that is largely unconstrained at $z\sim 5$ (Delvecchio et al., 2021). Further, while radio excess is a clear indicator of AGN activity, X-ray AGN are not necessarily radio loud.

To determine whether XS55 is intrinsically cold, or optically thick dust is making it appear cold, we apply three methods of diagnosing optically thick dust from Jin et al. (2022): (1) We compare the molecular gas masses derived from [CI](1-0) emission with gas mass inferred by both thin and thick $M_{\rm dust}$. We find that the lower gas mass from the thick dust model is in good agreement with the [CI]-derived gas mass, assuming $\delta_{\rm gdr}=100$ and $\alpha_{\rm[CI]}=17\pm 3\,{\rm M_{\odot}\,K^{-1}\,km^{-1}\,s\,pc^{-2}}$, while the gas mass from the thin dust model is likely overestimated. (2) We estimate the dust opacity at 100 $\mu$m, $\tau_{\rm 100\mu m}=\kappa\rho R_{\rm e}$, where $\kappa$ is adopted from Jones et al. (2013), and $\rho$ is the volumetric dust density assuming spherical symmetry, using $R_{\rm e}$ from dust continuum ($R_{\rm e}<1.76\,{\rm kpc}$). This yields $\tau_{\rm 100\mu m,thin}>1.6$, and $\tau_{\rm 100\mu m,thick}>1.2$. Furthermore, we calculate the surface SFR density and find $\Sigma_{\rm SFR}>37\,{\rm M_{\odot}\,yr^{-1}\,kpc^{-2}}$. With both $\tau_{\rm 100\mu m}>1$ and $\Sigma_{\rm SFR}>20\,{\rm M_{\odot}\,yr^{-1}\,kpc^{-2}}$, the dust is hence optically thick (Jin et al., 2022). (3) We derive a lower limit of IR luminosity surface density $\Sigma_{\rm IR}>3.7\times 10^{11}\,{\rm L_{\odot}\,kpc^{-2}}$, and place it on the $\Sigma_{\rm IR}-T_{\rm dust}$ diagram. As seen in Fig. 5-right, the thin case is clearly violating the blackbody Stefan-Boltzmann law, and the optically thick case is more favourable. Therefore, these pieces of evidence together imply that the dust of XS55 is optically thick in FIR.

Interestingly, the recovered $T_{\rm dust}$ assuming optically thick dust remains cold. Comparing to the redshift evolution of $T_{\rm dust}$ in main sequence galaxies (Schreiber et al., 2018; Jin et al., 2022), the thick dust solution is 0.13 dex ($\sim 4\sigma$) below the relation (Fig. 5). This makes XS55 one of the coldest DSFGs at $z>4$ discovered to date (e.g., Faisst et al., 2020; Jin et al., 2022; Algera et al., 2024a). Since dust temperature is proportional to the radiation field (i.e., $\langle U\rangle=(T_{\rm d}/18.9)^{6.04}$, Magdis et al., 2012), and the radiation field is proportional to the ratio of star formation efficiency to metallicity (i.e., $\langle U\rangle\propto{\rm SFE}/Z$, Magdis et al., 2012), and given that XS55 has typical SFE of DSFGs, the intrinsically low dust temperature could suggest a high metallicity in XS55, where the cooling is more efficient.

Interestingly, the CO(5-4) is more extended than the dust continuum size and shows a tentative velocity gradient. These properties could suggest a rotating molecular disk (e.g., Rizzo et al., 2023) with compact star-formation (e.g., Cochrane et al., 2019). However, given the low spatial resolution of the CO data, high resolution data are needed to confirm whether it is a rotating disk (e.g., Rowland et al. 2024). As suggested by the complex JWST morphology, galaxy interaction or merger can also be accounted for the extended CO(5-4) emission. Furthermore, given the possible AGN nature of XS55 and that the CO(5-4) is only resolved in one direction, a molecular outflow driven by a central AGN (e.g., Lutz et al., 2020) is also a potential scenario. To disentangle the above scenarios, a high resolution [CII] follow-up would be ideal to reveal the kinematics of XS55.

In nearly all terms, XS55 is a normal but massive main sequence star-forming galaxy at $z=4.8214$, within 2$\sigma$ of the Schreiber et al. (2015) main-sequence relation (${\rm\Delta MS}({\rm SFR/SFR_{MS}})=2.4\pm 0.9$). Its star formation efficiency, and thereby depletion time of $300\pm 100\,{\rm Myr}$, is comparable to other optically faint galaxies at similar redshift (e.g., Jin et al., 2019, 2022). However, with the compact stellar ($R_{e}=1.72\pm 0.13\,{\rm kpc}$) and dust continuum ($R_{e}<1.76\,{\rm kpc}$) sizes, XS55 falls 2.6$\times$ below the mass-size relation of main sequence galaxies (Ward et al., 2024). The compact size, together with a massive amount of dust, could explain the optically faint nature of this source, in agreement with observations from Gómez-Guijarro et al. (2023) and simulations from Cochrane et al. (2024).