ALMA resolves the first strongly-lensed Optical/NIR-dark galaxy

We present high-resolution ($\lesssim0.1$arcsec) ALMA observations of the strongly-lensed galaxy HATLASJ113526.2-01460 at redshift $z\sim3.1$ discovered in the Gama 12$^{\rm th}$ field of the Herschel-ATLAS survey. The gravitationally lensed system is remarkably peculiar in that neither the background source nor the foreground lens show a clearly detected optical/NIR emission. We perform accurate lens modeling and source morphology reconstruction in three different (sub-)mm continuum bands, and in the C[II] and CO(8-7) spectral lines. The modeling indicates a foreground lensing (likely elliptical) galaxy with mass $\gtrsim10^{11}\, M_\odot$ at $z\gtrsim1.5$, while the source (sub-)mm continuum and line emissions are amplified by factors $\mu\sim6-13$. We estimate extremely compact sizes $\lesssim0.5$ kpc for the star-forming region and $\lesssim 1$ kpc for the gas component, with no clear evidence of rotation or of ongoing merging events. We perform broadband SED-fitting and retrieve the intrinsic de-magnified physical properties of the source, which is found to feature a very high star-formation rate $\gtrsim10^3\, M_\odot$ yr$^{-1}$, that given the compact sizes is on the verge of the Eddington limit for starbursts; the radio luminosity at 6 cm from available EVLA observations is consistent with the star-formation activity. The galaxy is found to be extremely rich in gas $\sim10^{11}\, M_\odot$ and dust $\gtrsim10^9\, M_\odot$. The stellar content $\lesssim10^{11}\, M_\odot$ places the source well above the main sequence of starforming galaxies, indicating that the starburst is rather young with estimated age $\sim10^8$ yr. Our results indicate that the overall properties of HATLASJ113526.2-01460 are consistently explained by in-situ galaxy formation and evolution scenarios.


INTRODUCTION
Sub-millimetre galaxies (SMGs) are the main protagonists of the star formation at early cosmic times (Blain 1996, Casey et al. 2014. It is well established, that a substantial contribution at the peak of the cosmic Star Formation Rate (SFR) density comes from these heavily dust-obscured objects, featuring a sub-millimeter (submm) flux density S 870µm 1 mJy and extremely high SFRs, up to ∼ 10 3 M yr −1 (e.g. Simpson et al. 2020). Because of their huge dust content these objects are heavily obscured in optical bands and extremely bright in far-infrared (FIR)/sub-mm bands where the light of newborn stars, reprocessed by dust, is re-emitted. Moreover, SMGs have been identified as the progenitors of massive quiescent early-type galaxies and constitute the ideal laboratories to test galaxy evolutionary models. For example, in in-situ coevolutionary scenarios (Lapi et al. 2014, 2018, Pantoni et al. 2019, the intense star formation activity is accompanied by an exponential growth of the active nucleus, whose feedback will eventually sweep away the interstellar medium. The star formation is thus stopped on a relatively short timescale while the nucleus shines as an optical quasar.
In the last years, an even more extreme population of heavily obscured SMGs has been discovered. These objects are missed in optical/near-IR (NIR) surveys and have been found up to very high redshifts (z∼ 6; Riechers et al. 2013, 2020, Marrone et al. 2018. These heavily obscured star-forming galaxies often lack of a counterpart even in deep NIR observed-frame Hubble Space Telescope (HST) (Wang et al. 2019, Gruppioni et al. 2020 or either show extreme red colors (H − 3.6 µm > 4; see e.g. Wang et al. 2016) and are visible only from observed-frame mid-IR(MIR) images performed e.g. with the Spitzer /Infrared Array Camera (IRAC). Samples of optical/NIR dark objects have been detected by observing deep CO line emission (Riechers et al. 2020), and have been efficiently selected in sensitive radio observations (Talia et al. 2021, Enia et al. 2022). These peculiar objects provide a significant and previously unknown contribution to the cosmic SFR density at z 3 estimated to be at least 10% up to 25−40% with respect to the one inferred from UV-selected populations (Wang et al. 2019, Williams et al. 2019, Gruppioni et al. 2020, Talia et al. 2021, Enia et al. 2022.
The studies conducted so far are however limited by the poor angular resolution and sensitivity in MIR/FIR and sub-mm bands, causing confusion problems and prohibiting a detailed investigation of the physical properties of optical/NIR dark galaxies and the conditions of their Interstellar Medium (ISM). In the last years, Atacama Large Millimeter/submillimeter Array (ALMA) deep field observations strongly improved the quality of the observations of high redshift dusty galaxies, detecting SMGs up to flux density limits of S 870µm ∼ 0.1 − 1 mJy , Dunlop et al. 2017, Franco et al. 2018, Hatsukade et al. 2018). However, even high angular resolution studies indicate that these objects are extremely compact, with typical intrinsic sizes of a few tenths of an arcsec (Pantoni et al. 2021, hence very hard to resolve. Gravitational lensing enables the observation of regions in the luminosity-redshift space of these sources, that would be otherwise unattainable with current instrumentation in reasonable integration times. The gravitational magnification of the foreground lens increases the apparent luminosity proportionally to the magnification µ and stretches the angular sizes by a factor √ µ. This behavior offers the unique possibility of studying down to sub-kpc scales the properties of objects otherwise not exceptionally bright, massive, or peculiar, and belonging to the dusty star-forming galaxy population bulk at the peak of cosmic star formation. Several works demonstrated the effectiveness of sub-mm surveys in selecting strong lensing events adopting a flux density threshold of 100 mJy at 500 µm, in correspondence of a steep drop in the number counts of dusty starforming galaxies at sub-mm wavelengths (Blain 1996, Negrello et al. 2010, Lapi et al. 2012 where, thanks to the magnification, they emerge as the bright tail of the population count distribution, thus minimizing the probability of possible contaminants, such as flat spectrum radio sources and low redshift spiral galaxies. Moreover, in the FIR/sub-mm regime, while the high−z lensed dusty star-forming galaxies are particularly bright, negligible signal comes from the foreground lenses, which are often massive ellipticals at z < 1 that dominate the signal in optical bands. Several surveys conducted in the last decade with the Herschel Space Observatory led to the discovery of numerous strong lensing events. The Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) identified 11 lensed galaxies over 95 deg 2 (Wardlow et al. 2013);Nayyeri et al. (2016) selected other 77 candidate lensed galaxies in the HerMES Large Mode survey (HeLMS; Oliver et al. 2012) and in the Herschel Stripe 82 Survey (HerS; Viero et al. 2014). In particular, the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) is the widest area (600 deg 2 ) extragalactic survey undertaken with Herschel and has provided a sample of more than a hundred thousands dusty star-forming galaxies at high redshift. Among the H-ATLAS survey, a sample of 80 candidate strongly lensed dusty star-forming galaxies has been selected in Negrello et al. (2017) by means of a simple flux density selection (S 500µm > 100 mJy). Only 21 of them have been confirmed to be lensed thus far. Recently, another sample of 11 candidates has been selected by Ward et al. (2022) in the H-ATLAS Third Data Release conducted in the South Galactic Pole (SGP). The recent work of Shu et al. (2022) exploited lensing effects generated from galaxy clusters in order to systematically search for optically dark galaxies. Follow-ups performed with JCMT/SCUBA (∼ 850 µm) and ALMA (∼ 870 µm) for their sample reach a flux limits ∼ 3 times deeper than blank fields, highlighting the capabilities of gravitational lensing in detecting even more hidden and dark objects.
In this work, we present the lens modeling, the source reconstruction, and Spectral Energy Distribution (SED) analysis of HATLASJ113526.3−014605 (J1135 hereafter), also called G12v2.43 or G12H43, an optical/NIR dark strongly lensed galaxy at z = 3.1276 belonging to the Negrello et al. (2017) lensed candidate sample, featuring a flux density at 500 µm amounting to 204 ± 8.6 mJy. The plan of the paper is the following: in Sect. 2 we present the target of our analysis, and describe the archival ALMA observations and the available ancillary data in other bands; Sect. 3 and 5 describe respectively the lens modeling and source reconstruction and the SED fitting analysis; finally, we discuss and summarise our results in Sect. 6 and 7. Throughout the work, we adopt the standard flat ΛCDM cosmology (Planck Collaboration et al. 2020) with rounded parameter values: matter density Ω M = 0.32, dark energy density Ω Λ = 0.63, baryon density Ω b = 0.05, Hubble constant H 0 = 100h kms −1 Mpc −1 with h = 0.67, and mass variance σ 8 = 0.81 on a scale of 8 h −1 Mpc. At the redshift of the source 1 arcsec corresponds to 7.8 kpc.
2. THE TARGET J1135 is part of the sample of 80 (candidate) strongly lensed galaxies ) located in the equatorial GAMA 12 th field (RA=11:35:26, dec=-01:46:07, J2000). The spectroscopic redshift of z = 3.1276 of the background lensed source was obtained from blind CO searches with the Zpectrometer ultrawideband spectrometer on the Green Bank Telescope (GBT) (Harris et al. 2012) and confirmed by the Northern Extended Millimeter Array (NOEMA) observations (Yang et al. 2017). So far, no redshift measurement is available for the foreground lens. Andreani et al. (2018) presented observations of high CO transition (J=7-6) obtained with the Atacama Pathfinder EXperiment (APEX)/SEPIA band 5 receiver for the background object. From the comparison of the CO(7-6) transition with the CO(1-0) and CI(2-1) the authors pointed out to the presence of a large excitation status in the ISM of J1135.
Moreover, Vishwas et al. (2018) reported bright [OIII] 88 µm emission for J1135 detected through the z Early Universe Spectrometer (ZEUS-2) on APEX attributed to ionized hydrogen regions around massive stars. From the SED-fitting of the multi-band photometry of J1135, the authors predicted J1135 to be a young, gas-rich starburst galaxy.
The object has also been targeted by Submillimeter Array (SMA) high spatial resolution (∼ 0.8 arcsec) observations described in Bussmann et al. (2013), but only marginally resolved. For this reason, its lensed nature has been debated in the works described above.

ALMA observations
The object is part of low ( 2 arcsec) resolution observations in band 3 (2017.1.01694.S, PI: Oteo) 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) whose main goal was to investigate outflows in high red-shift star-forming galaxies by tracing OH+ and CO(9-8) lines.
In the following, we describe the calibration, imaging and analysis of further data sets with the highest angular resolution available in the ALMA Science Archive for J1135. These spatially resolved ALMA follow-ups reveal an almost complete Einstein ring, confirming out of any doubt the lensed nature of this system. The object has been target of ALMA Cycle 4 highresolution follow-up in band 8 (2016.1.01371.S, PI: Amit) aimed at resolving the lensed morphology of the source and tracing the continuum at ∼ 0.7 mm and the C[II] 158 µm FIR line. The continuum was observed exploiting four base-bands of width 1.98 GHz,centred at 472.284,470.451,460.409,458.534 GHz and composed by 128 channels each.
We re-calibrated the raw data using the Common Astronomy Software Applications (CASA) package version 4.7.2 and running the provided calibration scripts. The continuum subtraction was manually done using the task uvcontsub. Imaging has been performed manually adopting a Briggs weighting scheme, which assumes a robustness factor of 0.5. The properties of the generated images are reported in Table 1, the continuum cleaned images are shown in Fig. 1, and Fig. 2 reports the C[II] channel maps rebinned in a 20 km s −1 interval.
The second data-set we examine is part of the ALMA Cycle 6 project (2018.1.00861.S, PI: Yang) carried out with the goal of tracing H 2 0 and CO (J=8−7) lines in candidate lensed galaxies at high redshift (z ∼ 2 − 4) in band 6 and 7. Both observations are performed with 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 2 0(J=2 0,2 -1 1,1 ) and CO(J=8-7) are targeted with two spectral windows respectively centered at 239.376 GHz and 223.583 GHz, while other two 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 2 O(J=3 2,1 -3 1,2 ) and H 2 O(J=4 2,2 -4 1,3 ) lines, while continuum is observed in two windows centered at 280.314 and 294.266 GHz Calibration is performed running the available pipeline scripts in CASA version 5.4.0-68. Imaging is performed manually adopting a Briggs weighting scheme in both band 6 and 7, with robustness parameter equal to 0.5. We image the CO(8-7) line performing an automatic continuum-subtraction. The main features of the ALMA data analysed in this work and the properties of the final images are reported in Table 1. Note that H 2 O data cubes included in Cycle 6 observations will not be analysed in this paper.

Data analysis
The flux densities derived for the continuum emission of the lensed source are reported in Tab. 3. We also include the flux density value measured from the archival image of the Band 3 continuum emission mentioned in Sect. 2. Flux density uncertainties are computed including a 5% estimation of the flux calibration accuracy: with σ image being the image noise. By fitting the resolved ALMA spectral lines with a single Gaussian profile we obtain the full-width half maximum (FWHM) values for both CO(8-7) and C[II] lines corresponding to 215. ± 4 and 181 ± 5 km s −1 respectively, in concordance to what is found by GBT and NOEMA CO and H 2 0 lines analysed in Harris et al. (2012) and Yang et al. (2017). The peak is detected at ν obs = 460.504 ± 0.003 GHz for C[II] and ν obs = 223.356 ± 0.0011 GHz for the CO(8-7), confirming the redshift estimate by Harris et al. (2012) of 3.127 whose associated uncertainties are δ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 Fig.4. Following Carilli & Walter (2013) we compute the observed magnified C[II] and CO(8-7) luminosities expressed in units of K km s −1 as: Where S line ∆v is the measured flux of the line profile (in units of Jy km s −1 ) and D L is the luminosity distance. The luminosities expressed in L are computed as L line = 3×10 −11 ν 3 rest L line . The final values computed the C[II] and the CO(8-7) lines are summarised in Table  2. 2.3. Other bands J1135 is covered by several surveys, such as the Kilo-Degree Survey (KiDS, de Jong et al. 2013) and the Hyper Suprime-Cam Subaru Strategic Program in the UV/optical bands (Aihara et al. 2018(Aihara et al. , 2022, the VISTA Kilo-degree Infrared Galaxy Public Survey (VIKING, Edge et al. 2013), and the UK Infrared Deep Sky Survey Large Area Survey (UKIDSS-LAS, Lawrence et al. 2007) surveys in the Near-IR (NIR), the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) in the MIR. PACS and SPIRE FIR observations are reported in HATLAS First and Second Data Release catalogues (Valiante et al. 2016, Maddox et al. 2017. Moreover, the source is covered by the VLA Faint Images of the Radio Sky at Twenty-Centimeters (FIRST, 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 (AO) 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 un-sufficient sensitivity and angular resolution it is not possible to unambiguously confirm whether it is associated to the foreground lens or the background source.
The object is also detected in MIR observations available in the Spitzer /IRAC Data Archive (PI: Cooray) and 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, in particular follow-ups in C-band centred ad ∼ 6 GHz (project code: 16A-240, PI: Smith). Data are processed by running the calibration scripts, cleaning is performed manually with CASA adopting an interactive mask. The final image reaches a mean rms of ∼ 0.013 mJy beam −1 and a restored beam ellipse of 1.13×0.84 arcsec (see Fig. 5 ).
The multi-band (optical-to-mid IR) image cut-outs of J1135 are reported in Figure 6. A faint emission at ∼ 4 σ emerges starting from the VIKING H-band and is detected in both IRAC channels with a S/N 6, but the angular resolution is not sufficient to resolve any lensingfeatures (e.g arcs) in the NIR/MIR regime. Flux densities are estimated by performing aperture photometry with an aperture diameter of 2 arcsec for NIR VIKING images and 6 arcsec for Spitzer /IRAC images. Table 3 summarises the photometry for J1135, we report upper limits for non-detections (i.e. emission with S/N 3 ) .

LENS MODELING AND SOURCE RECONSTRUCTION
In order to reconstruct the intrinsic background source morphology we perform lens modeling analysis with the open source Python 3.6+ code PyAutoLens (Nightingale et al. 2018(Nightingale et al. , 2021, which implements the Regularized Semi-Linear Inversion (SLI) Method described in Warren & Dye (2003) together with the adaptive source  plane pixelization scheme described in  adapted to interferometric data as done in Dye et al. (2018Dye et al. ( , 2022, Enia et al. (2018), Massardi et al. (2017), Maresca et al. (2022) and detailed in Appendix A.

Lens model
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 modelled as a Singular Isothermal Ellipsoid (SIE; Kormann et al. 1994), i.e. an elliptical power-law density distribution which goes as ρ ∝ 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 , 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 counter-clockwise from the positive x-axis, and the factor f = (1 − q)/(1 + q) where q is the ratio be-tween the semi-major and semi-minor axis. The final expressions for the elliptical components are: (2φ), (3) 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 in a first non-linear search aimed at setting priors for the lens model, this allows us to exclude regions in the parameter space corresponding to un-physical solutions where the code could get stuck. The best-fit lens model parameters are then used as priors for a second search aimed at initialising the inversion, improving the computational process. A final search is then performed to fully optimise the lensmodel parameters. The fit is performed on a number of pixels delimited by a circular mask, where the ra-  dius changes according to the resolution of the cleaned ALMA image, in order to obtain a satisfactory fit without exceeding in terms of computational cost. The output best-fit parameters and their uncertainties are reported in Table 4. Fig. 7 and Fig. 8 shows the original lens-plane image, the model image, the residual map and the reconstructed source for the three ALMA continuum bands and the CO(8-7) and C[II] emission lines respectively. Differences in the retrieved physical scales values reflect the heterogeneity of the data adopted in this work, which are the product of different array configurations and angular resolutions.
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 indication of rotation or outflow (see Fig. 9). The resulting source physical properties are reported in Table 5. The magnification factor is computed as the ratio between the lensed and unlensed surface brightness. The effective radius is computed from the area enclosing all the pixels with signal-to-noise ratio 3 and 5 in the reconstructed source plane (SP) as A SP = πr 2 eff . The resulting reconstructed source contains only pixels excluded from the masked lensed image area. This key information allows us to retrieve the intrinsic properties of the lensed background object.

The lens
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 an insufficient sensitivity and/or angular resolution. As showed 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 those 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 HST WFC3/F110 data and it is not possible to establish a priori whether it is originated by the lens or by the lensed object.
We therefore assume the lens to be 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 constrain its luminosity by means of the Einstein (total) mass resulting from the lens modeling (M E ∼ 1.3 × 10 11 M ). Specifically, we adopt the template for an elliptical galaxy with 2 Gyr age from the SWIRE library (Polletta et al. 2007). The resulting SED of this template overlapped with the photometry reported in Table 3 is showed in Fig. 10. We find the contribution from the lens to be negligible for the flux densities from the H and Ks VIKING bands up to the higher wavelengths, hence 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. We compute the intrinsic Far Infrared Luminosity (FIR), defined in the wavelength range of 8-1000 µm, by de-magnifying and fitting the FIR-to-sub-mm flux densities available for G12H43. We consider the Herschel /PACS and Herschel /SPIRE photometry from the Negrello et al. (2017) sample, we include the SCUBA-2 880µm integrated flux density reported in Bakx et al. (2018), and the flux density value measured in the 0.64, 1.04, and 1.3 mm continuum ALMA image. We use a single temperature modified black body under the optically-thin approximation, where the dust emissivity index is fixed at β = 1.5 (Nayyeri et al. 2016;Negrello et al. 2017), while the spectrum normalisation and the dust temperature (T dust ) are kept as free parameters. The model (S ν,best ) which minimises the χ 2 is then integrated over the wavelength range 8-1000 µm as follows: The best-fit spectral energy distribution is represented in Figure 11, corresponding to a dust temperature T dust = 41.1 ± 2.9 K and to a resulting far infrared luminosity of L FIR = log(L/L ) = 12.91 ± 0.01. By assuming a power law spectrum S ν = ν α with average radio spectral index α = −0.7 ± 0.14, we compute the rest-frame radio luminosity L 1.4GHz at 1.4 GHz as: where S ν,o ∝ ν α is observed monochromatic flux density at 6 GHz corrected for a magnification factor computed as the mean of the magnification factors in output from the lens modeling of the ALMA continuum emission. ν e and ν o are the emitted and the observed frequency and D L is the luminosity distance. We obtain Figure 7. Results of the lens modeling and source reconstruction procedure for continuum data. From the first column to the right: the original ALMA image, the best-fit lensed model image, the fit residuals, and the reconstructed source. The colour bar indicates the surface brightness in units of Jy arcsec −1 . From the first raw: continuum emission in bands 8, 7, and 6. Note that the y axis is inverted. log L 1.4GHz = 24.2 ± 0.2 W Hz −1 . Finally, we derive the far-infrared/radio correlation as: inferring a value of q FIR = 2.69 ± 0.41.

SED FITTING
By correcting the available photometric information for the magnification factor we can retrieve the intrinsic physical properties of J1135. To achieve this goal, we perform Spectral Energy Distribution (SED) fitting with the e Code Investigating GAlaxy Emission (CIGALE, Boquien et al. 2019). CIGALE is a Python SED fitting code able to reproduce broad-band uv-to-radio photometric data according to the energy balance (i.e the energy coming from the stellar uv-NIR emission is the same as the one re-emitted by the dust in the MIR and FIR regime). The main physical properties are estimated by comparing the observed galaxy SED with the modelled one by means of a χ 2 and bayesian statistics. We exploit the available broad-band photometry described in Sect. 2.3 and the continuum ALMA emission, including a 3σ upper limit for non-detections. For low-resolution data, we correct the flux density values for the average magnification described in the previous sectiion. As described in Sect. 3.2, we adopt the assumption that the observed photometry belongs only to the lensed source. In the following, we describe the modules adopted for the SED-fitting procedure.
The stellar emission is computed following the Bruzual & Charlot (2003)   0.05. We assume a delayed star formation history, which predicts a nearly linear increase of the SFR: where t 0 being the age of the onset of star formation, and τ the time at which the SFR peaks.  Table 3 from HCS/g to SCUBA/850 µm are represented as red points. Upper limits at 3σ are showed as arrows.
In order to model the effect of the dust attenuation on FUV-optical light we adopt the modified Charlot & Fall (2000) prescriptions, where the attenuation is agedependent and described by two different power-laws, one for the ISM and one for the Birth Clouds (BC).  Wright et al. (2010) 2 From the HATLAS Data Release 1 catalogue described in Valiante et al. (2016) 3 From the HATLAS Data Release 2 catalogue described in Maddox et al. (2017) 4 From the Herschel bright sources (HerBS) sample (Bakx et al. 2018) The attenuation slopes are assumed to be -0.7 and the V-band attenuation is computed as: In our analysis we assume A ISM V spanning from 0.3 to 5.0 and a µ spanning from 0.3 to 0.6.
Following Draine & Li (2007), dust emission is modelled as two separated components: a diffuse one, illuminated with a single radiation field (U min ) originated by a general stellar population; and a second component is closely associated to regions in which the starformation occurs, heated by a variable radiation field Figure 11. Best-fit FIR to sub-mm rest-frame SED of J1135. Red points are the observed flux densities and errors and the black line is the best-fitting modified black body spectrum. The grey shaded area represents the 68% confidence interval for the best-fit model. described with a power-law profile with index α and defined between two values U min and U max . In particular, we use the most recent and refined version of this model which accounts also for dust-mass renormalisation (Draine et al. 2014).
The resulting value of the FIRRC parameter q FIR ∼ 2.7 computed in Sect. 4, is used as a prior for the CIGALE synchrotron module to fit radio flux density at 6 GHz assuming a fixed slope α =-0.7 as in Eq. 5.
The best-fit model is presented in Fig. 12 and the resulting best physical properties are summarised in Table  6. 6. DISCUSSION Taking advantage from the SED-fitting results, we are able to investigate the ISM conditions of J1135 and its evolutionary state.

Stellar and gas masses
The bunch of available data allow us to estimate the gas content by adopting several empirical calibrators. First, we directly estimate the gas mass from the C[II] following Zanella et al. (2018), we assume α C[II] ≡ M gas /L C[II] = 22M /L , which is calibrated on starburst galaxies spanning a redshift range z ∼ 2−6. Secondly, in order to estimate the molecular gas content (M H2 ) we derive the L CO(1−0) from the de-magnified L CO(8−7) luminosity. We then follow Fujimoto et al. (2022) adopting a conversion factor of L CO(1−0) = 1.5L CO(7−6) estimated for high redshift starburst galax- Table 4. Parameters for the best-fit lens model. θE is the Einstein radius measured in arcsec, the lens positions are measured in arcsec and are referred to the centre of the ALMA observation. q and φ are the axis ratio and the positional lens angle respectively, derived from the elliptical components as described in section 3. φ is defined from the positive horizontal axis.   Figure 12. Best-fit UV to radio observed-frame SED of J1135. Green arrows are 3σ upper limits, purple circles are the observed flux densities and errors. The black line is the best-fitting modified black body spectrum. ies in literature (e.g. Riechers et al. 2013). This conversion factor is referred to a different transition, corresponding to higher luminosity values of the CO-SLED (Yang et al. 2017), for this reason the resulting value of L CO(1−0) = 1.2 × 10 10 K km s −1 pc 2 is considered as an upper limit. This estimate is consistent with the value of L CO(1−0) ∼ 1.5 × 10 10 found by Harris et al. (2012) adopting an indicative magnification factor of 10. We then compute the molecular gas mass assuming two different values of α CO = 0.8 − 4.6. The value of α = 0.8 is calibrated from local ULIRGs with super-solar metallicity (Downes & Solomon 1998), while the higher value is calibrated in the Milky Way (Solomon & Barrett 1991). The molecular gas ISM 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 19 erg s −1 Hz −1 M −1 .
Finally, we convert the dust content resulting from the SED fitting into gas assuming a variable gas-to-dust ratio of δ GDR = 30 − 92 referred to typical solar and super solar metallicity following Magdis et al. (2012) and Fujimoto et al. (2022). The values obtained for the molecular gas masses are reported in Table 7.
The stellar mass estimate in output from the SED fitting must be considered as an upper limit. Indeed given the lack of a clear detection in NIR images it is not possible to correctly estimate the contribution coming from the lens (see Section 3.2 for a further discussion). Moreover, the dark-nature of this object hinders a complete sampling of the optical and NIR part of the SED. Aside from the value reported in Table 6, we compute the stellar mass assuming a typical stellar-to-dust mass ratio of δ SDR ≈ 100, obtaining a value of M STD * ∼ 2 × 10 11 M , in agreement with the SED fitting estimate.

ISM properties
From the gas mass values reported in Sect. 6.1, we estimate a depletion timescale of τ depl ∼ 10 8 yr. Moreover, the inferred stellar mass implies τ SFR ∼ 10 8 yr, 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 to the presence of young and massive stars ionising the surrounding HII regions. The same authors found no significant AGN contribution from the SED analysis, consistent with what we infer from the FIRRC (q FIR ≈ 2.7), which is indicative of a star-formation dominated object.
The hypothesis of J1135 being a compact starburst is also supported by the source reconstruction of the highest angular resolution ALMA continuum emission at 640 µm shown in Figure 7, where the effective radius reaches values of ∼ 400 pc. Similar physical scales are reached by the C[II] line emission (see Table 5 and Fig.  8). The C[II] is a fine structure line predominantly originated from high−z photon-dominated regions (PDR) and is typically used as a cool interstellar gas tracer and as a 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, Díaz-Santos et al. 2017, Smith et al. 2017) and highredshift star-forming galaxies (Stacey et al. 2010, Gullberg et al. 2015. This drop is found to reach very low values (L C[II] /L FIR ≈ 10 −4 ) in spatially resolved studies (e.g. Lagache et al. 2018, Gullberg et al. 2015, Rybak et al. 2019. For J1135, we infer a C[II]/FIR ratio of L C[II] /L FIR ≈ 5.4 × 10 −4 . Similar values are found for other strongly lensed galaxies among the HATLAS sample. For example, Rybak et al. (2020) reported a deficit down to ∼ 3 × 10 −4 for spatially resolved ALMA data of SDP.81 (Partnership et al. 2015, Rybak et al. 2015a,b, Dye et al. 2015, Hatsukade et al. 2015, Hezaveh et al. 2016) at z = 3.042. Lamarche et al. (2018) found similar values (∼ 2 × 10 −4 ) 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. From R eff,640µm we infer a star-formation surface density of ρ SFR ∼ 1600 M yr −1 kpc −2 , which is on the verge of the Eddington limit for a radiation pressure supported starburst (Andrews & Thompson 2011, Simpson et al. 2015. This value, is compatible with the possible explanation of the deficit to be attributed to a lower increase of the C[II] emission with respect to the FIR.

Evolutionary interpretation
By inspecting the HST/WFC3 image we find no evidence for galaxy companions of J1135 within a radius of at least ∼ 5 arcsec, corresponding to ∼ 40 kpc, so that we can exclude a merger-induced origin of the starburst. Thus the ISM conditions and the physical properties discussed so far can be interpreted in the light of in-situ galaxy formation scenarios (Lapi et al. 2014, 2018, Mancuso et al. 2017, Pantoni et al. 2019. In particular, the properties of J1135 are consistent with a compaction phase (see Fig. 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 to be an off-main sequence object in a early evolutionary stage, which will eventually move towards the main-sequence locus as the stellar mass content increases. Finally, the star formation will either progressively decrease as the galaxy exhaust its gas reservoir or will be abruptly stopped by the action of the feedback from an AGN (Mancuso et al. 2017).

SUMMARY AND CONCLUSIONS
In this work we have investigated the nature of the strongly-lensed galaxy HATLASJ113526.2-01460 (namely, J1135) at redshift z ≈ 3.1, discovered by the Herschel satellite in the GAMA 12 th field of the Herschel -ATLAS survey. We have performed detailed lens modeling and have reconstructed the source morphology in three different (sub-)mm continuum bands, and in the spectral emission of the C[II] and CO(8-7) lines. We have also exploited a wealth of photometric ancillary data to perform broadband SED-fitting and to retrieve intrinsic (i.e., corrected for magnification) physical properties. Our main findings are summarized below: • The lens modeling indicates that the foreground lens is constituted by a (likely elliptical) galaxy with mass 10 11 M at z 1.5, while the source is found to be an optical/NIR dark, dusty starforming galaxy whose (sub-)mm continuum and line emissions are amplified by factors µ ∼ 6 − 13.
• The emission of J1135 is extremely compact, with sizes 0.5 kpc for the star-forming region and 1 kpc for the gas component, with no clear evidence of rotation or of ongoing merging events.
• J1135 features a very high star-formation rate 10 3 M yr −1 , that given the compact sizes is on the verge of the Eddington limit for starbursts. The radio luminosity at 6 cm from available EVLA observations is consistent with the star-formation activity, so that no significant contribution from a central AGN is emerging (see also Vishwas et al. 2018).
• J1135 is found to be extremely rich in gas ∼ 10 11 M and dust 10 9 M . The stellar content 10 11 M places J1135 well above the main sequence of starforming galaxies, indicating that the starburst is rather young with an estimated age ∼ 10 8 yr, and that the stellar mass should at least double before star formation is quenched.
• The properties of J1135 can be consistently explained in terms of in-situ galaxy formation and evolution scenarios as typical of a rather young dusty starforming galaxy caught in the compaction phase.
In the next future, observations coming from the James Webb Space Telescope (JWST) will be crucial to shed further light on the nature of this obscured object and its foreground lens in the near-and mid-IR regime. Moreover, X-Ray follow-up, coupled with the available ALMA data, are required to establish the presence of the dust-enshrouded AGN and to better investigate the interplay between star-formation and the nuclear activity . The SLI formalism can be extended also to interferometry , Enia et al. 2018, Maresca et al. 2022, modeling a set of visibility data, i.e. the result of the correlation of signals coming from an astrophysical source and collected by the antennae array, whose Fourier transform gives the source surface brightness distribution. Performing an inversion directly on the Fourier space (or uv -plane) circumvents the issue of dealing with artifacts and noise correlation arising in the image as a consequence of a poor sampling of the uv -plane.
Following a similar formalism with respect to the one used in Dye et al. (2018), we introduce the rectangular matrix f ij containing the fluxes of the i-th pixel in the source plane and the respective j-th image-plane pixel. Analogously, complex visibilities from the lensed image are collected rectangular matrix g ij , which are the Fourier transform of the i source pixels in unit surface brightness computed at the j-th visibility point in the uv -plane. For each j-th visibility corresponding to the source pixel surface brightnesses s i , the model visibility set can be described as i s i g ij .
Given a set of observed visibilities V obs , the merit function can be described as: computed over a total of I Delaunay pixels and J visibilities. σ j are the 1σ uncertainties on the observed visibilities rescaled adopting the CASA task statwt to match their absolute value. The last term in the expression describes the regularization, where λ is a constant determining the strength of the regularization, and H the regularization matrix. The values s i , represented by the vector S which best reproduces the observed image-plane visibilities, can therefore be derived minimizing the merit function G. The solution to this linear problem is given by: where F and D are respectively the matrices F ij = J n=1 (g R in g R jn + g I in g I jn /σ 2 n ) and D i = J n=1 (g R in V R n + g I in V I n /σ 2 n ). For a fixed mass model, the image plane pixels are traced back to the source plane and grouped together by means of a k-clustering algorithm, comparing each source-pixel with the neighbors sharing a direct vertex. This procedure results in new source plane's centres, used to trace a Delaunay grid. When dealing with a large number of visibilities, the computational efficiency and the memory costs are greatly improved by performing non-uniform Fast Fourier Transform (NUFFT) algorithm, implemented in PyAutoLens exploiting the PyNUFFT (Lin 2018) library and the linear algebra package PyLops (Ravasi & Vasconcelos 2020).