# ALMA 1.3 mm Survey of Lensed Submillimeter Galaxies (SMGs) Selected by
Herschel: 111Herschel is an ESA space observatory with science instruments
provided by European-led Principal Investigator consortia and with important
participation from NASA.
Discovery of Spatially Extended SMGs and Implications

Fengwu Sun  Steward Observatory, University of Arizona, 933 N. Cherry Avenue,
Tucson, 85721, USA  Eiichi Egami  Steward Observatory, University of Arizona,
933 N. Cherry Avenue, Tucson, 85721, USA  Timothy D. Rawle  ESA, Science
Operations Department, STScI, Baltimore, MD 21218, USA  Gregory L. Walth  The
Observatories of the Carnegie Institution for Science, 813 Santa Barbara
Street, Pasadena, CA 91101, USA  Steward Observatory, University of Arizona,
933 N. Cherry Avenue, Tucson, 85721, USA  Ian Smail  Centre for Extragalactic
Astronomy, Department of Physics, Durham University, South Road, Durham, DH1
3LE, UK  Miroslava Dessauges-Zavadsky  Observatoire de Genève, Université de
Genève, 51, Ch. des Maillettes, 1290 Versoix, Switzerland  Pablo G. Pérez-
González  Centro de Astrobiología, Departamento de Astrofísica, CSIC-INTA,
Cra. de Ajalvir km.4 E-28850–Torrejón de Ardoz, Madrid, Spain  Johan Richard
Université Lyon, Université Lyon1, ENS de Lyon, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, Saint-Genis-Laval, France.  Francoise Combes
Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, 75014
Paris, France  Collège de France, 11 Place Marcelin Berthelot, 75231 Paris,
France  Harald Ebeling  Institute for Astronomy, University of Hawaii, 2680
Woodlawn Drive, Honolulu, HI 96822, USA  Roser Pelló  Aix Marseille
Université, CNRS, CNES, LAM (Laboratoire d’Astrophysique de Marseille), UMR
7326, 13388, Marseille, France  Paul Van der Werf  Leiden Observatory, Leiden
University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands  Bruno Altieri
European Space Astronomy Centre, ESA, Villanueva de la Cañada, 28691, Madrid,
Spain  Frédéric Boone  CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028,
Toulouse Cedex 4, France  Antonio Cava  Department of Astronomy, University of
Geneva, 51 Ch. des Maillettes, 1290 Versoix, Switzerland  Scott C. Chapman
Department of Physics and Astronomy, University of British Columbia, 6225
Agricultural Road, Vancouver, BC V6T 1Z1, Canada  National Research Council,
Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E
2E7, Canada  Department of Physics and Atmospheric Science, Dalhousie
University, Halifax, NS B3H 4R2, Canada  Benjamin Clément  Institute of
Physics, Laboratory of Astrophysics, Ecole Polytechnique Fdrale de Lausanne
(EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland  Alexis Finoguenov
Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2,
00014 Helsinki, Finland  Kimihiko Nakajima  National Astronomical Observatory
of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Wiphu Rujopakarn
Department of Physics, Faculty of Science, Chulalongkorn University, 254
Phayathai Road, Pathumwan, Bangkok 10330, Thailand  National Astronomical
Research Institute of Thailand (Public Organization), Don Kaeo, Mae Rim,
Chiang Mai 50180, Thailand  Kavli Institute for the Physics and Mathematics of
the Universe (WPI),The University of Tokyo Institutes for Advanced Study, The
University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Daniel Schaerer
Observatoire de Genève, Université de Genève, 51, Ch. des Maillettes, 1290
Versoix, Switzerland  Ivan Valtchanov  European Space Astronomy Centre, ESA,
Villanueva de la Cañada, 28691, Madrid, Spain

(Received September 18, 2020; Revised December 9, 2020; Accepted December 24,
2020)

###### Abstract

We present an ALMA 1.3 mm (Band 6) continuum survey of lensed submillimeter
galaxies (SMGs) at $z=1.0\sim 3.2$ with an angular resolution of $\sim 0\farcs
2$. These galaxies were uncovered by the Herschel Lensing Survey (HLS), and
feature exceptionally bright far-infrared continuum emission
($S_{\mathrm{peak}}\gtrsim 90$ mJy) owing to their lensing magnification. We
detect 29 sources in 20 fields of massive galaxy clusters with ALMA. Using
both the Spitzer/IRAC (3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m) and ALMA data,
we have successfully modeled the surface brightness profiles of 26 sources in
the rest-frame near- and far-infrared. Similar to previous studies, we find
the median dust-to-stellar continuum size ratio to be small
($R_{\mathrm{e,dust}}$/$R_{\mathrm{e,star}}$ = 0.38$\pm$0.14) for the observed
SMGs, indicating that star formation is centrally concentrated. This is,
however, not the case for two spatially extended main-sequence SMGs with a low
surface brightness at 1.3 mm ($\lesssim 0.1$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$),
in which the star formation is distributed over the entire galaxy
($R_{\mathrm{e,dust}}$/$R_{\mathrm{e,star}}$ $>$ 1). As a whole, our SMG
sample shows a tight anti-correlation between
($R_{\mathrm{e,dust}}$/$R_{\mathrm{e,star}}$) and far-infrared surface
brightness ($\Sigma_{\mathrm{IR}}$) over a factor of $\simeq$ 1000 in
$\Sigma_{\mathrm{IR}}$. This indicates that SMGs with less vigorous star
formation (i.e., lower $\Sigma_{\mathrm{IR}}$) lack central starburst and are
likely to retain a broader spatial distribution of star formation over the
whole galaxies (i.e., larger $R_{\mathrm{e,dust}}$/$R_{\mathrm{e,star}}$). The
same trend can be reproduced with cosmological simulations as a result of
central starburst and potentially subsequent “inside-out” quenching, which
likely accounts for the emergence of compact quiescent galaxies at $z\sim 2$.

submillimeter: galaxies — galaxies: evolution — galaxies: high-redshift

††journal: ApJ††facilities: ALMA, Spitzer (IRAC), Herschel (PACS and SPIRE),
HST (ACS and WFC3-IR)††software:  astropy (Astropy Collaboration et al., 2013;
Price-Whelan et al., 2018), SExtractor (Bertin & Arnouts, 1996), CASA
(McMullin et al., 2007), emcee (Foreman-Mackey et al., 2013), MAGPHYS (da
Cunha et al., 2008; Battisti et al., 2019), GALFIT (Peng et al., 2010)

tablenum SIXtablenum

##  1 Introduction

As the most vigorous stellar nursery in the Universe, submillimeter galaxies
(SMGs222Observed 1.3 mm flux density at $\gtrsim 1\,$mJy in this work. See
discussion of definition in Hodge & da Cunha (2020).) are discovered in
abundance at $z>1$, contributing $\sim 20\%$ of the cosmic star formation rate
density up to $z\sim 4$ (e.g., Swinbank et al., 2014; Casey et al., 2014a).
Due to the critical role played by dust grains in the interstellar medium
(ISM), commonly produced by asymptotic giant branch (AGB) stars and supernovae
(SNe), these galaxies are observed to be highly dust-obscured in the rest-
frame UV/optical bands (e.g., Whitaker et al., 2017; Dudzevičiūtė et al.,
2020). In the far-infrared (FIR), dust grains, heated up by the intense star
formation, emit thermal continuum radiation accessible through
submillimeter/millimeter observations. These galaxies are found to host
massive gas reservoirs with relatively short gas depletion time scales
($t_{\mathrm{dep}}\sim 10^{2}$ Myr; e.g., Tacconi et al. 2008; Bothwell et al.
2013; Miettinen et al. 2017). After the truncation of sufficient gas supply,
SMGs are believed to evolve towards compact quiescent galaxies (cQs) seen at
lower redshift (e.g., Toft et al., 2014; Simpson et al., 2014), which will
eventually become massive elliptical galaxies in the local Universe
potentially through additional gas-poor mergers (e.g., van Dokkum, 2005; Oogi
& Habe, 2013).

The trigger mechanism of SMGs remains a subject of debate. As the local
analogs of SMGs, ultra-luminous infrared galaxies (ULIRGs) are observed to
result from major mergers, with compact and prominent star-forming regions in
their nuclei (e.g., Sanders et al., 1988). At $z>1$, similar scenarios have
been proposed by certain galaxy evolution theories (e.g., Narayanan et al.,
2010; McAlpine et al., 2019), but minor mergers (e.g., Gómez-Guijarro et al.,
2018) and secular starburst (e.g., Davé et al., 2010) are other possible
physical explanations supported by either observational evidence or
theoretical frameworks.

Most recently, powerful ground-based interferometers like ALMA have started to
reveal the compact dust continua of SMGs by high-resolution imaging (e.g.,
Simpson et al., 2015; Hodge et al., 2016, 2019; Elbaz et al., 2018; Puglisi et
al., 2019; Gullberg et al., 2019; Tadaki et al., 2020). With a typical half-
light radius of $\sim$1–2 kpc, these intense star-forming regions are still
larger than those in the local (U)LIRGs, while their sizes do match with cQs
at slightly lower redshift (e.g., van Dokkum et al., 2008; van der Wel et al.,
2014). In addition, the number density and clustering properties of SMGs also
coincide with those of cQs, indicating an underlying evolutionary connection
(Hickox et al., 2012; Simpson et al., 2014; An et al., 2019; Dudzevičiūtė et
al., 2020). Therefore, understanding how star formation ceases in SMGs (i.e.,
quenching) holds important clues to connecting actively star-forming galaxies
and their (possible) red and dead descendants.

Spatially resolved studies have suggested that the massive star-forming
galaxies at $z\simeq 1-2$ exhibit a rising specific star-formation rate (sSFR;
SFR per unit stellar mass) profile from their center to the outskirts,
indicating that the fade-out of star formation commences from the galactic
center (e.g., Tacchella et al., 2015, 2018; Nelson et al., 2016; Spilker et
al., 2019). Such an “inside-out” process of quenching can also be reproduced
within the context of cosmological simulation (Tacchella et al., 2016).
Various physical mechanisms have been proposed to interpret the quenching
process, including the gas consumption by star formation, gas outflows driven
by stellar and supermassive black hole feedbacks (Dekel & Silk, 1986; Di
Matteo et al., 2005), as well as suppression of exterior gas supply through
shocking heating due to gravitationally infalling gas (Dekel & Birnboim, 2006,
2008).

Although emerging observational clues suggest an inside-out fashion of stellar
mass assembly and gas depletion, it is yet to be confirmed as the standard
process in the evolution of SMGs with compact and powerful star-forming
regions. The morphological modeling of SFR and stellar mass profile requires
high-resolution imaging of SFR tracers (e.g., H$\mathrm{\alpha}$, UV/FIR
continuum) and stellar component (the rest-frame optical/near-infrared
continuum). However, the accuracy of measurements at shorter wavelength is
clearly subject to the strong dust extinction in the center of SMGs (e.g.,
Simpson et al., 2017; Lang et al., 2019). Meanwhile, H$\mathrm{\alpha}$-based
SFR profile is also sensitive to the active galactic nuclei (AGN)
contribution, since the AGN fraction is higher in compact star-forming
galaxies compared to more extended ones at similar redshifts (Barro et al.,
2013). In a nutshell, it is necessary to develop novel modeling techniques for
stellar mass and SFR that are less sensitive to dust extinction and AGN
contribution.

To address this issue, one possible solution is to observe at longer
wavelength. Since most of the star formation in SMGs is obscured by dust, FIR
surface luminosity can represent the surface SFR with a sufficient accuracy.
In order to avoid the heavy dust obscuration of stellar continuum in SMGs, it
is also better to sample the rest-frame near-infrared (NIR) bands rather than
the optical ones (e.g., HST/WFC3-IR F160W samples the rest-frame $\sim$ $V$
band for $z=2$ galaxies). However, current sensitive mid-infrared imaging
instrument (e.g., Spitzer/IRAC) cannot allow such a study with ALMA-like
angular resolution as required.

In this regard, cluster-lensed SMGs (e.g., Smail et al., 1997; Swinbank et
al., 2010) can be useful targets to provide the morphological evidence.
Magnified by gravitational lensing provided by a foreground cluster, these
targets are sufficiently bright at multiple wavelengths, and their angular
sizes are also stretched significantly (note, however, that the lensing effect
conserves surface brightness). Compared to galaxy-lensed SMGs with a bright
lensing galaxy always in the front, cluster-lensed ones are often free from
blending with foreground galaxies, ensuring simplicity of the morphological
modeling on their stellar component. Strong magnification gradients on the
$\simeq$kpc scales relevant for resolving galaxies are also much less a
concern than in galaxy-lensed cases. During the course of our Herschel Lensing
Survey (HLS, Egami et al. 2010; Egami et al., in prep.), we uncovered a
substantial number of lensed SMGs with exceptionally bright FIR continuum
($S_{\mathrm{peak}}\gtrsim 90\,$mJy). We then carried out observations of
their dust continua in the ALMA Band 6 at 1.3 mm, as well as stellar continua
with Spitzer/IRAC at 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m. In this work, we
present the observations and analyses of the stellar and dust components in
cluster-lensed SMGs at both integrated and spatially resolved scales.

This paper is arranged as follows: Section 2 introduces all the obtained ALMA
and Spitzer/IRAC data and corresponding data reduction techniques, with
several ancillary data from various sources. Section 3 presents the
fundamental analysis of our data, including detection, photometry and surface
brightness profile modeling. In Section 4 we perform spectral energy
distribution (SED) modeling and show the statistical results of galaxy
properties from both the integrated and spatially-resolved analyses. We
discuss the underlying physics and make necessary comparison with both
observational evidence and theoretical predictions in Section 5. The
conclusions and broader implication can be found in Section 6. Throughout this
paper, we assume a flat $\Lambda$CDM cosmology with $h=0.7$ and
$\Omega_{m}=0.3$. The AB magnitude system (Oke & Gunn, 1983) is used to
express source brightnesses in the optical and NIR.

##  2 Observations and Data

###  2.1 The Sample

To discover and study a significant sample of gravitationally lensed SMGs, we
have conducted an extensive imaging survey of massive galaxy clusters in the
FIR using the Herschel Space Observatory (Pilbratt et al., 2010), known as the
HLS (Egami et al. 2010). The target clusters were selected mainly from the
samples produced by the following three surveys: (1) the ROSAT All-Sky Survey
(RASS) with the X-ray-luminous cluster sample tabulated by H. Ebeling (private
communication), (2) the COnstrain Dark Energy with X-ray (CODEX) survey, which
utilizes the combination of the RASS X-ray and Sloan Digital Sky Survey (SDSS)
optical data (Finoguenov et al., 2020), and (3) the South Pole Telescope (SPT)
survey, which selected clusters via the Sunyaev-Zel’dovich (SZ) effect (the
SPT-SZ survey; Bleem et al. 2015). The full HLS cluster sample will be
presented and described by the forthcoming survey paper (Egami et al., in
prep.).

With a substantial number of SMG detections at $z_{\mathrm{phot}}\gtrsim 1$,
we specifically selected a subset of exceptionally bright sources and obtained
ALMA follow-up observations. The selection criteria used were: (1) SPIRE color
$S_{500}/S_{250}>0.4$ to ensure the selection of $z\gtrsim 1$ sources, (2) FIR
continuum peak ($S_{\mathrm{peak}}$) brighter than 90 mJy if the source is
within 1′ from the cluster center, (3) $S_{\mathrm{peak}}>$150 mJy or with a
spectroscopic redshift if the source is beyond 1′ from the cluster center, and
(4) observable with ALMA (Dec. $<+30^{\circ}{}$). The second criterion is
designed to select bright sources at $L_{\mathrm{IR}}\gtrsim 10^{13.1}$
$\mu^{-1}\,\mathrm{L}_{\odot}$ at $z\sim 2$. The third criterion ensures that
the resultant sample includes the brightest sources in the HLS data even if
some of them may be boosted by a galaxy component on the line of sight. Due to
the coarse resolution of SPIRE, the $S_{\mathrm{peak}}$ quoted here is the sum
of all sub-mm sources within a radius of $\sim 15$′′. Multiple source systems
will be decomposed later using the ALMA data, and these individual sources
will not necessarily satisfy the same Herschel selection criteria.

We eventually constructed a sample of 20 sources based on the criteria listed
above.333One source in the cluster field SPT J0345-6419 was observed with ALMA
but was later identified as a low-redshift IR-bright galaxy rather than a
$z\gtrsim 1$ SMG. We therefore removed it from our sample. There are also
several other HLS sources satisfying the same brightness criteria, but we do
not incorporate them here since they were not observed by ALMA or lie at
different redshift range ($z\sim 5$; e.g., HLS0918, Combes et al. 2012, Rawle
et al. 2014; HLS0257, Sun et al. in prep.; HLS2043, Zavala et al. 2015; Walth
et al. in prep.).

###  2.2 Herschel

The HLS has performed far-infrared imaging observations of 581 massive galaxy
clusters with a total observing time of 418.7 hours at two typical depths. The
HLS-deep survey imaged 54 clusters deeply with PACS (Poglitsch et al., 2010)
and SPIRE (Griffin et al., 2010) through an open-time key program (44 targets;
PI: Egami; PID: KPOT_eegami_1) and an open-time Cycle 2 program (10 targets;
PI: Egami; PID: OT2_eegami_5). Three of the clusters studied here (MACS
J1115.8+0129, Abell 2813 and Abell 3088) were observed by HLS-deep, and
therefore have a five-band coverage with both PACS (100/160
$\mathrm{\SIUnitSymbolMicro}$m) and SPIRE (250/350/500
$\mathrm{\SIUnitSymbolMicro}$m). All of the PACS 100 and 160
$\mathrm{\SIUnitSymbolMicro}$m observations consist of two orthogonal scan
maps, each comprising 18–22 repetitions of 13 parallel 4-arcmin scan legs. The
SPIRE observations for the two Abell clusters were performed with 20
repetitions in the large scan map mode, each with two 4′ scans and cross-scans
(1.6 h scan for three bands simultaneously). MACS1115 was observed through
11-repetition small scan maps, and each repetition consisted of one scan and
one cross-scan of 4′ length (0.4 h scan).

The HLS-snapshot survey obtained shallower SPIRE-only data for 527 clusters
through two open-time programs during Cycle 1 and 2 (PI: Egami; PID:
OT1_eegami_4, OT2_eegami_6), providing SPIRE data for the remaining 17
clusters studied here. With shallow observations ($3-8$ min scan) in the small
scan map mode, HLS-snapshot provides nearly confusion-limited images in all
three SPIRE bands (typical RMS noise $\sim$ 10
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}^{-1}$
at 250, 350 and 500 $\mathrm{\SIUnitSymbolMicro}$m, compared with the
confusion noise levels of 5.8, 6.3 and 6.8
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}^{-1}$
measured by Nguyen et al. 2010).

Our SPIRE images were produced via the standard reduction pipeline in HIPE
v12.2 (Ott, 2010), and the processing routine was detailed in Rawle et al.
(2016) for Herschel coverage of the HST Frontier Fields (HFF, Lotz et al.
2017). The observation ID (OBSID) and total scan time of each obtained
Herschel/SPIRE observation are summarized in Table 1. HLS-deep PACS images
were generated with UniMap (Piazzo et al., 2015) with a pixel scale of
1$\farcs$0 at 100 $\mathrm{\SIUnitSymbolMicro}$m and 2$\farcs$0 at 160
$\mathrm{\SIUnitSymbolMicro}$m, also detailed in Rawle et al. (2016).

###  2.3 ALMA

ALMA Band 6 observations were carried out through project 2015.1.01548,
2016.1.00372 and 2017.1.01658 (PI: Egami) between April 30, 2016 and September
30, 2018. Since four sources exhibit extended structures in the SPIRE images
(FWHM$>$20′′ at 250 $\mathrm{\SIUnitSymbolMicro}$m), we requested multi-
pointing observations for these special cases. We observed all of our 20
targets in one of two spectral window settings. For 15 sources without
previous spectroscopic redshift determination, we performed continuum-only
observations with a central frequency at 233 GHz (corresponding to 1.287 mm).
For 5 sources with prior spectroscopic redshift information, we acquired both
dust continuum and at least one CO line spectrum, and thus the final effective
frequencies of these continuum products range from 224 to 238 GHz. The
diameter of the ALMA field of view (FoV) at the requested frequencies is 25′′.
A brief summary of our ALMA observations is also presented in Table 1.

The data were taken in various weather conditions with a median precipitable
water vapor (PWV) of 0.77 mm, with the 16th to 84th percentile ranging from
0.56 to 1.84 mm. A median angular resolution of 0$\farcs$26 was achieved, and
the median maximum recoverable scale was 1$\farcs$7.

All the ALMA data were reduced with CASA (McMullin et al., 2007) with the
pipelines v4.7.2 and v5.4.1 for observations obtained in different cycles.
Before the formal reduction work, we first checked the combined continuum
image and spectral cube of each target, delivered by the ALMA archive. If any
obvious source was detected above a 4$\sigma$ significance, we would examine
the spectral cubes, searching for possible spectral line features, and both
line and continuum would be imaged in the natural/Briggs weighting and uv-
tapered modes separately. If undetected, we would only produce the continuum
images. We performed continuum imaging at four different levels of synthesized
beam size: Briggs weighting (robust=0.5), natural weighting (robust=2), 1′′
uv-tapering and 2′′ uv-tapering. These settings were used for visualizing both
compact and extended emission structures in the SMGs. Interactive cleaning was
performed during each imaging process. The noise level of the final continuum
products is $0.11^{+0.04}_{-0.05}$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}^{-1}$
with a 1′′-tapered beam, which we used to obtain photometry for most of the
targets.

###  2.4 Spitzer/IRAC

We obtained Spitzer/IRAC Channel 1 (3.6 $\mathrm{\SIUnitSymbolMicro}$m) and
Channel 2 (4.5 $\mathrm{\SIUnitSymbolMicro}$m) images through various
programs. The majority of our data were from Program 12095 and 90218 (both PI:
Egami), covering 13 of 20 clusters. We also included other data with public
access on the Spitzer Heritage Archive444https://sha.ipac.caltech.edu/. This
includes Programs 80168, 90213 (PI: Bouwens), 12005, 14281 (PI: Bradac),
60099, 80012 (PI: Brodwin), 60034 (PI: Egami), 30344 (PI: Jarvis), 80162,
90233 (PI: Lawrence), 70149 (PI: Menanteau), 80066 (PI: Rawle), 61061 (PI:
Sheth), 12123 (PI: Soifer), 40370, 80096 (PI: Stanford), 14061, 60194 (PI:
Vieira). These programs provide 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m
coverages for all the 20 clusters. All of these fields are observed with
cycling sub-pixel dithering patterns with four or five dithering points at
least.

We started our IRAC data processing from the archive-delivered level 1 (BCD)
products. A standard and automatic MOPEX reduction routine was applied with an
output pixel size of 0$\farcs$6 pixel-1. We registered the output frames with
the GAIA DR2 (Gaia Collaboration et al., 2018) using SExtractor (for catalog
extraction; Bertin & Arnouts 1996) and SCAMP (for astrometric computation;
Bertin 2006). This achieved a final astrometric error of $\lesssim 0\farcs 1$
with the produced IRAC images. The median $5\sigma$ point-source depth in each
field, estimated from the variance of sky background, is presented in Table 1.

###  2.5 Other ancillary data

For a number of clusters, we also included other ancillary data to improve the
quality of analysis, mainly for the optical SED fitting. This provides more
accurate extinction and stellar mass estimates, compared with IRAC-only
analyses:

MACS J0553.4-3342— this cluster was observed with the HST treasury program
RELICS (Coe et al., 2019), and therefore we used its 7-band HST data for
photometry and SED fitting (ACS/F435W, F606W, F814W and WFC3-IR/F110W, F125W,
F140W, F160W).

MACS J1115.8+0129— this cluster was observed with the HST treasury program
CLASH (Postman et al., 2012), and therefore its 9-band HST data were utilized
(ACS/F435W, F606W, F775W, F814W, WFC3-IR/F105W, F110W, F125W, F140W, F160W).
We directly used the processed data of MACSJ0553 and MACSJ1115 available on
MAST555Mikulski Archive for Space Telescopes (MAST),
https://archive.stsci.edu/.

RXC J2332.4-5358— this cluster was observed through the HST SNAP program 12884
(PI: Ebeling). The WFC3-IR/F110W and F140W images were obtained with an
integration time of 706 s per filter. We reduced the data with Drizzlepac
v2.1.17 (Gonzaga & et al., 2012) under the PyRAF environment with an output
pixel size of 0$\farcs$06 pixel-1.

Aperture photometry of SMGs in HST images is conducted with SExtractor. We do
not obtain flux measurement of HLS0553-C (blended with a star $\sim
5^{\mathrm{m}}$ brighter than the SMG in the F814W band) and HLS2332-C
(blended with an irregular galaxy).

RXC J1314.3-2515— This cluster was observed with WFCAM on UKIRT in both the
$J$ and $K$ bands (PI: Walth), and here we only use the NIR photometry of the
lensed SMG in this cluster field.

No additional optical/NIR data was included for the analysis of remaining
sources. The Herschel sources in RXC J1314.3-2515, MACS J0455.2+0657 and MACS
J0600.1-2008 were observed with JCMT/SCUBA-2 at 850
$\mathrm{\SIUnitSymbolMicro}$m (Cheale et al., 2019), and here we quote the
850 $\mathrm{\SIUnitSymbolMicro}$m flux densities to improve the quality of
far-IR SED modeling.

##  3 Results

###  3.1 ALMA detection and photometry

To obtain reliable and complete ($\gtrsim 1$ mJy) ALMA detections of lensed
SMGs in all 20 observed cluster fields, we used SExtractor v2.19.5 (Bertin &
Arnouts, 1996) to derive uniform source extraction in all ALMA continuum image
products (without primary beam correction). Based on a quick visual inspection
of natural-weighted image products without uv-tapering (median beam FWHM
$=0\farcs 28$), we found that all of the obvious sources detected in our data
were spatially resolved. Since uv-tapering can increase the detectability of
extended structures, though at the expense of resolution, we performed the
source extraction with the 1′′-tapered image products.

We ran SExtractor for source detection and automatic photometry with Kron-like
elliptical apertures (PHOT_AUTOPARAMS values of 1.8 and 2.5) in the primary-
beam-uncorrected maps where the primary beam response is greater than 0.2. We
estimated the photometric errors ($\sigma_{\mathrm{aper}}$) based on aperture-
to-beam size ratio ($\Omega_{\mathrm{aper}}/\Omega_{\mathrm{beam}}$) and
continuum RMS ($\sigma_{\mathrm{RMS}}$) according to the following equation:

$\sigma_{\mathrm{aper}}={\sigma_{\mathrm{RMS}}}\cdot\sqrt{\frac{(\Omega_{\mathrm{aper}}/\Omega_{\mathrm{beam}})^{2}}{1+(\Omega_{\mathrm{aper}}/\Omega_{\mathrm{beam}})/2}}$
(1)

which is derived from a simulation of applying random apertures (enclosed area
as $\Omega_{\mathrm{aper}}$) on Gaussian-blurred (kernel area as
$\Omega_{\mathrm{beam}}$) Gaussian white noise maps. We also evaluated the
photometric uncertainty by directly applying random apertures on source-free
regions in ALMA primary-beam-uncorrected maps. We measured the standard
deviation of flux densities enclosed within the apertures of identical size,
and the results are consistent with the prediction by Equation 1. We then
corrected the flux densities and their errors for the gain of primary beam.

Based on this method, we detected 77 sources at S/N$>$3.0 in all ALMA images.
Here the S/N is defined as the ratio between the aperture-photometry flux
density and its uncertainty ($f/\sigma_{\mathrm{aper}}$). To eliminate
possible false detections, we studied the false detection rate in Appendix A,
and found that an S/N cut at 4.0 would ensure the total false detection number
to be $\lesssim 1$. We therefore detected 28 sources at S/N$>$4.0 based on the
1′′-tapered images. All of the sources were detected in the area where the
primary beam response is greater than 0.5, except for HLS2155-A.

We also included another S/N$>$4.0 detection in 2′′-tapered images, namely
HLS0840. This source is exceptionally extended with an ALMA-measured 1.3 mm
effective radius $R_{\mathrm{e,ALMA}}$ of $3\farcs 7\pm 1\farcs 7$ with a
relatively low surface brightness. It is split into multiple S/N$\sim$3
components in our 1′′-tapered image but remains as a single S/N=4.5 source in
the 2′′-tapered map. We also detected a point-like $0.5\pm 0.1$ mJy source at
its center in the 0$\farcs$2-resolution image. This might represent the core
of the galaxy, although it only contributes $13\pm 4$% of the total flux
density in the 2′′-tapered map. We included this source because of the
robustness of the detection in the 2′′-tapered images, and no other similar
example was found in our data.

We also studied the completeness of our source extraction in Appendix A. We
conclude that the completeness of point-like sources at S/N=5.0 in 1′′-tapered
map (0.79 mJy, under the assumption of a median continuum RMS and a primary
beam correction of $\sqrt{2}$) is 80$\pm$4%.

In one cluster field, CODEX 52909, we did not detect any significant source.
The SPIRE source in CODEX 52909 is extended, and it can be a composite of
several ALMA sources at S/N$\simeq$3–3.5 with reddened IRAC counterparts. To
avoid any confusion of $\sim 3\sigma$ fake sources in this mosaic ALMA FoV, we
did not analyze this field any further.

Because of the noise fluctuation in the ALMA maps, the fluxes of sources at
low S/N tend to be overestimated, known as flux boosting effect (e.g., Geach
et al., 2017; Stach et al., 2019). We examined this effect for sources
detected at relatively low S/N (HLS0546, HLS0840, HLS0043-B and HLS2155-A;
Table 2). Assuming the surface brightness profiles measured in Section 3.3, we
simulated the visibility data for these sources 10–15 times per each with
CASA, and the RMS noise of mock data was controlled to match with that of our
observations. We then applied the same imaging and source extraction routine,
measured the median output flux densities and compared them with those of
input models. Based on our simulations, no conspicuous flux boosting effect
was identified.

###  3.2 Multi-wavelength photometry

Figure 1: Postage stamp images of six cluster-lensed SMGs in this work, and
the remaining ones are shown in Figure 17. The galaxies are at the center of
each image. For each source, we show their ALMA 1.3 mm Briggs-weighted map in
the first column (resolution $\Delta\theta\sim 0\farcs 2$), and uv-tapered map
in the second column ($\Delta\theta\sim 1$′′; $2$′′-tapered for HLS0546 and
HLS0840). The synthesized ALMA beam is shown as a hatched ellipse at the
lower-left corner, and the source sizes and fluxes are noted in the tapered
images. IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m maps after the
neighborhood subtraction with GALFIT are shown in the third/forth columns. The
best-fit models, convolved with the PSFs, are overlaid as white contours. The
source sizes and magnitudes are noted in the upper-right and lower-right
corners of and IRAC images.

####  3.2.1 Herschel

Herschel/PACS 100 and 160 $\mathrm{\SIUnitSymbolMicro}$m flux densities were
measured with an aperture radius of 5′′, and the uncertainty was inferred from
the variance of sky background. We adopted the aperture correction factors
suggested in PACS data handbook. HLS1115 was observed to be 20.7$\pm$3.1 mJy
at 100 $\mathrm{\SIUnitSymbolMicro}$m and 50.6$\pm$3.8 mJy at 160
$\mathrm{\SIUnitSymbolMicro}$m. HLS0307-28-A was observed to be 50.7$\pm$2.0
mJy at 100 $\mathrm{\SIUnitSymbolMicro}$m and 110.1$\pm$4.7 mJy at 160
$\mathrm{\SIUnitSymbolMicro}$m. HLS0307-28-B was blended with a low-redshift
source and remained undetected in PACS images, and the two sources in A2813
fell outside of the PACS footprint.

We measured source flux densities using PSF photometry in Herschel/SPIRE 250,
350 and 500 $\mathrm{\SIUnitSymbolMicro}$m images. We assumed point-source
models for the majority of ALMA-detected SMGs, except for two extended sources
($R_{\mathrm{e}}>1\farcs 6$), HLS0840 and HLS0546, which were fit with a 2D
Gaussian model. PSF photometry was conducted using GALFIT (Peng et al., 2010).
We used ALMA source coordinates as prior source positions in the SPIRE maps.
For single-source cases, we floated the source position in the SPIRE images
due to the coarse spatial sampling. If multiple sources exist in a field, the
relative positions of their models are fixed. In several cluster fields, we
also included other PACS (A2813, A3088) or WISE-detected (RXCJ2155) sources to
optimize the fitting through multi-component decomposition.

28 out of 29 ALMA-detected sources were successfully extracted from all three
SPIRE bands. HLS0043-B was undetected in the SPIRE 350 and 500
$\mathrm{\SIUnitSymbolMicro}$m images, and therefore we only present its
$3\sigma$ upper limit ($\sim$11 mJy).

####  3.2.2 Spitzer/IRAC

IRAC photometry was performed with two main methods. Since our targets are
located in cluster fields, some of them may be blended with foreground cluster
members, decreasing the accuracy of aperture photometry. We use GALFIT to
model the source brightness and morphology of 24 cluster-lensed SMGs. The IRAC
warm-mission PSFs released by Hora et al. (2012) for Channel 1 and 2 were used
for model convolution. The orientation angle of spacecraft was also considered
in our modeling routine.

Sérsic and PSF source models were assumed for different sources in IRAC maps.
We adopted a Sérsic model for the majority of SMGs, and a range of Sérsic
index ($n$) between 0.2 and 4.0 was allowed. If the best-fit $n$ is beyond
this range, we fixed it at 0.2, 0.5 (as Gaussian), 1.0 (as exponential) or
1.5, depending on the goodness of fit. 15 sources have well constrained Sérsic
indices in at least one IRAC band, and nine sources were fit with fixed $n$.
We also applied a 2D Gaussian model for all the SMGs. Due to the large PSF
size (FWHM$\sim$1$\farcs$8) and relative small source size (the median half-
light radius along the semi-major axis is $R_{\mathrm{e}}=0\farcs 8$), the
degree of freedom for GALFIT modeling was limited, and thus Gaussian fits do
not show significant deviations in magnitude or source size from Sérsic fits,
as also mentioned by Puglisi et al. (2019).

HLS0043-B and HLS0307-B were fit with PSF models, since their faintness
($\sim$ 21 mag) and heavy blending with bright foreground sources ($\sim$ 17
mag) led to a divergence in the Gaussian/Sérsic modeling. The upper limit of
their $R_{\mathrm{e}}$ is estimated by the minimum measurable deviation from
PSF size ($\sim$0$\farcs$4). HLS1314 is observed as a lensed arc in the NIR
images (UKIRT/WFCAM) and poorly modeled with a single Gaussian/Sérsic profile.
Therefore, we performed aperture photometry of this source after subtracting
nearby sources. Since three clumps were seen in HLS1314 in the 1′′-tapered
ALMA map, we also modeled its morphology with triple Sérsic profiles at the
positions of the ALMA clumps. When the photometric measurements of the three
clumps were combined, the total magnitudes derived in the IRAC images were
recovered.

Another five sources, namely HLS1124-A/B, HLS0455, HLS0505 and HLS0600,
exhibited irregular morphologies in the ALMA continuum maps, revealing galaxy-
lensed rings or multiple components at a resolution of $\sim 0\farcs 2$. All
of these sources are blended with nearby sources, and their morphologies
cannot be well quantified through Gaussian or Sérsic models.

To perform reliable IRAC photometry on these five sources, we adopted their
ALMA continuum images as their morphological models in the IRAC bands. We
clipped their Briggs-weighted, native Band-6 continua at $3\sigma$ and then
convolved them with the corresponding PSFs as the source models. Based on the
$i$-band optical images from the Pan-STARRS DR1 and DES DR1 (Flewelling et
al., 2016; Abbott et al., 2018), we set up 2D Gaussian models for their nearby
sources, which were also convolved with the IRAC PSFs. We then used an MCMC
routine (emcee; Foreman-Mackey et al., 2013) to fit in the IRAC images the
brightnesses of the ALMA sources, the effective radii and brightnesses of the
optical sources, and the sky background. Photometry of the ALMA sources was
then performed with the residual maps after the best-fit models for the nearby
sources were subtracted. Figure 16 shows the ALMA images of the SMGs, optical
images of nearby sources, and IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m
images, before and after this MCMC neighborhood subtraction routine.

Figure 2: IRAC 3.6 $\mathrm{\SIUnitSymbolMicro}$m (upper panel) and 4.5
$\mathrm{\SIUnitSymbolMicro}$m (lower panel) radial surface brightness
profiles of three representative sources in this study (postage stamp images
shown in Figure 1). Source names with their $R_{\mathrm{e}}$, obtained with
GALFIT, are labeled in the legends. Radial profiles of IRAC PSFs (Hora et al.,
2012) are shown as grey dashed lines.

###  3.3 Quantitative morphological modeling

Figure 3: Visibility profiles of three representative sources in this study
(same as the ones in Figure 2). Source names with their $R_{\mathrm{e}}$,
obtained through a CASA uvmodelfit routine, are labeled in the upper-right
legend. Solid lines present the best-fit Gaussian profiles, and dashed lines
indicate the corresponding exponential profile with the same $R_{\mathrm{e}}$.
All points and lines are normalized according to their peak flux density.

Morphological modeling of our sources in the IRAC images was performed along
with photometry using GALFIT, as detailed in Section 3.2. For 24 of our
sources with successful IRAC surface brightness profile modeling, we measured
a median Sérsic index of 0.9 and 1.0 at 3.6 and 4.5
$\mathrm{\SIUnitSymbolMicro}$m, which is close to $n_{\mathrm{med}}=1.2\pm
0.3$ reported by Chen et al. (2015) for $z\simeq 1-3$ SMGs in the HST/WFC3-IR
F160W band. Among the six sources which we display their postage stamp images
in Figure 1, the IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m radial surface
brightness profiles of three representative sources are shown in Figure 2. The
median IRAC effective radius is around $0\farcs 8$, corresponding to 6.5 kpc
at $z=2$ in physical scale before lensing correction. For the five sources
with heavy blending and complex morphology (Figure 16), we did not measure
their structural profiles in either the IRAC or ALMA images. Here, we
concentrate on the modeling of structural profile of the remaining 24 lensed
SMGs at 1.3 mm. As already described, we continue to decompose HLS1314 as
three sources, so 26 sources are studied in this subsection.

We obtained structural parameters of our sources in the ALMA data, modeling
their visibility data using uvmodelfit (for single-source case) and uvmultifit
(for multiple-source case; Martí-Vidal et al. 2014), both under CASA
environment. We assumed Gaussian models for the surface brightness profile of
all sources, consistent with the modeling procedure applied by Puglisi et al.
(2019) and Tadaki et al. (2020). We note that Hodge et al. (2016) measured a
median Sérsic index of $n=0.9\pm 0.2$ for $z\sim 2.5$ SMGs, and Gullberg et
al. (2019) measured $n=1.0\pm 0.1$ at 870 $\mathrm{\SIUnitSymbolMicro}$m.
These Sérsic indices are closer to exponential disk profiles ($n=1$) rather
than Gaussian ones ($n=0.5$). Since the $R_{\mathrm{e}}$ modeled by a Gaussian
profile is smaller than exponential one by only $\sim 0.045$ dex in Hodge et
al. (2016), less than the median uncertainty of $R_{\mathrm{e}}$ in this work
($\sim 0.063$ dex), we keep this Gaussian profile assumption for SMGs in the
ALMA data. We used ALMA source positions, flux densities and morphological
parameters, obtained through SExtractor photometry, as the initial guess of
source models in the uv-plane. Figure 3 shows three representative visibility
profiles of lensed SMGs with different values of $R_{\mathrm{e}}$ (length of
semi-major axis) shown in Figure 2. There is no clear difference between the
$\chi^{2}$ of the best-fit circular Gaussian and exponential models for
HLS0546 and HLS0553-C, and the Gaussian profile fits HLS0612 better.

There were two cases that required special attention for multiple-source
modeling on the uv-plane. HLS0612 was observed to have two minor nearby
components in the ALMA high-resolution image (Figure 1), while uvmultifit
favored a single-source model. Therefore, we discarded the three-source
fitting results. On the contrary, the goodness of a three-source fit with
HLS1314 was better than that of a single-source one, so we continued to
decompose it into HLS1314-A, B and C, consistent with the structural fitting
in the IRAC bands.

The effective radii and their uncertainty of all the 26 sources (24 SMGs and
HLS1314 with three components) are presented in Table 3. We also compared the
IRAC (average of 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m) and ALMA source
centroids, and did not find any detectable offset (Appendix B).

All of our measured effective radii at 1.3 mm are smaller than their
corresponding maximum recoverable angular scales for given ALMA antenna
configurations, except for HLS0840 (Figure 4). As a result, the Re measurement
of HLS0840 has a large error bar (1$\farcs$7). To examine the quality of uv-
profile modeling for sources with extended $R_{\mathrm{e}}$ or at relatively
low S/N, we also simulated and modeled the visibility data for the same
sources as described in Section 3.1 for flux boosting evaluation. The median
$R_{\mathrm{e}}$ measured from the 10–15 sets of simulated visibility data for
each source is consistent with that of the input model at the $1\sigma$
confidence level, demonstrating the overall validity of our source-size
measurements with ALMA.

Figure 4: Visibility profile of HLS0840 (blue filled circles), the source with
the most extended 1.3 mm continuum in our sample. Best-fit morphological
parameters obtained with uvmodelfit are noted in the upper-right corner. The
best-fit Gaussian profile is indicated as the orange solid curve.

###  3.4 Summary of ALMA and IRAC counterparts

We successfully identified the IRAC counterparts of all the 29 sources
discovered with ALMA at $\mathrm{S/N}>4.0$, indicating that our sample is free
from the contamination by any false detection. Our Monte-Carlo simulation also
suggests a low probability ($<2$%) of random association between ALMA and IRAC
sources. We display all the ALMA and IRAC images of sources in our sample in
Figure 1, 16 and 17. We notice that HLS0840 and HLS0546 are two spatially
extended sources ($R_{\mathrm{e,ALMA}}>R_{\mathrm{e,IRAC}}\sim 1$′′) at 1.3 mm
with exceptionally low surface brightness ($\lesssim 0.1$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$)
and only detectable in uv-tapered ALMA maps instead of those at native
resolution (Figure 1). The visibility profiles of these two sources are
displayed in Figure 3 and 4.

Although all of the sources were discovered in lensing cluster fields, we
identified five sources that are subject to additional boosting by foreground
galaxies, i.e., HLS1124-A/B, HLS0455, HLS0505 and HLS0600 as shown in Figure
16. Among the remaining cluster-lensed cases, HLS0553 and HLS2332 are two
spectroscopically confirmed triply-imaged systems behind the cluster fields of
MACSJ0553 and RXCJ2332, respectively. HLS0553-A/B/C have been reported at
$z=1.14$ by Ebeling et al. (2017), and HLS2332-A/B/C were reported at $z=2.73$
by Greve et al. (2012). We do not find any other multiply-imaged SMGs in our
sample.

For the sources with known spectroscopic redshifts (e.g., from a CO search
with the IRAM 30-m telescope; Dessauges-Zavadsky et al., in prep.), our ALMA
Band-6 observations also targeted CO lines when possible. We confirmed the
redshifts of HLS0553 and HLS2332, and obtained/confirmed redshifts for
additional seven sources (Sun et al., in prep): HLS0455 ($z=2.93$; Zavala et
al. 2015), HLS1314 ($z=1.45$), HLS1124-A/B ($z=1.80$), HLS0011-A/B ($z=2.27$)
and HLS0600 ($z=2.89$). We did not detect any spectral line feature in the
ALMA data of the remaining 16 sources, and thus their redshifts remain
unknown.

Among the full sample of 29 ALMA sources, we found 16 are potentially isolated
sources with no detectable companion brighter than 0.6 mJy at 1.3 mm. We also
identified six SMGs that exhibit close companions at similar redshifts, namely
HLS1124 (observed as 4 components and grouped as 2 in Table 2; hereby noted as
4/2 and same later), HLS0111 (2/2), HLS0600 (2/1), HLS0612 (3/1), HLS1314
(3/1) and HLS2104 (2/2). Such grouping is determined by the angular separation
and SExtractor deblending threshold ($\sim$3′′) on 1′′-tapered maps. Source
groups HLS0111, HLS0600, HLS1124 and HLS1314 are spectroscopically confirmed
within a maximum velocity separation of 800
$\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}$ of each other.

##  4 Physical Properties of the Detected Lensed SMGs

###  4.1 SED fitting and photometric redshift

We perform SED modeling with the high-$z$ extension of magphys (da Cunha et
al., 2008, 2015). We also use the photo-$z$ extension of magphys (Battisti et
al., 2019) to estimate the photometric redshift ($z_{\mathrm{phot}}$) and
physical properties simultaneously when the spectroscopic redshift
($z_{\mathrm{spec}}$) is unknown. magphys assumes a Chabrier IMF (Chabrier,
2003), a continuous delayed exponential star-formation history (SFH), a two-
component dust absorption law (Charlot & Fall, 2000) and energy balance
between dust absorption in the UV and re-emission in the infrared. Other key
components of the model assumptions in magphys have been detailed in Martis et
al. (2019).

We adopt our observed photometric measurements to feed the SED fitting routine
without applying any lensing magnification correction. This is due to a lack
of detailed lens models for several of the clusters in this study. Note,
therefore, that all the derived physical quantities are those including the
effects of lensing. We also assume that there is no differential magnification
effect (i.e., effective magnification factors change as function of the source
size and therefore as function of wavelength in general).

For the majority of our sample, their SEDs are modeled using six-band
photometry from 3.6 $\mathrm{\SIUnitSymbolMicro}$m to 1.3 mm. We also utilize
the ancillary data described in Section 2.5 to improve the fitting, especially
to constrain the amount of rest-frame optical dust extinction and therefore
stellar mass. All the individual best-fit SEDs are displayed in Figure 18, and
a summary of the best-fit galaxy properties is presented in Table 4.

The 16–50–84th percentile of the redshift distribution for our sample is
1.23–1.93–2.73, and the highest and lowest photometric redshifts are
$z_{\mathrm{phot}}=0.96^{+0.17}_{-0.45}$ (HLS1623) and
3.23${}^{+0.48}_{-0.59}$ (HLS0043-B). Here the $z_{\mathrm{phot}}$ is the
median of likelihood distribution of redshift for each source. We assess the
uncertainty of photometric redshift estimate based on the likelihood
distribution obtained with magphys. The typical uncertainty of derived
$z_{\mathrm{phot}}$ is $\Delta z=0.16^{+0.05}_{-0.04}\,(1+z_{\mathrm{phot}})$.
HLS1115 exhibits a small $\Delta z$ of 0.03 because of the existence of 9-band
HST data. We also evaluate the far-IR $z_{\mathrm{phot}}$ by matching with the
LIRG templates in Rieke et al. (2009), and the derived redshifts are
consistent with those by magphys within $\sim 1\sigma$ confidence interval.

We also study the dust temperature with far-IR data over a rest-frame
wavelength of 50 $\mathrm{\SIUnitSymbolMicro}$m, following Greve et al.
(2012). We fit the dust continua of all SMGs with modified black-body (MBB).
The dust absorption coefficient is assumed to be
$\kappa=0.040\times(\nu/250\,\mathrm{GHz})^{\beta}$ in the unit of
$\mathrm{m}^{2}\,\mathrm{k}\mathrm{g}^{-1}$, where $\nu$ is the frequency in
GHz in the rest-frame. We assume a fixed dust emissivity of $\beta=1.8$, which
is widely adopted in previous studies (e.g., Díaz-Santos et al., 2017;
Dudzevičiūtė et al., 2020). We also incorporate the uncertainty of redshift in
that of the dust temperature if the $z_{\mathrm{spec}}$ is unknown. The dust
masses derived from this MBB fitting are consistent with those from magphys,
listed in Table 4, with a $1\sigma$ dispersion of 0.08 dex.

HLS1623 was observed to be much fainter at 1.3 mm compared to the prediction
from the SPIRE SED. This may indicate that the SPIRE fluxes are contributed by
certain ALMA-undetected components in the same FoV. Here, we use its SPIRE-
only SED to assess the dust temperature, but we only use its IRAC and ALMA
flux densities for magphys SED fitting.

###  4.2 Multi-wavelength color and redshift

Figure 5: Left: FIR color-color diagram of all lensed SMGs in this work,
color-coded with their redshifts. The squares denote
$z_{\mathrm{spec}}$-confirmed sample and the circles denote
$z_{\mathrm{phot}}$ ones. Unlensed SMGs presented in Ikarashi et al. (2015)
are shown as black dots. A linear fitting through MCMC is plotted as black
solid line with grey line groups indicating the uncertainty. Right: FIR color
($f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$) versus redshift of all our sources,
color-coded with their IRAC color. Markers stand for the same in left panel.
Similar MCMC linear fitting is plotted as black solid line with grey line
groups for its uncertainty. HLS1623 shows exceptionally blue
$f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$ color ($\sim 170$), indicating the 1.3
mm source detection in this cluster field may be incomplete.

We first compare the FIR colors, namely $f_{250}/f_{500}$ and
$f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$, of all galaxies detected in our
survey (Figure 5; left). Here we do not use the color of $f_{250}/f_{350}$ or
$f_{350}/f_{500}$ due to their narrow ranges of distributions. We find the
majority of our galaxy exhibit a positive linear correlation between these two
color indices. An MCMC fitting through emcee suggests the power-law relation
as $\log(f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$})=(1.507\pm
0.247)\times\log(f_{250}/f_{500})+(1.13\pm 0.05)$. This is consistent with the
unlensed SMG sample presented in Ikarashi et al. (2015), for which we apply a
conversion from 1.1 mm flux densities to 1.3 mm by a factor of 0.53. This
indicates red $f_{250}/f_{500}$ and red
$f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$ colors will occur simultaneously,
basically controlled by source redshift and dust temperature.

We then compare $f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$ versus source
redshifts, and an exponential relation can be identified against both the
spectroscopic and photometric redshift samples (Figure 5; right). Similar MCMC
routine implies an underlying relation as
$\log(f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$})=-(0.376\pm 0.044)\times
z+(2.102\pm 0.095)$, with a standard dispersion of 0.16 dex for the measured
$f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$ ratio. This trend is also consistent
with the unlensed sample in Ikarashi et al. (2015) at similar redshift range.
Such a correlation suggests that the $f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$
color of SMGs, at least in this study, is only weakly affected by dust
temperature and mainly reflects the redshift.

On the other hand, the IRAC
$[3.6\,\mbox{$\mathrm{\SIUnitSymbolMicro}${m}}]-[4.5\,\mbox{$\mathrm{\SIUnitSymbolMicro}${m}}]$
color of our SMG sample seems to show a larger dispersion at any given
redshift (see the color coding in the right panel of Figure 5), and thus no
substantial correlation can be derived between the IRAC color and redshift.
This reflects the complexity of stellar age and dust absorption among our SMG
sample. However, we find that at a given redshift, a red IRAC color is likely
to occur simultaneously with a blue $f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$
color. The red IRAC color can be a signature of high dust extinction, and thus
a high dust column density and high IR surface brightness
($\Sigma_{\mathrm{IR}}$). Since $\Sigma_{\mathrm{IR}}$ is observed to be
correlated with dust temperature at various redshift ranges (e.g., Díaz-Santos
et al., 2017; Spilker et al., 2016), this would result in a bluer
$f_{350}/f_{1.3$\mathrm{m}\mathrm{m}$}$ color as we have shown.

Several $z\sim 1$ galaxies were observed with a red IRAC color ($\sim$0.5),
namely HLS0307-28-A, HLS1623, which suggests high dust extinction in these
systems ($A_{V}\gtrsim 7$) through magphys SED modeling. Such a high $A_{V}$
is not seen in AS2UDS sources with secure optical/NIR detection at
$\lambda_{\mathrm{obs}}\lesssim 2.2$ $\mathrm{\SIUnitSymbolMicro}$m (but are
founded in sources with IRAC-only detections; Dudzevičiūtė et al., 2020). Due
to the lack of deep rest-frame UV and optical data, it is not clear whether
the determination of such high $A_{V}$ values is reliable. In the case of
HLS0553 and HLS2332, where HST photometry exists for two of the three triplet
lensed images, the HST photometry obtained at $\lambda\leq 1.6$
$\mathrm{\SIUnitSymbolMicro}$m slightly reduces the estimated $A_{V}$ (by
$\sim 0.4$) and thus the best-fit stellar mass (by $\sim 0.2$ dex; note that
the median uncertainty is $0.32$ dex if HST data is excluded). This
underscores the difficulty of deriving dust absorption in galaxies with only
two-band IRAC observations of their stellar continua.

###  4.3 Star formation in lensed SMGs

Figure 6: sSFR (SFR$\cdot M_{\ast}^{-1}$) versus redshift of all lensed SMGs
in this work. Squares stand for $z_{\mathrm{spec}}$ confirmed sample and
circles stand for $z_{\mathrm{phot}}$ ones. We compare the distribution of our
sources with AS2UDS sample (black dots for sources brighter than 3.6 mJy at
850 $\mathrm{\SIUnitSymbolMicro}$m, i.e., the completeness limit of the
SCUBA-2, and grey dots for fainter ones which are below that limit;
Dudzevičiūtė et al. 2020), the main sequence of star-forming galaxies with a
stellar mass of $10^{11}$ $\mathrm{M}_{\ast}$ (grey solid line; Speagle et
al., 2014), and Hubble parameter $H(z)$ (black dashed line).

The dust-obscured fraction of SFR in a galaxy has been claimed to be
correlated with the stellar mass with no conspicuous evolution found with this
correlation from redshift 2.5 to 0 (Whitaker et al., 2017). For a galaxy with
$M_{\ast}=10^{10}$ $\mathrm{M}_{\odot}$, 80% of the total star formation is
obscured, and for typical $z\sim 2$ SMGs with a stellar mass of $\sim 10^{11}$
$\mathrm{M}_{\odot}$ (e.g., Hainline et al. 2011; Dudzevičiūtė et al. 2020),
this fraction is higher than 95%. Without appropriate lensing correction, we
cannot accurately determine the intrinsic stellar mass of our lensed SMG
sample. However, the median value is
$10^{11.8}\,\mu^{-1}$$\mathrm{M}_{\odot}$, still well above
$10^{10}\,$$\mathrm{M}_{\odot}$ if a $\mu=5$ magnification factor correction
is applied. Therefore, it is reasonable to assume that the majority of star
formation in our sample is dust-obscured, and that the observed IR luminosity
(and emitting regions) adequately represent the total SFR (and star-forming
regions).

We measure a median total IR luminosity of $L_{\mathrm{IR}}=10^{12.92\pm
0.07}\,\mu^{-1}$$\mathrm{L}_{\odot}$ and thus a median SFR of 552$\pm$93
$\mu^{-1}$$\mathrm{M}_{\odot}\,\mathrm{y}\mathrm{r}^{-1}$ without lensing
corrections through SED fitting. Figure 6 displays the sSFR versus redshift of
all the sources in our sample. The distribution of sSFR in redshift space is
consistent with those of galaxies on the so-called star-forming “main
sequence” (MS) with a stellar mass of $10^{11}$ $\mathrm{M}_{\odot}$ (Speagle
et al., 2014), which have a median sSFR of 1.1
$\mathrm{G}\mathrm{y}\mathrm{r}^{-1}$ and a $1\sigma$ dispersion of 0.46 dex.
If the differential magnification is negligible, sSFR will conserve by
lensing, and hence the distribution of sSFR should be the same among lensed
and unlensed SMGs.

We also compare the stellar mass doubling time scale ($t_{\ast}=$1/sSFR) with
the Hubble time ($t_{H}(z)=1/H(z)$) at the redshifts of SMGs in this work. The
ratio between these two quantities ($t_{\ast}/t_{H}$) can in principle act as
an indicator of whether a galaxy is undergoing a significant star-formation
event (e.g., Tacchella et al., 2018). We find a median ratio of 0.18$\pm$0.02,
suggesting that the majority of our sample are vigorously star forming
galaxies. This value is consistent with the median of AS2UDS SMGs at $1<z<3$
(Dudzevičiūtė et al., 2020), regardless of whether the unlensed SMG is
brighter than the single-dish SCUBA-2 detection limit at 850
$\mathrm{\SIUnitSymbolMicro}$m (3.6 mJy) or not. Except for HLS1314 (the
square point at the bottom of Figure 6), no other galaxy in this study shows
$t_{\ast}/t_{H}\geq 1$.

###  4.4 Dust to stellar mass ratio

Figure 7: Dust-to-stellar mass ratio versus sSFR, color-coded with dust
temperature, of all 29 SMGs in this work. Squares represent sources with
confirmed $z_{\mathrm{spec}}$ while circles for $z_{\mathrm{phot}}$-only ones.
Grey-shaded region represents $1\sigma$ distribution range for 24 SMGs in
Santini et al. (2010) at $z_{\mathrm{med}}=2$ (dust mass recomputed by Calura
et al. 2017), with grey dashed line for the median
$\log(M_{\mathrm{dust}}/M_{\ast})$. Orange-shaded region represents $1\sigma$
range for dusty ETGs at $z<0.06$ (Agius et al., 2013; Rowlands et al., 2015).
The two grey solid lines indicate the cases in which the gas depletion time
scale is identical to 1 Gyr or 100 Myr, assuming a constant gas-to-dust-ratio
of 100.

Similar to sSFR, the dust-to-stellar mass ratio ($M_{\mathrm{dust}}/M_{\ast}$;
also referred to as the specific dust mass) and dust temperature are key
observables of the properties of star-forming galaxies that are conserved in
lensing. Since dust is produced through the process of star formation, the
ratio between dust and stellar mass should be closely related to the sSFR.
Figure 7 plots the dust-to-stellar mass ratio versus sSFR in all 29 SMGs,
color-coded with dust temperature. A positive linear relation can be found
between these two quantities, although the normalization is subject to
$T_{\mathrm{dust}}$.

Although such a correlation is not a surprise, the wide range of the dust-to-
stellar mass ratio and sSFR of this 29-SMG sample is remarkable. Distributed
between $-3.4$ and $-1.5$, the $\log(M_{\mathrm{dust}}/M_{\ast})$ span of this
sample is similar to that of low-redshift galaxies, whose SFRs are distributed
between $10^{-1.5}$ and $10^{2}$
$\mathrm{M}_{\odot}\,\mathrm{y}\mathrm{r}^{-1}$ (da Cunha et al., 2010). The
upper end of our $M_{\mathrm{dust}}/M_{\ast}$ distribution matches with
previous SMG literature (e.g., Santini et al., 2010; Dudzevičiūtė et al.,
2020), while our sample does have a significant excess at
$\log(M_{\mathrm{dust}}/M_{\ast})\lesssim-3$. We perform a Kolmogorov–Smirnov
(K-S) test of $\log(M_{\mathrm{dust}}/M_{\ast})$ between our sample and AS2UDS
SMGs in Dudzevičiūtė et al. (2020, 454 SMGs at $1<z<3$), and a null hypothesis
that SMGs in these two samples share the same $M_{\mathrm{dust}}/M_{\ast}$
distribution can be rejected ($p$-value less than 0.01). The lower end of our
$M_{\mathrm{dust}}/M_{\ast}$ distribution matches with that of low-redshift
dusty early-type galaxies (ETGs; Agius et al., 2013), while still higher than
that of dust-poor early-type galaxies (Rowlands et al., 2012, 2015).

The number excess at the lower end of the $M_{\mathrm{dust}}/M_{\ast}$
distribution can actually be a signature of evolved systems with lower gas
fraction and dust destruction in the post-starburst stage (e.g., Rowlands et
al., 2015; Li et al., 2019). Applying a Scoville et al. (2016) conversion from
luminosity at a rest frame of 850 $\mathrm{\SIUnitSymbolMicro}$m
($L_{\nu,850}$) to the molecular gas mass ($M_{\mathrm{gas}}$), we find a low
gas fraction
($f_{\mathrm{gas}}=M_{\mathrm{gas}}/(M_{\mathrm{gas}}+M_{\mathrm{star}})$) of
$0.17\pm 0.08$ for four sources at $M_{\mathrm{dust}}/M_{\ast}<10^{-3}$. We
further discuss this issue of dust-to-mass ratio and its implication in
Section 5.4.2.

Assuming a gas-to-dust ratio of 100, we calculate the gas depletion time scale
($\tau_{\mathrm{dep}}=M_{\mathrm{gas}}/\mathrm{SFR}$) for our sample, and the
median value with 1$\sigma$ dispersion is $226_{-73}^{+196}$ Myr. This is
consistent with the median value of AS2UDS sample ($\sim 320$ Myr) recomputed
with similar method (Dudzevičiūtė et al., 2020).

###  4.5 Dust-to-stellar continuum size ratio

Figure 8: Effective radius ($R_{\mathrm{e}}$) ratio between ALMA and IRAC
counterparts, versus $R_{\mathrm{e}}$ measured in the IRAC bands (unit:
arcsec; before lensing correction). The circles present sources modeled with
Sérsic profiles in the IRAC bands, while stars present the two sources modeled
with PSF profiles. The symbols are color-coded with their $R_{\mathrm{e}}$
ratio between 3.6 and 4.5 $\mathrm{\SIUnitSymbolMicro}$m, usually close to 1.
The grey dashed line indicates the case in which the two effective radii are
identical, while the dotted line shows the median value of $R_{\mathrm{e}}$
ratios (0.38).

SMGs are generally found to host dust continua that are more compact than the
stellar ones (e.g., Hodge et al., 2016, 2019; Lang et al., 2019; Gullberg et
al., 2019). This can be interpreted as an evolutionary connection from SMGs to
cQs at slightly lower redshift: after star formation ceases in a $\sim$1 kpc-
scale region at the galaxy center, a cusp of stellar component will remain in
this region, which can be observed to be compact and quiescent (e.g., Toft et
al., 2014; Simpson et al., 2015; Barro et al., 2016a; Lang et al., 2019).

We adopt the half-light radii measured at 1.3 mm ALMA uv-plane
($R_{\mathrm{e,ALMA}}$) as the effective radius of dust continuum (noted as
$R_{\mathrm{e,dust}}$) and thus the star-forming region. We note that this
ignores any radial gradients of dust temperature or opacity that could alter
the measured size of dust continuum. Based on the FIRE-2 simulation, Cochrane
et al. (2019) predicted that the effective radius of dust emission goes up
with the observed wavelength in the submillimeter/millimeter. This is because
observations at longer wavelengths are more sensitive to the cooler gas and
dust components in the outer region. With an effective rest-frame wavelength
of $440_{-95}^{+137}\,$$\mathrm{\SIUnitSymbolMicro}$m for our sample, the
measured $R_{\mathrm{e,dust}}$ is not expect to vary by more than 10% due to
the difference in the sampled wavelength.

With the IRAC data, we define the effective radius as the geometric mean of
the effective radii measured at 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m
($R_{\mathrm{e,IRAC}}=R_{\mathrm{e,[CH1]}}^{1/2}\cdot
R_{\mathrm{e,[CH2]}}^{1/2}$). At $z_{\mathrm{med}}=1.9$, IRAC Channel 1/2
samples the rest-frame $J/H$ bands with a similar angular resolution. Adopting
the Calzetti et al. (2000) extinction law, dust extinction in the rest-frame
$J/H$ bands ($A_{J}$ and $A_{H}$) will only be 0.3/0.2 of $A_{V}$. For a
median $A_{V,\mathrm{med}}=2.9$ in this study, the dust extinction in the IRAC
3.6/4.5$\mathrm{\SIUnitSymbolMicro}$m bands is 0.9/0.6, which is only $\sim$8%
of the extinction seen in HST $J_{125}/H_{160}$ bands for $z\simeq 2$ SMGs in
Lang et al. (2019). Therefore, the intrinsic stellar distribution should be
close to the light profile in the IRAC bands, without significant overestimate
of the effective radius of stellar mass distribution due to the concentration
of dust extinction at the galaxy center. We measure the effective radius ratio
between 3.6 and 4.5 $\mathrm{\SIUnitSymbolMicro}$m as $0.99^{+0.14}_{-0.19}$
for our sample, and such a consistency between the radii seen in the two bands
also suggests the weak influence of dust extinction on light profiles.
Therefore, we directly adopt the $R_{\mathrm{e,IRAC}}$ as a representation of
the effective radius of stellar mass (noted as $R_{\mathrm{e,star}}$).

One caveat is that the measured $A_{V}$ does not account for the stars which
are fully dust-obscured, and thus the $A_{V}$ derived from an energy-balanced
SED fitting code does not necessarily represent the extinction to the full
stellar mass component (e.g., Casey et al., 2014b). As pointed out in Lang et
al. (2019), it is possible that a compact and obscured stellar component
remains undetected at the ALMA/IRAC continuum centroid (e.g., Simpson et al.,
2017), resulting in an even more compact configuration of stellar mass than
the IRAC light profile.

The comparison of the effective radii of dust and stellar continua is shown in
Figure 8. Except for three sources (HLS1623, HLS0546 and HLS0840), 23 out of
the 26 sources have more compact ALMA dust continua relative to their stellar
components. The median dust-to-stellar continuum size ratio
($R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$) is 0.38, with a 1$\sigma$
distribution from 0.28 to 0.73. Lang et al. (2019) presented a
$R_{\mathrm{e,870\mbox{$\mathrm{\SIUnitSymbolMicro}${m}}}}/R_{\mathrm{e,star}}$
of $0.6\pm 0.2$, similar to our measurement. This is much smaller than the
dust-to-stellar size ratio of local spiral galaxies ($\sim 1.0$, e.g., Hunt et
al. 2015).

###  4.6 IR surface luminosity and dust temperature

Figure 9: Dust temperature versus IR surface luminosity of SMGs in this work
(red squares for sources with confirmed $z_{\mathrm{spec}}$ and circles for
$z_{\mathrm{phot}}$-only ones). We also compare our data with galaxy sample
from GOALS (local LIRG/ULIRGs in orange; Díaz-Santos et al. 2017), AS2UDS
($z_{\mathrm{med}}=2.7$ SMGs in black; Dudzevičiūtė et al. 2020; Gullberg et
al. 2019), SPT ($z_{\mathrm{med}}=4.3$ SMGs in green; Spilker et al. 2016;
Strandet et al. 2016) and KINGFISH (nearby galaxies in blue; Skibba et al.
2011; Hunt et al. 2015).

Since all of the 26 sources (except for heavily-blended galaxy-lensed cases)
are resolved in the ALMA 1.3 mm maps, we are able to obtain the surface
luminosity of their dust continua. This should also be a conserved quantity
independent of lensing magnification. We hereby define the IR surface
luminosity as $\Sigma_{\mathrm{IR}}=L_{\mathrm{IR}}/(2\pi
R_{\mathrm{e,ALMA}}^{2})$, and thus $\Sigma_{\mathrm{IR}}$ is the average IR
surface luminosity within a radius of $R_{\mathrm{e,ALMA}}$. Note that at a
rest-frame wavelength of $\sim$440 $\mathrm{\SIUnitSymbolMicro}$m, the
$R_{\mathrm{e,ALMA}}$ is expected to trace the size of cold dust component and
star-forming region better than the far-IR emission (e.g., 70
$\mathrm{\SIUnitSymbolMicro}$m continuum size used in Díaz-Santos et al.,
2017). Due to the effect of radial $T_{\mathrm{dust}}$ gradient or dust
optical depth, the $R_{\mathrm{e}}$ at the wavelength of FIR SED peak would be
smaller than the $R_{\mathrm{e,ALMA}}$ (e.g., Cochrane et al., 2019), and thus
the $\Sigma_{\mathrm{IR}}$ might be underestimated. The conversion factor
between IR and 1.3 mm surface brightness
($\Sigma_{\mathrm{IR}}/\Sigma_{\mathrm{1.3mm}}$) is $10^{10.65\pm 0.20}$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}\,\mathrm{L}_{\odot}^{-1}\,\mathrm{k}\mathrm{p}\mathrm{c}^{2}$
for sources in our sample, and the factor between surface SFR density and IR
luminosity ($\Sigma_{\mathrm{SFR}}/\Sigma_{\mathrm{IR}}$) is $10^{-10.17\pm
0.07}$
$\mathrm{M}_{\odot}\,\mathrm{y}\mathrm{r}^{-1}\,\mathrm{L}_{\odot}^{-1}$.

We plot the dust temperature versus IR surface luminosity in Figure 9,
comparing the distribution with those of the GOALS (local LIRG/ULIRGs; Díaz-
Santos et al. 2017), AS2UDS ($z_{\mathrm{med}}\simeq 2.7$ SMGs; Gullberg et
al. 2019; Dudzevičiūtė et al. 2020), SPT ($z_{\mathrm{med}}\simeq 4.3$ SMGs;
Spilker et al. 2016; Strandet et al. 2016) and KINGFISH (nearby galaxies;
Skibba et al. 2011; Hunt et al. 2015) galaxies. Galaxies in this work exhibit
a wide 3 dex range of $\Sigma_{\mathrm{IR}}$ from $10^{9.4}$ to $10^{12.4}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$. Such a range of
$\Sigma_{\mathrm{IR}}$ does resemble that of the local (U)LIRGs in the GOALS
sample.

We find that the dust temperature of SMGs in this work is barely correlated
with the IR surface luminosity. At $\Sigma_{\mathrm{IR}}<10^{11.5}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, the lensed SMGs
exhibits comparable dust temperature as local LIRGs (Díaz-Santos et al., 2017,
also assumed $\beta=1.8$). This indicates that the
$\Sigma_{\mathrm{IR}}-T_{\mathrm{dust}}$ relation for LIRGs may not evolve
with redshift up to $z\sim 2$. Symeonidis et al. (2013) analyzed Herschel-
selected LIRGs at $z\lesssim 1$ and suggested no significant evolution of dust
temperature at a constant IR luminosity ($L_{\mathrm{IR}}\lesssim
10^{11.5}\,$$\mathrm{L}_{\odot}$) across $z=0\sim 1$. However, ULIRGs
($L_{\mathrm{IR}}\geq 10^{12}\,$$\mathrm{L}_{\odot}$) at $z\gtrsim 1$ are much
cooler than their local analogs or more optically thick. Assuming that
$\Sigma_{\mathrm{IR}}$ is well correlated with the intrinsic $L_{\mathrm{IR}}$
(reported at various redshifts, e.g., Rujopakarn et al., 2011; Lutz et al.,
2016; Fujimoto et al., 2017), our result is consistent with Symeonidis et al.
(2013) but extending out to $z\sim 2$.

At $\Sigma_{\mathrm{IR}}>10^{11.5}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, the lensed SMGs
in our sample seem to show lower dust temperature than both local ULIRGs and
$z\sim 4$ SPT sources (biased towards galaxy-lensed cases), but consistent
with $z_{\mathrm{med}}=2.7$ SMGs in AS2UDS sample ($\Sigma_{\mathrm{IR}}$
based on 870 $\mathrm{\SIUnitSymbolMicro}$m continuum size; Gullberg et al.
2019; Dudzevičiūtė et al. 2020). However, such a comparison is limited to the
sample size since only two independent sources in our sample (HLS0553 and
HLS2332) are at $\Sigma_{\mathrm{IR}}>10^{12}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$.

###  4.7 Dust-to-stellar size ratio versus IR surface luminosity

Figure 10: Left: Effective radius ratio between dust and star, versus FIR
surface luminosity density of 26 sources in this work (circles with bold black
edges), 14 SMGs from Lang et al. (2019, stars with slim black edges), 6 cSFGs
from Barro et al. (2016a, rightward triangles). Symbols are color-coded with
their surface stellar-mass densities. The best-fit linear relation for 26
sources in this work is shown as the black dashed line, and the best-fit
relation of the combined sample is shown as the black solid line, with the
grey lines indicating the uncertainty. The solid magenta line shows the track
of the two-component model discussed in Section 5.2. Right: The deviation of
effective radius ratio from the best-fit
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}-\Sigma_{\mathrm{IR}}$ relation (i.e.,
the black solid line in the left panel), versus surface stellar-mass density
$\Sigma_{\mathrm{star}}$. The symbols are the same as in the left panel but
color-coded with galaxy-wide sSFR obtained through SED fitting. The best-fit
linear regression (the black dotted line) suggests a positive correlation
between the two parameters, although the intercept might be a function of
sSFR. Median $\Sigma_{\mathrm{star}}$ of $M_{\mathrm{star}}=10^{11}$
$\mathrm{M}_{\odot}$ early and late-type galaxies (ETGs/LTGs) at $z=1.75$ are
plotted as the orange and green solid vertical lines (van der Wel et al.,
2014).

We then investigate how the dust-to-stellar continuum size ratio is correlated
with other physical quantities. Similar to Lang et al. (2019), we do not find
any obvious correlation between the size ratio and sSFR on integrated-galaxy
scales. However, we do find that this continuum size ratio is well correlated
with FIR surface luminosity, plotted in the left panel of Figure 10. We find
that with the increase of $\Sigma_{\mathrm{IR}}$, the dust-to-stellar
continuum size ratio decreases slowly. We fit a linear relation through emcee,
and the best-fit relation for our 26 sources is shown as follows:

$\log(\frac{R_{\mathrm{e,dust}}}{R_{\mathrm{e,star}}})=(-0.23\pm
0.03)\cdot\log(\frac{\Sigma_{\mathrm{IR}}}{10^{10}})-(0.10\pm 0.04)$ (2)

where $\Sigma_{\mathrm{IR}}$ is in the unit of
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$.

We also incorporate six compact star-forming galaxies (cSFGs) presented in
Barro et al. (2016a) and 14 SMGs presented in Lang et al. (2019) into a
combined dataset of 46 sources. All of these sources are around $z=2.2\pm 0.3$
with accurate ALMA image-plane morphology modeling at 870
$\mathrm{\SIUnitSymbolMicro}$m. The $R_{\mathrm{e,star}}$ of the six cSFGs in
Barro et al. (2016a) was measured using HST/WFC3 F160W images without any
correction for dust extinction ($A_{V,\mathrm{SED}}$ is as low as $1.3\sim
1.6$). $R_{\mathrm{e,star}}$ of 14 SMGs in Lang et al. (2019) was also based
on F160W imaging but corrected with a pixel-to-pixel $A_{V}$ map. A similar
MCMC fitting to these 46 sources shows the following
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}-\Sigma_{\mathrm{IR}}$ relation:

$\log(\frac{R_{\mathrm{e,dust}}}{R_{\mathrm{e,star}}})=(-0.19\pm
0.04)\cdot\log(\frac{\Sigma_{\mathrm{IR}}}{10^{10}})-(0.09\pm 0.04)$ (3)

which is consistent with the fitting of 26 HLS sources within $1\sigma$. Such
a linear relation is relatively tight over the 3.1 dex range of
$\Sigma_{\mathrm{IR}}$ distribution. We measure a standard deviation of the 46
sources from the best-fit correlation as 0.21 dex, and the typical uncertainty
of effective radius ratio is $0.11\pm 0.06$ dex.

We also fit the data with least-squares method and bootstrapping. The derived
relation and its uncertainty are consistent with those by MCMC (slope is
$-0.27\pm 0.04$). If $L_{\mathrm{IR}}$ and $R_{\mathrm{e,star}}$ are (i)
random variables with a narrow distribution and (ii) independent of
$R_{\mathrm{e,dust}}$ distributed over a wide range, then one should expect to
derive a linear relation between
$\log(R_{\mathrm{e,dust}}/R_{\mathrm{e,star}})$ and
$\log(\Sigma_{\mathrm{IR}})$ with a slope of $-0.5$. We also generate mock
dataset of $L_{\mathrm{IR}}$, $R_{\mathrm{e,dust}}$ and $R_{\mathrm{e,star}}$
that satisfies our selection criteria (i.e., above certain thresholds of
$L_{\mathrm{IR}}$ and $\Sigma_{\mathrm{IR}}$) and matches our measurements
with respect to the variance, and the resultant slope is $-0.46$. However, our
fittings suggest that such a slope can be ruled out at a confidence of
$>5\sigma$, indicating that the observed relation is physical and not a direct
consequence of the relatively wider distribution range of
$R_{\mathrm{e,dust}}$.

###  4.8 Stellar surface density

One of the remarkable properties of $z\sim 2$ compact quiescent galaxies (cQs)
is their inferred ultra-high stellar surface density ($\gtrsim 10^{10}$
$\mathrm{M}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$ within effective
radius; e.g., van Dokkum et al. 2008; van der Wel et al. 2014). Such a compact
configuration of the stellar component indicates a previous star-formation
history at a similarly compact scale. This could be achieved through a nuclear
starburst, since galaxy interaction/merger and disk instability can cause gas
inflows and thus trigger the central stellar density enhancement (e.g.,
Hopkins et al., 2008).

We define $\Sigma_{\mathrm{star}}$ as the average stellar mass within
$R_{\mathrm{e,star}}$, namely $\Sigma_{\mathrm{star}}=M_{\ast}/(2\pi
R_{\mathrm{e,star}}^{2})$, similar to the definition of
$\Sigma_{\mathrm{IR}}$. The $\Sigma_{\mathrm{star}}$ of our SMG sample ranges
between $10^{8.6}$ and $10^{11.0}$, with a median value of $10^{9.4\pm 0.1}$
(unit: $\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$), similar to
that seen for SMGs by Hodge et al. (2016). Just like $\Sigma_{\mathrm{IR}}$,
$\Sigma_{\mathrm{star}}$ is a conserved quantity with respect to gravitational
lensing.

We first compare the $\Sigma_{\mathrm{star}}$ of our sample with the general
galaxy population in similar redshift and mass ranges. van der Wel et al.
(2014) showed that for early- and late-type galaxies at $z\sim 1.75$ with a
stellar mass of $10^{11}$ $\mathrm{M}_{\odot}$, which is close to that of our
SMG sample, the median $\Sigma_{\mathrm{star}}$ should be $10^{9.6\pm 0.3}$
and $10^{8.8\pm 0.4}$
$\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, respectively.
Therefore, the majority of galaxies in our sample do match the stellar surface
density of early-type galaxies rather than the late-type, although a test on
the $M_{\ast}-R_{\mathrm{e}}$ plane (e.g., Hodge et al., 2016) cannot be
performed due to the incompleteness of lensing magnification.

We then color-code the
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}-\Sigma_{\mathrm{IR}}$ plot with
$\Sigma_{\mathrm{star}}$ as shown in the left panel of Figure 10. Although a
tight linear relation can be found between the two quantities, we notice that
the deviation of continuum size ratios from the best-fit relation seems to be
correlated with $\Sigma_{\mathrm{star}}$. We then plot this residual of
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ versus $\Sigma_{\mathrm{star}}$ in
the right panel of Figure 10. An MCMC linear fitting under the assumption of
equal weighting of all data points suggests a positive relation between the
two quantities as:

$\Delta[\log(\frac{R_{\mathrm{e,dust}}}{R_{\mathrm{e,star}}})]=(0.24\pm
0.07)\cdot\log(\frac{\Sigma_{\mathrm{star}}}{10^{10}})+(0.17\pm 0.05)$ (4)

where $\Sigma_{\mathrm{star}}$ is in the unit of
$\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$.

This indicates that the relation shown as Equation 3 holds for the majority of
$z\sim 2$ SMGs with a stellar mass surface density of roughly $10^{9}\sim
10^{10}$ $\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$. With the
enhancement of the central stellar mass surface density above $\sim 10^{10}$
$\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, SMGs will show
larger $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ ratio than the regular
relation. This may reflect not only a newly-formed cuspy stellar profile in
the galaxy center, but also potentially the quenching of concentrated star
formation and dissipation of dust remnant through multiple physical process
(e.g., stellar or SMBH feedback), which we discuss further in Section 5.3.3.

We also note that the vertical scatter in the right panel of Figure 10 is
related with the sSFR at the galaxy-integrated scale, as color-coded in the
diagram. At a given $\Sigma_{\mathrm{star}}$, the increase of galaxy-wide sSFR
will lead to the growth of $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$, which
can be a trivial result of the fact that since sSFR is proportional to FIR
luminosity and thus $R_{\mathrm{e,dust}}$ when $\Sigma_{\mathrm{IR}}$
conserves.

##  5 Discussion

###  5.1 Discovery of spatially extended SMGs with low surface brightness

One major discovery of this work is the existence of spatially extended SMGs
with a low IR surface luminosity ($\Sigma_{\mathrm{IR}}<10^{10}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$), e.g., HLS0546
and HLS0840. These galaxies exhibit extended dust continua compared with the
stellar ones, and their FIR surface luminosities are $\gtrsim$1 dex lower than
those of typical SMGs at $z=1\sim 3$ ($\Sigma_{\mathrm{IR}}\gtrsim 10^{11}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$; e.g., Hodge et
al. 2016; Gullberg et al. 2019 and this work).

A galaxy with a similar $\Sigma_{\mathrm{IR}}$ has been reported by the ALMA
Frontier Fields Survey (e.g., A2744-ID05, $\Sigma_{\mathrm{IR}}=10^{9.9}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$; González-López
et al. 2017; Laporte et al. 2017). However, the dust continuum of this galaxy
is relatively compact ($R_{\mathrm{e,dust}}=1.3$ kpc) even compared with the
size of stellar continuum measured with HST/WFC3-IR ($R_{\mathrm{e,star}}=4.0$
kpc). Such a small dust-to-stellar size ratio is different from the extended
nature of the SMGs reported here. HLS0546 and HLS0840 are also different from
the spatially extended SMG at $z=2.8$, SMM J02399-0136
($\Sigma_{\mathrm{IR}}=10^{10.2}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$; Genzel et al.
2003; Ivison et al. 2010) because of even lower $\Sigma_{\mathrm{IR}}$ and a
lack of merger feature for our sample. Most recently, Tadaki et al. (2020)
reported the discoveries of HST-selected massive ($M_{\mathrm{star}}>10^{11}$
$\mathrm{M}_{\odot}$) star-forming galaxies at $z\sim 2$ with
$R_{\mathrm{e,dust}}\simeq 5$ kpc. The dust-to-stellar size ratio and IR
surface luminosity of these sources are generally consistent with those of
spatially extended SMGs reported here.

It should be noted that currently spatially extended $z\sim 2$ SMGs with a low
IR surface luminosity density ($\Sigma_{\mathrm{IR}}\simeq 10^{9.5}$
$\mathrm{L}_{\odot}\cdot\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$) can only be
detected and resolved through ALMA observations with deep integration (e.g.,
pointed observations such as Tadaki et al. 2020; blank-field surveys such as
1.2 mm ASPECS, González-López et al. 2020; Aravena et al. 2020) or of lensing-
cluster fields (this work or ALMA Cycle 6 large program ALCS, Kohno et al., in
prep.). Previous ALMA Band-6 surveys in cosmological deep fields (e.g., HUDF,
Dunlop et al. 2017; GOOD-S, Franco et al. 2018; Hatsukade et al. 2018) can
only reach a 4$\sigma$ depth of $\gtrsim 0.24$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$
at 1.3 mm, which is not sufficient to detect these low-surface-brightness
galaxies ($\lesssim 0.1$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$).
Without lensing, these extended SMGs would remain barely resolved with a
$\sim$1′′ beam, and therefore further uv-tapering of the data would hardly
improve the detectability of these galaxies. Galaxy-lensing usually results in
difficulty of source-plane reconstruction, especially in the rest-frame
optical due to the existence of a bright foreground lensing galaxy,
introducing large uncertainties into multi-wavelength comparison at spatially-
resolved scales. In contrast, cluster-lensed SMGs, often with a magnification
factor of $\gtrsim 5$ and reduced contamination from foreground objects, are
significantly stretched spatially ($R_{\mathrm{e}}\gtrsim 1$′′). Therefore,
the detectability of cluster-lensed intrinsically extended SMGs can be
substantially improved by uv-tapering (as shown in Figure 1) or high-
sensitivity facilities with a larger beam size and recoverable angular scale.

###  5.2 Central starburst vs. galaxy-wide star formation

Qualitatively, one simple way to explain the observed correlation between
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ and $\Sigma_{\mathrm{FIR}}$ seen in
Figures 10 (the left panel) is to assume a model in which a compact dust
component ($R_{\mathrm{e,dust}}\sim 1$ kpc) with a varying surface luminosity
($\Sigma_{\mathrm{IR\ }}$ $\simeq$ $0-10^{12.5}$
$\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$) is superposed on an
extended component ($R_{\mathrm{e,dust}}\sim 5$ kpc, $\Sigma_{\mathrm{IR}}\sim
10^{9.5}$ $\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$). In this
model, the former corresponds to central starburst while the latter
corresponds to galaxy-wide star formation (e.g., a star-forming disk). Such a
two-component model was also suggested by Gullberg et al. (2019) based on the
morphological evidence from an SMG stacking analysis. In this model, as the
compact component (i.e., the central starburst) increases its brightness,
$\Sigma_{\mathrm{IR}}$ increases while
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ decreases, and the trend will
reverse when the compact component fades.

In the left panel of Figure 10, we overplot the behavior of this two-component
model (the solid magenta line), using a set of representative values for
various parameters. More specifically, the effective radii of the compact and
extended dust components were fixed to 1 and 5 kpc. The infrared luminosity of
the extended component was fixed to $10^{12.1}$ L☉ while that of the compact
component was varied between 0 and $10^{12.9}$ L☉ to mimic the rise/decline of
the central starburst. $R_{\rm e,dust}$ was then measured for the combined
source. For the stellar component, $R_{\rm e,star}$ and $M_{*}$ were assumed
to be 3.5 kpc and $10^{11}$ M☉. As Figure 10 shows, this simple two-component
model reproduces the observed trend well.

###  5.3 Structural evolution of SMGs

####  5.3.1 Main sequence offset

To evaluate the star-forming properties of the observed SMGs further, we
calculate the offset between the observed SFR and that of the expected on the
star-forming main sequence, which is a function of redshift and stellar mass
(e.g., Speagle et al., 2014; Schreiber et al., 2015). This quantity is
commonly referred as the main-sequence offset
($\Delta\mathrm{MS}=\log[\mathrm{SFR}/\mathrm{SFR}_{\mathrm{MS}}(M_{\mathrm{star}},z)]$),
and here the SFR on the main sequence ($\mathrm{SFR}_{\mathrm{MS}}$) is
calculated using the formula given by Speagle et al. (2014). Following Aravena
et al. (2020), we also define the boundaries of the main sequence as
$\Delta\mathrm{MS}=\pm 0.4$, and classify sources above/below the main
sequence as starburst/passive galaxies.

The accurate derivation of $\Delta\mathrm{MS}$ requires the knowledge of
intrinsic SFR and stellar mass and thus the magnification factor ($\mu$).
Based on published cluster mass models of MACS1115 and MACS0553 (Oguri, 2010;
Zitrin et al., 2015; Ebeling et al., 2017), we derive a median magnification
factor of $\sim 5$ for HLS1115 and HLS0553-A/B/C. Therefore, we assume a
lensing magnification of $\mu=5$ for all the 26 sources uniformly with
morphological measurements. We show that for a $z\sim 2$ lensed SMGs with an
intrinsic stellar mass of $10^{11}$ $\mathrm{M}_{\odot}$, an uncertainty of
$0.5$ dex (i.e., a factor of $\sim 3$) with a lensing magnification factor
would lead to an error of only 0.12 dex with $\Delta\mathrm{MS}$. This is
smaller than the uncertainty of observed, lensing-boosted $M_{\mathrm{star}}$
and SFR. Therefore, $\Delta\mathrm{MS}$ is a relatively robust quantity
against the uncertainty of magnification factor.

Figure 11 displays the distribution of main-sequence offsets and stellar
masses for the joint 46-source sample. A general agreement between the lensed
and unlensed sample is clear, justifying the use of $\mu=5$ as a
representative magnification factor for the lensed sample. One may notice an
anti-correlation between the $\Delta\mathrm{MS}$ and $M_{\mathrm{star}}$.
However, this could be the consequence of a selection bias against relatively
low-mass ($M_{\mathrm{star}}\lesssim 10^{10.8}\,$$\mathrm{M}_{\odot}$)
galaxies without starburst.

Among the full sample, six (nine) sources can be classified as starburst
(passive) galaxies. The remaining 31 sources are therefore galaxies on the
main sequence. K-S tests suggest no significant difference among the redshift
distribution of the three subsamples.

Figure 11: Main sequence offset ($\Delta\mathrm{MS}$, defined as
$\log[\mathrm{SFR}/\mathrm{SFR}_{\mathrm{MS}}(M_{\mathrm{star}},z)]$) versus
stellar mass ($M_{\mathrm{star}}$) of the joint 46-source sample. The symbols
are the same as Figure 10 though color-coded with their redshifts. We assumed
a uniform magnification factor of $\mu=5$ for 26 sources presented in this
work. $\Delta\mathrm{MS}=\pm 0.4$ are specially noted with dashed lines which
differentiate the regions for passive, main-sequence and starburst galaxies.

Figure 12 (the top panel) shows the relation between the IR surface brightness
and main-sequence offset. There is a positive correlation between these two
quantities, and an MCMC fitting to the joint 46-source sample suggests
$\log(\Sigma_{\mathrm{IR}})=(1.10\pm 0.19)\Delta\mathrm{MS}+(11.42\pm 0.10)$,
although the dispersion is considerable (0.62 dex) for the measured IR surface
luminosity. Such a large dispersion was also seen in Elbaz et al. (2018),
indicating that central starbursts may or may not be present in galaxies which
apparently fall near the main sequence. This relation demonstrates that
sources with larger $\Delta\mathrm{MS}$ are generally galaxies with more
vigorous central star-forming activities and thus higher surface density of IR
luminosity. Note, however, that the spatially extended SMGs reported in
Section 5.1 are on the star-forming MS, indicating that extended dust continua
reflect active star formation over the whole galaxy.

####  5.3.2 MS offset versus dust-to-stellar size ratio

Figure 12: Top: Infrared surface luminosity ($\Sigma_{\mathrm{IR}}$) versus
main-sequence offset ($\Delta$MS). The empty circles denote the 26 sources in
this work, and the black dots denote the 20 sources reported in Lang et al.
(2019) and Barro et al. (2016a). Typical uncertainty of each source is shown
as the grey error bars at the lower-left corner. The median value of all the
sources in the three bins of $\Delta$MS (i.e., considered as passive, main-
sequence and starburst galaxies) are shown as black squares with error bars
denoting the 16-84th percentile of distribution. Best-fit linear relation is
shown as the grey solid line. Bottom: Dust-to-stellar source size ratio
($R_{\mathrm{dust}}/R_{\mathrm{star}}$) versus main-sequence offset. Symbols
are the same as the top panels. The black empty arrows denote the theoretical
evolutionary path of a $M_{\mathrm{star}}\sim 10^{10}$ $\mathrm{M}_{\odot}$
galaxy, computed from the models of Tacchella et al. (2016). Four proposed
evolutionary phases, namely pre-blue nugget (pre-BN), blue nugget (BN), post-
blue nugget (post-BN) and quenching phase (Q), are labelled in colored bold
texts although we stress that the term “blue” does not apply to these highly
obscured systems. Typical $R_{\mathrm{dust}}/R_{\mathrm{star}}$ of local
spiral galaxies is shown as the horizontal dashed line (Hunt et al., 2015). In
both panels, the solid magenta line shows the track of the two-component model
discussed in Section 5.2.

According to Figure 10 and Section 4.7, one should further expect an anti-
correlation between the $\Delta\mathrm{MS}$ and
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$. Figure 12 seems to show a weak
anti-correlation between the dust-to-stellar size ratio and main-sequence
offset. A least-squares linear fitting with bootstrapping suggests the slope
is less than zero but only at a $2.4\sigma$ significance. Similar conclusion
has been made in Lang et al. (2019) where the authors claimed no conspicuous
relation was found between sSFR and $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$.

In Figure 12, we also plot the two-component model discussed in Section 5.2
(the solid magenta line). Unlike the left panel of Figure 10, this simple
model does not produce a good fit in either plot. This is because in these
plots, the observed SMGs exhibit a large spread in $\Sigma_{\mathrm{IR}}$ and
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$, reflecting the diversity of
apparently main-sequence SMGs in the sense that some of them are dominated by
central starburst while others are dominated by galaxy-wide star formation (as
noted by Puglisi et al. 2019). Since the two-component model, as defined in
Section 5.2, assumes a transition of SMGs from those dominated by central
starburst to those dominated by galaxy-wide star formation as sSFR decreases
from the starburst to main-sequence range, it fails to reproduce the main-
sequence SMGs dominated by central starburst. In order to reproduce this
population of compact main-sequence SMGs with this type of two-component
model, it would be necessary to decrease the luminosity of galaxy-wide star
formation further.

####  5.3.3 Comparison with theoretical evolutionary tracks

To better understand the observed trends and underlying properties of the
observed SMGs, we further investigate the distribution of
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}-\Delta\mathrm{MS}$ in the lower panel
of Figure 12 by comparing with theoretical predictions from cosmological
simulations. By stacking the simulated galaxy profiles of different
evolutionary phases, Tacchella et al. (2016) suggested that the radial profile
of sSFR declines from center to outskirts during the central starburst phase,
resulting in a concentrated star-forming region with a smaller radius when
compared with that of evolved stellar continuum. However, after the gas
compaction and central enhancement of SFR, quenching then starts from the
center as a combined effect of gas depletion due to star formation, feedback
and truncation of further gas inflow. Compared with the galaxy center, the
outskirts can still retain a ring-like star-forming region, increasing the
ratio of $R_{\mathrm{e,SFR}}/R_{\mathrm{e,star}}$.

Through this so-called compaction and “inside-out” quenching scenario, the
main-sequence offset (and thus IR surface luminosity) of an SMG can increase
and subsequently decline, and the dust-to-stellar size ratio will decline
first and then rise. Such behaviors of a galaxy evolution model are generally
consistent with the trends as we see in Figure 10 and 12.

To conduct a quantitative comparison with the theoretical predictions, we
utilize the radial profiles of $\Sigma_{\mathrm{SFR}}$ and
$\Sigma_{\mathrm{star}}$ for 26 simulated galaxies presented in Tacchella et
al. (2016). The median stellar mass of these galaxies is $\sim 10^{10}$
$\mathrm{M}_{\odot}$ at $z=2$, and Tacchella et al. (2016) stacked the radial
profiles into four evolutionary phases, namely pre-blue nugget (pre-BN), blue
nugget (BN), post-blue nugget (post-BN), and quenching (Q) phase with
increasing time. Here, the BN refers to a massive compact star-forming galaxy
of high central density in stellar mass, following Zolotov et al. (2015),
although the blue stellar population in our systems is heavily dust reddened
and therefore would appear red, making “blue” nugget somewhat a misnomer.
Using the stacked $\Sigma_{\mathrm{SFR}}(r)$ and $\Sigma_{\mathrm{star}}(r)$
profiles, we compute the main-sequence offset and the ratio between the
effective radii of star-forming region and stellar component in these four
phases, and overlay the evolutionary trend with arrows in the bottom panel of
Figure 12. We assume the BN phase occurs at $t_{z=2}$, i.e., the age of the
universe at $z=2$ in the unit of Gyr. The remaining three phases are assumed
to occur at $t=t_{z=2}-0.4$, $+0.3$ and $+1.0$ Gyr, consistent with the time
range reported in Tacchella et al. (2016).

We find that the observed distribution of
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ and $\Delta\mathrm{MS}$ generally
matches the theoretical evolutionary tracks. The slopes of these arrows are
less than zero, consistent with the tentative negative slope reported in
Section 5.3.2, and the vertical offset between the pre-BN $\rightarrow$ BN and
post-BN$\rightarrow$Q trend can also explain the large dispersion of
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ seen around the main sequence.

We find that the five main-sequence SMGs with more extended dust continua
(i.e., $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}>1$; including the two extended
SMGs reported in Section 5.1) hold a higher stellar mass surface density (mean
$\Sigma_{\mathrm{star}}=10^{9.9\pm 0.4}$
$\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$) than the 26 main-
sequence compact sources at $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}<1$ (mean
$\Sigma_{\mathrm{star}}=10^{9.4\pm 0.1}$
$\mathrm{M}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$). This indicates
that SMGs with more extended dust continua are likely at a later evolutionary
phase than those compact ones.

One minor mismatch is that the compact/extended SMGs exhibit a even
lower/higher $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$ than the theoretical
predictions. This result is tentative because the SMG sample is biased to a
higher stellar mass (median
$M_{\mathrm{star}}=10^{11.1}\,$$\mathrm{M}_{\odot}$). By dividing the
simulated galaxies at $M_{\mathrm{star}}=10^{10.2}\,$$\mathrm{M}_{\odot}$ into
two mass bins, Tacchella et al. (2016) reported that galaxies with higher
stellar masses show a more noticeable evolutionary pattern of the
$\mathrm{SFR}(r)$ profile. This will likely result in a broader distribution
of $R_{\mathrm{e,SFR}}/R_{\mathrm{e,star}}$ at various evolutionary phases.
Therefore, we suggest that our observations are consistent with the
theoretical evolutionary track of galaxies in a gas compaction and potentially
subsequent “inside-out” quenching process, which can be the driver of the
tight $R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}-\Sigma_{\mathrm{IR}}$
correlation seen in Figure 10.

####  5.3.4 Further test of the evolutionary sequence

Future observations of spatially extended SMGs may testify whether they are
undergoing a starburst and subsequent quenching process. This can be performed
through NIR spectroscopy by providing tighter constraints on instantaneous SFR
(e.g., through H$\mathrm{\alpha}$ spectroscopy, sensitive to star formation in
previous $\sim 3-10$ Myr; Kennicutt & Evans 2012; although this is subject to
dust obscuration) and comparing with the dust-based one which traces the star
formation at a longer duration ($\sim 100$ Myr). If the properly dust-
corrected H$\mathrm{\alpha}$ SFR is significantly smaller than the dust-based
one, then it is very likely that the galaxy is in a quenching phase, since the
H$\mathrm{\alpha}$ emission can also be contributed by AGN activity. In
addition to this, NIR spectroscopy can also probe the stellar age (e.g.,
through 4000 $\AA$ break; Barro et al. 2016b; Newman et al. 2015, 2018) and
thus the evolutionary stage of the galaxy (but only for the less obscured
stellar populations).

Previous works also reported a stronger [N II] $\lambda$6585 flux compared
with that of H$\mathrm{\alpha}$ in a $z=1.67$ fast-quenching galaxy (Barro et
al., 2016b) and $z\gtrsim 2$ cluster-lensed cQs (Newman et al., 2018). A high
value of the [N II]/H$\mathrm{\alpha}$ ratio is typical of the LINER galaxies
in the BPT diagram (Baldwin et al., 1981), suggesting a suppressed star
formation, which has been observed in low-redshift red-sequence or post-
starburst galaxies (e.g., Yan et al., 2006), although possible AGN
contamination will be an issue with LINER-like spectra in general. For most of
the targets in this work, such an analysis is inaccessible due to the general
lack of $J/H$-band photometry and spectroscopy, and thus the true SFH can only
be determined accurately with more observational data.

###  5.4 The evolutionary picture of SMGs

####  5.4.1 Trigger mechanism of SMGs

Our ALMA observations suggest that late-phase wet-wet major merger could be
the triggering mechanism for $\sim$27% (6 out of 22 cases) of the SMGs studied
in this work because of the existence of companions (angular separation of
$2\farcs 5_{-1.8}^{+5.7}$ in the image plane) at similar redshift (Section
3.4; e.g., Bournaud et al. 2011). Our derived major-merger pair ratio
($\gtrsim$27%) is consistent with Fujimoto et al. (2017, as 27%) and An et al.
(2019, as 22%), but note that we do not have a uniform survey of companions
due to the cluster lensing effect, and therefore our measurement should be a
lower limit.

The remaining 16 sources in this work are potentially isolated SMGs with no
companion brighter than $\sim 0.6$ mJy at 1.3 mm, a quarter of the median ALMA
flux density we measured for the primary sources. Since a stellar mass ratio
as $4:1$ between the major and minor components is widely adopted to
distinguish major and minor merger (e.g., Lotz et al., 2011; Man et al.,
2016), under the assumption that FIR flux ratio equals to stellar mass ratio
for galaxies in a merging system, our observations suggest no evidence of
major-merger companion for $\lesssim$73% of the SMGs in this lensed sample.
Note that we cannot rule out the possibility that some SMGs are in the late
phase of major merger (e.g., separation is $\lesssim 1$ kpc between various
components), which should be further tested through high-resolution
observations of gas kinematics (e.g., Litke et al., 2019; Neeleman et al.,
2019).

Based on the companion search, the co-centered distribution of dust/stellar
continua and the evolutionary trends presented in this work, we conclude that
SMGs could be triggered by a variety of mechanisms including major merger,
minor merger and secular burst, consistent with the conclusions in Fujimoto et
al. (2017), Rujopakarn et al. (2019), Lang et al. (2019) and Jiménez-Andrade
et al. (2020). Such a diversity of triggering mechanism is different from what
we have seen in galaxies of comparable infrared luminosities in the local
Universe, i.e., ULIRGs, in which the major merger is the prevalent mode of
igniting a circumnuclear starburst (e.g., Sanders et al., 1988).

We find no clear difference of the central and total sSFR between the 27% SMGs
with close FIR-bright companions and the rest of our sample. This suggests
that other mechanisms like minor mergers (e.g., Gómez-Guijarro et al., 2018),
mergers with gas-poor companions and secular inflow of gas (e.g., Dekel &
Burkert, 2014) may trigger the star formation in the core of SMGs at
comparable intensity as the major merger. However, we also stress that there
are a number of caveats in our analysis of the merging fraction, most notably
being the need for companions to be FIR-bright and to lie close enough in the
source plane to suffer comparable amplification, which may weaken these
conclusions.

####  5.4.2 The connection between SMGs and cQs

Previous studies have suggested that SMGs are linked to the cQs seen at
slightly lower redshift. The evidence includes the matched number densities
(e.g., Simpson et al., 2014; Dudzevičiūtė et al., 2020), clustering properties
(e.g., Hickox et al., 2012; An et al., 2019) and size/mass similarities (e.g.,
Barro et al., 2016a; Lang et al., 2019). In this subsection, we discuss about
the possible evolutionary connection between SMGs and cQs using those physical
quantities that are independent of lensing effects.

We show that a subset of SMGs in our sample exhibits similar dust-to-stellar
mass ratios (four at $M_{\mathrm{dust}}/M_{\mathrm{star}}<10^{-3}$) as early-
type galaxies at low redshift (Figure 7). Such a low dust-to-stellar mass
ratio is consistent with the stacked $z\sim 1.8$ quiescent galaxies (Gobat et
al., 2018). Assuming a gas-to-dust ratio of 100, these galaxies also match
with the $z\sim 0.7$ post-starburst (PSB) galaxies reported in Suess et al.
(2017). One caveat is that the lack of NIR data may lead to an overestimate of
dust extinction and thus stellar mass through an energy-balance approach in
the SED fitting (e.g., HLS1623), but we show that for HLS1314 with accurate
$J/K$-band photometry, its $A_{V}$ is tightly constrained and thus its
$\log(M_{\mathrm{dust}}/M_{\mathrm{star}})$ can be determined as $-3.25\pm
0.10$, similar to the values for PSBs and early-type galaxies quoted above.

We also show that the stellar surface density ($\Sigma_{\mathrm{star}}$) of
SMGs in this work match better with early-type than late-type galaxies at
similar redshift (van der Wel et al., 2014), in the right panel of Figure 10.
Since the spatial distribution of star-forming regions in SMGs is typically
more compact than those of stars, the ongoing intense star formation will lead
to an even higher $\Sigma_{\mathrm{star}}$ in later phases, increasing the
difference from the typical value of late-type galaxies at that cosmic age.
Such a comparison has been conducted by Barro et al. (2013), Hodge et al.
(2016) and Lang et al. (2019) via HST imaging and interpreted as structural
consistency between the stellar components of SMGs and quiescent galaxies,
suggesting a possible evolutionary link after the cessation of star formation
in SMGs.

We further show that the IR surface luminosity and spatial extent of extended
SMGs in this work may match with SMGs in a transitional phase to quiescent
galaxies. Gullberg et al. (2019) showed the existence of an extended dust
component ($R_{\mathrm{e,dust}}\sim 4$ kpc) of typical $z\sim 3$ SMGs through
stacking 153 of them in the ALMA Band 7. Such an extended component is
reported to contribute to $\sim 13\%$ of the total emission at 870
$\mathrm{\SIUnitSymbolMicro}$m, and the corresponding surface luminosity is
$\Sigma_{\mathrm{IR}}\sim 10^{9.9}$
$\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, under the
assumption of typical SED of an SMG at $z\sim 3$. Furthermore, in the sample
of Gullberg et al. (2019), the FIR surface brightness of the compact dust
component ($R_{\mathrm{e,dust}}\sim 1$ kpc) decreases with the declining of
total $L_{\mathrm{IR}}$ while it remains the same for the extended component
(Figure 13).

In this work, we do discover SMGs with an IR surface brightness comparable to
the faint and extended component of SMGs reported by Gullberg et al. (2019),
and such a comparison is presented in Figure 13. Taking HLS0546
($\Sigma_{\mathrm{IR}}=10^{9.8\pm 0.3}$
$\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$) as an example,
assuming a reasonable lensing magnification of $\mu\simeq 5$, the physical
effective radius of the dust continuum will be $6.3\pm 1.0$ kpc. Consider the
radial gradient of dust temperature, this result could be consistent with the
size of the extended dust continuum of a typical SMG observed at rest-frame
220 $\mathrm{\SIUnitSymbolMicro}$m (Gullberg et al., 2019). This magnification
assumption will also lead to a de-magnified stellar mass of
$M_{\ast}=10^{11.0\pm 0.4}$ $\mathrm{M}_{\odot}$ and
$R_{\mathrm{e,star}}=3.5\pm 0.6$ kpc, within a 1$\sigma$ distribution of
early-type galaxies seen at the given redshift (van der Wel et al., 2014).
Therefore, it is possible that after the dissipation of central star formation
on a compact physical scale ($1-2$ kpc), SMGs will maintain a low-level IR
surface brightness over a more extended galaxy structure ($\gtrsim 4$ kpc), as
suggested by the evolutionary tracks in Section 5.3.3. At the same time, the
previous concentrated star formation would lead to the formation of compact
($\sim 2$ kpc) quiescent and spheroidal galaxies within a timescale of $\sim
300$ Myr. Therefore, spatially extended SMGs with large dust-to-stellar radii
ratio in this study may provide possible evidence of evolutionary connection
between typical SMGs and compact quiescent galaxies at slightly lower
redshift.

Figure 13: 1.3 mm surface brightness ($\Sigma_{1.3\,$\mathrm{m}\mathrm{m}$}$)
versus intrinsic SFR (assuming a lensing magnification of $\mu=5$) of resolved
SMGs in our sample (grey circles with black edges). Gullberg et al. (2019)
reported the existence of both compact (shown as triangles) and extended dust
components (squares) of $z\sim 3$ SMGs through stacking. Assuming typical SED
of a SMG at $z=3$, we convert the 870 $\mathrm{\SIUnitSymbolMicro}$m surface
brightness in Gullberg et al. (2019) to that at 1.3 mm by dividing a factor of
2.7, which is $\sim 0.2$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$
for extended components. This value matches those of low-surface-brightness
SMGs in our sample ($\Sigma_{1.3\,$\mathrm{m}\mathrm{m}$}\sim 0.1$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$).

##  6 Summary

We have obtained and analysed ALMA 1.3 mm, Herschel/SPIRE 250/350/500
$\mathrm{\SIUnitSymbolMicro}$m and Spitzer/IRAC 3.6/4.5
$\mathrm{\SIUnitSymbolMicro}$m data of 29 lensed SMGs in 20 cluster fields.
These SMGs were discovered by the Herschel Lensing Survey as exceptionally
bright sources in the far infrared ($S_{\mathrm{peak}}\gtrsim 90$ mJy). We
have carried out modeling of their structural profiles in both the IRAC and
ALMA bands, as well as their SEDs from NIR to millimeter wavelengths. Since
accurate lens models are not yet available for many of the observed SMGs, we
focus our discussion on quantities that are independent of lensing effects,
such as surface brightness and size ratios. When necessary, we also assumed a
canonical magnification factor of $\mu\simeq 5$. The main results of this
study are the following:

  1. 1.

29 sources were detected in our uv-tapered ALMA 1.3 mm maps at S/N$>$4.0. Five
sources are identified as galaxy-lensed or highly-blended cases, and their
IRAC fluxes were successfully decomposed using ALMA and optical priors. The
remaining 24 sources are cluster-lensed SMGs, and all of their dust continua
were spatially resolved (HLS1314 is resolved as three components).

  2. 2.

Because of gravitational lensing, we were able to resolve the structural
profile of stellar continuum in 24 out of 26 SMGs in Spitzer/IRAC bands. 88%
of the SMGs in this study show smaller half-light radii in dust continua
compared with the stellar ones. The medium
$R_{\mathrm{e.dust}}/R_{\mathrm{e,star}}$ was found to be 0.38$\pm$0.14. Two
sources (HLS0840 and HLS0546) are discovered as spatially extended main-
sequence SMGs with low surface brightness ($\lesssim 0.1$
$\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{-2}$
at 1.3 mm) and only detectable in uv-tapered ALMA maps.

  3. 3.

We fit the SED of all the observed SMGs with MAGPHYS, deriving their physical
properties such as SFR and stellar/dust mass. The dust-to-stellar mass ratio
of SMGs is correlated with sSFR (at galaxy scale; subjected to dust
temperature $T_{\mathrm{dust}}$), and four of them exhibit
$\log(M_{\mathrm{dust}}/M_{\mathrm{star}})<-3$, resembling dusty early-type
galaxies rather than typical SMGs at $z\sim 2$.

  4. 4.

We find that the IR surface luminosity ($\Sigma_{\mathrm{IR}}$) of SMGs in our
sample spans over a wide range of 3 dex, similar to that covered by local
(U)LIRGs in the GOALS sample (Díaz-Santos et al., 2017). At a given
$\Sigma_{\mathrm{IR}}$ below $10^{11.5}$
$\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, our SMGs show
consistent dust temperature as local LIRGs, indicating no significant
evolution of $\Sigma_{\mathrm{IR}}-T_{\mathrm{dust}}$ relation for LIRG-like
galaxies from $z\sim 2$ to the present Universe. At
$\Sigma_{\mathrm{IR}}>10^{11.5}$
$\mathrm{L}_{\odot}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-2}$, the
$T_{\mathrm{dust}}$ of ULIRGs at $z\sim 2$ are lower than those of the local
ones as reported previously (e.g., Symeonidis et al., 2013; Dudzevičiūtė et
al., 2020).

  5. 5.

We find that the IR surface luminosity of SMGs is anti-correlated with the
dust-to-stellar source size ratio ($R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$).
Compared with a simple analytic model and cosmological simulations (Tacchella
et al., 2016), this relation could be interpreted by the morphological
evolution of SMGs through central starburst and potentially subsequent inside-
out quenching: the central starburst leads to an increase in
$\Sigma_{\mathrm{IR}}$ and compaction of star formation (i.e, smaller
$R_{\mathrm{e,dust}}$), and in a following quenching phase, SFR quickly
declines from the center (smaller $\Sigma_{\mathrm{IR}}$) but retains a
considerable level at the outskirts, resulting in a larger
$R_{\mathrm{e,dust}}/R_{\mathrm{e,star}}$.

  6. 6.

The distribution of spatial offset between the dust and stellar continua in
SMGs does not support any intrinsic offset rather than astrometric
uncertainty. $\gtrsim 27\%$ of the SMGs can be visually identified as late-
stage wet-wet major-merger pair (consistent with Fujimoto et al., 2017; Lang
et al., 2019). Therefore, it is possible that SMGs could be triggered through
a variety of mechanisms, including major/minor merger and secular evolution,
although lensed samples are not well suited to identify early-stage mergers
due to the spatially varying magnification.

  7. 7.

Our observations are consistent with the hypothesis that SMGs are progenitors
of the compact quiescent galaxies at slightly lower redshift. This is because
(1) low-sSFR SMGs show dust-to-stellar mass ratios
($M_{\mathrm{dust}}/M_{\mathrm{star}}$) comparable to those of compact early-
type or post-starburst galaxies, (2) the high stellar surface densities of
SMGs match those of early-type galaxies rather than late-type ones at similar
redshift, and finally (3) spatially extended SMGs in this work exhibit low IR
surface luminosity, matching some expected properties of SMGs in the
transition phase to cQs.

The discovery of spatially extended SMGs in this study expands the population
of submillimeter-selected galaxies to a lower FIR surface luminosity limit,
and therefore provides direct observational constraints on the evolution of
dusty starburst galaxies which may be in the quenching phase. Further studies
of these (possibly) transitional systems will help to establish the
evolutionary picture of massive galaxy at the epoch when the cosmic SFR
density peaks (Madau & Dickinson, 2014).

E.E. would like to thank Observatoire de Lyon and European Space Astronomy
Centre (ESAC) Faculty Council for hosting and supporting his long-term visits,
during which this project was originally started. I.R.S. acknowledges support
from the Science and Technology Facilities Council (grant number
ST/T000244/1). H.E. gratefully acknowledges financial support from Space
Telescope Science Institute for program SNAP-12884. W.R. is supported by the
Thailand Research Fund/Office of the Higher Education Commission Grant Number
MRG6280259 and Chulalongkorn University’s CUniverse. This paper makes use of
the following ALMA data: ADS/JAO.ALMA#2015.1.01548, #2016.1.00372 and
#2017.1.01658.S. ALMA is a partnership of ESO (representing its member
states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and
ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic
of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.
This work is based on observations made with the Spitzer Space Telescope,
obtained from the NASA/ IPAC Infrared Science Archive, both of which are
operated by the Jet Propulsion Laboratory, California Institute of Technology
under a contract with the National Aeronautics and Space Administration.
Herschel is an ESA space observatory with science instruments provided by
European-led Principal Investigator consortia and with important participation
from NASA. SPIRE has been developed by a consortium of institutes led by
Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China);
CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm
Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ.
Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development
has been supported by national funding agencies: CSA (Canada); NAOC (China);
CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC,
UKSA (UK); and NASA (USA). This research is based on observations made with
the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are
associated with program SNAP 12884, CLASH Treasury Program (GO 12065) and
RELICS Treasury Program (GO 14096).

##  Appendix A False detection rate and Completeness of ALMA Detections

Figure 14:  Left: Cumulative false detections per simulated single-pointing
FoV, as a function of measured S/N through SExtractor aperture photometry. We
then cut our detection at 4$\sigma$ (red dashed line) to eliminate fake
detections. Right: Source detection completeness for simulated point sources,
as a function of source flux density assuming a continuum RMS of 0.105
mJy/beam and primary beam correction factor of $\times\sqrt{2}$). We plot the
completeness curves at various setting of S/N cut, and highlight S/Ncut=4.0
curve in red solid steps.

We investigated the false detection rate and completeness of our ALMA
detections. We generated Gaussian white noise maps and smoothed them with a 2D
Gaussian kernel to simulate the noise distribution of ALMA continuum images
obtained with a 1′′ beam. We then run SExtractor on the mock noise maps with
the same settings used for the scientific maps, and treated all the extracted
sources as false detections. This experiment demonstrates that for a typical
ALMA Band 6 single-pointing FoV at this uv-taper, we expect to detect $\sim
1.69$ fake sources at S/N$>$3, $\sim 0.84$ fake sources at S/N$>$3.5, but only
$\sim 0.06$ at S/N$>$4 (Figure 14, left). Therefore, for the purpose of
eliminating false detections, we applied an S/N cut of our detections at
S/N=4. This ensures the total false detection number to be $\lesssim 1$ in all
of our ALMA maps.

We also studied the point-source detection completeness through similar
simulation. We generate 100 sources with a given flux density on a Gaussian
white noise map, comparable to the combined area of 20 single-pointing FoVs.
We then blur the image with a 2D Gaussian kernel of FWHM=1′′, simulating the
1′′-tapered continuum images. This mock image was fed into the same SExtractor
routine for source detection, and we counted the number of recovered
detections. Such an experiment is repeated at various input source strength,
and the completeness as a function of source flux density is shown in the
right panel of Figure 14. Based on our simulations and a S/N cut of 4.0, our
detection is $47\pm 5\%$ complete for point sources at S/N=4.0, corresponding
to a flux density of 0.63 mJy for median continuum RMS (0.105 mJy) and half-
radius primary beam correction ($\times\sqrt{2}$). At S/N=5.0 (0.79 mJy for
median continuum RMS), the point-source completeness is $80\pm 4\%$.
Therefore, our ALMA survey is fairly complete for any point source brighter
than 0.9 mJy.

Note that Franco et al. (2018) showed the completeness of ALMA detection
decreases dramatically for larger galaxy sizes. With 0$\farcs$6-tapering,
their 1.1 mm observation is 94% complete for point sources at 1.2 mJy but only
9% for FWHM=0$\farcs$6 sources. We notice that 7 out of 26 sources with our
FIR morphology modeling exhibit effective radii $R_{\mathrm{e}}>0\farcs 5$.
Assuming similar completeness curves as Franco et al. (2018) while scaled to
0.1 mJy continuum RMS and 1′′-tapering, our completeness can be lower than 50%
at 1.4 mJy for these extended sources.

##  Appendix B Spatial offset between ALMA and IRAC counterparts

Figure 15: Left: Spatial offsets between the dust and stellar continua of
SMGs. Lensed SMGs in this work are plotted as red circles, and 73 SMGs in
CANDELS field are plotted as black dots (ALMA positions by Fujimoto et al.
2017; matched with HST counterparts in van der Wel et al. 2012). The only case
where offset is larger than 1′′, namely HLS0546, are specially noted. Middle:
Histograms of dust/stellar offset distribution in this work (red) and Fujimoto
et al. (2017, black). Best-fit normal distribution profiles are shown as solid
lines, with their standard deviations ($\sigma$) labeled. Right: K-S test of
the dust/stellar offset distribution in this work (red steps), compared with
cumulative distribution function (CDF) of normal distribution (black dashed
line). Assuming a typical IRAC position error at 0$\farcs$20 ($\sim$10% of PSF
FWHM), the hypothesis that the dust and stellar continuum of SMG is co-
centered cannot be ruled out ($p$-value as 0.909).

Through our IRAC morphological fitting in Section 3.2, the sky coordinates of
19 out of 26 ALMA sources were allowed to vary, so we compare these IRAC
positions (average of 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m) with ALMA
positions obtained from the uv-plane fitting. We find that the mean offset
between ALMA and IRAC sources is $0\farcs 05\pm 0\farcs 05$ in RA and $0\farcs
00\pm 0\farcs 04$ in DEC. Assuming an uncertainty of 0$\farcs$2 with IRAC
images (i.e., $\sim 10\%$ of PSF FWHM), the expected standard error of the
dust/stellar offset is $\sim$0$\farcs$05 if the stellar and dust continua are
co-centrally distributed, matching our result.

We compare our lensed SMGs with unlensed sample compiled in Fujimoto et al.
(2017). In their work, the authors conducted a large-sample analysis of
$\sim$1 mm ALMA sources on the uv-plane, and a significant number of sources
are also covered in the HST CANDELS fields (van der Wel et al., 2012),
enabling precise measurements of the offset between the dust and stellar
continua in SMGs at $z\sim 2.5\pm 0.5$. We cross-matched these two catalogs
and identified the HST counterparts of 73 SMGs with a maximum allowed
separation of 1′′. The spatial offset between the dust and stellar continuum
center in each source is plotted in the left panel of Figure 15. The
histograms of the dust/stellar offsets in these two SMG samples are shown in
the middle panel of Figure 15, and the Gaussian fittings of these histograms
suggest a standard deviation between dust and stellar continuum centroids as
$0\farcs 23\pm 0\farcs 07$ in our lensed sample, and $0\farcs 27\pm 0\farcs
03$ in the sample of Fujimoto et al. (2017).

Such a comparison suggests an intrinsically small spatial offset between the
centroids of stellar/dust continua in lensed SMGs, since our offset
measurements are amplified by lensing. Though HST has a higher angular
resolution and thus a higher precision in determining the centroid of SMG in
the optical/NIR bands, strong dust obscuration could affect the observed
morphology of stellar continuum in HST (e.g., Hodge et al. 2016). Such an
effect can probably bias the direct spatial mapping between observed stellar
light in the near-infrared and the intrinsic stellar mass (also suggested by
Lang et al., 2019), result in an underestimate of the uncertainty for stellar
centroid determination in HST images (e.g., Chen et al., 2015).

The only SMG showing offset larger than 1′′ is HLS0546, which is an extended
SMG ($R_{\mathrm{e,dust}}=1\farcs 7\pm 0\farcs 3$) with high Sérsic index
($n\sim 2.3$) for stellar component. These features may suggest that HLS0546
has already entered a late phase of SMG evolution, showing cuspy stellar
component at its center and ongoing star formation at outskirts as discussed
in Section 5.3.3. This can cause the observed offset if the star-forming
regions are not symmetrically distributed in the galaxy disk plane.

We further show that such a distribution of dust/stellar offsets may be
primarily due to astrometric errors through a K-S test. We keep the assumption
of the stellar position uncertainty as 0$\farcs$2 for IRAC images, and the
uncertainties of dust continua are determined during the uv-plane fitting
(median value is 0$\farcs$04). We cannot rule out the null hypothesis that the
sample of spatial offsets relative to astrometric errors are drawn from a
Gaussian distribution (p-value as 0.909). Assuming a slightly lower IRAC
astrometric error (e.g., 0$\farcs$15) will not change the conclusion. Since
the co-centered distribution of stellar and dust continuum will result in a
Gaussian distribution of spatial offsets relative to their uncertainty, this
K–S test demonstrates that we cannot rule out the co-centered distributing
scenario between the ongoing star formation and the evolved stellar population
in $z\sim 2$ SMGs.

##  Appendix C ALMA, IRAC images and SED plots of all sources

Here we show the multi-wavelength images, i.e., ALMA 1.3 mm (native and
tapered) and Spitzer/IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m images, of
five highly blended SMGs in Figure 16 and 20 cluster-lensed SMGs in Figure 17.
We also present magphys best-fit SED of all sources in Figure 18.

Figure 16: Illustration of the IRAC photometry of five highly-blended SMGs
with complex morphology. Note we define two sources in CODEX 39326 (first
row), namely HLS1124-A (the upper-right faint one) and HLS1124-B (the lower
bright one with three sub-components in ALMA map), both at $z=1.80$ with
CO(6-5) detection (Sun et al., in prep.). We display their Briggs-weighted
ALMA 1.3 mm continuum images with their flux densities at first column,
optical images with ALMA contours (levels: 4, 10$\sigma$) at second column.
Original IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m images are displayed at
Column 3 and 5, with the best-fit optical source models shown in contours and
magnitudes shown in text. Residual images after nearby source subtraction are
shown at Column 4 and 6, with ALMA-model contours and measured magnitudes in
text.

Figure 17: Postage stamp images of the remaining 20 cluster-lensed SMGs in
this work. The layouts of images are the same as those in Figure 1. Note that
we do not subtract other sources when we display the ALMA images of
HLS1314-A/B/C, and therefore three components can be identified in their
1′′-tapered ALMA images.

Figure 18: magphys best-fit SEDs of all the 29 sources in this work. Source
names with their best redshifts are noted above the plots. The photometric
data are plotted as open red squares, and the best-fit SEDs are shown as the
black solid lines. Key parameters derived from the SED fitting are presented
in Table 4. The measured Herschel flux densities of HLS1623 are shown as upper
limits, since they are likely partially contributed by other sources not
detected by ALMA.

Table 1: Summary of Herschel/SPIRE, ALMA and Spitzer/IRAC observations  Cluster Name | CoordinatesaaCoordinates of the ALMA pointing centers, not the exact positions of the ALMA-detected sources. |  Herschel/SPIRE | ALMA Band 6 |  Spitzer/IRAC
---|---|---|---|---
| RA | Dec | OBSID |  $t_{\mathrm{obs}}$bbTotal observation time. For the ALMA observations obtained in the “single” pointing mode, this is the full on-source integration time. For SPIRE and the ALMA multiple-pointing observation, this indicates the full scan time for the scientific targets. | PointingccMode of ALMA observations (single or multiple pointings). |  $t_{\mathrm{obs}}$bbTotal observation time. For the ALMA observations obtained in the “single” pointing mode, this is the full on-source integration time. For SPIRE and the ALMA multiple-pointing observation, this indicates the full scan time for the scientific targets. | RMSddContinuum RMS noise with a 1′′ uv-tapered beam. | Program ID |  $5\sigma$ DeptheeMedian $5\sigma$ point-source depth in the final IRAC 3.6/4.5 $\mathrm{\SIUnitSymbolMicro}$m image products.
|  |  |  | (s) |  | (s) | ($\mathrm{m}\mathrm{J}\mathrm{y}\,\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}^{-1}$) |  | (mag, mag)
SPT J0114–4123 | 1:14:39.38 | –41:24:03.1 | 1342247867 | 169 | single | 907 | 0.038 | 12095 | 22.27, 22.37
SPT J0307–5042 | 3:07:49.38 | –50:41:39.1 | 1342240119 | 169 | single | 302 | 0.111 | 12095 | 22.28, 22.39
SPT J0505–6145 | 5:05:33.14 | –61:43:55.0 | 1342229231 | 169 | single | 195 | 0.116 | 80162 | 21.23, 21.60
SPT J0546–5345 | 5:46:39.76 | –53:45:18.2 | 1342240055 | 445 | single | 302 | 0.104 | 60099, 70149 | 22.78, 22.84
SPT J0612–4317 | 6:12:03.68 | –43:17:10.8 | 1342240063 | 169 | single | 302 | 0.148 | 80012, 90233 | 21.42, 21.78
Abell 2813 | 0:43:35.41 | –20:42:00.4 | 1342188582 | 5803 | single | 121 | 0.169 | 60034 | 23.62, 23.65
Abell 3088 | 3:07:08.09 | –28:40:17.4 | 1342188659 | 5803 | single | 604 | 0.067 | 80066 | 22.69, 22.65
CODEX 35646 | 16:23:46.51 | 26:34:11.7 | 1342239984 | 169 | single | 302 | 0.126 | 90218 | 22.28, 22.26
CODEX 39326 | 11:24:02.22 | 24:04:38.4 | 1342256839 | 169 | single | 532 | 0.101 | 12095 | 22.25, 22.17
CODEX 52909 | 11:53:20.37 | 07:56:00.3 | 1342247962 | 169 | multiple | 171 | 0.150 | 90218 | 22.38, 22.29