ALMA Resolves the First Strongly Lensed Optical/Near-IR-dark Galaxy

The object is part of the low-resolution (≲2'') observations in band 3 (2017.1.01694.S; PI: Oteo), which were aimed at tracing dense molecular gas through J = 4–3 transitions of HCN, HCO+, and HNC molecules. J1135 was also included in a project (2019.1.00663.S; PI: Butler) with the main goal of investigating the outflows in high-redshift star-forming galaxies, by tracing the OH+ and CO(9–8) lines.

In the following, we describe the calibration, imaging, and analysis of further data sets with the highest available angular resolution in the ALMA Science Archive for J1135. These spatially resolved ALMA follow-ups reveal an almost complete Einstein ring, confirming with no doubt the lensed nature of this system.

The object has been the target of an ALMA Cycle 4 high-resolution follow-up in band 8 (2016.1.01371.S; PI: Vishwas), which aimed to resolve the lensed morphology of the source as well as trace the continuum at ∼0.7 mm and the C[ii] 158 μm FIR line. The continuum was observed to exploit four basebands of width 1.98 GHz, centered at 472.284, 470.451, 460.409, and 458.534 GHz, and composed of 128 channels each.

We recalibrated the raw data by using the Common Astronomy Software Applications (CASA; McMullin et al. 2007) package, version 4.7.2, and running the provided calibration scripts. The continuum subtraction was done manually, using the task uvcontsub. The imaging was also performed manually, by adopting a Briggs weighting scheme, which assumes a robustness factor of 0.5. The properties of the generated images are reported in Table 1, while the continuum-cleaned images are shown in Figure 1. Figure 2 reports the C[ii] spectral line profile and the moment maps corresponding to the integrated brightness, the velocity distribution, and the velocity dispersion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second data set that we examine is part of the ALMA Cycle 6 project (2018.1.00861.S; PI: Yang), which was carried out with the goal of tracing H₂O and CO (J = 8−7) lines in candidate lensed galaxies at high redshift (z ∼ 2–4) in bands 6 and 7. Both observations were performed using the same configuration, with a maximum baseline of 1397 m and four spectral windows of 1.875 GHz bandwidth and 240 × 7.8 MHz channels each. In band 6, the H₂O(J = 2_0,2−1_1,1) and CO(J = 8−7) are targeted with two spectral windows, centered at 239.376 GHz and 223.583 GHz, respectively, while the two other windows, centered at 235.940 and 221.705 GHz, are dedicated to continuum observations. In Band 7, two spectral windows, centered at 281.766 and 292.621 GHz, target the H₂O(J = 3_2,1−3_1,2) and H₂O(J = 4_2,2−4_1,3) lines, while the continuum is observed in two windows centered at 280.314 and 294.266 GHz.

Calibration is performed by running the available pipeline scripts in CASA version 5.4.0–68. Imaging is performed manually, adopting a Briggs weighting scheme in both bands 6 and 7, with a robustness parameter equal to 0.5. We image the CO(8−7) line by performing an automatic continuum subtraction. Figure 3 shows the CO(8−7) spectral line profile and moment maps.

Zoom In Zoom Out Reset image size

The main features of the ALMA data analyzed in this work and the properties of the final images are reported in Table 1. Note that the H₂O data cubes included in the Cycle 6 observations will not be analyzed in this paper.

The flux densities derived for the continuum emission of the lensed source are reported in Table 2. We also include the flux density values measured from the archival image of the Band 3 continuum emission mentioned in Section 2. The flux density uncertainties are computed by including at least a 5% estimation of the flux calibration accuracy (ALMA proposer's Guide): ⁸

$$\begin{eqnarray}&&\delta {S}_{\mathrm{image}}=\sqrt{N{({\sigma }_{\mathrm{image}})}^{2}+{(0.05\times {S}_{\mathrm{image}})}^{2}},\end{eqnarray} \tag{ 1 }$$

with σ_image being the rms and N the number of resolution elements inside the aperture adopted to extract the flux density.

By fitting the resolved ALMA spectral lines with a single Gaussian profile, we obtain the FWHM values for both the CO(8−7) and C[ii] lines, corresponding to 215. ± 4 and 181 ± 5 km s⁻¹, respectively, in concordance with what is found for the GBT and NOEMA CO and H₂O lines analyzed in Harris et al. (2012) and Yang et al. (2017). The peaks are detected at ν_obs = 460.504 ± 0.003 GHz for C[ii] and ν_obs = 223.356 ± 0.0011 GHz for CO(8−7), confirming the redshift estimate by Harris et al. (2012) of 3.127, with associated uncertainties of δ z_C[ii] = ±0.005 and δ z_CO(8–7) = ±0.003, respectively. The observed magnified line profiles measured within a region containing the whole source emission are shown in Figures 2 and 3. Following Carilli & Walter (2013), we compute the observed magnified C[ii] and CO(8−7) luminosities, expressed in units of K km s⁻¹, as:

$$\begin{eqnarray}&&L{{\prime} }_{\mathrm{line}}=3.25\times {10}^{7}\times {S}_{\mathrm{line}}{\rm{\Delta }}v\displaystyle \frac{{D}_{L}^{2}}{{(1+z)}^{3}{\nu }_{\mathrm{obs}}^{2}},\end{eqnarray} \tag{ 2 }$$

where S_line Δ v is the measured flux of the line profile (in units of Jy km s⁻¹) and D_L is the luminosity distance. The luminosities expressed in L_⊙ are computed as L_line = 3 × 10⁻¹¹ ν³ _rest L'_line. The final values computed for the C[ii] and CO(8−7) lines are summarized in Table 3.

J1135 is covered by several surveys, such as the Kilo-Degree Survey (de Jong et al. 2013) and the Hyper Suprime-Cam Subaru Strategic Program in the UV/optical bands (Aihara et al. 2018, 2022), the VISTA Kilo-degree Infrared Galaxy Public Survey (VIKING; Edge et al. 2013) and the UK Infrared Deep Sky Survey Large Area Survey (Lawrence et al. 2007) in the NIR, and the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) in the MIR. PACS and SPIRE FIR observations are reported in the H-ATLAS first and second data release catalogs (Valiante et al. 2016; Maddox et al. 2018). Moreover, the source is also covered by the VLA Faint Images of the Radio Sky at Twenty-Centimeters (Becker et al. 1995) survey in the radio band, where no emission is detected.

High-resolution NIR follow-up observations are available for J1135. The target was observed as part of the Cycle 19 HST/WFC3 snapshot program (PI: Negrello), at a wavelength of λ = 1.15 μm (see Negrello et al. 2014 for further details of the observations), and with the Keck telescope in adaptive optics in the Ks band (Calanog et al. 2014). No successful detection has been found in the Ks image, while a marginal emission (≲3σ) is present in the HST image; however, given the insufficient sensitivity and angular resolution, it is not possible to unambiguously confirm whether it is associated with the foreground lens or the background source.

The object has also been detected in MIR observations that are available in the Spitzer/IRAC Data Archive (PI: Cooray; ID: 80156) and are described in Ma et al. (2015), covering IRAC channel 1 and channel 2, at 3.6 μm and 4.8 μm, respectively.

In addition, we find EVLA radio data in the NRAO Archive, particularly follow-ups in the C band centered at ∼6 GHz (project code: 16A-240; PI: Smith). The data are processed by running the calibration scripts, while cleaning is performed manually with CASA, adopting an interactive mask. The final image reaches a mean rms of ∼0.013 mJy beam⁻¹ and a restored beam ellipse of 113 × 084 (see Figure 4 ).

Zoom In Zoom Out Reset image size

The multiband (optical-to-MIR) image cutouts of J1135 are reported in Figure 5. A faint emission at ∼4σ emerges, starting in the VIKING H band and being detected in both IRAC channels, with a signal-to-noise ratio (S/N) ≳6, but the angular resolution is not sufficient to resolve any lensing features (e.g., arcs) in the NIR/MIR regime. The flux densities are estimated by performing aperture photometry, with aperture diameters of 2'' for NIR VIKING images and 6'' for Spitzer/IRAC images. Table 2 summarizes the photometry for J1135; we report the upper limits for nondetections (i.e., emission with S/N ≲ 3).

Zoom In Zoom Out Reset image size

In reconstructing the source's light profile, we first need to assume a density profile for the mass of the foreground object. The lens is modeled as a singular isothermal ellipsoid (Kormann et al. 1994), i.e., an elliptical power-law density distribution that follows ρ ∝ r^−γ , with r being the elliptical radius, and with a fixed slope value γ = 2. The profile is described by five parameters: the Einstein radius θ_E, the lens centroid positions x_c , y_c , and the first and the second ellipticity components of the elliptical coordinate system (e_x , e_y ). The latter originate from two quantities: the positional angle (ϕ), defined counterclockwise from the positive x-axis, and the factor f = (1–q)/(1 + q), where q is the ratio between the semimajor and semiminor axis. The final expressions for the elliptical components are:

$$\begin{eqnarray}&&\begin{array}{l}{e}_{y}=f\times \sin \,(2\phi ),\\ {e}_{x}=f\times \cos \,(2\phi ).\end{array}\end{eqnarray} \tag{ 3 }$$

From several tests, we have verified that the shear component does not improve the model, but rather worsens the fit results, and for this reason it has been omitted.

PyAutoLens performs lens fitting through the nested sampling algorithm Dynesty (Speagle 2020), which samples the parameter space and computes the posterior probability distributions for the parameters of a given lens model.

Our searching chain consists of two steps. (1) We search for the best-fit lens model through nonlinear parametric fitting. We assume the source light to be described by a Sérsic profile and simultaneously fit the ALMA data in all three bands (including spectral lines). (2) Keeping the best-fit lens model parameters obtained in the first step fixed, we then perform the inversion. The fit is performed on a number of pixels delimited by a circular mask, where the radius changes according to the resolution of the cleaned ALMA image, in order to obtain a satisfactory fit that is not excessive in terms of the computational cost.

The best-fit lens model is described by the following parameters: ${\theta }_{\mathrm{Ein}}^{\mathrm{mass}}={0.4241}_{-0.0005}^{+0.0005}$, ${y}^{\mathrm{mass}}=-{0.2329}_{-0.0005}^{+0.0005}$, ${x}^{\mathrm{mass}}={0.1494}_{-0.0009}^{+0.0011}$, $q={0.637}_{-0.001}^{+0.001}$, and $\phi =-{35.31}_{-0.04}^{+0.04}$. The best-fit parameters and their uncertainties in the outputs from the parametric nonlinear search and inversion procedure are reported in Table 4. This key information allows us to retrieve the intrinsic properties of the lensed background object, which will be discussed in detail in Section 5. The resulting reconstructed source contains only pixels that have been excluded from the masked lensed image area. The magnification factors are computed as μ = A_IP,5σ /A_SP,5σ , where A_IP,5σ and A_SP,5σ are the areas enclosing significant (i.e., >5σ) pixels in the reconstructed image plane (IP) and reconstructed SP, respectively. The noise is estimated as the rms in the reconstructed source map. From the area enclosing all the pixels with S/Ns ≳3 and ≳5 in the reconstructed SP, the effective radius can be computed as ${r}_{\mathrm{eff}}={({A}_{\mathrm{SP}}/\pi )}^{0.5}$. In order to determine the uncertainties of these parameters, we exploit the set of samples provided by the nonlinear search that was performed during the inversion. Each sample corresponds to a set of inversion parameters (i.e., the regularization coefficient and the pixelization shape) that were evaluated and accepted by the nonlinear search. The uncertainties are then retrieved as the 16th and 84th quantiles of the parameter distribution drawn from ∼200 accepted samples.

Table 4. Output Properties of the Lens Modeling and Source Reconstruction Analysis

	μrec	θrec	Reff,3σ	Reff,5σ	Reff,par
		(arcsec)	(pc)	(pc)	(pc)
Band 8	${12.80}_{-0.18}^{+0.32}$	0.03	${541}_{-10}^{+18}$	${392}_{-13}^{+5}$	${571}_{-142}^{+160}$
Band 7	${8.34}_{-0.05}^{+0.2}$	0.07	${518}_{-70}^{+35}$	${314}_{-24}^{+64}$	${995}_{-34}^{+38}$
Band 6	${6.51}_{-1.56}^{+1.48}$	0.1	${1146}_{-158}^{+149}$	${770}_{-183}^{+157}$	${789}_{-65}^{+75}$
C[ii]	${6.01}_{-0.38}^{+0.44}$	0.04	${596}_{-42}^{+43}$	${525}_{-50}^{+40}$	${764}_{-50}^{+110}$
CO(8−7)	${6.22}_{-0.87}^{+0.79}$	0.1	${1200}_{-89}^{+106}$	${970}_{-105}^{+84}$	${936}_{-270}^{+248}$

Note. From the left: the magnification factor, the demagnified angular resolution, and the effective radius for 3σ and 5σ emission. The last column reports the best-fit value for the effective radius of the Sérsic profile adopted in the first step of our analysis.

Download table as:  ASCIITypeset image

It is crucial to point out that the robustness of these quantities strongly depends on the noise covariance that, in the case of interferometric source reconstruction, is difficult to quantify (Rizzo et al. 2021; Stacey et al. 2021). In order to understand how much the estimation of the noise impacts on our results, we compare the effective radii inferred from the source reconstruction with the values that we report for the best-fit effective radii of the Sérsic profile (R_eff,par) used in the nonlinear parametric fit described in Section 3.1. The values are shown in Table 4. We find that the effective radii inferred from the two different methods are broadly consistent with one another. The only exception is for Band 7, where R_eff,par is higher with respect to the values computed from the source reconstruction, likely due to a higher noise level. Therefore, the parameter uncertainties from the inversion are likely to be underestimated. Figures 6 and 7 show the original lens plane image, the model image, the normalized residual map, and the reconstructed source for the three ALMA continuum bands and the CO(8−7) and C[ii] emission lines, respectively. The differences in the retrieved physical scale values reflect the heterogeneity of the data adopted in this work—the product of different array configurations and angular resolutions. The central feature that emerges in the normalized residuals of the continuum in bands 8, 7, and 6 could originate from the faint foreground lens, the nature of which is discussed in detail in the following section. Our attempt to remove the emission by light profile subtraction was unsuccessful, mainly because of its extreme compact size and faintness compared to the background lensed source. However, the presence of the lens in the continuum bands does not have a considerable impact on the lens model in bands 8 and 7, where we estimate its contribution to be no higher than ∼4% with respect to the lensed source flux density (i.e., computed by including only the light originating by the arcs in the lens plane). One caveat concerns the band 6 continuum, where the contribution of the central emission is estimated to be up to ∼28% of the lensed source flux density, and where the lens modeling results are more difficult to interpret. In this case, the uncertainties have an impact on the estimates of the magnification factors and effective radii, and are included in the errors of these quantities. Since the band 6 continuum is only used in the SED fitting (see Section 4), and given that it is not part of the general discussion, its overall impact on our analysis is negligible. We have decided to keep the lens modeling and the estimated quantities in this paper, since they could be informative for future works regarding J1135.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, we reconstruct the velocity map for the CO(8−7) line by dividing and modeling the emission in three different velocity bins. As there is no significant difference in the reconstructed emission in the bins, we cannot claim any indications of rotation or outflow (see Figure 8). This is also noticeable in the first moment maps shown in Figure 3, where no strong velocity gradients are visible along the widths of the arcs. The velocity dispersion peak is co-spatial with the integrated brightness peak, located in the southern region of the arc, and corresponding to the northern clump in the reconstructed velocity map. Our attempt to reconstruct the C[ii] velocity field was unsuccessful, because of the poor S/N of the individual velocity bins, even though a modest velocity gradient is visible in Figure 2, peaking at the same position as the velocity dispersion.

Zoom In Zoom Out Reset image size

One peculiar aspect of the J1135 gravitational lensed system is the faintness of the foreground object. Indeed, no redshift estimate is available for the lens galaxy, and no clear detection is measured from the photometric images, likely due to insufficient sensitivity and/or angular resolution. A faint central emission is detected in the residuals of the lens modeling, and visible in Figure 6, but it is not clear whether it should be unambiguously associated with the lens or with some spurious emission related to the noise. As shown in analogous studies, and as revealed by HST/NIR high-resolution images (e.g., Negrello et al. 2014), the foreground object usually dominates the emission in these bands, with a progressively higher contribution coming from the background source at higher wavelengths. For this reason, in order to achieve reliable results from the SED fitting procedure, it is essential to fit and subtract the light profile of the foreground galaxy. In this case, however, only a marginal emission (≲3σ) comes from the HST WFC3/F110 data, and it is not possible to establish a priori whether it originates from the lens or from the lensed object. We therefore explore the assumption of the lens as a massive elliptical, and attribute its faintness to its relatively high redshift (e.g., z ≳ 1.5). We model the SED of the foreground object according to this assumption, and we constrain its luminosity by means of the Einstein (total) mass resulting from the lens modeling (M_E ∼ 1.15 × 10¹¹ M_⊙). Specifically, we adopt the template for an elliptical galaxy with a 2 Gyr age from the SWIRE library (Polletta et al. 2007). The resulting SED of this template, overlapped with the photometry reported in Table 2, is shown in Figure 9. We find the contribution from the lens to be negligible for the flux densities from the H and Ks VIKING bands up to higher wavelengths, so no lens subtraction is needed. The situation is less clear for the marginal HST WFC3/F110 detection, and, for this reason, we consider this value as an upper limit.

Zoom In Zoom Out Reset image size

Table 5 reports the best-fit dust mass estimated by CIGALE. The main advantage of choosing the Draine et al. (2014) multiparameter library is due to its physical motivation, as the dust is indeed described as a mixture of carbonaceous and amorphous silicate grains. This results in a more robust estimate of the dust mass (M_dust) with respect to a single-temperature modified blackbody (MBB) fit, which tends to underestimate M_dust by a factor of ∼2 (Magdis et al. 2012; Berta et al. 2016), resulting in a wrong gas mass derivation. Another issue has to do with the dust temperature estimation of high-redshift star-forming galaxies. Different results in the literature (e.g., Riechers et al. 2013; Spilker et al. 2016; Scoville et al. 2017; Simpson et al. 2017; Jin et al. 2019; Cortzen et al. 2020; Jin et al 2022) suggest that the widely adopted optically thin MBB approximation may not be sufficient for inferring the highest dust temperatures for dusty and optically thick galaxies. For example, Jin et al. (2019) and Jin et al (2022) reported the presence of a population of compact high-redshift (z ∼ 4) starbursts, selected in the FIR with Herschel and detected with ALMA and NOEMA, showing abnormally cold dust temperatures. This behavior can either be associated with low star formation efficiency (SFE), accompanied by the rapid enrichment of metals, or with a dust continuum in an optically thick FIR regime (Cortzen et al. 2020).

Table 5. Best-fit Output Parameters from CIGALE—from the First Row: Dust Luminosity, Dust Mass, SFR, Stellar Mass, and the IRRC Parameter

SED-fitting Results
$\mathrm{log}{L}_{\mathrm{dust}}$ (L⊙)	13.03 ± 0.06
$\mathrm{log}{M}_{\mathrm{dust}}$ (M⊙)	9.06 ± 0.04
logSFR (M⊙ yr−1)	2.97 ± 0.08
$\mathrm{log}{M}_{\star }$ (M⊙)	≲11.75
qIR	2.84 ± 0.09

Download table as:  ASCIITypeset image

Following Cortzen et al. (2020), we compute the dust temperatures using three different approaches: from the CIGALE results, we adopt the physically motivated Draine et al. (2014) model, we assume an optically thin MBB, and we assume an optically thick MBB. The Draine & Li (2007) and Draine et al. (2014) model assumes an optically thin dust emission and does not provide a luminosity weighted temperature. However, following Draine (2011), we retrieve the dust temperature as ${T}_{\mathrm{dust},\mathrm{DL}14}=20{U}_{\min }^{1/6}$ K, where ${U}_{\min }$ is the value of the minimum intensity of the radiation field, as inferred from CIGALE. From the best-fit value of ${U}_{\min }=44.3\pm 10.6$, we obtain a dust temperature of T_dust,DL14 = 37.7 ± 1.5 K.

By exploiting the available FIR-to-submm photometry, we then fit a single-temperature MBB under the optically thin approximation, described as

$$\begin{eqnarray}&&{S}_{\nu }\propto \displaystyle \frac{{\nu }^{3+\beta }}{{e}^{(h\nu /{{kT}}_{d,\mathrm{thin}})}-1},\end{eqnarray} \tag{ 6 }$$

where k is the Boltzmann constant and β is the dust emissivity index, which is here assumed to be β = 2. Similarly, for the optically thick regime, we fit the single-temperature MBB defined as

$$\begin{eqnarray}&&{S}_{\nu }\propto (1-{e}^{-{\tau }_{\nu }})\displaystyle \frac{{\nu }^{3}}{{e}^{(h\nu /{{kT}}_{d,\mathrm{thick}})}-1},\end{eqnarray} \tag{ 7 }$$

with ${\tau }_{\nu }={(\nu /{\nu }_{0})}^{\beta }$ and ν₀ being the rest-frame frequency corresponding to a dust opacity equal to unity.

The best-fit results for both the optically thin and optically thick approximations are represented in Figure 11. The resulting dust temperatures are T_dust,thin = 38.6 ± 1.1 K and T_dust,thick = 64.8 ± 2.8 K, respectively, where the former is consistent with T_dust,DL14.

Zoom In Zoom Out Reset image size

From our analysis, we derive a significant discrepancy between T_dust,thin and T_dust,thick, which is apparently in agreement with what has been observed for similar objects in the literature. In particular, the recent work of Jin et al (2022) has presented useful diagnostics that can be adopted to identify whether an optically thick model is more appropriate for describing dust emission (e.g., the SFE and the position relative to the IR luminosity surface density diagram). By taking advantage of these indicators, we investigate the nature of the dust continuum emission of J1135. The star formation surface density that we retrieve in Section 5.3 is well above the threshold suggested by the authors for the optically thick regime (i.e., Σ_SFR ≳ 20 M_⊙ yr⁻¹ kpc⁻²). Moreover, the dust temperature that is inferred from the optically thin approximation implies ${{\rm{\Sigma }}}_{\mathrm{IR}}\gtrsim 1.5\times {10}^{5}{T}_{\mathrm{dust},\mathrm{thin}}^{4.21}$, which is inconsistent with the Stefan–Boltzmann law, representing the upper boundary for an object to emit as a blackbody. In Figure 12, we show the dust temperature against the intrinsic IR luminosity (∝SFR), color coded in relation to the effective radius for J1135. We compare our source to other samples of DSFGs, according to its position with respect to the modified Stefan–Boltzmann law inferred by Yan & Ma (2016), for different effective radii. As expected, temperatures of the order of T_dust,thin and T_dust,DL14 would imply effective radii ≳1.5 kpc at the inferred L_dust (hereafter, L_IR), which are well above what we measure from the reconstructed dust continuum morphology. In this plot, we therefore assume the case T_d = T_dust,thick, obtaining a more consistent result. This is indeed in concordance with what Jin et al (2022) predicted for apparently cold starbursts, disfavoring the low-efficiency star formation mechanism accompanied by fast metal enrichment.

Zoom In Zoom Out Reset image size

We then compare J1135 with other samples of DSFGs: we include the sample of apparently cold FIR-selected starbursts from Jin et al (2022) in the redshift range 3 ≲ z ≲ 6, the sample of FIR-to-submm-selected lensed quasars (1.5 ≲ z ≲ 2.5) from Stacey et al. (2021), and, finally, the sample of lensed DSFGs selected by Herschel (Nayyeri et al. 2016; Negrello et al. 2017), and analyzed in works by Nayyeri et al. (2016) and Dye et al. (2018), distributed over the range 1 ≲ z ≲ 4 (the values are reported in Table A.1. of Stacey et al. (2021). J1135 shows smaller sizes compared to other lensed DSFGs at similar luminosities (L_IR ≳10¹³ L_⊙) and temperatures that are comparable with some of the warmer lensed quasars.

The available data allow us to estimate the gas content by adopting several empirical calibrators. First, we directly estimate the gas mass from C[ii]: following Zanella et al. (2018), we assume α_C[II] ≡ M_gas/L_C[II] = 22 M_⊙/L_⊙, which is calibrated from starburst galaxies spanning a redshift range z ∼ 2–6. Second, in order to estimate the molecular gas content (${M}_{{{\rm{H}}}_{2}}$), we derive the $L{{\prime} }_{\mathrm{CO}(1-0)}$ from the demagnified $L{{\prime} }_{\mathrm{CO}(8-7)}$ luminosity. We then follow Fujimoto et al. (2022), by adopting a conversion factor of $L{{\prime} }_{\mathrm{CO}(1-0)}=1.5L{{\prime} }_{\mathrm{CO}(7-6)}$, as estimated for high-redshift starburst galaxies in the literature (e.g., Riechers et al. 2013). This conversion factor is referred to a different transition, corresponding to the higher luminosity values of the CO-SLED (Yang et al. 2017), and for this reason the resulting value of $L{{\prime} }_{\mathrm{CO}(1-0)}\sim 1.6\times {10}^{10}$ K km s⁻¹ pc² is considered an upper limit. This estimate is consistent with the value of $L{{\prime} }_{\mathrm{CO}(1-0)}\sim 1.5\times {10}^{10}$ that was found by Harris et al. (2012), by adopting an indicative magnification factor of 10. We then compute the molecular gas mass by assuming two different values of α_CO = 0.8–4.6. The value of α = 0.8 is calibrated from local ULIRGs with supersolar metallicity (Downes & Solomon 1998), while the higher value is calibrated from the Milky Way (Solomon & Barrett 1991).

The molecular gas content can also be estimated by means of the empirical calibration (Scoville et al. 2017) as α ≡ < L_ν850μm/M_gas > = 6.7 ± 1.7 × 10¹⁹ erg s⁻¹ Hz⁻¹ ${M}_{\odot }^{-1}$.

Finally, we convert the dust content into gas, by assuming a variable gas-to-dust ratio of δ_GDR = 30–92, referring to the typical solar and supersolar metallicities, following Magdis et al. (2012) and Fujimoto et al. (2022). The values obtained for the molecular gas masses are reported in Table 6. For consistency, we adopt the dust mass value inferred from CIGALE, but it is worth checking how the dust opacity could affect our results. However, we caution the reader against deriving the dust mass through the effective dust temperatures inferred from a single-temperature MBB, as this estimate does not reflect the actual dust content of the ISM (see, e.g., Scoville et al. 2016). We compute the temperature-weighted dust masses in the optically thick regime, by adopting the value of T_dust,thick obtained in Section 5.1. We infer $\mathrm{log}{M}_{\mathrm{dust},\mathrm{thick}}\,=8.7\pm 0.05$ M_⊙, which is ∼2.3 times lower than the CIGALE estimate, and would result in a gas mass range of $\mathrm{log}{M}_{\mathrm{gas},\mathrm{GDR}}\sim (10.2\mbox{--}10.7)$ M_⊙.

The stellar mass estimate that was output from the SED fitting must be considered an upper limit. Indeed, given the lack of a clear detection in the NIR images, it is not possible to correctly estimate the contribution coming from the lens (see Section 3.3 for a further discussion). Moreover, the dark nature of this object hinders a complete sampling of the optical and NIR parts of the SED. Aside from the value reported in Table 5, we compute the stellar mass by assuming a typical stellar-to-dust mass ratio of δ_SDR ≈ 100, obtaining a value of ${M}_{* }^{\mathrm{STD}}\sim 1.1\times {10}^{11}$ M_⊙, which is in agreement with the SED fitting estimate.

From the reconstructed continuum originating from the dust (Figure 6), a clear "clumpy" double-peaked structure is visible. First, we compare the effective radii inferred from the dust continuum at different wavelengths. Excluding the lower–angular resolution reconstruction at 1.3 mm (band 6), the most resolved 0.64 and 1.0 mm emissions do not show particular differences in their spatial extensions that are ascribable to a dust temperature gradient, even though further observations are needed in order to explore this possibility. In Figure 13, we then overlap the SP maps reconstructed from our lens modeling for three different tracers. We show the C[ii] and CO(8−7) line emission and the dust continuum at 640 μm, the latter corresponding to the data set with the highest angular resolution. The first visible difference between the three emissions concerns their spatial extensions: while the C[ii] and dust continuum occupy similar areas, the CO(8−7) is more extended. The size discrepancy is most likely associated with the differences in angular resolution of the respective data sets (see also the R_eff and θ values reported in Table 4), rather than being an intrinsic morphological difference. We also note that this difference is present in the dust continuum for the same data set as CO(8−7), i.e., band 6, where the estimated effective radii reach up to 1.1 kpc at 3σ. This comparison shows a co-spatial emission for the dust continuum and C[ii] line, even though the peaks are located at different positions. The clumpy structures also seem to be present in the reconstructed C[ii] line emission, while they are not likely to be resolved for the CO(8−7) line, which extends in a more ellipsoidal profile.

Zoom In Zoom Out Reset image size

From the average of the gas mass values reported in Section 5.2, we estimate the depletion timescale as τ_depl ≃ M_gas/SFR ≃ 10⁸ yr, which translates to a high SFE, reaching up to ∼10⁻⁸ yr⁻¹. Moreover, the inferred stellar mass implies a star formation burst duration of τ_SFR ≃ M_*/SFR ≃ 10⁸ yr, which is indicative of a young galaxy offset from the main-sequence locus of star-forming galaxies at z ∼ 3 (Speagle et al. 2014).

Our results are consistent with the expectations reported in Vishwas et al. (2018), where the analysis of the Lyman continuum photons required to sustain the luminosity of the O[iii] 88 μm line pointed out the presence of young and massive stars ionizing the surrounding H ii regions. The same authors found no significant AGN contribution from the SED analysis, consistent with what we infer from the IRRC (q_IR ≈ 2.8), which is indicative of a star formation–dominated object. However, it should be pointed out that in the absence of a good sampling of the MIR part of the SED, the presence of the AGN in J1135 is still arguable. Moreover, the CO(8−7) line, associated with high transitions, can point toward the presence of large amounts of energy that are linked with activity from a heavily dust-embedded central nucleus, as well as a strong star formation activity.

The hypothesis of J1135 being a compact starburst is also supported by the source reconstruction of the highest–angular resolution ALMA dust continuum emission at 640 μm and 1.1 mm, shown in Figure 6, where the effective radius reaches an average value of ∼350 pc.

We now focus on the properties of the C[ii] line of J1135. The C[ii] line is a fine-structure line that predominantly originates from high-z photon-dominated regions and is typically used as a cool interstellar gas tracer as well as an SFR estimator (see Casey et al. 2014 for a review). A well-known deficit in the C[ii]/FIR ratio is observed in both nearby (e.g., Luhman et al. 2003; Diaz-Santos et al. 2017; Smith et al. 2017) and high-redshift star-forming galaxies (Stacey et al. 2010; Gullberg et al. 2015). This drop is found to reach very low values (L_C[II]/L_IR ≈ 10⁻⁴) in spatially resolved studies (e.g., Gullberg et al. 2015; Lagache et al. 2018; Rybak et al. 2019). For J1135, we infer a C[ii]/IR ratio of L_C[II]/L_IR ≈ 5.4 × 10⁻⁴. Similar values are found for other strongly lensed galaxies among the H-ATLAS sample. For example, Rybak et al. (2020) reported a deficit down to ∼3 × 10⁻⁴ for spatially resolved ALMA data of SDP.81 at z = 3.042 (Dye et al. 2015; Hatsukade et al. 2015; Partnership et al. 2015; Rybak et al. 2015a, 2015b; Swinbank et al. 2015; Tamura et al. 2015; Hezaveh et al. 2016). Lamarche et al. (2018) found similar values (∼2×10⁻⁴) for SDP.11 at z = 1.7, even though our galaxy shows a more compact morphology in the C[ii] emission with respect to other objects. We compute the star formation surface density by assuming ${\rho }_{\mathrm{SFR}}=0.5\times \tfrac{\mathrm{SFR}}{\pi {R}_{\mathrm{eff},5\sigma }^{2}}$ (Stacey et al. 2021), where R_eff,5σ is the average of the effective radii of the 5σ continuum dust emission at 640 μm and 1.1 mm reported in Table 4. We infer a value of ρ_SFR ∼ 1200 M_⊙ yr⁻¹ kpc⁻², which is consistent with a galaxy being on the verge of the Eddington limit for a radiation pressure–supported starburst (Andrews & Thompson 2011; Simpson et al. 2015). This result is compatible with the possible explanation of the deficit being attributed to a lower increase of the C[ii] emission with respect to the IR.

By inspecting the HST/WFC3 image, we find no evidence for galaxy companions of J1135 within a radius of at least ∼5'', corresponding to ∼40 kpc, although the detection of possible closer and fainter sources is hindered by the current data sensitivities and angular resolutions. Moreover, from the CO(8−7) reconstructed image and the C[ii] first moment map, we find no clear evidence of a complex kinematic structure possibly associated with a merging event. With no further hints pointing toward this scenario, we are led to interpret the ISM conditions and physical properties discussed so far in the light of in situ galaxy formation scenarios (Lapi et al. 2014; Mancuso et al. 2016a, 2016b, 2017; Lapi et al. 2018; Pantoni et al. 2019). In particular, the properties of J1135 are consistent with a compaction phase (Barro et al. 2014; Ikarashi et al. 2015, 2017; Kocevski et al. 2017; Silverman et al. 2019; Valentino et al. 2020; Stacey et al. 2021; and see also Figure 3 in Lapi et al. 2018), in which the dust-enshrouded star formation activity increases at an almost constant rate in the inner regions of the galaxy where the stellar mass is being accumulated. At this stage, the in situ scenario envisages the galaxy as an off-main-sequence object at an early evolutionary stage, which will eventually move toward the main-sequence locus as the stellar mass content increases. Finally, the star formation will either progressively decrease, as the galaxy exhausts its gas reservoir, or it will be abruptly stopped by the action of the feedback from an AGN (Mancuso et al. 2016a, 2017).