# Planck's dusty GEMS. V. Molecular wind and clump stability in a strongly lensed star-forming galaxy at $z = 2.2$

R. Canameras, N. P. H. Nesvadba, M. Limousin, H. Dole, R. Kneissl, S.

Koenig, E. Le Floc'h, G. Petitpas, D. Scott

► To cite this version:

R. Canameras, N. P. H. Nesvadba, M. Limousin, H. Dole, R. Kneissl, et al., Planck's dusty GEMS. V. Molecular wind and clump stability in a strongly lensed star-forming galaxy at $z = 2.2$. Astronomy and Astrophysics - A&A, 2018, 620, pp.A60. 10.1051/0004-6361/201833679. hal-02374144

HAL Id: hal-02374144

https://hal.science/hal-02374144

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Planck's dusty GEMS

V. Molecular wind and clump stability in a strongly lensed star-forming galaxy
at z = 2.2*

R. Cañamas<fnref n="1" />,<fnref n="2" />, N. P. H. Nesvadba<fnref n="2" />, M. Limousin<fnref n="3" />, H. Dole<fnref n="2" />, R. Kneissl<fnref n="4" />,<fnref n="5" />, S. Koenig<fnref n="6" />, E. Le Floc’h<fnref n="7" />,
G. Petitpas<fnref n="8" />, and D. Scott<fnref n="9" />

Received 19 June 2018 / Accepted 20 September 2018

ABSTRACT

We report the discovery of a molecular wind signature from a massive intensely star-forming clump of a few $10^9 M_\odot$, in the strongly gravitationally lensed submillimeter galaxy “the Emerald” (PLCK_G165.7+49.0) at $z = 2.236$. The Emerald is amongst the brightest high-redshift galaxies on the submillimeter sky, and was initially discovered with the *Planck* satellite. The system contains two magni- ficient structures with projected lengths of 28.5” and 21” formed by multiple, near-infrared arcs, falling behind a massive galaxy cluster at $z = 0.35$, as well as an adjacent filament that has so far escaped discovery in other wavebands. We used HST/WFC3 and CFHT optical and near-infrared imaging together with IRAM and SMA interferometry of the CO(4–3) line and 850 µm dust emission to characterize the foreground lensing mass distribution, construct a lens model with LENS TOOL, and calculate gravitational magnification factors between 20 and 50 in most of the source. The majority of the star formation takes place within two massive star-forming clumps which are marginally gravitationally bound and embedded in a $9 \times 10^{10} M_\odot$, fragmented disk with 20% gas fraction. The stellar continuum morphology is much smoother and also well resolved perpendicular to the magnification axis. One of the clumps shows a pronounced blue wing in the CO(4–3) line profile, which we interpret as a wind signature. The mass outflow rates are high enough for us to suspect that the clump might become unbound within a few tens of Myr, unless the outflowing gas can be replenished by gas accretion from the surrounding disk. The velocity offset of –200 km s$^{-1}$ is above the escape velocity of the clump, but not that of the galaxy overall, suggesting that much of this material might ultimately rain back onto the galaxy and contribute to fueling subsequent star formation.

**Key words.** galaxies: high-redshift – galaxies: evolution – galaxies: star formation – galaxies: ISM – infrared: galaxies – submillimeter: galaxies

**1. Introduction**

Rapid, intense star formation that occurs in dusty star-forming galaxies at $z \sim 1-4$ (Casey et al. 2014) is expected to dominate the cosmic star-formation rate density at these epochs (e.g., Dole et al. 2006) and corresponds to the early growth phase of giant ellipticals seen in high-density regions of the local Universe (e.g., Lilly et al. 1999; Swinbank et al. 2006). Intense star formation is sustained for timescales up to a few hundred Myr in these high-redshift galaxies. Many recent studies propose that the global properties of the molecular gas reservoirs, including gas fractions, determine whether these galaxies will fall on the high-redshift main-sequence

of star formation (Daddi et al. 2010; Genzel et al. 2010), or be in the starburst mode or the transition regime (e.g., Tacconi et al. 2013; Dessauges-Zavadsky et al. 2015; Lee et al. 2017; Cañamas et al. 2017a). Variations in star-formation efficiency might also play a role in this regard (e.g., Hodge et al. 2015; Genzel et al. 2015; Usero et al. 2015).

The role of local mechanisms such as star-formation feed-
back, winds and turbulence in shaping the interstellar medium
of these galaxies, in regulating their star-formation activity and
in triggering their major growth phase is still a matter of active
debate. For instance, large-scale outflows are a major compo-
nent of galaxy evolution models (e.g., Hopkins et al. 2006), since
they can affect and even quench star formation within the hosts
by expelling the gas to the circumgalactic medium. Molecular
outflows are ubiquitous in nearby ultraluminous infrared galaxies
(ULIRGs; e.g., Weiß et al. 1999; Sturm et al. 2011; Cicone et al.

<footnote n="1">¹ Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, 2100 Copenhagen, Denmark
e-mail: canameras@dark-cosmology.dk</footnote>

<footnote n="2">² Institut d'Astrophysique Spatiale, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Bât. 121, 91405 Orsay, France</footnote>

<footnote n="3">³ Aix Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="4">⁴ European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Casilla 19001, Santiago, Vitacura, Chile</footnote>

<footnote n="5">⁵ Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central Offices, Alonso de Cordova 3107, Casilla 763-0355, Santiago, Vitacura, Chile</footnote>

<footnote n="6">⁶ Chalmers University of Technology, Onsala Space Observatory, Onsala, Sweden</footnote>

<footnote n="7">⁷ Laboratoire AIM, CEA/DSM/IRFU, CNRS, Université Paris-Diderot, Bât. 709, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="8">⁸ Harvard-Smithsonian Center for Astrophysics, 02138 Cambridge, MA, USA</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, 6658 Vancouver, British Columbia, Canada</footnote>

<footnote>* Based on data obtained with the following programs: IRAM Plateau de Bure Interferometer program ID: X0AE. Canada-France-Hawaii Telescope program ID: 14AF06. Submillimeter Array program ID: 2013B-S050. Spitzer Space Telescope program ID: 10010.</footnote>
2014; Veilleux et al. 2013) and are attributed to either feedback from star formation, or from a central AGN or both. At high redshift, outflows have been almost exclusively detected in ionized gas (e.g., Barger et al. 1999; Nesvadba et al. 2007; Newman et al. 2012), so it remains unclear how intense winds affect the molecular gas reservoirs. This emphasizes the need to increase the number of high-redshift dusty star-forming galaxies with measurements of the local gas kinematics and of the stellar mass, gas mass and star-formation surface densities (e.g., Hatsukade et al. 2015). This must be done down to typical disk-fragmentation scales (Toomre 1964; Escala & Larson 2008), in order to probe local energy injection from a range of feedback processes and to characterize the resolved Schmidt-Kennicutt law (e.g., Swinbank et al. 2011).

Star formation in more than half of the high-redshift dust-obscured galaxies appears to occur within massive giant star-forming clumps of $10^7-10^9 M_\odot$ and size of about 1 kpc or less (e.g., Elmegreen & Elmegreen 2005), embedded within more diffuse disks (Swinbank et al. 2011). These clumpy structures were originally identified in rest-frame UV and optical studies, and their properties play a central role in the overall evolution of the host galaxies (e.g., Tacconi et al. 2013; Mayer et al. 2016; Dessauges-Zavadsky et al. 2017; Cava et al. 2018). Clumps with sufficiently long lifetimes of a few 100 Myr survive the feedback from young stellar populations and could migrate inward to form the central bulge of galaxies (e.g., Ceverino et al. 2010; Bournaud et al. 2014). However, processes such as clump mergers, gas accretion, dynamical interactions within the disk and star-formation feedback could lead to their dissolution on much shorter timescales (e.g., Tamburello et al. 2015). A better understanding of high-redshift star-formation process therefore requires us to resolve the gas, dust and stellar properties of these clumps. This is best achieved in very strongly gravitationally lensed galaxies for which the lensing magnifications extend the resolution limits of current facilities and boost the apparent source brightness.

Here we present optical and near-infrared imaging and submillimeter and millimeter interferometry of PLCK_G165.7+49.0, a strongly gravitationally lensed dusty star-forming galaxy at $z = 2.236$. This source was identified as part of our Planck's Dusty Gravitationally Enhanced subMillimeter Sources (GEMS) follow-up program of 11 of the brightest high-redshift galaxies on the submillimeter sky discovered with the Planck all-sky survey and Herschel space observatory (Planck Collaboration Int. XXVII 2015; Cañamas et al. 2015; Planck Collaboration Int. XXXIX 2016). PLCK_G165.7+49.0 comprises a bright submillimeter arc, which we refer to as the “Emerald”, near several other extended arcs falling behind a rich foreground environment, as recently discovered with CFHT imaging (Cañamas et al. 2015; Cañamas 2016).

The luminous dusty starburst galaxy PLCK_G165.7+49.0 has apparent far-infrared (FIR) luminosity of $\mu L_{FIR} = (1.0 \pm 0.1) \times 10^{14} L_\odot$ originating from $\mu M_d = (5.1 \pm 0.1) \times 10^9 M_\odot$ of dust heated to a temperature of $T_d = 42.5 \pm 0.3$ K (Cañamas et al. 2015, hereafter C15), where $\mu$ indicates the gravitational magnification factor. The far-infrared radio correlation does not suggest the presence of a radio-loud AGN in this system, and photometric constraints from WISE at 22 $\mu$m, IRAS at 60 and 100 $\mu$m, and SPIRE at 250 $\mu$m do not suggest more than at most a few percent AGN contamination to the overall FIR luminosity (C15). We also detected luminous CO(3-2) line emission with the wide-band heterodyne receiver EMIR on the 30-m telescope of IRAM,

with an integrated flux of $\mu I_{CO} = 25.4 \pm 0.3$ Jy km s$^{-1}$ and a line full-width-at-half-maximum (FWHM) of 580 km s$^{-1}$ (C15).

To further characterize this source, we obtained CFHT and Spitzer optical and near-infrared imaging, as well as subarcsecond submillimeter and millimeter interferometry of the dust and CO(4-3) line emission with the Submillimeter Array (SMA) and the IRAM Plateau de Bure Interferometer (PdBI). We used these data to characterize the foreground lensing potential, a hitherto unknown massive galaxy cluster at $z = 0.348$ with an adjacent filament, to calculate a strong lensing model with LENSTOOL, and to characterize the gas kinematics and spatially resolved dust and star formation properties of the background source. To constrain the lensing model and characterize the stellar components in the rest-frame UV, including their morphologies, we also used recently obtained HST/WFC3 imaging through the F110W and F160W filters, which are described in more detail in Frye et al. (2018).

We present our analysis as follows. In Sect. 2 we describe our observations, data reduction, and the construction of the photometric catalogs and spectral imaging maps. In Sect. 3, we characterize the foreground structure with a three-way approach, by quantifying the local overdensity using an adaptive kernel density estimator, by identifying the red sequence of passively evolving member galaxies of the foreground structure, and by estimating photometric redshifts. In Sect. 4 we compute a strong lensing model based on the results of Sect. 3 and the positions and brightnesses of foreground and background galaxies. In Sect. 5 we characterize the intrinsic stellar, dust, gas and star formation properties of the Emerald. In Sect. 6, we discuss the stability of the star-forming clumps, present the first detection of a molecular wind at high-redshift, and investigate whether it results from the kinetic energy and momentum injection from star formation. We then conclude with a summary in Sect. 7.

Throughout the paper we have adopted the flat $\Lambda$CDM cosmology from Planck Collaboration XIII (2016), with $H_0 = 67.81$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.308$, and $\Omega_\Lambda = 1 - \Omega_M$. At the redshift $z = 2.236$ of PLCK_G165.7+49.0, this corresponds to a luminosity distance $d_{L,bg} = 18.25$ Gpc, with a projected scale of $8.40$ kpc arcsec$^{-1}$. With the same cosmology, the luminosity distance of a foreground source at $z = 0.348$ is $d_{L,fg} = 1.90$ Gpc. At that redshift, $5.07$ kpc corresponds to $1''$ on the sky. All magnitudes are in the AB system.

## 2. Observations and data reduction

### 2.1. Optical and near-infrared imaging

We obtained optical and near-infrared wide-field imaging of PLCK_G165.7+49.0 and surrounding sky with MEGACAM and WIRCAM on the Canada-France-Hawaii Telescope (CFHT) through the r-, z-, J-, and Ks-band filters, and with the IRAC camera through the 3.6 $\mu$m and 4.5 $\mu$m filters on the Spitzer Space Telescope.

At the CFHT, PLCK_G165.7+49.0 was observed during several nights between March and May 2014, as part of program 14AF06 (PI: Nesvadba). We obtained a total of 40 min and 49 min of on-source observing time with MEGACAM through the r- and z-band filters, respectively, and 93 and 52 min with WIRCAM through the J- and Ks-band filters, respectively. The seeing was between 0.8 and 1.1" in the optical and the J-band, and 0.7" in the Ks-band. The near-infrared detectors were read out once every 10s in J, and once every 15 s in Ks.
Optical and near-infrared (NIR) images were bias and dark-frame subtracted, respectively, and all frames were flat-fielded before being released to the principal investigator, as is customary at the CFHT. We used these preprocessed frames, and subtracted the sky from the near-infrared images by averaging over the ten frames that had been taken most closely in time to a given frame, then subtracted the average from this science frame. Individual frames were aligned relative to each other and to the world coordinate system with the astrometric tools SWARP and SCAMP (Bertin 2010a,b), resampled to 0.3" pixel scale, cropped to 5' × 5', and flux calibrated relative to the Sloan Digital Sky Survey (SDSS, Alam et al. 2015) and the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2006) for the optical and near-infrared images, respectively. We obtained a relative calibration and zero-point uncertainties below 0.05 mag in each band, by fitting the spectral energy distribution of nearby non-saturated stars with blackbody spectral energy distributions.

The IRAC images were obtained as part of program 10010 (PI: Nesvadba) during *Spitzer* observing cycle 10 on 8 July 2014, and were observed through the 3.6 µm and 4.5 µm filters as part of the warm mission. The total observing time per band was 1200 s, composed of individual exposures of 100 s duration. Basic calibrated data were released after a preliminary processing conducted by the *Spitzer* Science Center standard pipeline. The dark currents and flat fields were automatically calibrated and subtracted during this stage. In both channels the frames were flux calibrated, combined into mosaics with 0.60" pix⁻<fnref n="1" /> sampling, corrected for cosmic ray artifacts and astrometrically calibrated with external 2MASS catalogs. After they were released to us, we again used SWARP and SCAMP to put these images onto a common reference frame with our ground-based data.

## 2.2. Aperture photometry

We selected our sources from the *Ks*-band imaging to approximate a mass selection at intermediate redshifts, and used SExTRACTOR (Bertin & Arnouts 1996) to measure aperture magnitudes of 3.9" diameter in all CFHT bands. We verified carefully that these apertures were large enough to minimize flux losses, while being small enough that sources were not blended. For extended sources, we used the corrected isophotal magnitudes from SExTRACTOR computed down to the 3σ isophotes. Positional uncertainties relative to the *Ks* band have rms ≈ 0.1", about 10–15% of the FWHM size of the point spread function (PSF).

The IRAC images have substantially larger PSFs, with FWHM sizes of around 1.7", so that blending becomes more important, in particular in the denser regions of our foreground structures. We measured magnitudes within the same apertures of 3.9" as for the ground-based data, and applied aperture correction factors of 1.4 and 1.5 in the 3.6 and 4.5 µm bands, respectively, following Barmby et al. (2008) and Martinache et al. (2018). Positional uncertainties are between 0.2" and 0.5" relative to the *Ks*-band image.

The final catalog includes 737 objects down to 3σ limiting AB magnitudes of 25.5, 23.7, 23.9 and 23.0 mag in the *r*, *z*, *J*, and *Ks* bands, respectively. In the 3.6 µm and 4.5 µm *Spitzer* channels, the limiting magnitudes are 23.5 and 23.4 mag.

## 2.3. SMA 850 µm interferometry

The 850 µm continuum (345 GHz) of PLCK_G165.7+49.0 was observed with the SMA on 12 December 2013, with a total integration time of seven hours as part of program 2013B-S050 (PI:

Nesvadba) in the compact (COM) configuration at 339.15 GHz, and with another seven hours as part of program 2016B-S005 (PI: Nesvadba) in the extended (EXT) configuration at the same frequency on 25 March 2017. Combining data from both runs gives a beam with FWHM size of 0.9" × 0.75" at PA = 39°. The data were taken under excellent conditions, with precipitable water vapor below 2 mm, and with individual scan durations of 30 s beam⁻<fnref n="1" />. In the COM run, we used Callisto as flux calibrator, and 1153+495 and 1159+292 were used as phase calibrators with the bandpass calibrated on 3C 279. In the EXT run, we used 1159+292 and 1146+399 as phase calibrators, Titan as flux calibrator, and 3C84 to calibrate the bandpass.

Data were reduced and calibrated with the MIR package, and Fourier transformed and deconvolved with MIRIAD. Images were created using “Briggs” weighting, with a parameter robust = 0.5 on the MIRIAD task invert, as is standard at the SMA. A comparison between the COM maps and the single-dish flux densities from SCUBA-2 suggests that we recover at least 80–90% of the total flux density at 850 µm (C15). We reached an rms noise of about 1.6 mJy beam⁻<fnref n="1" /> using the SMA.

We also observed PLCK_G165.7+49.0 with the very extended (VEXT) configuration of SMA on 20 January 2015, with a total integration time of seven hours and 0.3" beam size. The rms reached with these observations was 1.75 mJy beam⁻<fnref n="1" />. PLCK_G165.7+49.0 was not detected in these observations, suggesting that the flux was either resolved out or that the surface brightness was too faint even at the center of the clumps to be measured with the small beam size. At any rate, the non-detection suggests that there are no high surface-brightness clumps of dust emission in PLCK_G165.7+49.0 that are more compact than about 0.3" in the image plane.

## 2.4. IRAM CO(4–3) interferometry and spectral line fitting

The CO(4–3) line from PLCK_G165.7+49.0 was observed with the PdBI with six antennas in the B configuration in Band 2 on March 19 2014. The shortest and longest baselines in the data set are 88 m and 452 m and thus, the observations are sensitive to scales smaller than approximately 5". The phase center of the observations was located at α=11h27m14.60s and δ = +42°28′25.0". The 2 mm receivers were tuned to a sky frequency of 142.47 GHz, corresponding to the rest-frame wavelength of CO(4–3) at a redshift z = 2.23606. The WideX correlator with its 3.6 GHz bandwidth at a spectral resolution of 1.95 MHz provided a velocity coverage of 7500 km s⁻<fnref n="1" /> with 4.1 km s⁻<fnref n="1" /> wide channels. The bright quasar 3C 84 was used for bandpass calibration, LkHa 101 was observed as primary flux calibrator, and we regularly observed the nearby quasars 1150+497 and 1128+385 for gain and phase calibration.

The data were calibrated and imaged within GILDAS/CLIC and MAPPING<fnref n="1" />. A few outliers in the visibilities of channel 42 were removed using the uvClip task in MAPPING. To image the data, we used the standard clean procedure together with a mask that was carefully adapted to each individual frequency channel, and applied “natural” weighting. This resulted in a beam size of 0.76" × 0.75" at a position angle of 92°. The data reach an rms of 0.6 mJy beam⁻<fnref n="1" /> per spectral channel, which are 42 km s⁻<fnref n="1" /> wide. After calibration and imaging the data cube was exported as a fits file for analysis.

We used MPFIT (Markwardt 2009) to construct maps of the line fluxes, local velocities relative to z = 2.23606, and FWHM

<footnote n="1">¹ http://www.iram.fr/IRAMFR/GILDAS</footnote>
line widths from our CO(4–3) data cube. These maps are shown in Fig. 8. In a first step, we fit a single Gaussian component to spectra extracted from apertures of 3 × 3 spatial pixels, or 0.6" × 0.6", slightly less than the beam size. This maximizes the signal-to-noise ratio while causing no loss of spatial information. We only fit pixels in which the line was detected at ≥3σ. In a small region near the critical line and extending to about 1.3" from it on either side (see Fig. 8), single-Gaussian fits lead to significant residuals (≥3σ). We adopted a two-component fit in this region to include this component in our analysis. The integrated spectrum and maps corresponding to this secondary component are shown in Fig. 9 and in the lower panel of Fig. 8, respectively. The total CO flux extracted from this region is μ I<sub>CO(4–3)</sub> = 19.5 ± 0.4 Jy km s<sup>-1</sup>, corresponding to a total luminosity, μ L'<sub>CO(4–3)</sub> = 2.9 × 10<sup>11</sup> K km s<sup>-1</sup> pc<sup>2</sup>, uncorrected for the gravitational magnification factor μ.

## 2.5. Ancillary data sets

We obtained HST/WFC3 imaging through the F110W and F160W filters for a subset of the GEMS through program 14223 (PI: Frye). For details of the observations and data reduction see Frye et al. (2018). A major justification of that observing program was to refine the strong lensing model obtained with LENSTOOL as presented in the present paper, which depends critically on measuring accurate positions, colors and morphologies of the multiply imaged sources from high-resolution imaging, as HST ideally provides. We therefore included these data in our photometric catalogs, and based our lens model on the colors and positions measured from these HST images, which have a point spread function with a FWHM size of about 0.15".

We also used spectra from the 12th data release of the Sloan Digital Sky Survey (Alam et al. 2015), available for four targets, to validate our photometric redshifts through spectroscopy, and to confirm the cluster redshift obtained through a red-sequence analysis. Three red-sequence galaxies also have SDSS spectra, including the two main galaxies of the groups directly adjacent to the gravitationally lensed arcs (see Fig. 2). Their spectral properties are listed in Table 1.

# 3. Multiple extended arcs behind a rich foreground environment

Since the GEMS were unknown targets previous to our follow-up observations, and do not fall into well characterized regions of the night sky, we must apply some care in characterizing the foreground mass distributions that are responsible for the gravitational lensing effect.

Figure 1 shows a wide view of the field around the line of sight toward PLCK_G165.7+49.0, revealing multiple red, strongly gravitationally lensed arcs. These arcs are mainly distributed along two extended curves with projected lengths of 28.5" and 21", respectively, which we will refer to as the western and eastern arcs, respectively, and which fall between two groups of early-type galaxies that have similar SDSS spectroscopic redshifts around z = 0.35 (Table 1) and might represent substructure within a galaxy cluster. This makes the field of PLCK_G165.7+49.0 the richest environment toward any of our GEMS. Both arcs are also resolved in the narrow direction, with widths of typically about 0.6".

Only subcomponents of the western arc are associated with dust continuum and CO line emission, probed with the SMA and PdBI, and will be discussed in more detail in Sect. 5. The stellar

<table><caption>**Table 1.** Characteristics of three massive galaxies within the foreground structure of PLCK_G165.7+49.0 taken from the 12th data release of the Sloan Digital Sky Survey.</caption><thead><tr><th>ID</th><th>RA<br>(J2000)</th><th>Dec<br>(J2000)</th><th>z<sub>spec</sub></th><th>σ<sub>abs</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>G1</td><td>11:27:16.59</td><td>+42:28:41.0</td><td>0.34788 ± 0.00007</td><td>323 ± 22</td></tr><tr><td>G2</td><td>11:27:16.69</td><td>+42:28:38.2</td><td>0.33767 ± 0.00004</td><td>307 ± 14</td></tr><tr><td>G3</td><td>11:27:13.75</td><td>+42:28:22.6</td><td>0.34770 ± 0.00012</td><td>271 ± 35</td></tr></tbody></table>

**Notes.** The columns are: source name; right ascension and declination; spectroscopic redshift; and velocity dispersion measured from stellar absorption lines.

continuum emission along the western arc is globally diffuse and filamentary, with one pronounced clump near the center of the component aligned with the submillimeter emission (see Fig. 1). One subcomponent of the western arc is particularly bright in the submillimeter, producing most of the dust continuum emission of PLCK_G165.7+49.0. Its integrated 850 µm flux density is 48 mJy, about two thirds of the total 850 µm flux density of 71.6 mJy measured with the SMA within the half-power beam width of SCUBA-2. In the following discussion we refer to this submillimeter arc as “the Emerald”.

The CO interferometry shows that this arc probes a single source at a redshift z = 2.236, with a velocity gradient consistent with the presence of two merging images, as further discussed in Sect. 4. Toward the southeast and northwest are two additional, compact and more moderately magnified CO emitters, which also fall within the 20"–30" beam of SPIRE and typical single-dish observations in the submillimeter and millimeter wavelength range. We will refer to these sources as “Co–S” and “Co–N” (see Fig. 1). They are also detected in CO(4–3) and are at similar redshifts to the Emerald.

We followed Song et al. (2012a,b) in adopting a three-way approach to characterize this environment in a more rigorous, quantitative way. Firstly, we calculated the projected density along the line of sight, using the adaptive kernel density estimation (AKDE) algorithm, which adopts a non-parametric, scale-independent smoothing technique to calculate the local density around each detected source (see Ferdosi et al. 2011; Pisani 1996; Valtchanov et al. 2013). We used the package already described by Valtchanov et al. (2013) to determine the position and significance of overdensity peaks. Secondly, we searched for the presence of a red sequence in optical and near-infrared color-magnitude diagrams, which would be the clearest signature of a dense region in three dimensions, and provide robust and accurate redshift constraints (Fassbender et al. 2008, 2011). Thirdly, we used the publicly available Bayesian photometric redshift package (BPZ, Benítez 2000) to estimate photometric redshifts and study the redshift distribution along the line of sight. We will start by describing the analysis of the local source density projected onto the sky around the sightline toward PLCK_G165.7+49.0, and come back to the second and third steps in the following two subsections.

## 3.1. Local projected source density

As a first step to quantify the nature of this structure, we estimated the density distribution of sources selected in the CFHT Ks-band within a field-of-view of 5' × 5' around PLCK_G165.7+49.0. We applied the AKDE on the 519 detections, finding a density peak of 6.4σ significance at about 10"
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Top panel: HST/WFC3 F110W and F160W band two-color image of the Planck's Dusty GEMS PLCK_G165.7+49.0 and surrounding field of view. The orange solid line shows the critical curve at z = 2.236, for the gravitational lensing potential from our best-fitting lens model, derived using the position and multiplicity of the lensed images identified with HST and the PdBI (see text and Table 3). Black diamonds indicate the members of the foreground structure that were included in the lens model, and white triangles mark the position of the two compact submm emitters. The left inset shows a 9" × 9" wide enlargement of the submm arc at z = 2.236, "the Emerald". In the right inset, we show the dust continuum at 850 µm from the SMA, with contours starting at 4σ and increasing in steps of 4σ; the beam size is shown in the lower left corner and the blue arrow marks the position of the stellar continuum clump further discussed in the text. The bar in the upper left corner shows a projected distance of 10", corresponding to 50 kpc at z = 0.348 and 84 kpc at z = 2.236. Bottom panel: 15" × 15" postage stamps centered on the Emerald (blue symbols) from our optical and infrared imaging obtained with CFHT and Spitzer.</figcaption></figure>

from the Emerald, and secondary peaks of 4.6σ and 4.1σ at sep-
arations of 1.3' and 2.9', respectively. We took advantage of the
exceptionally large field-of-view of WIRCAM of 20' × 20' to
estimate the background standard deviation of the AKDE, σ,
at the same depth and observing conditions. We also applied
the AKDE to 5' × 5' wide areas randomly distributed over the
WIRCAM field-of-view, but far away from the GEMS and

surrounding galaxy overdensities, and adopted the median den-
sity in these fields as a conservative approximation of the source
density in the field. Comparison with the density contours
obtained in WIRCAM's entire field-of-view indicated that the
size of the field neither affects the position of the AKDE peaks
nor the morphology of the overdensity contours, as expected for
a scale-free method.
<figure><img src="image_3.png" /></figure>

**Fig. 2.** Distribution of near-infrared sources toward PLCK_G165.7+49.0, including candidate members of the foreground structure. The grayscale shows the CFHT $K_s$-band image and blue contours indicate local overdensities from the adaptive kernel density estimate, from 1 to 6σ. A main peak at 6.4σ near the Emerald and two secondary peaks at >4σ are found (blue triangles), indicating a particularly rich environment. The position of the Emerald (the submm arc) is shown as a white diamond. Red and yellow circles indicate galaxies that fall onto the red sequence, or have photometric redshifts consistent at the 2σ level with z = 0.35, respectively. The three galaxies marked G1, G2 and G3 have spectroscopic redshifts around z = 0.35 available in the SDSS, as reported in Table 1. We show the 5' × 5' field-of-view used to characterize the foreground mass distribution. North is up and east is to the left.

The density contours are shown in Fig. 2. In the following subsections, we further characterize the most significant overdensity associated with PLCK_G165.7+49.0. Its 4σ contours define a region with a projected major axis length of 73".

## 3.2. Red sequence analysis and cluster redshift

Massive galaxy populations in dense environments fall within a tight, well defined region in optical and near-infrared color-magnitude diagrams. This “red sequence” is a signature of their early, rapid, and very uniform formation history for a short epoch at high redshift, followed by passive evolution for most of cosmic history. The tightness and uniformity of the red sequence makes it an excellent tool for determining redshifts of the overall structure, often with uncertainties better than $\Delta z = 0.05$ out to redshifts $z \approx 1$ (e.g., Fassbender et al. 2011). In the case of the GEMS, we can use the red sequence to search for the presence of one or multiple massive structures along the line of sight, which could contribute to the gravitational magnification of the GEMS, and to help determine their redshifts.

The $r - K_s$ versus $K_s$ color-magnitude diagram of the field-of-view around PLCK_G165.7+49.0 in Fig. 3 shows a clear red sequence, with colors $r - K_s \approx 2.0$, which is associated with the massive early-type galaxies seen in Fig. 2, and which we can use to determine a redshift of the overall structure, to measure its extent, and to identify the member galaxies.

We followed Song et al. (2012a,b) in modeling synthetic spectra of early-type galaxies using the stellar population synthesis tool of Bruzual & Charlot (2003), and to populate an artificial red sequence matched to that observed in the Coma cluster at $z = 0.0023$ (Eisenhardt et al. 2007). More explicitly, we adopted short starbursts with e-folding time $\tau = 50$ Myr and a Chabrier initial mass function, starting at $z = 3$, and followed by passive

<figure><img src="image_4.png" /></figure>

**Fig. 3.** Color-magnitude diagram of a 5' × 5' field-of-view surrounding PLCK_G165.7+49.0, showing the $r - K_s$ color versus $K_s$-band magnitude. Gray and yellow lines show the expected position of the red sequence for a range of redshifts between $z = 0$ and $z = 1.0$ and between $z = 0.25$ and $z = 0.45$, respectively. The red line indicates the best-fit sequence at $z = 0.35$. Red stars indicate galaxies which have spectroscopic redshifts from the SDSS or Frye et al. (2018) falling in the range $z = 0.350 \pm 0.025$.

evolution until today. We used six different templates from the Padova library covering a wide range of metallicities from 0.05 to 2.5 times the solar metallicity, and extracted the spectral energy distributions for several ages corresponding to redshifts $z = 0.0$ to $z = 1.0$ in steps of $\Delta z = 0.05$. These SEDs were then rescaled to best reproduce the color and slope of the red sequence measured on the Coma cluster, using the eight color-magnitude diagrams from Eisenhardt et al. (2007) with the same set of filters.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Difference between the photometric redshift distribution of galaxies within 1' of PLCK_G165.7+49.0, and the distribution obtained in the rest of the 5' × 5' field, both normalized to the same area (black histogram). This shows a strong excess of sources in the inner field, which is consistent with $z \approx 0.35$ (dotted vertical line). We used only the most robust photometric redshifts, and show the sum of their redshift probability distribution functions obtained with BPZ (in orange).</figcaption></figure>

We redshifted these SEDs to our grid of redshifts out to $z = 1.0$, and convolved them with the transmission curves of the CFHT filters to obtain a range of red sequence models in each color-magnitude diagram. For PLCK_G165.7+49.0, which shows a pronounced red sequence at $r - K_s \approx 2$, we find a best-fitting model at $z = 0.35$, in excellent agreement with the spectroscopic redshifts given in the SDSS for the three galaxies within the groups adjacent to the submm components (Table 1) and with those presented in Frye et al. (2018, see Fig. 3). The galaxies with follow-up spectroscopy in Frye et al. (2018) were selected to fall onto the red sequence shown in Fig. 3. Their redshifts, which all fall into the range expected from a massive cluster, are thus a direct confirmation of our photometry.

We selected galaxies that fall within the limits in the color-magnitude diagram set by our models for $z = 0.3$ and $z = 0.4$, with an additional magnitude cut to only include galaxies with $K_s < 21$. To reduce the number of interlopers, we performed a similar selection in the $J - K_s$ versus $K_s$ color-magnitude diagram, which also exhibits a narrow red sequence at $J - K_s \approx 0.5$, and considered only galaxies that fulfill both criteria as robust members of the cluster red sequence. We find 76 members within the 5' × 5' field. Many of these members fall along a diagonal axis from the northeast to west (as seen in Fig. 2). About half of the red-sequence galaxies lie in a 2.5' × 1.5' wide region within the 3$\sigma$ AKDE contours, further indicating that the overdensity of NIR sources is not due to chance alignment along the line-of-sight. This suggests that PLCK_G165.7+49.0 falls behind a massive, extended galaxy cluster with adjacent filament, which has so far not been identified by other surveys.

### 3.3. Photometric redshifts

We used the BPZ package of Benítez (2000) to estimate photometric redshifts from our $r$, $z$, $J$, and $K_s$-band photometry, which probes the 4000 Å break for a redshift range of the member galaxies of the lensing structure of $z \approx 0.2-1.0$, which appeared most likely prior to our analysis. We also included the shallower and bluer, publicly available SDSS *ugriz* photometry to improve the robustness of our redshift estimates for the brightest targets, or upper limits in the bluer bands.

The BPZ algorithm returns redshift probability distribution functions based on fitting a set of template SEDs. We used the standard set of templates provided by BPZ and, since we are targeting fields with known bright FIR emitters, we also added two strongly reddened SEDs with $A_V = 3$ and 5 mag, respectively. Including the IRAC photometry did not improve the robustness of our estimates, due to the greater photometric uncertainties, and because these wavebands probe relatively flat, featureless spectral regions of the SEDs for galaxies at low and intermediate redshifts, where most of our sources lie. Because of this, we did not include these two bands in our analysis.

In the 5' × 5' field surrounding PLCK_G165.7+49.0, we identified 187 galaxies with reliable photometric redshifts (parameters ODDS $> 0.9$ and $\chi^2 < 10$ in BPZ). Comparison between these photometric redshifts and spectroscopic redshifts of 30 sources, taken from the SDSS and follow-up spectroscopy with MMT/HECTOSPEC and Gemini/GMOS (Frye et al. 2018), showed that our estimates are robust, with an average scatter $|z_{spec} - z_{BPZ}|/(1 + z_{spec}) \approx 0.09$.

As shown in Fig. 2, about 80% of the red-sequence galaxies have photometric redshifts consistent at the $2\sigma$ level with $z = 0.35$, which makes them good candidates for being members of the foreground structure. We computed the difference between the redshift distribution of sources within 1' of PLCK_G165.7+49.0 and associated with the AKDE peak, and the redshift distribution of sources in the rest of the field, after normalizing to the same area (Fig. 4). This shows a strong excess of sources at $z = 0.3-0.4$ toward PLCK_G165.7+49.0 and the presence of a massive structure in this redshift range, which is consistent with the redshift of $z = 0.35$ found from the red sequence analysis and the spectroscopic redshifts of galaxies G1, G2 and G3 from SDSS.

## 4. Gravitational lens modeling

The Planck's Dusty GEMS PLCK_G165.7+49.0 falls behind a very rich, so far unexplored galaxy environment, as discussed in the previous section. The western and eastern extended arcs seen in the $K_s$-band, including the long-wavelength emission from the Emerald, Co-S and Co-N, fall between two compact groups of early-type galaxies at a common redshift $z = 0.35$, which likely probe the inner region of a massive galaxy cluster. The small redshift offset in the SDSS spectra between the three galaxies of $\Delta z = 0.01$ shows that both are part of a single bound structure. Our red sequence and photometric redshift analysis also shows that this structure extends from north-east to west over about 3', perhaps representing a massive filament (Fig. 2). The structure is not detected in the Planck catalog of clusters identified with the Sunyaev-Zel'dovich effect (Planck Collaboration XXXII 2015), but has a faint X-ray counterpart in the Rosat All-Sky Survey (Voges et al. 2000). All of this suggests that PLCK_G165.7+49.0 is magnified by a massive dark-matter halo underlying the cluster at $z = 0.35$, in addition to several member galaxies of this cluster at similar redshifts forming a bimodal mass distribution, consistent with our spectroscopic and photometric constraints.

We performed a strong lensing analysis with LENSTOOL (Jullo et al. 2007), by modeling the mass distribution toward the Emerald and neighboring western and eastern arcs. LENSTOOL is a publicly available Bayesian lens modeling package, which uses the number of arclets detected in the image plane, their association in multiply imaged systems of the same regions in the source plane, and their positions relative to the critical line, in
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left panel: As Fig. 1, with labels showing the multiple images identified for our gravitational lens modeling, either using IRAM CO(4–3) interferometry (red squares), or HST imaging (white squares). The yellow triangle indicates the position of the additional counter-image of system #8 predicted by the best-fit model. Right panel: Magnification map obtained for the best-fitting lens model, by taking the median of the maps generated by LENSTOOL for each MCMC realization. The inset shows the position of the SMA continuum emission from the Emerald (white contours, starting at +4σ and increasing in steps of +4σ). Cyan contours indicate areas where the average of the relative difference in magnification between the best-fit and each of the four alternative models is higher than 30%. The systematic errors induced by the different mass distributions in each model are much lower at the position of the Emerald thanks to the spectroscopic identification of the submm images.</figcaption></figure>

order to derive a best-fitting lensing potential that is responsible for the gravitational amplification.

## 4.1. Identification of multiply imaged systems

The LENSTOOL software relies on two sets of observational parameters, which can be obtained from the imaging – the number and position of multiple lensed images, and the position, brightness, and structural parameters of the foreground lensing galaxies. The parametrization of the dark-matter halo underlying the galaxy cluster must also be provided.

In the following analysis, we use the HST/WFC3 F110W and F160W images also presented by Frye et al. (2018) to measure the position of the faint lensed arcs and foreground galaxies, and our IRAM CO(4–3) interferometry to constrain the positions of the images for the long-wavelength emitters. All images are magnified by the same lens regardless of the wavelength in which they are studied, so combining both sets of constraints in a joint analysis provides the most robust model of the underlying mass distribution. LENSTOOL determines the best-fit model by minimizing the positional offsets between the measured and reconstructed image positions. We have considered a model adequate when the rms of all offsets, rms<sub>img</sub, is of order 0.1"–0.2", the PSF of the HST/WFC3 imaging.

The eastern and western extended arcs seen in the $K_s$-band are resolved into several fainter arclets in the HST imaging (see Figs. 1 and 2). We identified compact clumps within each arclet and use their F110W–F160W color, their morphology, and spatial distribution to combine them into seven multiply imaged systems with unambiguous associations. Given the lack of spectroscopic redshifts for these near-infrared-selected systems, we

only used those with the most secure image configuration and colors. As shown in Fig. 5, the eastern arc comprises systems #5, #6 and #7, with another nearby system #8, and we identify systems #2, #3, and #4 in the western arc. The LENSTOOL modeling constrains their redshifts, and we require that $z < 4.5$, as suggested by the non-detection of the $Ly\alpha$ emission line in Gemini/GMOS spectra (Frye et al. 2018). The similar colors and low angular separation between systems #2, #3 and #4 suggest that they are subcomponents of a single background galaxy and we therefore assume a common redshift.

In addition, we used the gas kinematics in PLCK_G165.7+49.0 from the CO(4–3) emission line to probe the lensing configuration of the submm components. We identified the NIR counterpart of the Emerald in the F160W image with system #1 at $z = 2.236$. Figure 8 shows that the direction of the velocity gradient is flipped in the northern and southern parts of the submm arc, and so are the distributions of the line widths shown in the same figure. This parity inversion is a clear signature that the Emerald contains two merging images of the same source. This is further shown by the small line widths in the center of the arc, which would be difficult to explain for two independent, and partially overlapping sources. Moreover, the line profiles and velocity offsets of the two fainter, smaller submillimeter images (north and south from the arc, which we label Co–N and Co–S), have different line profiles and velocity offsets (see Table 6 and Fig. 7). While the overall properties of Co–N suggest that this is another, less strongly magnified image of the same galaxy that is also seen in the submm arc, Co–S seems to be another galaxy at 18 kpc projected distance in the source plane. The properties of these galaxies are discussed in more detail in Sect. 5.4.
<table><caption>**Table 2.** Parameters of the foreground mass distribution inferred by LENSSTOOL, for our best lensing model and the four alternative models.</caption><thead><tr><th>Model</th><th>rms<sub>img</sub><br>(")</th><th>Component</th><th>ΔRA<br>(")</th><th>ΔDec<br>(")</th><th>ε</th><th>θ<br>(deg)</th><th>r<sub>core</sub><br>(kpc)</th><th>r<sub>cut</sub><br>(kpc)</th><th>σ<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">Best</td><td rowspan="3">0.21</td><td>Large scale<br>G4</td><td>39.9 ± 4.0<br>[7.0]</td><td>-20.3 ± 2.2<br>[-1.5]</td><td>0.69 ± 0.07<br>[0.1]</td><td>-33 ± 1<br>[-43]</td><td>142 ± 10<br>[0.25]</td><td>[500]<br>[70]</td><td>1081 ± 112<br>233 ± 14</td></tr><tr><td>L* galaxy</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[0.25]</td><td>59 ± 17</td><td>267 ± 17</td></tr><tr><td>Fixed center</td><td>0.64</td><td>Large scale<br>G4</td><td>[31.2]<br>[7.0]</td><td>-18.8<br>[-1.5]</td><td>0.73 ± 0.03<br>[0.1]</td><td>-31 ± 1<br>[-43]</td><td>59 ± 5<br>[0.25]</td><td>[500]<br>[70]</td><td>759 ± 25<br>293 ± 9</td></tr><tr><td rowspan="2">Faint arcs</td><td rowspan="2">0.43</td><td>L* galaxy</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[0.25]</td><td>40 ± 19</td><td>213 ± 30</td></tr><tr><td>Large scale<br>G4</td><td>28.6 ± 1.0<br>[7.0]</td><td>-14.4 ± 0.7<br>[-1.5]</td><td>0.70 ± 0.09<br>[0.1]</td><td>32 ± 1<br>[-43]</td><td>82 ± 12<br>[0.25]</td><td>[500]<br>[70]</td><td>851 ± 58<br>243 ± 9</td></tr><tr><td rowspan="2">Non-cored</td><td rowspan="2">0.34</td><td>L* galaxy</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[0.25]</td><td>63 ± 23</td><td>296 ± 28</td></tr><tr><td>Large scale<br>G4</td><td>24.9 ± 0.5<br>[7.0]</td><td>-11.6 ± 0.5<br>[-1.5]</td><td>0.64 ± 0.09<br>[0.1]</td><td>31 ± 2<br>[-43]</td><td>[20]<br>[0.25]</td><td>[500]<br>[70]</td><td>573 ± 26<br>284 ± 11</td></tr><tr><td rowspan="3">NFW</td><td rowspan="3">0.30</td><td>L* galaxy</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[0.25]</td><td>133 ± 29</td><td>301 ± 16</td></tr><tr><td>Large scale<br>G4</td><td>37.8 ± 1.6<br>[7.0]</td><td>-20.1 ± 0.9<br>[-1.5]</td><td>0.68 ± 0.05<br>[0.1]</td><td>-32.5 ± 0.6<br>[-43]</td><td>3.4 ± 0.7<sup>(1)</sup><br>[0.25]</td><td>377 ± 53<sup>(2)</sup><br>[70]</td><td>...</td></tr><tr><td>L* galaxy</td><td>...</td><td>...</td><td>...</td><td>...</td><td>[0.25]</td><td>74 ± 21</td><td>296 ± 24</td></tr></tbody></table>

**Notes.** “Large scale” refers to the extended dark-matter halo associated with the lensing structure. Positional offsets are given in arcseconds relative to $\alpha = 11h27m16.6s$ and $\delta = +42^\circ 28'38.8''$. The ellipticity of the mass distribution, $\epsilon$, is given as $(a^2 - b^2)/(a^2 + b^2)$. Parameters in brackets are fixed, and errors correspond to $1\sigma$ confidence intervals for the best-fit parameters. For the NFW model, (1) and (2) refer to the concentration parameter, $c$, and the scale radius, $r_s$, respectively.

In total, the arc and the other seven systems at unknown redshifts provide us with 28 constraints, which we can use to infer the properties of the lensing potential. The positions of all images included in the analysis are listed in Table 3 and plotted in Fig. 5.

## 4.2. Foreground mass distribution modeling

The LENSSTOOL software inverts the lensing equation for the above input parameters, and derives the optimal set of parameters with Markov chain Monte Carlo (MCMC) simulations. We used 3000 MCMC iterations, and a large-scale dark-matter halo underlying the overall structure as the main deflector, together with the perturbations from smaller halos associated with the cluster members. We describe each mass component with a pseudo-isothermal elliptical mass distribution (PIEMD, Elíasdóttir et al. 2007), which has a radial profile characterized by a central mass surface density, a central velocity dispersion, and core and cutoff radii, $r_{core}$ and $r_{cut}$ (Jullo et al. 2007). We note that the velocity dispersion of the PIEMD is not numerically identical to the physical velocity dispersion as measured from stellar absorption lines (see Limousin et al. 2007a, for further details).

The large number of available constraints allows us to find the best solutions for the projected position on the sky, the ellipticity and the position angle of the main underlying dark-matter halo, as well as its velocity dispersion and the core and cutoff radii. Dark-matter halos are also assigned to individual galaxies on the red sequence, considering only galaxies within 1' of PLCK_G165.7+49.0, and rejecting those with photometric redshifts that are inconsistent with $z = 0.35$ at $\ge 2\sigma$. The selection includes four galaxies with spectroscopic confirmation, three from the SDSS (G1, G2 and G3) and one from MMT/HECTOSPEC (Frye et al. 2018). In total, we include 21 galaxies as perturbers in our model (shown with yellow and red circles in Fig. 2).

To find the best-fitting foreground mass distribution we favor the simplest parametrization that provides a rms consistent with our requirements. We were unable to reproduce the position of lensed images and arcs without including a massive dark-matter

component centered near the western group that underlies the overall structure. Given the lack of spectroscopic redshifts for systems #5, #6, #7 and #8, we did not find a solution where adding a second potential associated with the eastern group would have improved the fit. Models using this second large-scale dark-matter halo are poorly constrained and result in mass components with unrealistically high ellipticities. We therefore describe the mass distribution toward PLCK_G165.7+49.0 using a single large-scale dark-matter halo and galaxy-scale perturbers.

The cutoff radius of the main halo is not well constrained by the lensing configuration, and we therefore set it to 500 kpc. We tested carefully that the outcome of the lensing model depends only very weakly on the precise value of this parameter. We varied the ellipticity of the mass distribution between 0 and 0.8, the core radius between 30 and 200 kpc, and the velocity dispersion between 400 and 2000 km s$^{-1}$. We also allowed the position of the potential to vary within $\pm 10''$ with respect to the center of the western group.

The position, ellipticity and position angle of individual galaxy halos are matched to the light profiles in the $K_s$-band, while we let their cutoff radii and velocity dispersions scale with the galaxy luminosity following the two relationships

$$r_{\text{cut}} = r_{\text{cut}}^{*} \left( \frac{L}{L^{*}} \right)^{1/2} \quad \text{and} \quad \sigma = \sigma^{*} \left( \frac{L}{L^{*}} \right)^{1/4} \qquad (1)$$

Here $L$ is the luminosity of individual galaxies, and $L^*$ the characteristic luminosity of a galaxy at $z = 0.35$. We adopted $K = 16.0$ for an $L^*$ galaxy at $z = 0.35$ (de Propris et al. 1999) and varied the associated characteristic cutoff radius, $r_{cut}^*$, and velocity dispersion, $\sigma^*$, between 50 and 150 kpc and between 150 and 300 km s$^{-1}$, respectively, following, for example, Limousin et al. (2007a) and Richard et al. (2014). We held their core radius fixed at $r_{core} = 0.25$ kpc, as usually done in comparable studies in the literature (e.g., Brainerd et al. 1996; Limousin et al. 2007b; Richard et al. 2014). The rms of systems #5 to #8 is dominated by the dark-matter halo mass profile of a single foreground galaxy, labeled “G4” in Fig. 5. Accordingly, we determined the velocity dispersion of this halo separately.
<table><caption>**Table 3.** Multiply imaged systems used to calculate the best-fit model of the foreground mass distribution.</caption><thead><tr><th rowspan="2">Image ID</th><th rowspan="2">RA</th><th rowspan="2">Dec</th><th rowspan="2">z<sub>spec</sub></th><th colspan="5">z<sub>opt</sub></th></tr><tr><th>Best</th><th>Fixed center</th><th>Faint arcs</th><th>Non-cored</th><th>NFW</th></tr></thead><tbody><tr><td>1.1</td><td>171.81167</td><td>42.472683</td><td>2.236</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>1.2</td><td>171.81128</td><td>42.473342</td><td>2.236</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>2.1</td><td>171.81197</td><td>42.471372</td><td>—</td><td>2.2 ± 0.1</td><td>3.4 ± 0.3 <sup>(1)</sup></td><td>2.2 ± 0.2 <sup>(2)</sup></td><td>1.1 ± 0.5</td><td>1.9 ± 0.2</td></tr><tr><td>2.2</td><td>171.80998</td><td>42.474375</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>2.3</td><td>171.80819</td><td>42.475792</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>3.1</td><td>171.81230</td><td>42.471500</td><td>—</td><td>2.3 ± 0.1</td><td>3.4 ± 0.3 <sup>(1)</sup></td><td>2.2 ± 0.2 <sup>(2)</sup></td><td>1.0 ± 0.4</td><td>2.0 ± 0.2</td></tr><tr><td>3.2</td><td>171.81004</td><td>42.474744</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>3.3</td><td>171.80872</td><td>42.475778</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>4.1</td><td>171.81252</td><td>42.471500</td><td>—</td><td>2.4 ± 0.1</td><td>3.4 ± 0.3 <sup>(1)</sup></td><td>2.2 ± 0.2 <sup>(2)</sup></td><td>1.0 ± 0.4</td><td>2.0 ± 0.2</td></tr><tr><td>4.2</td><td>171.80974</td><td>42.475153</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>4.3</td><td>171.80925</td><td>42.475558</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>5.1</td><td>171.81668</td><td>42.474678</td><td>—</td><td>2.2 ± 0.4</td><td>3.6 ± 0.4</td><td>2.0 ± 0.3</td><td>1.2 ± 0.6</td><td>1.8 ± 0.4</td></tr><tr><td>5.2</td><td>171.81513</td><td>42.476181</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>5.3</td><td>171.81402</td><td>42.478100</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>6.1</td><td>171.81725</td><td>42.474458</td><td>—</td><td>2.1 ± 0.3</td><td>3.5 ± 0.3</td><td>2.0 ± 0.2</td><td>1.2 ± 0.5</td><td>1.7 ± 0.4</td></tr><tr><td>6.2</td><td>171.81499</td><td>42.476697</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>6.3</td><td>171.81440</td><td>42.477931</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>7.1</td><td>171.81637</td><td>42.474692</td><td>—</td><td>2.3 ± 0.3</td><td>3.0 ± 0.3</td><td>2.0 ± 0.3</td><td>1.2 ± 0.5</td><td>1.8 ± 0.5</td></tr><tr><td>7.2</td><td>171.81503</td><td>42.475978</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>8.1</td><td>171.81590</td><td>42.478017</td><td>—</td><td>3.5 ± 0.3</td><td>1.0 ± 0.6</td><td>3.5 ± 0.4</td><td>1.8 ± 0.5</td><td>3.5 ± 0.3</td></tr><tr><td>8.2</td><td>171.81584</td><td>42.478325</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>8.3</td><td>171.81591</td><td>42.477315</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>9.1</td><td>171.82636</td><td>42.481387</td><td>—</td><td>—</td><td>—</td><td>1.5 ± 0.8</td><td>—</td><td>—</td></tr><tr><td>9.2</td><td>171.82567</td><td>42.481899</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

**Notes.** Colors and positions of most images are measured from the HST imaging, except for system #1 which represents the PdBI CO(4–3) morphology of the Emerald, the bright submm arc. Systems #1 to #8 are the most securely identified from their position, color and morphologies and were used in each of the five models, while system #9 was only considered in one alternative model and ignored in the remaining four (Ø symbols). Here $z_{opt}$ refers to the redshifts derived from our lens models, and their $1\sigma$ confidence intervals, with indices (i) indicating systems that were fitted with a common redshift. Measuring spectroscopic redshifts for other multiply imaged systems would allow us to discriminate between models and to better constrain the foreground dark-matter distribution.

The best-fit mass model reproduces the image positions with an rms<sub>img</sub> of 0.21". The modeled parameters of the PIEMDs are summarized in Table 2. The main dark-matter halo broadly follows the light distribution, as it is oriented on the same axis as the filamentary structure. Its ellipticity is close to the upper range expected from cosmological simulations (Despali et al. 2017). The main halo has a large core, with $r_{core} > 100$ kpc, and is offset by about 8" toward the east from the center of the western group of galaxies. This is most likely due to the lack of constraints on the opposite side of the potential, as previously encountered in other studies of strong lensing clusters. The best-fit potential therefore corresponds to a bimodal mass distribution induced by the two groups of cluster members, with the large-scale dark-matter halo producing an additional convergence term.

In Fig. 6 we show the critical line and image-plane morphology of the best-fit model reconstructed by LENSTOOL, as well as the internal and external caustic lines and morphology of the Emerald in the source plane. In addition to images previously identified to constrain the mass distribution, the model predicts a fourth counter-image of system #8, within 1" from a faint near-infrared source detected with HST (see Fig. 5). It also predicts a third counter-image to the Emerald (system #1), at a position consistent with that of Co-N in the dust and CO(4–3) maps.

We explored the systematic errors on the magnification factors and source plane properties of the Emerald induced by the mass parametrisation and identification of multiple images,

through deriving a grid of alternative models. We successively included and excluded some multiple image systems without spectroscopic redshifts, and some galaxies from the scaling relations. We also fitted and restricted different parameters of the large-scale PIEMD, and we tested different mass profiles (following, e.g., Limousin et al. 2016). After excluding models resulting in unphysical mass distributions and/or predicting a greater number of bright images than observed in HST bands, we obtained four reasonable models with rms<sub>img</sub> in the range 0.2"–0.6", with the following features:

1. we used a large-scale PIEMD with position fixed to the center of the western group;

2. we tested the effect of adding candidate counter-images on the opposite side of the foreground groups (see Table 3);

3. we derived a non-cored model with $r_{core}$ fixed to 20 kpc for the main halo;

4. we determined the impact of the mass-density slope degeneracy by replacing the main PIEMD with an NFW profile (while still describing individual member galaxies with PIEMDs).

All models produce similar numbers of images at the correct positions, although they have different underlying mass distributions (Table 2) and predicted redshifts for the multiple images (Table 3). Our current photometric identification of multiply imaged system thus induces some degeneracies in the model which prevent us from deriving robust constraints on
the underlying dark-matter distribution. However, this is not the main focus of this paper. Instead, we estimated the corresponding systematic uncertainties in magnification, by computing the absolute value of the difference between the magnification map of the best-fit model and the four alternatives. Although the differences in the predicted magnification factors are significant for some systems of multiple images in Table 3, due to the different properties of the main lensing potential (see also Limousin et al. 2016), they are not a major concern for the resulting gravitational magnification of submm emitters, which is our focus in this paper. Figure 5 illustrates that the average difference in magnification per pixel remains below 30% toward all components of PLCK_G165.7+49.0 including the Emerald; this is because their spectroscopic identification in the submm forces all models to converge locally.

The F110W and F160W WFC3 images presented in Fig. 1 illustrate that near-infrared emission from the Emerald is extended and diffuse, with fairly uniform surface brightness that varies by not more than about a factor of two or three, and over small scales comparable to the size of the PSF. The only exception is a brighter and very circular clump near the critical line, which has a significantly bluer color than the surrounding stellar emission (blue arrow on Fig. 1). The FWHM size of this clump is 0.45" × 0.33" and it is therefore well resolved with WFC3 at 0.15" × 0.15" PSF. While the diffuse continuum is clearly part of the high-redshift galaxy, the clump could in principle also be an intervening dwarf galaxy, perhaps a member of the main lensing cluster. The hypothesis of an intervening galaxy would in particular explain its very round, symmetric morphology. To evaluate the potential impact of this source for the lensing geometry, we ran LENSTOOL for both hypotheses, finding that this clump, which might be a dwarf galaxy at z = 0.35, would introduce a brightness difference of a factor of 0.8 between images #1.1 and #1.2 from the arc (as observed), but also distort image #1.2 in ways that are not compatible with the stellar morphology seen with HST. The impact of an intervening source at most other redshifts would be even smaller. While the observational evidence is thus not conclusive, neither scenario introduces major systematic uncertainties into our analysis.

## 4.3. Magnification factors

To obtain the magnification factors corresponding to our best-fit model in each pixel in the image plane, we calculated a magnification map for each MCMC realization at the redshift of the source, $z = 2.236$, and computed the median of these maps. To determine the intrinsic properties of the stellar continuum, dust continuum at 850 µm, and molecular gas in the Emerald, and to account for possible differential magnification effects, we computed the luminosity-weighted magnification factors, $\mu$, for each component separately, using the same pixel scale as in each image. The three factors are deduced from the best-fitting lensing model combining gas and stellar constraints, by using pixels above a $3\sigma$ threshold in the dust continuum and CO(4–3) line flux maps, and by using pixels included in the SExtractor segmentation map for the F160W band.

LENSTOOL also computes distributions for the magnification factors at the position of each image. For the submm arc, the resulting luminosity-weighted average values and $1\sigma$ confidence intervals are $\mu_{\text{dust}} = 29.4 \pm 5.9$, $\mu_{\text{gas}} = 24.1 \pm 4.8$ and $\mu_{\text{stars}} = 34.1 \pm 6.8$, for the dust, gas and stellar component, respectively. The variations are caused by small morphological differences between the couterparts and small positional offsets between the dust and gas peaks of a few tenths of an arcsecond, less than

the beam size. This shows that the impact of differential lensing between the gas and the dust is not larger than other systematic effects when deriving spatially-integrated results for the Emerald. The difference in magnification between the dust and stellar components are also minor compared to the multiwavelength configurations of other high-redshift SMGs strongly lensed by galaxy clusters (e.g., MacKenzie et al. 2014; Timmons et al. 2016). Magnification factors change by up to 30% in our alternative models and are consistent with those measured from light-traces-mass models in Frye et al. (2018), suggesting that the remaining model degeneracies are not a major concern for the analysis of the Emerald.

We find $\mu = 3.8 \pm 0.5$ and $\mu = 6.1 \pm 0.9$ for the dust continuum emission in the more compact submm sources south (Co–S) and north (Co–N) from the arc, respectively. These values become about 20% lower when considering instead their gas emission.

## 4.4. Source-plane reconstruction

In Fig. 6 we show the reconstructed source plane morphologies of the gas and dust in the Emerald. We used the best-fitting lens model of the foreground cluster to infer the intrinsic size and position of the source seen as an extended arc. The source falls very near the critical line, so that the magnification varies by a factor of at least ten between the center and the two extreme ends (see inset in Fig. 5). For such a configuration, reconstructing correctly the source plane morphology of the dust and the molecular gas requires a specific procedure, such as those described in more detail in Johnson et al. (2017) and Fu et al. (2012). We followed a similar approach and described the source-plane profile of the gas and dust components with simple 2D elliptical Gaussian models, computed the associated brightness intensity in the image plane, convolved with the beam, and compared with the PdBI and SMA images.

The CO(4–3) emission is modeled with six parameters in the source plane, its position, semi-major and semi-minor axes, orientation and total flux density. We drew these parameters from Gaussian distributions centered on the values derived from the MCMC calculations, and used LENSTOOL to ray-trace the associated source profile to the image plane through the best-fitting potential derived in Sect. 4.2. We then convolved with the PdBI beam and computed the residual between the modeled and observed images. The $\chi^2$ was derived from the residual map, by considering all pixels that fall within the $3\sigma$ contours of the spectrally-integrated CO line emission. We explored the parameter space by iterating 1000 times, and adopted the source profile associated with the lowest $\chi^2$ in the image plane as the best fit. We performed a similar reconstruction of the HST/F110W emission. For the dust continuum, we assumed a circular Gaussian profile given that the arc is not resolved in the tangential direction by the SMA beam. The resulting maps and source plane models are shown in Fig. 6.

The best-fit position of the stellar continuum emission is consistent with that of the gas and dust centroids within the uncertainties. All components also have comparable projected spatial extents, resulting in a size of 2.7 kpc × 1.7 kpc in the source plane, elongated along position angle $PA = -73^\circ$ (east from north). Systematic uncertainties on these size estimates are up to about 25%, according to our alternative lens models. The stellar continuum is spatially resolved along both the major and minor axes, so that these sizes are both intrinsic; however, the major axis of the gas and dust measurement are dominated by the beam size. It is nonetheless encouraging to find similar sizes, and we will in the following discussion assume that the dust, gas,
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Source-plane reconstruction of the gas, dust and stellar emission from the Emerald, using the best-fitting lens model of the massive foreground cluster. Top row: (a) PdBI map summed over the frequency channels of the CO(4–3) line. The black lines show contours at –2σ (dashed lines), +2 and +5σ (solid lines). The beam of the PdBI at the redshifted frequency of CO(4–3) is plotted in the lower-left corner. (b) Best-fit source model ray-traced to the image plane and convolved with the PdBI beam (see details in the text). Orange and blue solid lines show the critical and caustic curves at z = 2.236, respectively, obtained from the best-fitting lens model. The box indicates the size and position of the enlarged region in the rightmost panel. (c) CO(4–3) residual and contours at –2, +2 and +5σ. (d) Best-fit model of the CO source, truncated at its FWHM (yellow ellipse). The bar in the lower-right corner illustrates the physical scales at z = 2.236. Center row: (e) SMA 850 µm image in the extended configuration and –2, +2 and +5σ contours (black lines). The SMA beam is shown in the lower-left corner. (f) Best-fit model convolved with the SMA beam in the image plane. (g) 850 µm residual and contours at –2, +2 and +5σ. (h) Best-fit azimuthally-averaged model of the 850 µm source, truncated at its FWHM (red circle). Bottom row: (i,j,k) Same as before, for the HST F110W imaging, within the field-of-view indicated in panel (e). (l) Best-fit model of the rest-frame optical stellar continuum in the source plane (blue ellipse), compared with those of CO(4–3) and 850 µm continuum. We also show the direct ray-tracing of the WFC3 F110W emission to the source plane (dark blue contours drawn at +2, +3 and +4σ), and the highly elongated PSF of the WFC3 image in the source plane (dark blue ellipse in the lower-left corner).</figcaption></figure>

and stellar components are distributed over similar regions in the source plane.

Since the HST imaging resolves the rest-frame optical con-
tinuum of the Emerald, we can reconstruct the stellar morphol-
ogy in the source plane by directly ray-tracing the observed
WFC3/F110W image pixel by pixel through the lensing mass
distribution with the *cleanlens* algorithm (Sharon et al. 2012);
*cleanlens* models the intrinsic morphology by allowing the
source-plane pixels to have arbitrary distortions and sizes
to match those of the un lensed source. The delensed stel-
lar continuum in the Emerald has a very elongated shape,
as expected given the reconstructed profile of the point-
spread function of the HST, which has an axis ratio of about
0.15.

Finally, the source-plane reconstruction allows us to infer
the relative positions of the arc and the two additional galax-
ies Co–N and Co–S. The reconstruction shows that the intrin-
sic positions of the arc and Co–N are consistent within 1σ,
for the best-fit and alternative models. This suggests that Co–
N is another image of the same galaxy as the one forming the
arc (i.e. the Emerald). For our best-fit model, Co–S falls at
a projected distance of about 18 kpc from the Emerald in the
source plane and is likely another, interacting or merging galaxy.
Alternative models suggest an offset of up to about 30 kpc.
Although these models give worse fits to the lensing configu-
ration, their results do provide a rough (and perhaps overly pes-
simistic) constraint on the systematic uncertainties that we might
expect.
<table><caption>**Table 4.** Intrinsic dust, gas, and star formation properties of the Emerald, spatially-integrated over the submm arc, as inferred from the single-dish observations and analysis presented in C15 and Cañamas et al. (2018).</caption><thead><tr><th>Quantity</th><th>Value</th><th>Unit</th><th>μ</th></tr></thead><tbody><tr><td>L<sub>FIR,gb</sub></td><td>(1.4 ± 0.3) × 10<sup>12</sup></td><td>L<sub>☉</sub></td><td>29.4 ± 5.9</td></tr><tr><td>L<sub>FIR,tpl</sub></td><td>(1.8 ± 0.4) × 10<sup>12</sup></td><td>L<sub>☉</sub></td><td>29.4 ± 5.9</td></tr><tr><td>SFR</td><td>176 ± 35</td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>29.4 ± 5.9</td></tr><tr><td>M<sub>d</sub></td><td>(7.4 ± 1.5) × 10<sup>7</sup></td><td>M<sub>☉</sub></td><td>29.4 ± 5.9</td></tr><tr><td>I<sub>CO(4-3)</sub></td><td>0.49 ± 0.10</td><td>Jy km s<sup>-1</sup></td><td>24.1 ± 4.8</td></tr><tr><td>L<sub>CO(4-3)</sub></td><td>(2.4 ± 0.5) × 10<sup>7</sup></td><td>L<sub>☉</sub></td><td>24.1 ± 4.8</td></tr><tr><td>L'<sub>CO(4-3)</sub></td><td>(7.6 ± 1.5) × 10<sup>9</sup></td><td>K km s<sup>-1</sup> pc<sup>2</sup></td><td>24.1 ± 4.8</td></tr><tr><td>M<sub>mol</sub></td><td>(1.1 ± 0.2) × 10<sup>10</sup></td><td>M<sub>☉</sub></td><td>24.1 ± 4.8</td></tr></tbody></table>

**Notes.** Each quantity is denensed using the relevant gravitational magnification factor, $\mu$, obtained in Sect. 4.3, corrected for the fraction of single-dish fluxes emitted by sources Co–N and Co–S, and divided by a factor of 2 to account for the two merging images in the arc. Errors include statistical uncertainties on $\mu$. We list the FIR luminosities integrated over the 8–1000 µm range, for a simple graybody function and for mid-infrared-to-millimeter templates (see C15).

# 5. Stellar population, molecular gas, dust, and star formation in the Emerald

## 5.1. Intrinsic integrated properties

After characterizing the foreground environment along the line of sight toward PLCK_G165.7+49.0 and the lensing configuration, we now discuss the properties of the Emerald, the main submm arc. The integrated properties are derived in C15, and Harrington et al. (2018) also present an integrated CO(1–0) spectrum of this source obtained with the Green Bank Telescope (GBT). However, neither analysis included a detailed lens model, so that the intrinsic properties of the Emerald could not be given. The interpretation of these results is also complicated by the presence of the two sources Co–S and Co–N within the 20′′–30′′ beams of Herschel/SPIRE and typical single-dish telescopes in the submillimeter and millimeter.

We used the lens modeling from Sect. 4.3 to derive intrinsic source properties from the observed integrated results. The total flux of the main arc and sources Co–N and Co–S measured on the SMA EXT map is (71.6 ± 0.6) mJy, corresponding to 80% of the flux measured with SCUBA-2. Assuming that all sources have similar dust temperatures, our 850 µm SMA map suggests that (81 ± 1)% of the total FIR luminosity detected with Herschel and the single-dish radio telescopes are emitted by the main arc, and (12 ± 1)% and (7 ± 1)% by sources Co–N and Co–S, respectively. Accounting for the magnification factors derived in Sect. 4.3, this results in intrinsic star-formation rates of SFR<sub>arc</sub> = (176 ± 35) M<sub>☉</sub> yr<sup>-1</sup> for the main arc, and SFR<sub>Co-N</sub> = (252 ± 39) M<sub>☉</sub> yr<sup>-1</sup> and SFR<sub>Co-S</sub> = (236 ± 28) M<sub>☉</sub> yr<sup>-1</sup> for the compact sources. We obtain a lower delensed star-formation rate for the Emerald compared to Co–N likely because, due to the parity inversion, none of the two merging images forming the submm arc are complete images of the intrinsic source. These star-formation rates are given for a Chabrier (2003) stellar initial mass function (IMF), and are therefore a factor 1.8 lower than when applying the popular prescription by Kennicutt (1989). Table 4 gives a summary of the intrinsic spatially-integrated dust, gas, and star formation properties of the Emerald inferred from the observational results of C15 and Cañamas et al. (2018), by

correcting for the strong lensing magnification and by a factor of two to account for the two images forming the arc.

The Emerald is a ULIRG with an intrinsic far-infrared luminosity, $L_{FIR} = 1.8 \times 10^{12} L_{☉}$, and intrinsic dust and molecular gas masses of $7.4 \times 10^7 M_{☉}$ and $1.1 \times 10^{10} M_{☉}$, respectively. The latter value assumes a “ULIRG” conversion factor between CO luminosity and molecular gas mass of $0.8 M_{☉}/(\text{K km s}^{-1} \text{ pc}^2)$, and a line luminosity ratio $R_{4,1} = L'_{CO(4-3)}/L'_{CO(1-0)} = 0.55$, as found empirically by, for example, Spilker et al. (2014) and Danielson et al. (2011). These masses are fully consistent with those measured for many submillimeter galaxies in the field (e.g., Tacconi et al. 2008), and the Emerald is thus a representative member of this class of galaxy.

A significantly lower ratio, $R_{4,1} = 0.25 \pm 0.10$, is found if we use the integrated CO(4–3) line luminosity we measured with the IRAM 30-m telescope, $\mu L'_{CO(4-3)} = (46.0 \pm 3.1) \times 10^{10} \text{ K km s}^{-1} \text{ pc}^2$ (Cañamas et al. 2018), and the CO(1–0) line luminosity from Harrington et al. (2018), who recently detected this line with the GBT for a few lens candidates drawn from the Planck all-sky survey, including five GEMS. Their CO(1–0) flux estimates are overall unusually high compared to other high-redshift galaxies, and highlight that this also results in exceptionally high gas masses, high gas-to-dust ratios, and low gas excitations, akin to local ULIRGs rather than other high-redshift galaxies (in spite of their high star-formation intensities).

In Fig. 7 we compare the Harrington et al. (2018) CO(1–0) line profile with the PdBI CO(4–3) line profile integrated over the arc. We find significant differences, with about 30% of the CO(1–0) flux emitted at velocities outside the range covered by CO(4–3); this is not expected if both lines trace the same clouds. The velocity range in the wings is well-matched with the velocity range of the sources Co–N and Co–S, which also fall within the FWHM ≈ 21′′ beam of the GBT at about 36 GHz, the redshifted frequency of CO(1–0) in the Emerald.

If we subtract the wings of the CO(1–0) line from the GBT spectrum, and correct for missing flux in our PdBI spectrum spatially-integrated over the arc, we find $R_{4,1} \approx 0.45$, comparable to the typical values found by Spilker et al. (2014) and Danielson et al. (2011). However, this higher value is also a lower limit to the intrinsic ratio in the Emerald, because of contamination with Co–N and Co–S. The source confusion within the beam of the GBT therefore prevents an in-detail comparison between the $J = 1-0$ and $J = 4-3$ CO line emissions and motivates our choice of adopting a fiducial value of $R_{4,1} = 0.55$.

## 5.2. Dust and gas morphology and gas kinematics

We show the SMA 850 µm morphology and the CO morphology and kinematics of the Emerald in Fig. 8. The long-wavelength emission is associated with the central regions of the western arc of 28.5" shown in Fig. 1. The gas and dust emissions are extended along an arc of about 5.3" length, which is not spatially resolved along the direction parallel to the critical line (and perpendicular to the magnification axis). The FWHM size along the minor axis is 0.84", compared to a beam size of 0.76" × 0.75". The Emerald is composed of four clumps in total, with the two southernmost clumps being strongly blended and difficult to separate. This is better seen at the somewhat higher resolution of the CO than the dust image. Within the beam size, the gas and dust morphologies are broadly consistent with each other, except for a small, roughly 0.1" positional offset for some clumps, which is much less than the beam size in either image.

About 60% of the dust emission originates from the clumps, the rest from a fainter component between. It is difficult to
<figure><img src="image_8.png" /></figure>

**Fig. 7.** PdBI CO(4–3) line profiles of the Emerald (the main arc of PLCK_G165.7+49.0, red line in the *left panel*) and the two nearby sources, Co–N and Co–S (red and blue lines in the *right panel*), compared to the GBT CO(1–0) line profile of Harrington et al. (2018, black line in both panels). The wings in the spatially-integrated CO(1–0) line profile are not detected toward the arc, but their velocity range is well-matched with the range covered by the CO(4–3) line emission from sources Co–N and Co–S, which are blended with the arc in the GBT beam.

<figure><img src="image_9.png" /></figure>

**Fig. 8.** Top panel: dust and gas morphology and gas kinematics of the Emerald, the main submm arc of PLCK_G165.7+49.0, from single-component Gaussian line fits. Left to right panel: maps of the SMA 850 µm dust and PdBI CO(4–3) morphologies, relative velocity, and Gaussian line FWHM. The arc is composed of two multiple images of the same source, #1.1 on the southern portion and #1.2 on the northern one, respectively (see also Fig. 5). We identify two star-forming clumps per image and focus most of the analysis on the two northern ones, called #1.2 S and #1.2 N. The SMA and PdBI beam FWHMs are 0.90" × 0.75" and 0.76" × 0.75", respectively. Velocity offsets and line widths are given in km s⁻¹. Contours on the two leftmost panels show the SMA 850 µm dust continuum, starting at +6σ and increasing in steps of +2σ, and those on other panels show the stellar continuum in the F110W band from HST/WFC3. Bottom, left to right panel: same maps for the secondary Gaussian component detected at ≥ 3σ over a region of about 2.5" × 1.5" near the center of the arc.

constrain whether this component is also clumpy or more dif-
fuse because the emission is faint and blended with the clumps,
and the intraclump regions are spatially not well resolved.

The velocity structure of the arc is seen in our PdBI CO(4–3)
maps (Fig. 8), which we constructed in the way described in
Sect. 2.4. Relative velocities increase from southeast to north-
west from –63±15 km s⁻¹ to +336±15 km s⁻¹, and then decrease
again to –70 ± 15 km s⁻¹. Gaussian line widths, σ, are rather
moderate and range from 93 ± 25 km s⁻¹ to 190 ± 25 km s⁻¹;
they are shown in the right panel of Fig. 8. It is interesting that
the line widths in the center of the arc, near the critical line of
the best-fitting lensing model, are more narrow than further out.

If several, partially overlapping images were present we would
expect the line widths to increase because of blending, and if the
velocity offsets across the major axis were indicating a merger
of two rotating gaseous disks, we would also expect a higher tur-
bulence and broader line profiles at this position. We therefore
consider the narrow lines as additional, supporting evidence that
the arc comprises two lensed images of the same region, with a
parity inversion on both sides of the critical line.

We can estimate a dynamical mass from the velocity gradi-
ent, v, of 380 km s⁻¹ and the intrinsic FWHM size of the Emer-
ald of 2.7 kpc (Sect. 4.4). If we assume, as is usually done, that
the gradient encompasses both sides of a disk, then we find a
dynamical mass of $M_{dyn} = (v/\sin i)^2 R/G = 1.1 \times 10^{10} M_{\odot}$, where R is the disk radius, and G the gravitational constant. This estimate ignores possible inclinations, i, or a mismatch between lensing and the rotational major axis, which may lead to incomplete sampling of the rotation curve. The resulting mass is much lower than the sum of the gas and stellar masses of $5.1 \times 10^{10} M_{\odot}$ (Sects. 5.1 and 5.3), implying that we are either observing a disk nearly seen face-on at $\le 20^\circ$ inclination angle, a disk magnified approximately along the kinematic minor axis, or a disk where only one side of the rotation curve is being magnified. In the latter case, we obtain $M_{dyn} = 9.1 \times 10^{10} M_{\odot}$, by setting velocity $v = 380 \text{ km s}^{-1}$ and $R = 2.7 \text{ kpc}$, about 40% larger than the baryonic mass. The discrepancy would be alleviated if we had used a Salpeter (1955) initial mass function instead of the typically adopted Chabrier (2003) IMF, and as favored by the lens of another GEMS, the Ruby (Cañamas et al. 2017a, C17a hereafter). Neither mass estimate, however, places the Emerald outside the typical mass range of massive, dusty, intensely star-forming galaxies at $z \sim 2$ (e.g., Casey et al. 2014).

### 5.2.1. Systematic effects of lensing on the kinematic measurements

Strong gravitational lensing near the critical line affects the brightness distribution of the gas and dust emission, and obviously plays a large role in deriving integrated continuum and line fluxes and related quantities. Integrated line profiles and velocity offsets, however, are a convolution of the intrinsic kinematics and emission-line surface-brightness distribution, and might therefore also be affected by the details of the lens reconstruction, and could potentially add important systematic uncertainties to estimates of intrinsic source properties. To investigate the impact of gravitational lensing on the emission-line parameters, we have extracted the spectra within three apertures: including the region showing a blue emission-line component in Fig. 8; a $0.6''$ wide annulus around that region; and the remainder of the Emerald, before and after our pixel-by-pixel correction for the gravitational magnification. We then compared the results of our Gaussian line fitting in each case.

The difference of the kinematic properties of each line component in the image and source plane are in fact very small. When extracting the spectra without the lensing correction, we find relative velocities of $318 \text{ km s}^{-1}$, $84 \text{ km s}^{-1}$, $161 \text{ km s}^{-1}$, and $-7 \text{ km s}^{-1}$ for the narrow and broad components associated with the clump, the spectra from the annulus and the remainder of the Emerald, respectively. With the exception of the annulus, these are within $10 \text{ km s}^{-1}$ of the values measured from the spectra after correcting for the local magnification factors. The FWHM line widths are $218 \text{ km s}^{-1}$, $418 \text{ km s}^{-1}$, $443 \text{ km s}^{-1}$, and $360 \text{ km s}^{-1}$, even closer to those measured after the correction. This suggests that the impact of gravitational lensing on the derived emission-line kinematics is negligible, at least for our analysis of the Emerald.

### 5.2.2. Systematics on the gas kinematics from clump identification

As a final test of our analysis of the resolved dust and CO line emission in the Emerald, we also used CLUMPFIND (Williams et al. 1994) to quantify the surface-brightness distribution and kinematic substructure in the arc of the Emerald in a more reproducible way than the one used when identifying these structures by eye in a rather heuristic way. CLUMPFIND is a publicly available IDL-routine that identifies contiguous

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Spatially-integrated CO(4–3) spectrum of the region shown in Fig. 8 where we detect a secondary emission-line component. The red line shows the double-Gaussian fit to the overall spectrum and blue lines show the individual blueshifted and systemic components.</figcaption></figure>

structures within an imaging spectroscopy data cube, starting with the brightest peak in the cube, and then lowering the flux threshold with a step size that can be selected by the user. We started at $11.8 \text{ mJy beam}^{-1}$ and decreased the flux in steps of $1.8 \text{ mJy beam}^{-1}$ ($3 \times$ the rms of our data cube), until we reached a threshold of $3 \text{ mJy beam}^{-1}$ (corresponding to $5\sigma$).

CLUMPFIND finds the same two extended components, which we also identified by eye, and in addition, the secondary, blueshifted component in the center of the source, which it identifies as a single structure. In addition to the structure within the arc, CLUMPFIND also identifies the two separate sources, Co–N and Co–S, each associated with a single clump. This confirms our by-eye analysis. All individual clumps are detected at $S/N = 12$ and 20 per beam in the CO line emission and dust continuum, respectively, well above the signal-to-noise ratios at which Hodge et al. (2016) found spurious clumpiness in high-resolution dust imaging of high-redshift galaxies with ALMA.

### 5.3. Resolved stellar component

The near-infrared HST/WFC3 emission from the western arc extends over a total length of $28.5''$, and consists of at least 12 multiple images. We focus again on the Emerald, the central region of the arc, also detected in CO and $850 \mu\text{m}$ dust emission. Its stellar continuum morphology is best seen in the bottom panel of Fig. 6. In the rest-frame optical this component consists mainly of diffuse emission, with a bright clump near the center of the arc, where the gravitational magnification is greatest. The intrinsic FWHM size of this region is $2.7 \text{ kpc} \times 1.7 \text{ kpc}$ along the major and minor axis, respectively, using our lensing model presented in Sect. 4.3. The size of the arc in the WFC3 images is comparable to that seen in the dust and gas (Sect. 5.2), but the morphologies are not strictly the same, since the stellar continuum reaches the highest surface-brightnesses where the gas and dust emissions are faintest. In turn, the brightest clumps seen at long wavelengths correspond to rather faint regions in the stellar continuum. It is entirely possible that this is mainly a sign of variations in dust and cloud cover.
<table><caption>**Table 5.** Intrinsic photometry of the Emerald obtained from our new optical/near-infrared imaging with the CFHT, HST and *Spitzer*, and by correcting the source-integrated (sub-)millimeter flux densities presented in C15 for the contribution of Co-N and Co-S.</caption><thead><tr><th>Band</th><th>Flux density</th><th>Unit</th></tr></thead><tbody><tr><td>CFHT/MEGACAM <i>r</i></td><td>&lt;8.5</td><td>nJy</td></tr><tr><td>CFHT/MEGACAM <i>z</i></td><td>&lt;29.3</td><td>nJy</td></tr><tr><td>HST/WFC3 F110W</td><td>47.8 ± 3.8</td><td>nJy</td></tr><tr><td>CFHT/WIRCAM <i>J</i></td><td>64.5 ± 5.3</td><td>nJy</td></tr><tr><td>HST/WFC3 F160W</td><td>110 ± 9</td><td>nJy</td></tr><tr><td>CFHT/WIRCAM <i>K<sub>s</sub></i></td><td>345 ± 13</td><td>nJy</td></tr><tr><td><i>Spitzer</i>/IRAC 3.6 µm</td><td>1.19 ± 0.01</td><td>µJy</td></tr><tr><td><i>Spitzer</i>/IRAC 4.5 µm</td><td>1.86 ± 0.02</td><td>µJy</td></tr><tr><td>WISE W3</td><td>9.5 ± 2.7</td><td>µJy</td></tr><tr><td>WISE W4</td><td>204 ± 19</td><td>µJy</td></tr><tr><td><i>Herschel</i>/SPIRE 250 µm</td><td>13.4 ± 0.3</td><td>mJy</td></tr><tr><td><i>Herschel</i>/SPIRE 350 µm</td><td>11.6 ± 0.2</td><td>mJy</td></tr><tr><td><i>Herschel</i>/SPIRE 500 µm</td><td>7.3 ± 0.2</td><td>mJy</td></tr><tr><td>JCMT/SCUBA-2 850 µm</td><td>1.4 ± 0.4</td><td>mJy</td></tr><tr><td>IRAM/GISMO 2 mm</td><td>123 ± 39</td><td>µJy</td></tr></tbody></table>

**Notes.** We demagnified the optical and near-infrared stellar continuum fluxes by the gravitational magnification factor $\mu_{\text{stars}} = 34.1 \pm 6.8$ and the dust continuum fluxes by $\mu_{\text{dust}} = 29.4 \pm 5.9$, both derived from our best-fitting lensing model. Since the stellar and dust components likely have similar contributions to the observed WISE emission, we used the average value of $\mu$ to correct the fluxes in the W3 and W4 bands (as done in Timmons et al. 2016). We also give the $3\sigma$ upper limits of non-detection in MEGACAM *r*- and *z*-bands.

We probed the intrinsic stellar continuum properties of the Emerald using the counterparts of the submm arc in our CFHT, WFC3 and IRAC imaging. To correct for some faint underlying continuum emission from member galaxies of the foreground cluster, we modeled the three lens galaxies in the western group with Sérsic profiles using GALFIT (Peng et al. 2010). The debленched fluxes of the high-redshift source are measured in the lens-subtracted residual images, with a prior in the $K_s$-band. These fluxes are then divided by a factor of two, to account for the two merging images producing the arc, and corrected for the magnification factor $\mu_{\text{stars}}$ reported above. We also corrected the SPIRE, SCUBA-2 and IRAM 30-m/GISMO single-dish flux densities for gravitational lensing using $\mu_{\text{dust}}$, divided by a factor of two to account for the two blended counter images, and corrected for the roughly 20% of the total flux from sources Co-N and Co-S. The resulting optical to submillimeter photometry is presented in Table 5.

Using the simple stellar population models (SSPs) from Bruzual & Charlot (2003), with solar metallicity, a Chabrier (2003) stellar IMF, and the Calzetti et al. (1994) extinction law, we obtain a best-fitting model for very young stellar populations of about 25 Myr and $A_V \approx 2.7$ mag, with a goodness-of-fit of $\chi^2 = 2.7$. We also obtain good fitting results when using exponentially declining star-formation histories instead, with ages below 50 Myr.

We then used MAGPHYS (da Cunha et al. 2008, 2015) to better constrain the amount of energy reprocessed by dust through our long-wavelength photometry, and to derive a robust stellar mass estimate for the Emerald. MAGPHYS has proven to give consistent results with other codes using different template libraries (e.g., Nayyeri et al. 2017). A possible AGN contribution to the dust heating is neglected, but this is not a concern here, since

we have already shown that AGN heating does not dominate the far-infrared spectral energy distribution (C15).

We fit the delensed SED of the arc with the values listed in Table 5, assuming a redshift $z = 2.236$, and obtain the best-fit SED shown in Fig. 10, with a goodness-of-fit of $\chi^2 = 3$. The resulting intrinsic stellar mass of the source that produces the submm arc is $M_* = (4.1 \pm 0.4) \times 10^{10} M_\odot$, for a Chabrier IMF, with $A_V = 3.9 \pm 0.3$ mag, and an intrinsic star-formation rate (SFR) of $(74 \pm 3) M_\odot \text{yr}^{-1}$ (corrected for a gravitational magnification factor of $\mu = 34.1$). The corresponding stellar mass surface density is $(11 \pm 5) \times 10^9 M_\odot \text{ kpc}^{-2}$ for an intrinsic size of $2.7 \text{ kpc} \times 1.7 \text{ kpc}$ along the (delensed) major and minor axis, respectively. The dust properties inferred with MAGPHYS are consistent with those listed in Table 2 of C15, after correcting for the gravitational magnification where necessary.

The blue clump very near to the critical line, which may either be part of the Emerald or be an intervening source (see Sect. 4.2 and Fig. 1), contributes only marginally to the integrated flux of the Emerald. We do not include it in the SED fitting presented here, but we did run an alternative model which includes the clump, and found that the fit results remain within the $1\sigma$ uncertainties. We also checked that the small, spatial offset between the rest-frame optical and submillimeter and millimeter position of the Emerald (Sect. 4.4) have no significant impact on the SED modeling with MAGPHYS, by reproducing the fit without the WFC3 and CFHT fluxes. The results are consistent with those of the best-fit model within the $1\sigma$ uncertainties, likely because the stellar mass is mainly constrained by IRAC photometry in the rest-frame NIR.

The total star-formation rate found from the stellar continuum of $SFR_{\text{opt}} = 74 M_\odot \text{yr}^{-1}$ is much lower than that obtained from the far-infrared spectral energy distribution, $SFR_{\text{FIR}} = 176 M_\odot \text{yr}^{-1}$, after correcting for a magnification factor of $\mu = 29.4$ (see Sect. 4.3). This may indicate that most of the star formation in the Emerald is hidden behind high dust and gas column densities. This is also consistent with the morphologies of the dust, gas, and stellar component. High-surface brightness dust emission extends over much larger radii than the bright stellar continuum (although faint continuum emission is also probed over larger scales). Given that the stellar component falls very near the caustic line, and that the overall distribution of gas and dust in intense star-forming regions is likely very clumpy (e.g., Genzel et al. 2012; Swinbank et al. 2015; Iono et al. 2016), it seems plausible that we are fortuitously seeing along a relatively low-extinction sight-line into the starburst.

## 5.4. Additional sources

As mentioned previously, the Emerald is surrounded by the two neighboring sources Co-S and Co-N at projected distances of 5" and 13" toward southeast and northwest, respectively (see Fig. 11). These components were first introduced in Sect. 4.4, where we also listed their luminosity-weighted gravitational magnification factors, $\mu_{\text{Co-N}} = 6.1 \pm 0.9$ and $\mu_{\text{Co-S}} = 3.8 \pm 0.5$. The two sources are almost compact, with projected major axis length of about 1.1" and 0.9", respectively, in the CO line image obtained with the IRAM interferometer, for a beam size of $0.76'' \times 0.75''$. Our lensing model suggests that source Co-N is an additional counter-image of the same galaxy that also gives rise to the Emerald (seen with a much lower magnification factor), whereas Co-S is a distinct component in the source plane. Our best-fitting lensing model favors a scenario where Co-S is a companion galaxy about 18 kpc away from the Emerald, although the relative projected distance is uncertain given
<figure><img src="image_11.png" /></figure>

**Fig. 10.** Top panel: intrinsic spectral energy distribution of the Emerald from the optical to the millimeter. Red circles show the photometry of the arc presented in Table 5, corrected for the gravitational magnification factors of the stellar and dust components, as described in the text. Flux uncertainties are smaller than the symbols and downward arrows show the $3\sigma$ upper limits on the MEGACAM fluxes. The best-fit SED to the full wavelength range obtained with MAGPHYS is plotted as a solid blue curve, and the best-fit stellar continuum without dust attenuation is shown in light blue. The dash-dotted orange line indicates the best-fit template of the local starburst galaxy Arp 220, shifted to $z = 2.236$ and normalized to the flux density of the Emerald in the 350 µm band of SPIRE. Bottom panel: residuals of the best-fitting MAGPHYS model.

<figure><img src="image_12.png" /></figure>

**Fig. 11.** Dust and gas morphology of the Emerald and the two neighboring compact sources Co-S and Co-N toward southeast and northwest, respectively. The two insets have sizes of 3" × 3". The color scale displays the PdBI CO(4–3) emission and the contours show the SMA 850 µm dust continuum.

the remaining degeneracies in the model. Major and minor axis sizes in the source plane are 5.9 kpc × 3.1 kpc for Co-S and 6.0 kpc × 2.7 kpc for Co-N.

We show the integrated spectra of both galaxies in Fig. 7. Their intrinsic properties are listed in Table 6, and were derived with the same methods and assumptions as previously used for the Emerald. Values found for Co-N are consistent with those for the Emerald itself, further suggesting that we may be seeing another image of the same galaxy with a smaller magnification factor.

## 5.5. Star-formation law in the Emerald

Star-formation surface densities (star-formation intensities) and gas-mass surface densities are closely related through the Schmidt-Kennicutt relationship (Schmidt 1959; Kennicutt 1989), and allow us to characterize the star-formation processes in galaxies over more than six orders of magnitude in star-formation rate and gas density. The position of a galaxy within this diagram is an indicator of the efficiency with which gas is being turned into stars. While the main physical driver (or drivers) of this relationship are still a matter of active debate, different scenarios of how star formation is being regulated in these galaxies make different predictions for where a galaxy should fall within this model (which might not, however, necessarily be unique).

### 5.5.1. Star-formation intensities

In order to investigate the star formation properties of the Emerald with the aid of the Schmidt-Kennicutt diagram, we need to translate our spatially resolved dust and CO surface brightnesses into star-formation intensities and gas-mass surface densities. This requires several assumptions. First, we obtain star-formation intensities by assuming that all the dust within the Emerald has a single temperature, which corresponds to the value $T_d = 42.5$ K measured in C15 from the spatially-integrated FIR-to-millimeter spectral energy distribution. This allows us to convert the local surface brightnesses of the 850 µm continuum into local star-formation intensities, by setting SFR = $2.5 \times 10^{-44} L_{\text{FIR}}$, where SFR is given in $M_{\odot} \text{yr}^{-1}$. The far-infrared luminosity, $L_{\text{FIR}}$, is in erg s$^{-1}$ and $L_{\text{FIR}}$ is integrated over 8–1000 µm (Kennicutt 1989). We remind the reader that this estimate is appropriate for a Chabrier stellar IMF and is a factor 1.8 (0.26 dex) lower than that originally adopted by Kennicutt (1989). The Chabrier (2003) parametrization is currently the most commonly adopted in high-redshift studies, and we adopted this IMF to be consistent with these studies. We thus correct all values in this section and those adopted from the literature to take this into account, while understanding that this might introduce a bias, since massive galaxies may be better characterized by a Salpeter (1955) IMF (e.g., Cañamas et al. 2017b).

This approach results in central star-formation rate surface densities between 38 and $45 M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}$ in the four clumps seen in the dust continuum, in images #1.1 and #1.2 of the Emerald (Table 7). Since lensing conserves surface brightness, we did not need to correct these quantities for the magnification factors.

### 5.5.2. Gas-mass surface densities

The second quantity relevant for our analysis is the molecular gas mass surface density. We used the measured CO(4–3) fluxes, and followed Solomon et al. (1997) to convert them into mass surface densities of molecular hydrogen. To do so, we also needed to adopt a CO-to-H$_2$ conversion factor, a quantity whose value is still controversial. We adopted the “standard” conversion factor appropriate for ULIRGs (Downes & Solomon 1998), $\alpha_{\text{CO}} = 0.8 M_{\odot}/(\text{K km s}^{-1} \text{pc}^2)$, which is commonly used for intense star-burst galaxies at high redshift and also consistent with our analysis of the gas-to-dust ratios in the GEMS (see C15).

Another complication is that we used a mid-$J$ transition of CO(4–3), whereas the Schmidt-Kennicutt law was calibrated on CO(1–0). However, emission from the ground rotational state of CO can be significantly contaminated by diffuse molecular gas
<table><caption>**Table 6.** Properties of compact components Co-N and Co-S.</caption><thead><tr><th>Source</th><th>RA<br>(J2000)</th><th>Dec<br>(J2000)</th><th>D<sub>maj</sub> × D<sub>min</sub><br>(arcsec<sup>2</sup>)</th><th>PA<br>(deg)</th><th>S<sub>850</sub><br>(mJy)</th><th>L<sub>FIR</sub><br>(10<sup>12</sup> L<sub>☉</sub>)</th><th>SFR<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>I<sub>CO(4-3)</sub><br>(Jy km s<sup>-1</sup>)</th><th>L'<sub>CO(4-3)</sub><br>(10<sup>10</sup> K km s<sup>-1</sup> pc<sup>2</sup>)</th><th>M<sub>mol</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>Co-N</td><td>11:27:13.85</td><td>+42:28:16.56</td><td>1.1 × 0.9</td><td>60 ± 2</td><td>1.8 ± 0.3</td><td>2.6 ± 0.4</td><td>252 ± 39</td><td>0.72 ± 0.16</td><td>1.1 ± 0.2</td><td>1.6 ± 0.3</td></tr><tr><td>Co-S</td><td>11:27:13.89</td><td>+42:28:35.30</td><td>0.9 × 0.9</td><td>44 ± 2</td><td>1.7 ± 0.2</td><td>2.5 ± 0.3</td><td>236 ± 28</td><td>0.67 ± 0.13</td><td>1.1 ± 0.2</td><td>1.5 ± 0.3</td></tr></tbody></table>

**Notes.** Columns are: source name; right ascension and declination; apparent major axis and minor axis FWHM, and position angle measured on the PdBI CO(4–3) map; intrinsic flux density at 850 µm; delensed FIR luminosity integrated in the range 8–1000 µm and star-formation rate given for a Chabrier IMF; intrinsic velocity-integrated flux and luminosity of CO(4–3) line; molecular gas mass assuming $R_{4,1} = 0.55$ and $\alpha_{CO} = 0.8 M_{\odot}/(\text{K km s}^{-1} \text{pc}^2)$. Intrinsic quantities are corrected for the local magnification factors, $\mu_{CO-N} = 6.1 \pm 0.9$ and $\mu_{CO-S} = 3.8 \pm 0.5$, inferred in Sect. 4.3, and errors include statistical uncertainties on $\mu$.

<table><caption>**Table 7.** Properties of individual star-forming clumps in the Emerald.</caption><thead><tr><th>Source</th><th>RA<br>(J2000)</th><th>Dec<br>(J2000)</th><th>μ S<sub>850</sub><br>(mJy)</th><th>SFR<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>μ</th><th>I<sub>CO(4-3)</sub><br>(Jy km s<sup>-1</sup>)</th><th>M<sub>gas</sub><br>(10<sup>9</sup> M<sub>☉</sub>)</th><th>v<br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>α<sub>vir</sub></th></tr></thead><tbody><tr><td>1.2 N</td><td>11:27:14.67</td><td>+42:28:23.9</td><td>7.8 ± 0.6</td><td>100 ± 8</td><td>14.2</td><td>(2.3 ± 0.2) × 10<sup>-1</sup></td><td>5.0 ± 0.5</td><td>44 ± 9</td><td>374 ± 38</td><td>2.9 ± 0.6</td></tr><tr><td>1.2 S (m)</td><td>11:27:14.70</td><td>+42:28:23.0</td><td>8.1 ± 0.6</td><td>41 ± 3</td><td>35.0</td><td>(5.0 ± 0.5) × 10<sup>-2</sup></td><td>1.0 ± 0.1</td><td>334 ± 28</td><td>199 ± 17</td><td>1.4 ± 0.3</td></tr><tr><td>1.2 S (b)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>(5.0 ± 0.6) × 10<sup>-2</sup></td><td>1.2 ± 0.1</td><td>135 ± 17</td><td>333 ± 43</td><td>5.3 ± 1.1</td></tr></tbody></table>

**Notes.** We list observed values for continuum flux densities, $S_{850}$, and intrinsic values for line fluxes, star-formation rates and gas masses. The local magnification factors are listed in column $\mu$. Error bars give the rms for the continuum and uncertainties in the fit for emission-line data. For the clumps where we identify a secondary line component, we also list the velocity offset and FWHM line width at the position of that clump. The virial parameters, $\alpha_{vir}$, quantify the relative contribution from gravitational binding and turbulent energy.

outside of star-forming clouds (Ivison et al. 2011). A possible caveat of directly using CO(1–0) measurements is therefore that the total line emission might overestimate the gas-mass surface densities within the star-forming regions themselves. As already discussed in Sect. 5.1, we used a ratio of CO(4–3) to CO(1–0) of 0.55, as typically found for strongly gravitationally lensed, dusty, far-infrared-selected high-redshift galaxies (Spilker et al. 2014; Danielson et al. 2011). The difference compared with the lower estimate based on the CO(1–0) line detection of Harrington et al. (2018) is below 0.1 dex, and indicated by the small horizontal arrow in Fig. 12. Resulting gas-mass surface densities are 2400–2900 $M_{\odot} \text{ pc}^{-2}$.

### 5.5.3. Star-formation law

In Fig. 12 we show that the Emerald falls just below the relation-ship for high-redshift starburst galaxies. Filled red and empty blue stars indicate the central surface brightnesses of the four bright clumps shown in Fig. 8, and the residual intraclump emis-sion, respectively (see Sect. 6.2). Star-formation and gas mass surface densities in the Emerald are comparable to submillime-ter galaxies in the field (Bothwell et al. 2010), and are signifi-cantly lower than in the Ruby (C17a), SDP.81 (Hatsukade et al. 2015), and the Eyelash (Swinbank et al. 2011). This includes the bright star-forming clumps, but also the extended diffuse emis-sion, which has molecular gas mass surface densities of only a few $100 M_{\odot} \text{ pc}^{-2}$.

It is important to test if the lower surface densities are a con-sequence of our larger beam sizes. The Ruby has clump sizes of 0.1"–0.3", and if we run a toy model with clumps of such sizes, we see indeed a decrease in gas mass and star-formation surface density of a few 0.1 dex compared to the 0.1" beam with which the Ruby was observed (C17a), although the position rel-ative to the ridge lines of starburst and main-sequence galaxies does not change. However, clumps as compact as those in the Ruby should have clearly been seen in our SMA VEXT imag-ing. We ran a suite of simple toy models to investigate how many

sources with 0.2" FWHM size, comparable to the clump sizes in the Ruby, would be necessary to explain the 6-mJy minimal brightness we observed with the EXT configuration of the SMA, without violating the $3\sigma$ upper surface-brightness limit imposed by the non-detection of such clumps in the VEXT configuration down to rms = 1.75 mJy. This was only possible by populating the 0.8" beam with a near-uniform distribution of more compact sources of about the same brightness. This strongly disfavors the presence of bright, compact clumps within the Emerald with sizes much less than the 0.6" minor axis size seen in stellar light. Beam-smearing effects in Fig. 12 must therefore be small, which implies star-formation intensities of a few tens of $M_{\odot} \text{ yr}^{-1} \text{ kpc}^{-2}$ and gas densities of a few thousand $M_{\odot} \text{ pc}^{-2}$.

## 6. Feedback and disk fragmentation

### 6.1. Global disk (Toomre) stability

Most galaxies observed at redshifts $z \ge 2$ with resolved data sets and discernable velocity gradients fall near the critical value for rotationally supported gas, meaning that they have Toomre parameters $Q \sim 1$. For thin, uniform disk models, this implies that their gas reservoirs are globally stable against gravitational collapse on kpc scales. Several authors have recently pointed out that gas-rich, clumpy, star-forming galaxies that are already frag-mented, may represent more complex environments where gas at $Q > 1$ can still form clumps. This pushes the critical stability parameter below which the gas becomes unstable to gravitational collapse to somewhat higher values, of order $Q \sim 2–3$ (e.g., Inoue et al. 2016).

The Emerald is no exception in this regard. Using $v_c = 380 \text{ km s}^{-1}$, $M_{\text{tot}} = 9.1 \times 10^{10} M_{\odot}$, and $M_{\text{gas}} = 1.1 \times 10^{10} M_{\odot}$, and setting $Q = \sigma_0/v_c \times a \times M_{\text{tot}}/M_{\text{gas}}$ (Genzel et al. 2011), we find $Q = 1.3$. For this estimate we adopted the lowest value of $\sigma_0$ that we could find within the Emerald, $\sigma_0 = 95 \text{ km s}^{-1}$, which is probably the best approximation of the gas kinematics outside the brightest star-forming knots, as suggested by Inoue et al. (2016).
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Spatially-resolved Schmidt-Kennicutt diagram in the Emerald, for the four star-forming clumps identified with CLUMPFIND (red stars) and the intraclump regions (open blue stars), the spatially resolved pixel-by-pixel analysis of the Ruby (light blue circles, Cañameras et al. 2017a) and its luminosity-weighted properties (black stars), and for other galaxy samples taken from the literature. These samples include spatially-resolved studies of the Eyelash (green triangles, Swinbank et al. 2011) and SDP.81 (yellow triangles, Hatsukade et al. 2015), submillimeter galaxies at $z \sim 2$ (green upside-down triangles, Bothwell et al. 2010), normal star-forming galaxies at $z = 1-2.3$ (blue squares, Tacconi et al. 2010), and local starbursts (black circles, Kennicutt 1998). The dashed lines labeled “SB” and “MS” show the ridge lines of “starburst” and “main-sequence” galaxies of Daddi et al. (2010), respectively (see also Genzel et al. 2010). Typical error bars are shown in the lower-right corner. The small horizontal and vertical arrows indicate the expected offsets when using $R_{4,1} = 0.45$ instead of 0.55, and the Salpeter stellar IMF instead of the Chabrier IMF, respectively (see text for details).</figcaption></figure>

Here *a* is a morphological parameter which is of order unity (see Genzel et al. 2011, for details). The value of *Q* = 1.3 is well within the range typically found for high-redshift galaxies. Even the combination of estimates of $v_c$, $M_{\text{gas}}$, $M_{\text{tot}}$, and $\sigma$ that would yield the highest possible value of *Q*, would give *Q* = 2.6, still in the range of marginally Toomre-stable disks expected from simulations of fragmenting galaxies.

## 6.2. Clump properties and stability

Many authors have already discussed the importance of clump survival for massive high-redshift galaxies (e.g., Cowie et al. 1996; Elmegreen 2007; Genzel et al. 2008; Bournaud & Elmegreen 2009; Förster Schreiber et al. 2009; Tacconi et al. 2013; Mayer et al. 2016; Dessauges-Zavadsky et al. 2017). Blue stellar continuum morphologies in about 50–60% of actively star-forming, UV/optically-selected galaxies at high redshift show that considerable fractions of the stellar mass in these galaxies are in giant clumps of 1 kpc or less in size and a few times $10^7$ up to $10^9 M_{\odot}$ in mass (e.g., Elmegreen & Elmegreen 2005; Förster Schreiber et al. 2011; Livermore et al. 2012). Moreover disk fragmentation and clump formation can also be important for dusty, far-infrared and submillimeter selected galaxies like the GEMS (Swinbank et al. 2010). If clumps are long-lived, they may sink toward the galaxy center within a few orbital times (few hundred Myr) and merge to form a bulge. This scenario has been put forward in particular in the context of early simulations of gas-rich, fragmented disk galaxies (e.g., Ceverino et al. 2010; Bournaud et al. 2014; Mandelker et al. 2014), but has recently been challenged by detailed hydrodynamic models (Tamburello et al. 2015; Mayer et al. 2016; Oklopčić et al. 2017) and observations of

blue, gravitationally lensed galaxies (Dessauges-Zavadsky et al. 2017; Tamburello et al. 2015). These studies favor a scenario where clumps are more marginally gravitationally bound, and may dissolve within a few tens of Myr if feedback becomes too strong, either through mass loss in the form of winds, or by producing turbulent velocities near or above the virial velocity, or a mixture of both (for the latter, see in particular Hayward & Hopkins 2017). Whether clumps survive seems to depend critically on the detail of how feedback is implemented in these simulations (Mayer et al. 2016).

In the Emerald, the degree of fragmentation between the stellar component on the one hand, and gas and dust on the other, is very different. Stars, dust, and gas in the Emerald are found over large and similar ranges in a disk of at least 2.7 kpc × 1.7 kpc in size, and with rather moderate gas and star-formation rate surface densities (Sect. 5.5). The stellar light distribution is very smooth. At most about 25% of the F110W emission from the arc is in high surface-brightness features, and less than a few percent in the F160W image. Overall, the surface brightness varies by less than factors of between two and three in the F160W image, and varies only on small scales comparable to the PSF size, which would be smoothed out if seen at a resolution comparable to those of our long-wavelength data.

In contrast, most of the gas and dust emission is concentrated in two bright clumps with radii $\le 100-250$ pc, each imaged to the north and south of the critical line (see Sect. 5.2 and Fig. 8). These are upper limits set by the beam size and local gravitational magnification factors listed in Table 7. The two clumps together comprise about 60% of the total dust emission in the Emerald. A fainter, diffuse component is also seen after modeling the clumps with a Gaussian beam, and removing their contribution to the emission. Moreover, the average stellar mass surface
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Spectra of individual clumps #1.2 S (**left panel**) and #1.2 N (**right panel**), extracted from our PdBI CO(4–3) data cube using 0.6" × 0.6" apertures, and corrected for gravitational magnification. Clump #1.2 N is well fitted with a single Gaussian function (red curve), while #1.2 S exhibits a secondary component, blueshifted with respect to the systemic component (blue curves), which we interpret as the signature of a stellar wind.</figcaption></figure>

density of $1.1 \times 10^{10} M_{\odot} \text{ kpc}^{-2}$ is a factor of a few greater than those of the molecular gas inferred in Sect. 5.5 (see Fig. 12). This results in a global gas-to-baryonic mass ratio in the Emerald of about 20%.

We show the spectra of the two clumps in Fig. 13, extracted over $3 \times 3$ pixel apertures in our CO line data (0.6" × 0.6"), and scaled to the flux density per beam of the central pixel, as is most appropriate for unresolved sources. We have focused on the two northern images of these clumps, #1.2 N and #1.2 S, which are easier to deblend. Their properties are derived following the same approach and assumptions used for the Schmidt-Kennicutt analysis in Sect. 5.5. Both clumps are massive and gas-rich, with $M_{\text{gas}} = 5.0 \times 10^9 M_{\odot}$ and $1.0 \times 10^9 M_{\odot}$, respectively, in the upper range of clump masses in other high-redshift galaxies with giant clumps (Genzel et al. 2011). Table 7 lists their individual properties, including star-formation rates of 40–100 $M_{\odot} \text{ yr}^{-1}$ (for a Chabrier IMF), and FWHM line widths of CO(4–3) of 200–375 km s$^{-1}$. Their sizes are smaller than those seen in field galaxies with very massive clumps (Förster Schreiber et al. 2011, but we note that these masses may be overestimated, Tamburello et al. 2015; Dessauges-Zavadsky et al. 2017). Gas mass and star-formation rate surface densities are also similar to both massive high-redshift clumps and submillimeter galaxies (Fig. 12).

We follow Oklopčić et al. (2017) in using the virial parameter, $\alpha_{\text{vir}}$, to investigate whether these clumps are gravitationally bound. The quantity $\alpha_{\text{vir}}$ was first introduced by Bertoldi & McKee (1992) and can be used as a measure to compare gravitational binding and turbulent energy in the gas, $\alpha_{\text{vir}} = 5\sigma^2 R/GM_{\text{gas}}$, where $\sigma$ is the velocity dispersion of the main emission line components of the clumps, R the clump radius, $M_{\text{gas}}$ the gas mass in this emission line, and G the gravitational constant. Values of $\alpha_{\text{vir}} \sim 1–5$ are typical for star-forming molecular gas at low redshift (Heyer et al. 2009), and have also been observed at very high star-formation intensities at high redshift (C17a). They are considered a signature of turbulence-regulated star formation (Krumholz & McKee 2005). Values below unity show that most of the gas within the structure is likely gravitationally bound.

With the radius set by the beam size and gravitational magnification at the cloud center, and cloud mass and velocity dispersion as listed in Table 7, we find values for $\alpha_{\text{vir}}$ of 1.4 and 2.9 for the two systemic components. $\alpha_{\text{vir}}$ would of course be lower if we had strongly overestimated R, as would seem possible for spatially unresolved clouds. However, this would be in contradiction with the non-detections of these clumps with the SMA in the VEXT configuration, which rule out very compact sizes with projected diameters much below about 0.6" (Sect. 5.5).

Our results suggest that the gas in these clumps is marginally gravitationally bound, as in many star-forming molecular clouds in the Milky Way, and perhaps also in intensely star-forming galaxies at high redshift out to the highest gas mass and star-formation surface densities (C17a). This is consistent with the simulations of Mayer et al. (2016) and Oklopčić et al. (2017) with detailed feedback descriptions. These authors argue that feedback in clumpy high-redshift galaxies with stellar masses and gas fractions comparable to the Emerald should maintain star-forming clumps near or above $\alpha_{\text{vir}} \sim 1$. It is possible that this is due to the relatively low global gas fraction in the Emerald of 20% (Bournaud et al. 2014; Renaud 2018), although the example of the Ruby, one of the most intensely star-forming galaxies at high redshift, also suggests $\alpha_{\text{vir}} \gtrsim 1$ (C17a). The two clumps we consider here have masses of about $10^9 M_{\odot}$, in the upper mass range in these simulations. We now analyze our evidence for feedback in the Emerald, before discussing its potential role for the clumps and disk in this galaxy.

### 6.3. Feedback from star formation

Several authors have recently discussed the importance of winds for cloud stability (Genel et al. 2012; Mayer et al. 2016; Oklopčić et al. 2017). If winds remove large fractions of the cloud mass in a few tens of Myr, they may lead to clump dissolution, even if the initial mass was large enough to form a self-gravitating and star-forming cloud. These timescales are shorter than the age of the young stellar population in the Emerald, and the smooth distribution of stellar light compared to the clumpy star formation may provide circumstantial evidence for this.

The spectra in Fig. 13 show that only one clump, #1.2 N, can be fitted well with a single Gaussian distribution, the other, #1.2 S, has a pronounced blue component with FWHM = 333 km s$^{-1}$ and a blueshift of -199 km s$^{-1}$ relative to the systemic line, which is also seen as a secondary component in the line maps in Fig. 8. Blueshifted emission-line components are generally viewed as a characteristics of winds, which may or may not escape from the host galaxy, depending on the outflow velocity and depth of the gravitational potential. While emission-line signatures of winds of warm ionized gas with velocities of a few hundred km s$^{-1}$ are common at high redshift (e.g., Le Tiran et al. 2011; Nesvadba et al. 2007; Genzel et al. 2011), this is, to our knowledge, the first such component seen in molecular gas in an intensely star-forming galaxy at high redshift. An example for a similar, but more extreme component in the CO line emission of a high-redshift quasar has been given by Nesvadba et al. (2011). At low redshift, multiple examples of wind components seen in molecular gas exist (e.g., Weiß et al. 1999; Sturm et al. 2011; Walter et al. 2017).
<table><caption>**Table 8.** Properties of the molecular wind toward star-forming clump #1.2 S.</caption><thead><tr><th>Source</th><th>Δν<br>(km s<sup>-1</sup>)</th><th>E<sub>Δν</sub><br>(10<sup>56</sup> erg)</th><th>P<sub>Δν</sub><br>(10<sup>49</sup> dyn s)</th><th>E<sub>σ</sub><br>(10<sup>56</sup> erg)</th><th>P<sub>σ</sub><br>(10<sup>49</sup> dyn s)</th><th>dE<sub>SF</sub>/dt<br>(10<sup>43</sup> erg s<sup>-1</sup>)</th><th>dP<sub>SF</sub>/dt<br>(10<sup>35</sup> dyn)</th></tr></thead><tbody><tr><td>1.2 S (b)</td><td>199 ± 28</td><td>4.7</td><td>2.4</td><td>7.1</td><td>5.0</td><td>2.6</td><td>2.0</td></tr></tbody></table>

**Notes.** The columns are: velocity offset between the wind and systemic component, Δν; energy (E<sub>Δν</sub>) and momentum (P<sub>Δν</sub>) derived from the observed velocity offsets; energy (E<sub>σ</sub>) and momentum (P<sub>σ</sub>) corresponding to the observed line widths; energy and momentum injection rates, dE<sub>SF</sub>/dt and dP<sub>SF</sub>/dt, from star formation. Momentum and energy injection rates are normalized to an outflow timescale of 50 Myr. We provide a measurement uncertainty only for the velocity offset. For the other estimates, typical measurement uncertainties are between 10 and 20%, and likely much smaller than the astrophysical uncertainties, which are difficult to quantify accurately.

Observing a blueshifted component may imply that we only see one side of the wind, which is lifted off an optically thick gas disk. Given that the region very near the critical line also has the most prominent stellar continuum emission, it is possible that the wind has lowered the overlying columns of dust and gas. We now analyze the physical properties of this outflow in more detail, and discuss its likely impact for the clump from which it originates, and the galaxy overall.

### 6.3.1. Kinetic energy and momentum estimates

With the same assumptions as in Sect. 5.5, we estimated an intrinsic molecular gas mass of $1.2 \times 10^9 M_\odot$ for the secondary component in clump #1.2 S. With this mass estimate, and following our earlier analysis of the Ruby (C17a), we were able to use the kinetic energy and momentum in this component as constraints to investigate whether the star formation is powerful enough to produce a wind with the observed properties.

We followed Heckman et al. (2015) to estimate the momentum, and C17a to estimate the energy injection rates from star formation into the gas. We discarded a contribution from an AGN, because the spectral energy distribution shows no evidence of one (C15). Heckman et al. (2015) showed that starbursts may inject $(4.8 \times 10^{33}) \times \text{SFR dyn of momentum flux into the gas per unit of stellar mass formed, where the star-formation rate is given in } M_\odot \text{ yr}^{-1}\text{. This estimate is matched to observed outflows in very vigorous low-redshift starbursts, and is also in the range of long-term winds proposed by Dekel & Krumholz (2013) for giant clumps. It includes contributions from radiation pressure as well as mechanical feedback and is valid for gas entrained in a hot wind, following the Chevalier & Clegg (1985) approach. For the } 41 M_\odot \text{ yr}^{-1} \text{ of star formation in #1.2 S, this corresponds to a momentum injection rate of } 2 \times 10^{35} \text{ dyn. Corresponding values for the other clump and the Emerald overall are } 5.0 \text{ and } 8.4 \times 10^{35} \text{ dyn, respectively, for star-formation rates of } 100 \text{ and } 176 M_\odot \text{ yr}^{-1}\text{.}$

As in C17a, we relied on Starburst99 (Leitherer et al. 1999) to estimate that continuous star formation over a few times $10^7 \text{ yr, solar metallicity and a Chabrier stellar IMF would produce a mechanical luminosity of } 10^{41.8} \text{ erg s}^{-1} \text{ for each newly formed solar mass. For SFR } = 41 M_\odot \text{ yr}^{-1}, \text{ this corresponds to } 2.6 \times 10^{43} \text{ erg s}^{-1} \text{ in clump #1.2 S, and to } 6 \text{ and } 11 \times 10^{43} \text{ erg s}^{-1} \text{ for star-formation rates of } 100 \text{ and } 176 M_\odot \text{ yr}^{-1} \text{ in clump #1.2 N and the Emerald, respectively. All results obtained from clump #1.2 S are also listed in Table 8.}$

Following our earlier analysis of the Ruby (C17a), we were able to constrain the kinetic energy and momentum of the gas. For the kinetic energy we set $E_{\Delta\nu} = 1/2 M v^2$ for bulk, and $E_\sigma = 3/2 M \sigma^2$ for unordered motion, and we derived the corresponding momentum by setting $P_{\text{wind}} = E_{\text{wind}}/v$ (see Table 8). With $\sigma_{\text{wind}} = 141 \text{ km s}^{-1}$ and $\Delta\nu = -199 \text{ km s}^{-1}$ velocity offset between blueshifted and systemic component, we find an

average combined kinetic energy of $11.8 \times 10^{56} \text{ erg in the outflowing gas, and a momentum of } 7.4 \times 10^{49} \text{ dyn s. These kinematics can be powered by the starburst in clump #1.2 S if the current momentum and energy injection rates are being maintained for at least 11.7 Myr and 1.5 Myr, respectively. This is less than the age of the starburst of about 50 Myr estimated in Sect. 5.3.}$

### 6.3.2. Discussion

Is the starburst in the Emerald powerful enough to unbind the gas from the clump, and from the galaxy overall? To test the first question, we assumed that the clump is a virialized gas sphere, and set $v_{\text{esc}} = \sqrt{2}v_{\text{vir}} = \sqrt{2MG/5R} = 97 \text{ km s}^{-1}$ for a clump with mass M and radius R as observed. $v_{\text{esc}}$ and $v_{\text{vir}}$ are the escape and virial velocity, respectively. The resulting escape velocity is significantly lower than the velocity offset of $199 \text{ km s}^{-1}$ found between the wind and systemic component (even while discarding projection effects), suggesting that much of the outflowing gas will escape from the clump.

However, the same does not necessarily hold for the galaxy overall. For a galaxy with at least $5 \times 10^{10} M_\odot$ of mass, as implied by our gas and stellar mass estimates, we would expect an escape velocity, $v_{\text{esc}}$, of at least $v_{\text{esc}} = \sqrt{2}v_c$, i.e., $380 \text{ km s}^{-1}$ for a disk with $R = 3 \text{ kpc, compared to a velocity offset of } 199 \text{ km s}^{-1}$. The same is suggested by a more detailed estimate following Ostriker & Shetty (2011) and Newman et al. (2012). We assumed that the wind is mainly momentum-driven, and set $(P_{\text{tot}}/dM_*/dt) = 7.2 \times f_{g,0.1} \times \Sigma_{d,1000}^2/\Sigma_{\text{SFR},100}$, where $(P_{\text{tot}}/dM_*/dt)$ is the characteristic momentum injection rate or, in other words, the total momentum, $P_{\text{tot}}$, injected by the star formation per unit mass formed, $dM_*/dt$. The quantity $f_{g,0.1}$ is the gas fraction in units of $0.1$, $\Sigma_{d,1000}^2$ is the disk mass surface density in units of $1000 M_\odot \text{ pc}^{-2}$, and $\Sigma_{\text{SFR},100}$ is the star-formation rate surface density in units of $100 M_\odot \text{ yr}^{-1} \text{ kpc}^{-2}$. With a gas-to-baryonic mass fraction of 0.2, disk mass surface density $1.3 \times 10^{10} M_\odot \text{ kpc}^{-2}$, and $\Sigma_{\text{SFR}} = 49 M_\odot \text{ yr}^{-1} \text{ kpc}^{-2}$ (for SFR = $176 M_\odot \text{ yr}^{-1}$ and a disk surface of $3.6 \text{ kpc}^2$), we find that each solar mass worth of stars formed must provide a characteristic momentum injection rate of $4990 \text{ km s}^{-1}$ to balance the hydrostatic mid-plane pressure of the galaxy. This value is much greater than expected for purely radiation-pressure driven winds ($1000 \text{ km s}^{-1}$, Murray et al. 2005), and also for winds including the momentum injection from supernovae and young stars ($2500-3000 \text{ km s}^{-1}$, Ostriker & Shetty 2011; Heckman et al. 2015). It is therefore unlikely that this wind is very efficient in removing gas from the Emerald.

Nonetheless, the wind has a potentially important impact on the future of the clump. Assuming that the gas in both components of clump #1.2 S can be modeled with the same CO-to-H$_2$ conversion factor, we find a mass-outflow rate of $dM/dt = 247 M_\odot \text{ yr}^{-1}$, for a dynamical time estimate
$t_{dyn} = R/v = 1.5$ Myr, using $R = 300$ pc and $v = 199$ km s$^{-1}$. This corresponds to a mass-loading factor of about five, not unusually high for a momentum-driven wind, and suggests that the clump may lose most of its mass in about 2 Myr. This is comparable to the free-fall time of gas with average density of a few times $100 \text{ cm}^{-3}$, in the range of the average density of a gas sphere with a mass of about $5 \times 10^8 M_{\odot}$ and radius of 100 pc. While this is a highly simplified toy model and neglects any contribution from gas accretion onto the clump (which can be considerable, e.g., Dekel & Krumholz 2013), it does highlight that these outflow times are very short, and that these clumps are likely to be transient structures whose survival depends sensitively on the dynamical equilibrium between gas accretion and outflows, as suggested previously by simulations, and as also found for giant molecular clouds in the Milky Way (e.g., Murray 2011).

Although much of the outflowing gas is likely to escape from the clump itself, it remains gravitationally bound to the host galaxy, and will therefore ultimately fall back and be available for star formation again. Galactic fountains are well studied in the nearby Universe (e.g., Marinacci et al. 2011; Sancisi et al. 2008), although their re-accretion times are very long compared to the relevant timescales in these clumps. However, given the much higher gas densities and accordingly shorter cooling times in our case, it is possible that the fall-back timescales are also much shorter. It may be that re-accretion of gas from such fountains contributes significantly to the accretion that is required to maintain the clumps marginally bound and star-forming over more extended periods of time, and potentially also to the turbulence observed in the Emerald. A detailed test of this scenario would require hydrodynamic modeling.

## 7. Summary and conclusions

We have presented a detailed study of the molecular gas, dust and stellar components in the Emerald (PLCK_G165.7+49.8), as observed with IRAM and SMA interferometry of the CO(4–3) line and 850 µm dust emission, and HST/WFC3 and CFHT optical and near-infrared imaging. The Emerald is a strongly gravitationally lensed dusty star-forming galaxy at $z = 2.236$, part of Planck’s Dusty GEMS and surrounded by two compact submillimeter sources and multiple extended, near-infrared gravitational arcs. It falls behind a previously unknown massive galaxy cluster and associated large-scale filament at $z = 0.348$. We characterized the foreground environment through the density distribution, color-magnitude diagram, and photometric and spectroscopic redshifts of member galaxies. We used the cluster members and eight multiply imaged systems of background sources to constrain the strong lensing mass distribution with LENSTOOL, and computed the magnification factors of the Emerald. The Emerald is composed of a main arc formed by two images of the same region in the source plane, with intrinsic size of about $2.7 \text{ kpc} \times 1.7 \text{ kpc}$, and magnification factors of $\mu_{\text{dust}} = 29.4 \pm 5.9$, $\mu_{\text{gas}} = 24.1 \pm 4.8$ and $\mu_{\text{stars}} = 34.1 \pm 6.8$ for the dust, gas, and stars. The two nearby compact submillimeter sources are magnified by more-moderate factors. One is likely an additional counter-image of the source producing the Emerald, the other is a separate galaxy at the same redshift. With 18 kpc projected distance it is probably a companion galaxy of the Emerald.

The Emerald has intrinsic properties of a typical ULIRG, with a far-infrared luminosity, $L_{\text{FIR}} = (1.8 \pm 0.4) \times 10^{12} L_{\odot}$, SFR of $(176 \pm 35) M_{\odot} \text{ yr}^{-1}$, and dust and molecular gas masses of $(7.4 \pm 1.5) \times 10^7 M_{\odot}$ and $(1.1 \pm 0.2) \times 10^{10} M_{\odot}$, respectively. Gas and dust morphologies are clumpy, and the gas kinematics exhibit strong velocity offsets, symmetric on both

sides of the arc due to the parity flip at the position of the critical curve. Interpreting the double-imaged velocity gradient of $380 \text{ km s}^{-1}$ over the arc as the partially-sampled rotation curve of a gaseous disk, we obtain a dynamical mass of $9.1 \times 10^{10} M_{\odot}$. Gaussian line widths are in the range 90–190 km s$^{-1}$, consistent with a marginally stable disk with a Toomre parameter $Q \sim 1.3$.

The stellar morphology seen with HST/WFC3 is diffuse and filamentary, in contrast to the dust and gas morphologies, which are dominated by two doubly imaged clumps. The total, intrinsic stellar mass and mass surface density are $(4.1 \pm 0.4) \times 10^{10} M_{\odot}$ and $(11 \pm 5) \times 10^9 M_{\odot} \text{ kpc}^{-2}$, respectively, similar to other high-redshift dusty starburst galaxies. We also find $A_V = 3.9 \pm 0.3$ mag and an average gas-to-baryon fraction of about 20%.

We probed the resolved Schmidt-Kennicutt law in the two clumps, finding star-formation intensities of about $40 M_{\odot} \text{ yr}^{-1} \text{ kpc}^{-2}$, placing them just below the relationship for high-redshift starburst galaxies, and near to the properties of other well-studied high-redshift clumps and dusty star-forming galaxies in the literature. The more diffuse intraclump emission with gas mass surface densities down to $600 M_{\odot} \text{ pc}^{-2}$ falls into a similar regime.

The two clumps are massive, with gas masses of 1.0 and $5.0 \times 10^9 M_{\odot}$, and star-formation rates of 40 and $100 M_{\odot} \text{ yr}^{-1}$, respectively, and upper limits on their FWHM sizes of 200–500 pc. They are marginally gravitationally bound. The integrated CO(4–3) line profile of one of these clumps shows a significant blueshifted wing, offset by about $-200 \text{ km s}^{-1}$ relative to the main emission-line component of the clump. We interpret this offset in the usual way, that is, as an outflow signature. To our knowledge, this is the first signature of a molecular wind arising from a massive clump in a dusty star-forming galaxy at high redshift. The kinetic energy and momentum injection rates from star formation in the nearby massive clump are sufficient to drive such a wind, and to make the gas escape from the clump, but probably not from the host galaxy itself. This molecular wind has a higher mass-loading factor than those measured for ionized winds in $z \sim 2$ galaxies, as expected in scenarios where clumps are short-lived and dissolve on timescales of a few tens of Myr, except if most of the gas is replenished by continuous accretion onto the clump. Unless we see the outflow at an extreme inclination angle, the gas will probably remain bound to the host galaxy and will form a galactic fountain, potentially contributing to the fueling of subsequent star formation in the Emerald.

**Acknowledgements.** The authors would like to thank the anonymous referee for providing useful comments that helped improve the paper. RC would like to thank Raphael Gavazzi, Claudio Grillo, Lise Christensen and Vera Patrício for useful discussions. We acknowledge the staff at IRAM and SMA for carrying out the observations on which this program is based. A major part of this work is based on observations carried out with the Plateau de Bure Interferometer of IRAM, and the Canada-France-Hawaii Telescope and the Submillimeter Array on top of Mauna Kea, Hawaii. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. The Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. This work is also based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. It also uses data from SDSS-III, funded by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The
SDSS-III web site is http://www.sdss3.org/. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. RC is supported by DFF – 4090-00079. ML acknowledges CNRS and CNES for support.

## References

Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12

Barger, A. J., Cowie, L. L., Smail, I., et al. 1999, AJ, 117, 2656

Barmby, P., Huang, J.-S., Ashby, M. L. N., et al. 2008, ApJS, 177, 431

Benítez, N. 2000, ApJ, 536, 571

Bertin, E. 2010a, Astrophysics Source Code Library [record ascl:1010.063]

Bertin, E. 2010b, Astrophysics Source Code Library [record ascl:1010.068]

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140

Bothwell, M. S., Chapman, S. C., Tacconi, L., et al. 2010, MNRAS, 405, 219

Bournaud, F., & Elmegreen, B. G. 2009, ApJ, 694, L158

Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57

Brainerd, T. G., Blandford, R. D., & Smail, I. 1996, ApJ, 466, 623

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Cañamas, R., 2016, PhD Thesis, University of Paris, Saclay, IAS – Institut d'astrophysique spatiale

Cañamas, R., Nesvadba, N. P. H., Guery, D., et al. 2015, A&A, 581, A105

Cañamas, R., Nesvadba, N., Kneissl, R., et al. 2017a, A&A, 604, A117

Cañamas, R., Nesvadba, N. P. H., Kneissl, R., et al. 2017b, A&A, 600, L3

Cañamas, R., Yang, C., Nesvadba, N. P. H., et al. 2018, A&A, 620, A61

Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Cava, A., Schäerer, D., Richard, J., et al. 2018, Nat. Astron., 2, 76

Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151

Chabrier, G. 2003, PASP, 115, 763

Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJ, 714, L118

Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, MNRAS, 410, 1687

de Propris, R., Stanford, S. A., Eisenhardt, P. R., Dickinson, M., & Elston, R. 1999, AJ, 118, 719

Dekel, A., & Krumholz, M. R. 2013, MNRAS, 432, 455

Despali, G., Giocoli, C., Bonamigo, M., Limousin, M., & Tormen, G. 2017, MNRAS, 466, 181

Dessauges-Zavadsky, M., Zamojski, M., Schäerer, D., et al. 2015, A&A, 577, A50

Dessauges-Zavadsky, M., Schäerer, D., Cava, A., Mayer, L., & Tamburello, V. 2017, ApJ, 836, L22

Dole, H., Lagache, G., Puget, J.-L., et al. 2006, A&A, 451, 417

Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

Eisenhardt, P. R., De Propris, R., Gonzalez, A. H., et al. 2007, ApJS, 169, 225

Elíasódtir, Á., Limousin, M., Richard, J., et al. 2007, ArXiv e-prints [arXiv:0710.5636]

Elmegreen, B. G. 2007, ApJ, 668, 1064

Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632

Escala, A., & Larson, R. B. 2008, ApJ, 685, L31

Fassbender, R., Böhringer, H., Lamer, G., et al. 2008, A&A, 481, L73

Fassbender, R., Böhringer, H., Nastasi, A., et al. 2011, New J. Phys., 13, 125014

Ferdosi, B. J., Buddelmeijer, H., Trager, S. C., Wilkinson, M. H. F., & Roerdink, J. B. T. M. 2011, A&A, 531, A114

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

Förster Schreiber, N. M., Shapley, A. E., Genzel, R., et al. 2011, ApJ, 739, 45

Frye, B. L., Pascale, M., Zitrin, A., et al. 2018, ApJ, submitted [arXiv:1805.04790]

Fu, H., Jullo, E., Cooray, A., et al. 2012, ApJ, 753, 134

Genel, S., Naab, T., Genzel, R., et al. 2012, ApJ, 745, 11

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101

Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

Harrington, K. C., Yun, M. S., Magnelli, B., et al. 2018, MNRAS, 474, 3866

Hatsukade, B., Tamura, Y., Iono, D., et al. 2015, PASJ, 67, 93

Hayward, C. C., & Hopkins, P. F. 2017, MNRAS, 465, 1682

Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147

Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092

Hodge, J. A., Riechers, D., Decarli, R., et al. 2015, ApJ, 798, L18

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1

Inoue, S., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 456, 2052

Iono, D., Yun, M. S., Aretxaga, I., et al. 2016, ApJ, 829, L10

Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913

Johnson, T. L., Sharon, K., Gladders, M. D., et al. 2017, ApJ, 843, 78

Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, New Journal of Physics, 9, 447

Kennicutt, Jr., R. C. 1989, ApJ, 344, 685

Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Le Tiran, L., Lehnert, M. D., van Driel, W., Nesvadba, N. P. H., & Di Matteo, P. 2011, A&A, 534, L4

Lee, N., Sheth, K., Scott, K. S., et al. 2017, MNRAS, 471, 2124

Leitherer, C., Schäerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Lilly, S. J., Eales, S. A., Gear, W. K. P., et al. 1999, ApJ, 518, 641

Limousin, M., Kneib, J. P., Bardeau, S., et al. 2007a, A&A, 461, 881

Limousin, M., Richard, J., Jullo, E., et al. 2007b, ApJ, 668, 643

Limousin, M., Richard, J., Jullo, E., et al. 2016, A&A, 588, A99

Livermore, F. C., Jones, T., Richard, J., et al. 2012, MNRAS, 427, 688

MacKenzie, T. P., Scott, D., Smail, I., et al. 2014, MNRAS, 445, 201

Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675

Marinacci, F., Fraternali, F., Nipoti, C., et al. 2011, MNRAS, 415, 1534

Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, Astronomical Society of the Pacific Conference Series, 411, 251

Martinache, C., Rettura, A., Dole, H., et al. 2018, A&A, in press, DOI: 10.1051/0004-6361/201833198

Mayer, L., Tamburello, V., Lupi, A., et al. 2016, ApJ, 830, L13

Murray, N. 2011, ApJ, 729, 133

Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569

Nayyeri, H., Cooray, A., Jullo, E., et al. 2017, ApJ, 844, 82

Nesvadba, N. P. H., Lehnert, M. D., Genzel, R., et al. 2007, ApJ, 657, 725

Nesvadba, N. P. H., Polletta, M., Lehnert, M. D., et al. 2011, MNRAS, 415, 2359

Newman, S. F., Genzel, R., Förster-Schreiber, N. M., et al. 2012, ApJ, 761, 43

Oklopčić, A., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 465, 952

Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Pisani, A. 1996, MNRAS, 278, 697

Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration Int. XXVII. 2015, A&A, 582, A30

Planck Collaboration XXXII. 2015, A&A, 581, A14

Planck Collaboration Int. XXXIX. 2016, A&A, 596, A100

Renaud, F. 2018, New A Rev., 81, 1

Richard, J., Jauzac, M., Limousin, M., et al. 2014, MNRAS, 444, 268

Salpeter, E. C. 1955, ApJ, 121, 161

Sancisi, R., Fraternali, F., Oosterloo, T., & van der Hulst, T. 2008, A&A Rev., 15, 189

Schmidt, M. 1959, ApJ, 129, 243

Sharon, K., Gladders, M. D., Rigby, J. R., et al. 2012, ApJ, 746, 161

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144

Song, J., Mohr, J. J., Barkhouse, W. A., et al. 2012a, ApJ, 747, 58

Song, J., Zenteno, A., Stalder, B., et al. 2012b, ApJ, 761, 22

Spilker, J. S., Marrone, D. P., Aguirre, J. E., et al. 2014, ApJ, 785, 149

Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJ, 733, L16

Swinbank, A. M., Chapman, S. C., Smail, I., et al. 2006, MNRAS, 371, 465

Swinbank, A. M., Smail, I., Longmore, S., et al. 2010, Nature, 464, 733

Swinbank, A. M., Papadopoulos, P. P., Cox, P., et al. 2011, ApJ, 742, 11

Swinbank, A. M., Dye, S., Nightingale, J. W., et al. 2015, ApJ, 806, L17

Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tamburello, V., Mayer, L., Shen, S., & Wadsley, J. 2015, MNRAS, 453, 2490

Timmons, N., Cooray, A., Riechers, D. A., et al. 2016, ApJ, 829, 21

Toomre, A. 1964, ApJ, 139, 1217

Usero, A., Leroy, A. K., Walter, F., et al. 2015, AJ, 150, 115

Valtchanov, I., Altieri, B., Berta, S., et al. 2013, MNRAS, 436, 2505

Veilleux, S., Meléndez, M., Sturm, E., et al. 2013, ApJ, 776, 27

Voges, W., Aschenbach, B., Boller, T., et al. 2000, IAU Circ., 7432

Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265

Weiß, A., Walter, F., Neininger, N., & Klein, U. 1999, A&A, 345, L23

Williams, J. P., de Geus, E. J., & Blitz, L. 1994, ApJ, 428, 693