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

A KECK ADAPTIVE OPTICS SURVEY OF A REPRESENTATIVE SAMPLE OF GRAVITATIONALLY LENSED STAR-FORMING GALAXIES: HIGH SPATIAL RESOLUTION STUDIES OF KINEMATICS AND METALLICITY GRADIENTS

NICHA LEETHOCHAWALIT¹, TUCKER A. JONES²,<fnref n="6" />, RICHARD S. ELLIS¹, DANIEL P. STARK³, JOHAN RICHARD⁴, ADI ZITRIN¹,<fnref n="6" />, AND MATTHEW AUGER⁵

¹ Cahill Center for Astronomy and Astrophysics, California Institute of Technology, MS 249-17, Pasadena CA 91125, USA; nichal@astro.caltech.edu

²Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu HI 96822, USA

³Department of Astronomy, University of Arizona, Tucson, AZ 85721, USA

⁴CRAL, Observatoire de Lyon, F-69561 Saint-Genis-Laval, France

⁵Institute of Astronomy, University of Cambridge, Cambridge CB3 0HA, UK

Received 2015 September 8; accepted 2016 January 28; published 2016 March 22

ABSTRACT

We discuss spatially resolved emission line spectroscopy secured for a total sample of 15 gravitationally lensed star-forming galaxies at a mean redshift of $z \simeq 2$ based on Keck laser-assisted adaptive optics observations undertaken with the recently improved OSIRIS integral field unit (IFU) spectrograph. By exploiting gravitationally lensed sources drawn primarily from the CASSOWARY survey, we sample these sub-L* galaxies with source-plane resolutions of a few hundred parsecs ensuring well-sampled 2D velocity data and resolved variations in the gas-phase metallicity. Such high spatial resolution data offer a critical check on the structural properties of larger samples derived with coarser sampling using multiple-IFU instruments. We demonstrate how kinematic complexities essential to understanding the maturity of an early star-forming galaxy can often only be revealed with better sampled data. Although we include four sources from our earlier work, the present study provides a more representative sample unbiased with respect to emission line strength. Contrary to earlier suggestions, our data indicate a more diverse range of kinematic and metal gradient behavior inconsistent with a simple picture of well-ordered rotation developing concurrently with established steep metal gradients in all but merging systems. Comparing our observations with the predictions of hydrodynamical simulations suggests that gas and metals have been mixed by outflows or other strong feedback processes, flattening the metal gradients in early star-forming galaxies.

Key words: galaxies: abundances – galaxies: evolution – galaxies: high-redshift – galaxies: kinematics and dynamics – galaxies: starburst

1. INTRODUCTION

Statistical surveys of star-forming galaxies at the peak era of activity, corresponding to a redshift interval of $z \simeq 1.5-3$, have been increasingly complemented over the past few years by spatially resolved spectroscopic observations. Early work focused on characterizing global trends such as the evolution of star formation rate (SFR), density with redshift (Madau & Dickinson 2014), and evolution of the “main sequence” of SFR as a function of galaxy stellar mass (Noeske et al. 2007; Peng et al. 2010). Spatially resolved data from integral field unit (IFU) spectrographs have provided complementary information on the ionized gas kinematics (Förster Schreiber et al. 2006, 2009; Stark et al. 2008; Law et al. 2009; Jones et al. 2010a; Wisnioski et al. 2011, 2015; Swinbank et al. 2012), the size and spatial distribution of giant star-forming regions (Genzel et al. 2008; Jones et al. 2010b; Livermore et al. 2012), and radial metal abundance gradients (Cresci et al. 2010; Jones et al. 2010b, 2013, 2015; Yuan et al. 2011; Swinbank et al. 2012; Stott et al. 2014; Troncoso et al. 2014). Two-dimensional (2D) kinematic data and the properties of clumpy star-forming regions provide valuable insight into the emergence of primitive disks, the role of instabilities, and the gradual assembly of central bulges possibly formed from inwardly migrating clumps (Genzel et al. 2011). The properties of these clumps provide important evidence of gravitational

instabilities in the primitive disks, whereas correlations between their ages and radial positions give insight into the morphological evolution of these early systems (Förster Schreiber et al. 2011).

Metallicity gradients also provide an important opportunity to study the amount of feedback, i.e., the energy returned to the surrounding medium by star formation and nuclear activity. Energy input can arise from various mechanisms including supernovae, radiation pressure, and cosmic rays. These effects combine to drive large-scale outflows of gas and metals. Outflows are ubiquitously observed from systems with high SFR densities (Heckman et al. 2001), including virtually all star-forming galaxies at high redshifts (e.g., Shapley et al. 2003). Cosmological simulations confirm that outflow mass loss rates must be comparable to and often larger than SFRs in order to explain measured stellar mass functions and the mass-metallicity relation (e.g., Vogelsberger et al. 2014; Crain et al. 2015). The rates of mass and metal loss are degenerate with rates of subsequent accretion and hence are poorly known. However, this “galactic fountain” cycle may be constrained via its imprint on metallicity gradients. To first order, high rates of outflow and subsequent accretion will redistribute heavy elements, resulting in a flatter gradient. Furthermore, the redistribution of interstellar gas will affect the radial profile of star formation and hence future metal production and feedback. Evidence for these effects has been recently observed in the form of flattened metal gradients in the extended disks of

<footnote n="6">⁶ Hubble Fellow.</footnote>
nearby galaxies (e.g., Bresolin 2011). Several groups have now
explored how various forms of feedback affect metal gradients
in simulations (e.g., Yang & Krumholz 2012; Gibson
et al. 2013; Anglés-Alcázar et al. 2014). Despite a range of
methods, there is a clear consensus that stronger feedback (that
is, a higher rate of energy injection from stellar winds,
supernovae, and other sources) results in flatter gradients due to
mixing of gas and metals over larger physical scales. Various
prescriptions are able to match observations of local galaxies
but predict different behavior at early times. We therefore seek
to constrain the degree of gas cycling via feedback by directly
measuring gradients at high redshifts.

Initial IFU surveys targeted modest numbers of galaxies
beyond *z* ≃ 2 (Förster Schreiber et al. 2006; Law et al. 2009;
Jones et al. 2010a), reflecting the challenges of securing
adequate signal to noise per spatial resolution element. The
advent of multi-IFU spectrographs (e.g., KMOS on the ESO
Very Large Telescope (VLT); Sharples et al. 2013) has recently
led to a significant improvement in survey capability. For
example, Wisnioski et al. (2015) report resolved kinematic data
for over 100 galaxies within the redshift range 2 < *z* < 2.7
using this impressive instrument. However, the multi-IFU
approach comes with a major limitation in terms of angular
resolution. Typical 0.′′6 seeing conditions correspond to a
physical scale of 5 kpc at *z* ≃ 2, whereas the half-light radius
of a typical L* galaxy at this redshift is only 2 kpc. Clearly only
the largest and most massive systems can be adequately probed
with seeing-limited data. In fact, even with adaptive optics, the
physical resolution is only ≈1 kpc (Genzel et al. 2006; Law
et al. 2009). Jones et al. (2010a, 2013) illustrate the difficulties
of correctly interpreting velocity fields and abundance
gradients from data with so few resolution elements across
each galaxy. As an example, early IFU surveys suggested that
compact galaxies at *z* ≃ 2 were characterized by dispersion-
dominated kinematics with little or no rotation (Förster
Schreiber et al. 2009). Deeper observations with adaptive
optics showed that most of these sources in fact harbor rotating
thick disks with circular velocities commensurate with their
lower masses compared to rotation-dominated systems (New-
man et al. 2013), confirming earlier results from lensed galaxies
with superior resolution (Jones et al. 2010a).

Studying gravitationally lensed systems with adaptive optics
represents a highly valuable route to addressing the physics of
galaxy formation, particularly for the less massive and more
abundant systems at *z* > 2. Lensing magnification enables
higher spatial resolution and better sampling than is otherwise
possible. In a pioneering study, Stark et al. (2008) illustrated
the potential of the lensing approach by securing resolved
kinematic data for MACSJ2135-0102, a *z* = 3.075 star-
forming galaxy magnified in angular size by a factor of eight
along its major axis by a foreground galaxy and galaxy cluster.
Using the Keck OSIRIS IFU spectrograph with adaptive optics,
the source-plane velocity field was sampled with a resolution of
120 pc leading to 20 independent points on its rotation curve.
Subsequently Jones et al. (2010a) studied a sample of six
lensed galaxies, demonstrating that a high fraction of such ≲L*
systems display well-ordered velocity fields, in contradiction to
the claims of earlier less well-sampled studies of unlensed
galaxies of similar luminosities (Law et al. 2009). Likewise,
Jones et al. (2010b) derived a radial metal abundance gradient
for a lensed *z* = 2.00 galaxy utilizing the spatial variation of
[N II]/Hα and [O III]/Hα on 500 pc scales for a system with a

half-light radius of 2.9 kpc. The necessity of mapping weaker
lines such as [N II] makes the metal gradient aspect particularly
challenging, but four systems were studied in Jones et al.
(2013) showing a diverse range of gradients, possibly
dependent upon the kinematic properties (see also Cresci
et al. 2010). Jones et al. (2013) introduced a simple model
which suggests that metal gradients and their evolution should
provide insight into the radial variation in the mass loading
factor governing the amount of outflowing gas. More recently,
a KMOS study of a larger sample of 21 galaxies at *z* ≃ 1
suggests the metal gradient tends to be less steep in more
intensely star-forming systems (Stott et al. 2014). In view of the
diverse results, clearly larger samples are required, particularly
at high redshift where only lensed galaxies provide the
necessary physical resolution.

This paper presents the results of a significantly enlarged
IFU survey of lensed galaxies. Two practical developments
have motivated us to extend the original sample presented in
Jones et al. (2013). First the number of lensed targets with
known spectroscopic redshifts has increased following the
CASSOWARY survey of lensed star-forming galaxies located
primarily in the Sloan Digital Sky Survey (Stark et al. 2013).
The second is a substantial improvement in the performance of
the OSIRIS IFU spectrograph following its transfer to Keck I.
This exploits the more powerful center launch laser which
offers a brighter and more compact LGS beacon, ensuring a
much improved Strehl. The installation of a new grating also
improved the instrumental throughput by a factor of 1.5–2.5
depending on wavelength.

The aim of the present survey is to exploit recent
improvements at the Keck Observatory to extend the original
campaign, obtaining IFU spectroscopy of a larger and more
representative sample of lensed galaxies. The improved
throughput has enabled us to discard the selection criterion
adopted by Jones et al. (2013) whereby targets were initially
pre-screened with a long-slit spectrograph (NIRSPEC) to
ensure adequately strong emission lines, possibly biasing the
sample to unusually active sources. Our goals are two-fold.
First, we aim to measure the radial metallicity gradient on sub-
kiloparsec scales for a representative sample of *z* ≃ 2 galaxies
and to examine further the origin of the diverse results obtained
by various groups. Jones et al. (2013) and Yuan et al. (2013)
argued that much of the discrepancy in metal gradient
measurements, e.g., Queyrel et al. (2012) might arise from
poorly sampled data, whereas Stott et al. (2014) suggested
correlations with the specific star formation may be the cause.
Second, we aim to characterize the kinematics and star
formation on sub-kiloparsec scales for a larger sample. Our
well-sampled kinematics and star formation morphologies
allow us to evaluate the utility and conclusions drawn from
complementary larger surveys being undertaken, e.g., with
KMOS (Stott et al. 2014; Wisnioski et al. 2015), at coarser ~5
kpc resolutions.

The plan of the paper follows. In Section 2 we discuss the
sample drawn primarily from the CASSOWARY survey and
the relevant selection criteria. We present the OSIRIS spectro-
scopic observations and their initial reduction. Section 3
discusses the reductions of the resolved spectroscopic data
into the source plane utilizing the available mass models for the
foreground lenses. We discuss the kinematic properties in
Section 4 and the metal gradients in Section 5. Section 6
overall (from ≃5 to 250 $M_\odot$ yr$^{-1}$), the bulk of our sample has values in the range 10 to 100, with a median of ≃40$M_\odot$ yr$^{-1}$. This is comparable to the rates observed in the unlensed surveys (Förster Schreiber et al. 2009; Law et al. 2009; Wisnioski et al. 2011, 2015). Near-infrared photometry, essential for deriving accurate stellar masses, is available for about half the objects in our sample (and hence are not listed in Table 1). However, they lie in the range log $M_*$ ≃ 9.0 – 9.6, which is significantly less than for those in the unlensed surveys that generally probe systems from 10$^{10}$ to 10$^{11} M_\odot$. Our survey therefore samples star-forming galaxies somewhat above the main sequence as defined by Behroozi et al. (2013).

We closely followed the observing technique described in Jones et al. (2010b). Observations were undertaken in one or more near-infrared passbands (J, H, or K) to secure spatially resolved data on the $H\alpha$ and [N II] emission lines. Where practical, we continued to observe $H\beta$, and the [O III] emission lines within a shorter wavelength band. The average spectral resolution is $R$ ≃ 3600 which corresponds to ≃3–6 Å across the J, H, K bands. We used the 100 mas pixel scale which gives a field of view of at least 1.6 × 6.4 arcsec. We took short exposures of each tip/tilt star to center the position before moving to the lensed galaxy. Each exposure comprised a number of 15 minute sub-exposures with an ABAB dither pattern ABAB of increment ~2–3 arcsec ensuring that the target was present in all frames. Exposure times varied from 1 to 4 hr. The seeing during the observations varied between 0.′′4 and 1.′′5. The median AO-corrected seeing was 0.′′17.

## 2.2. Data Reduction

We used the latest OSIRIS Data Reduction Pipeline (Larkin et al. 2006)<fnref n="7" /> to perform dark subtraction and cosmic-ray rejection, followed by a direct or scaled sky subtraction prior to spectral extraction, wavelength calibration, and telluric correction using faint standard stars. Adjacent exposures were used as sky reference frames. The data cubes from each exposure were finally combined using a $\sigma$-clipped mean.

## 2.3. Emission Line Fitting

We fit a Gaussian to the various emission lines to determine the line flux, velocity, and velocity dispersion. In all cases, $H\alpha$ is the most prominent line and a key indicator of the velocity field and SFR distribution. A four-parameter Gaussian curve is fit to the $H\alpha$ line for each spatial pixel using a weighted $\chi^2$ minimization procedure. In each data cube, we select a region devoid of emission lines and calculate the weight for each wavelength from a variance over this region, $w(\lambda) = V^{-1}(\lambda)$, for use in the Gaussian fitting. We require that each $H\alpha$ detection must be above $5\sigma$ and, where necessary, we spatially smooth the data cube with a Gaussian kernel of FWHM 3 pixels to increase the signal to noise. Given a $H\alpha$ detection, we fit a Gaussian of identical width and velocity to the [N II] emission. Where observed, we fit the [O III] and $H\beta$ emission lines independently from the $H\alpha$ and [N II] emission lines but with a similar procedure. The velocities and line widths obtained from the [O III] emission lines are consistent with those obtained from $H\alpha$ to within the $1\sigma$ uncertainties. The intrinsic velocity dispersion is calculated by subtracting the

instrumental resolution measured from OH sky lines $\sigma_{inst}$ ≃ 50 km s$^{-1}$ in quadrature from the best-fit line width.

## 2.4. Flux Calibration, Extinction, and SFR

For each galaxy observed in the $H\alpha$ observation band, we use the tip/tilt reference star to calibrate the absolute flux. We fit a PSF to the image of an integrated flux over the wavelength of the star and obtain a 1D star spectrum from the derived PSF. The flux calibration is then calculated at $H\alpha$ wavelength of the galaxy. The only exception is CSWA15 whose tip/tilt reference star is an extended galaxy and its flux calibration is measured from a UKIRT infrared standard star FS26 observed in the same night. We cross checked all the $H\alpha$ flux calibrations with the flux calibrations derived from standard stars (FS9, FS11, FS19, or FS26) taken at the end of the same nights. Flux calibration derived from these standard stars generally agrees to within 25%.

Ideally we would do the same process to the $H\beta$ observation band observations. However, there are only three galaxies with tip/tilt reference stars taken in the $H\beta$ bands (CSWA139, CSWA159, and Abell773). Hence, we use the flux calibration from the infrared standard stars to calibrate the $H\beta$ fluxes when the information from tip/tilt stars is not available. The flux calibrations derived from the three tip/tilt reference stars are in agreement with those from standard stars within 15%.

We follow Jones et al. (2013) in the calculation of $H\alpha$ extinction and SFR. In short, we used the Balmer line ratios $H\alpha/H\beta$ with a Calzetti et al. (2000) reddening curve to calculate dust extinction $E(B-V)$. We then calculate the extinction of $H\alpha$, $A_{H\alpha} = 3.33E(B-V)$. The SFR is computed from total $H\alpha$ flux corrected for extinction and lensing magnification with the $H\alpha$ SFR relation from Kennicutt (1998). For the three galaxies with no observations in $H\beta$ bands (CSWA11, CSWA28, and CSWA31), we use the weighted mean and the standard deviation of $A_{H\alpha}$ of the other 12 galaxies in our sample, $\bar{A}_{H\alpha} = 0.8 \pm 1.0$, to calculate the SFRs. The resulting $A_{H\alpha}$ and SFR are listed in Table 1.

# 3. SOURCE PLANE RECONSTRUCTION

By observing gravitationally lensed systems we can secure much higher spatial resolution and sampling for our targets than would otherwise be the case. Nonetheless, this gain in resolution is only possible by using accurate mass models for the lens system via which our observations in the image plane can be transferred into the (unlensed) source plane. The key to developing appropriate mass models is the correct identification of multiply imaged systems ideally with spectroscopic redshifts. Given the variety of lenses surveyed in this study, from SDSS galaxies to Abell clusters, it has not been possible to adopt a uniform approach to constructing these mass models across our sample.

When *Hubble Space Telescope (HST)* images or Gemini Science Archive images were not available, we took photometric BRI images in good seeing with the Echellette Spectrograph and Imager (ESI) on the Keck II telescopes. In some cases, e.g., CSWA159, mass models were already available in the literature (Dahle et al. 2013). Depending on the circumstances, three distinct methods were used in constructing new mass models as detailed below. The methods and imaging data employed to develop the mass models are summarized in Table 2.

<footnote n="7">⁷ Data Reduction Pipeline Version 3.2.</footnote>
<table><caption>**Table 2**
Mass Models</caption><thead><tr><th>ID</th><th>Method</th><th>Image</th><th>Photometric bands</th></tr></thead><tbody><tr><td>cswa11</td><td>LTM</td><td>Gemini</td><td>gri</td></tr><tr><td>cswa15</td><td>LTM</td><td>Keck ESI</td><td>BRI</td></tr><tr><td>cswa19</td><td>LTM</td><td>HST</td><td>GVIH</td></tr><tr><td>cswa20</td><td>MA</td><td>SDSS</td><td>gri</td></tr><tr><td>cswa28</td><td>LTM</td><td>HST</td><td>GVI</td></tr><tr><td>cswa31</td><td>LTM</td><td>Gemini</td><td>gri</td></tr><tr><td>cswa128</td><td>MA</td><td>SDSS</td><td>gri</td></tr><tr><td>cswa139</td><td>LTM</td><td>Keck ESI</td><td>VRI</td></tr><tr><td>cswa159</td><td>LTM</td><td>SDSS</td><td>gri</td></tr><tr><td>cswa165</td><td>MA</td><td>SDSS</td><td>gri</td></tr><tr><td>a773</td><td>JR</td><td>HST</td><td>YJH</td></tr></tbody></table>

We used the Light-Traces-Mass (LTM) method (Zitrin et al. 2015) to develop mass models for CSWA11, CSWA15, CSWA19, CSWA28, CSWA31, CSWA139, and CSWA159. Generally these lenses are dominated by a single galaxy with a few associated companions as they represent lenses with angular Einstein diameters in excess of the size of an SDSS fiber (>3 arcsec). The method assumes that the lensing mass distribution of lensing galaxies is described by an elliptical power law, with a Gaussian core for the brightest central galaxies, scaled according to their luminosities. The member galaxies are identified spectroscopically and/or via a color cut. To represent the dark matter, the mass distribution of the galaxies is smoothed by a 2D Gaussian kernel. Multiple images are then identified from their similar colors and/or spectroscopic redshifts, if available, in conjunction with a preliminary guess for the mass model. In many cases, spectroscopic redshifts are available in the literature (Dahle et al. 2013; Stark et al. 2013; Brewer et al. 2011) and from our OSIRIS observations. In a few cases where redshifts are not available for identified multiple images, we leave the redshift as a free parameter. In constructing the final best-fit mass model, we employ a Markov Chain Monte Carlo (MCMC) to minimize the image position $\chi^2$, i.e., the total distance between the predicted and actual position of the multiple images. The model consists of six basic parameters: the exponent of the power law, the smoothing Gaussian width, the relative weight between the galaxy component and the dark matter component, the overall normalization of the mass distribution, and the two external shear parameters. The mass-to-light ratio is fixed for all galaxies except that for the brightest central galaxies that we left as a free parameter in each model.

Some mass models for CASSOWARY sources, CSWA20, CSWA128, and CSWA165 were already available following modeling developed by Matt Auger (M.A.) as part of that survey (Stark et al. 2013). In these cases, gri imaging from SDSS was used to fit each lensing component with a singular isothermal ellipsoid following the procedure discussed in Auger et al. (2011). The same mass model is fitted to each of the filters simultaneously, although the amplitudes of the surface brightness profiles can vary.

To determine how uncertainties in the lensing models affect our analyses of both the kinematic state of each galaxy and its chemical gradient, we empirically analyzed a marginal case of the lensing models. Since CSWA28 has a relatively high magnification factor and contains only one set of multiple images in constraining the lensing model, i.e., the magnification is the least well-constrained, we selected CSWA28 as the

appropriate case. We constructed 20 different mass models following the parameters randomly drawn from the MCMC and performed the same kinematic property analysis (simple disk model fitting) as in Section 4 below. We found that in most of the cases, different lens models primarily affect the position of the source relative to the image but the morphology and size of the galaxy in the source plane is largely the same. From the kinematic model fitting in each source plane, the derived galaxy inclinations, position angles, and centers are in agreement given the uncertainties we quote for the best model. The resulting uncertainty in radius at each pixel is $\lesssim 10\%$. In a few rare cases, the result from kinematic model fittings yields a radius that is 30% higher. For CSWA28, if we have propagated the expected maximum uncertainty of 10% in radii into the metal gradient calculation, the final uncertainty in the derived $N_{2_{\text{pp04}}}$ metal gradient would have increased by only 2% from $\pm 27\%$ to $\pm 29\%$. Since CSWA28 is one of the least certain lens models, we expect that the uncertainties due to lensing models are not the dominance source of uncertainty and did not propagate this uncertainty into the subsequent kinetic modeling.

Finally, for the rich cluster, Abell773, the mass model is available from a detailed study conducted by Johan Richard (J. R.) (Richard et al. 2011, 2010; Limousin et al. 2012; Livermore et al. 2012). In this case, a parametric mass model of the central region of the cluster was developed using the LENSTOOL package (Jullo et al. 2007). Briefly, we assumed the cluster mass distribution follows a double pseudo-isothermal elliptical profile (Elíasdóttir et al. 2007) and added one or more central cluster members as smaller scale perturbations to the mass distribution. As above, optimization is conducted by minimizing the predicted and observed positions of the many multiple images in these well-studied clusters. Table 1 summarizes the overall 2D magnification $\mu$ appropriate for each source although, in the following analysis, full account is taken of the spatial dependence in transforming our data to the source plane.

# 4. KINEMATIC PROPERTIES

## 4.1. Methods and Motivation

The improved spatial sampling made possible by studying gravitationally lensed sources offers the key to resolving several uncertainties of interpretation of the kinematic and morphological properties of high redshift star-forming galaxies. These include discriminating between single rotating systems and close merging pairs as well as the validity of well-ordered disk models claimed to fit a high proportion of objects studied in the literature. For this study we also seek to examine possible correlations between the ratio of ordered and random motions ($v/\sigma$—typically used as a measure of the extent to which a system is or is not dispersion-dominated) with the presence of a metallicity gradient (see Section 5) as this might indicate important diagnostics of how disks assemble and mature over a key period in cosmic history. For these and other applications, clearly determining the reliability of the derived values of the rotational velocity and the dispersion is crucially important.

Early LGSAO-based data on lensed sources already highlighted the importance of securing well-sampled data (Jones et al. 2010a, 2010b). Prior to these studies it was commonly believed that systemic rotation was confined to the galaxies with larger stellar masses (Genzel et al. 2006; Law et al. 2009).
However, due to the small number of lensed sources available, limited statistics have remained a problem. Jones et al. (2010b) presented spatially resolved dynamics for six lensed systems claiming four of the six showed well-ordered velocity fields—a fraction consistent with that determined for larger, more luminous sources. Recently, Livermore et al. (2015) collated the available data on 17 lensed sources including 10 studied with natural guide-star-assisted AO using SINFONI on the VLT and 6 with LGSAO and Keck's OSIRIS, mostly from Jones et al. (2010b, 2013), similarly concluding ≈60% of these systems with stellar masses in the range $4 \times 10^8 < M/M_\odot < 6 \times 10^{10}$ are consistent with rotating disks. This fraction is lower than that deduced from the KMOS$^{3D}$ survey (Wisnioski et al. 2015) based on coarser-sampled data of a larger (≈30) sample of $1.9 < z < 2.7$ galaxies where, over their entire sample, it is claimed 83% of the sample is rotation-dominated with $v/\sigma > 1$ and at least 70% can be considered “disklike” systems.

Our newly enlarged set of well-sampled resolved spectroscopic data provides a valuable test of the above claims and we now examine both the validity of a simple disk model as a fit to our source plane velocity fields as well as the extent to which our systems are rotation- or dispersion-dominated. Our Hα emission maps and their associated radial velocity fields are presented in Figure 1. While our Hα-based radial velocity fields are in many cases indicative of velocity gradients, the associated Hα surface brightness distributions have irregular, asymmetric, and non-disklike shapes. Such irregular emission line morphologies are common at these redshifts and SFRs and contrast with more regular distributions in broad-band filters sampling the stellar distributions. To verify this, for each galaxy, we attempted to extract the total continuum spectrum adjacent to the Hα line. Although the resulting signal to noise of the continuum is limited beyond a scale of ≈1 kpc, it provides a valuable indicator of the center of each galaxy.

First we assess the validity of a rotating disk as a representation of the observed Hα radial velocity field following the method adopted in Jones et al. (2010b). Those authors used an arctangent function as the simplest model, vis-à-vis:

$$V(R) = V_0 + \frac{2}{\pi} V_c \arctan \frac{R}{R_t} \quad (1)$$

where R is the radius from the disk center,

$$R^2 = (x \cos \theta - y \sin \theta)^2 + \left( \frac{x \sin \theta + y \cos \theta}{\sin i} \right)^2.$$

The simple model has seven unknown parameters: the inclination i, position angle θ, the position of the disk center (αc, δc), the scale radius Rt, an asymptotic velocity Vc, and an overall systemic velocity V0. The variables x and y represent the distance of each pixel from the disk center. The resulting asymptotic velocity is corrected for the effect of inclination.

We use an MCMC method to find the best fitting disk model parameters. We fit the observed Hα radial velocity map with the model convolved with an elliptical point-spread function based on the reconstructed tip/tilt star image. In the MCMC fitting, we set flat priors within the parameter ranges as follows. The disk center is set to lie within the peak of the stellar continuum light distribution (which does not necessarily coincide with the Hα emission peak). The scale radius is confined to lie within the range of the observed Hα map. The asymptotic velocity is $\lesssim 3 \times \Delta V$ where $\Delta V$ is the maximum

velocity difference along the slit that passes through the disk center for all position angles. We allow the full range for inclination and position angles, i.e., $0^\circ \le i \le 90^\circ$ and $0^\circ \le \theta \le 360^\circ$. The uncertainty in each derived parameter is calculated from its posterior probability distribution. The results from these simple disk model fits are summarized in Table 3.

## 4.2. Results

Our data reveal many cases where there are significant kinematic deviations from this simple rotating disk model. Because our spatial resolution is <500 pc, we find the velocity field for the Hα emission line samples the local bulk motion of gas and is not smeared over a large area as would be the case in coarser-sampled data. Typically our resolution elements have velocity uncertainties of only a few km s$^{-1}$ which often results in large reduced $\chi^2$ despite similar rms velocity residuals of $\langle v_{\text{model}} - v_{\text{data}} \rangle \sim 30 \text{ km s}^{-1}$ compared to studies sampling $\sim$ kiloparsec scales (e.g., Swinbank et al. 2012).

These significant velocity deviations suggest that several of our systems are undergoing various stages of a merging process. One-dimensional velocity profiles in Figure 1 show many galaxies in our sample can be considered as a perturbed disk where the large-scale velocity map exhibits a clear gradient but the velocity dispersion peak is offset from our adopted center. In several cases (e.g., CSWA19 and CSWA139), two distinct Hα sources of different velocities can be discerned, consistent with early stages of a merging pair. There is good evidence from either the velocity profiles or the velocity deviations from our disk model, that CSWA19, 31, 128, 139, and 165 represent merging systems.

At this point, it is interesting to conjecture whether such complex kinematic patterns would have been discernible with a lower spatial resolution. As an experiment we re-analyzed two systems, CSWA19 and CSWA31, both clearly poor candidates for pure rotation-dominated systems, as if they had been observed at lower angular resolution. The source plane velocity fields were first smoothed with a typical PSF of 0.6, rebinned to a coarser resolution of 0.2, comparable to that achievable with the KMOS$^{3D}$ survey (Wisnioski et al. 2015), and re-analyzed with our disk model. The results of this experiment are shown in Figure 2. In both cases the $\chi^2$ fit is significantly improved and there is little evidence of any departure from a pure disk model. This accentuates the importance of optimally sampled data for a robust interpretation of the velocity field.

The second question posed above relates to the precision of the derived rotational velocity and, particularly, the ratio of systemic to random motions, which is such an important diagnostic of the maturity of early disk systems. To examine this, we first adopt the inclination and position angle values obtained from the disk model fitting. Despite the velocity deviations from the disk model discussed above, the posterior probability distributions of inclination and position angle are generally Gaussian in form suggesting reasonably good estimates of these parameters. Other methods for determining the inclination, e.g., analyzing the elliptical shape of the Hα or near-infrared light distributions (e.g., Newman et al. 2013; Wisnioski et al. 2015) are less applicable to our less luminous systems with irregular forms. As for the position angle, our disk model fits are generally consistent (within $2\sigma$) with values obtained from a “peak-to-peak” velocity method where a pseudo slit placed at the center is rotated to find to yield the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** From left to right: source plane Hα emission map, two-dimensional velocity field, and one-dimensional velocity (red diamonds) and dispersion (green squares) for each of our targets. Ellipses represent the source-plane PSFs. Solid lines show the slits used to extract the one-dimensional velocity and dispersion with position angles determined from our simple disk model fits (see the text in Section 4 for details). Black crosses mark the adopted centers of each galaxy.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 1. (Continued.)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 1. (Continued.)</figcaption></figure>

<table><caption>Table 3
Kinematic Properties of the Samples</caption><thead><tr><th rowspan="2">ID</th><th colspan="5">Simple-disk Model Results</th><th colspan="5">Peak-to-peak Method Results</th></tr><tr><th>χ<sup>2</sup><sub>red</sub></th><th>⟨V<sub>model</sub> − V<sub>data</sub>⟩<sub>rms</sub></th><th>i(°)</th><th>θ(°)</th><th>V<sub>c</sub> (km s<sup>−1</sup>)</th><th>V<sub>c</sub>/σ</th><th>ΔV (km s<sup>−1</sup>)</th><th>σ (km s<sup>−1</sup>)</th><th>θ'(°)</th><th>ΔV/2σ</th></tr></thead><tbody><tr><td>cswa11</td><td>17</td><td>40.2</td><td>64<sup>69</sup><sub>56</sub></td><td>16 ± 12</td><td>173<sup>216</sup><sub>125</sub></td><td>1.8<sup>2.5</sup><sub>1.1</sub></td><td>167 ± 53</td><td>94 ± 26</td><td>48</td><td>0.89 ± 0.19</td></tr><tr><td>cswa15</td><td>18</td><td>12.9</td><td>25<sup>29</sup><sub>19</sub></td><td>194 ± 10</td><td>29<sup>40</sup><sub>17</sub></td><td>0.8<sup>2.2</sup><sub>0.5</sub></td><td>48 ± 3</td><td>34 ± 5</td><td>205</td><td>0.71 ± 0.06</td></tr><tr><td>cswa19<sup>a</sup></td><td>9.6</td><td>13.6</td><td>49<sup>56</sup><sub>43</sub></td><td>155 ± 4</td><td>126<sup>159</sup><sub>83</sub></td><td>1.6<sup>2.1</sup><sub>1.0</sub></td><td>89 ± 2</td><td>78 ± 14</td><td>153</td><td>0.57 ± 0.05</td></tr><tr><td>cswa19<sup>b</sup></td><td>246</td><td>42.3</td><td>75<sup>77</sup><sub>65</sub></td><td>156 ± 4</td><td>377<sup>442</sup><sub>248</sub></td><td>3.6<sup>5.8</sup><sub>2.8</sub></td><td>120 ± 7</td><td>78 ± 14</td><td>150</td><td>1.54 ± 0.30</td></tr><tr><td>cswa20</td><td>35</td><td>21.8</td><td>40<sup>47</sup><sub>32</sub></td><td>162 ± 2</td><td>133<sup>446</sup><sub>99</sub></td><td>2.2<sup>3.0</sup><sub>1.1</sub></td><td>46 ± 14</td><td>62 ± 25</td><td>168</td><td>0.37 ± 0.09</td></tr><tr><td>cswa28</td><td>6.6</td><td>25.3</td><td>23<sup>47</sup><sub>18</sub></td><td>129 ± 10</td><td>91<sup>146</sup><sub>49</sub></td><td>1.6<sup>2.2</sup><sub>0.8</sub></td><td>29 ± 25</td><td>56 ± 12</td><td>127</td><td>0.26 ± 0.12</td></tr><tr><td>cswa31</td><td>93</td><td>20.8</td><td>22<sup>27</sup><sub>19</sub></td><td>249 ± 3</td><td>122<sup>135</sup><sub>110</sub></td><td>2.2<sup>2.6</sup><sub>1.8</sub></td><td>102 ± 12</td><td>56 ± 8</td><td>244</td><td>0.91 ± 0.08</td></tr><tr><td>cswa128</td><td>990</td><td>69.9</td><td>57 ± 3</td><td>210 ± 2</td><td>570<sup>578</sup><sub>554</sub></td><td>8.6<sup>2.4</sup><sub>4.8</sub></td><td>262 ± 8</td><td>66 ± 29</td><td>208</td><td>1.98 ± 0.44</td></tr><tr><td>cswa139</td><td>38</td><td>25.6</td><td>24 ± 4</td><td>297 ± 15</td><td>265<sup>290</sup><sub>211</sub></td><td>4.3<sup>6.5</sup><sub>1.9</sub></td><td>166 ± 27</td><td>62 ± 32</td><td>294</td><td>1.34 ± 0.36</td></tr><tr><td>cswa159</td><td>7.6</td><td>12.7</td><td>13<sup>19</sup><sub>10</sub></td><td>175 ± 9</td><td>141<sup>191</sup><sub>85</sub></td><td>2.3<sup>3.5</sup><sub>1.1</sub></td><td>79 ± 29</td><td>61 ± 22</td><td>164</td><td>0.65 ± 0.17</td></tr><tr><td>cswa165</td><td>141</td><td>17.3</td><td>57 ± 5</td><td>164 ± 5</td><td>359<sup>384</sup><sub>307</sub></td><td>5.9<sup>8.5</sup><sub>3.1</sub></td><td>83 ± 28</td><td>44 ± 19</td><td>218</td><td>0.94 ± 0.26</td></tr><tr><td>a773</td><td>17</td><td>53.7</td><td>22<sup>37</sup><sub>17</sub></td><td>144 ± 12</td><td>176<sup>225</sup><sub>108</sub></td><td>3.9<sup>3.3</sup><sub>2.2</sub></td><td>54 ± 22</td><td>45 ± 10</td><td>158</td><td>0.60 ± 0.28</td></tr></tbody></table>

**Notes.**

a Results from fitting to the main galaxy only.

b Results from fitting to both mergers together as one system.

largest velocity shear. With the exception in CSWA165 where the discrepancy between the two methods is large and the disk model fit is particularly poor, we adopt the position angle from the peak-to-peak method ($\theta'$ in Table 3).

Our sample can then be split into two subsets. For five galaxies (CSWA11, 15, 28, A773, 159) where the disk model represents an acceptable fit ($\chi^2_{red} < 20$), the ratio $v_c/\sigma$ is probably a valid indicator of whether the system is dynamically supported. Adopting the standard often-used threshold $v_c/\sigma \ge 1$ and calculating the intrinsic velocity dispersion $\sigma$ from the weighted velocity dispersions of each spaxel along the major axis, all but CSWA15 could be considered to be rotationally supported systems (see Table 3). However, the $v_c/\sigma$ ratio is clearly inappropriate for the remaining galaxies where the disk model is a much poorer fit. An alternative statistic sometimes used is the observed velocity ratio $\Delta V/2\sigma$ where $\Delta V$ is the maximum velocity shear obtained from the peak-to-peak method. Applying the threshold criterion of $\Delta V/2\sigma \ge 0.4$ for rotationally supported systems adopted by Förster Schreiber et al. (2009) to the entire sample, 9 out of 11 of our sources would then be classified as rotation-dominated systems. Figure 3 shows the correspondence between these two measures ($v_c/\sigma$, $\Delta V/2\sigma$) and the goodness of the disk model ($\chi^2_{red}$) for our sample. Broadly speaking, both above criteria are reasonably consistent in selecting rotationally supported systems. However, without the highly sampled data that allow

an adequate test of the goodness of fit, a much larger number would be incorrectly placed in this category.

We perform the same kinematic fitting method to the four galaxies published in Jones et al. (2013) and find consistent categorization using the $\chi^2_{red}$ criterion above. Jones et al. (2010a) already performed kinetic model fitting for J0744 and J1206 and we find very similar results. J1148, which Jones et al. (2013) considered a rotating system using a 1D velocity profile, has a $\chi^2_{red} \approx 6.5$ with a very low velocity residual rms of $7 \text{ km s}^{-1}$ suggesting a well-ordered disk. Finally, J1038, which is visibly a merger, has a $\chi^2_{red} \approx 380$ and is categorized as a perturbed disk.

It is illustrative, therefore, to compare our conclusions to the first results of the KMOS$^{3D}$ Survey which reports kinematic results of an IFU survey of ~200 galaxies at $0.7 < z < 2.7$ (Wisnioski et al. 2015). Our galaxy sample has an average $\Delta V/\text{sini}$ of $148 \text{ km s}^{-1}$ and average velocity dispersion of $60 \text{ km s}^{-1}$, both of which are broadly consistent with the range of $\Delta V/\text{sini} \sim 100-200$ and $\sigma \sim 30-80 \text{ km s}^{-1}$ in the KMOS$^{3D}$ data. However, according to Figure 3 only 4 out of 11 members (~36%) of our sample would be classified as rotationally supported galaxies via the combination of both criteria in Figure 3, whereas in the KMOS$^{3D}$ 74% of $z \sim 2$ galaxies were classified as being rotationally supported. The higher percentages in the KMOS$^{3D}$ analysis may arise at least in part from a combination of the lower spatial resolution where late merging systems are mistaken as a regular rotating system
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Potential dangers of interpreting kinematic data with poor sampling. Left: source plane Hα intensity map and velocity field in this study for CSWA19 (top row) and CSWA31 (bottom row). Right: respective simulated Hα intensity map and velocity field indicating the deterioration in resolution equivalent to that for non-lensed sources studied with KMOS. In both cases, the complex morphology and velocity structure is lost in the poorer sampled data leading to the erroneous conclusion of kinematically well-ordered systems.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of two methods for determining the degree of rotational support in our sample plotted against the goodness of fit in the disk model. $v_c/\sigma$ from our disk fits and $\Delta V/2\sigma$ from the peak-to-peak method is considered in the context of the reduced $\chi^2$. Horizontal dashed lines represent the usually adopted thresholds, $\Delta V/2\sigma \ge 0.4$ and $v_c/\sigma > 1$, for rotationally supported systems and the vertical dashed line represents $\chi^2 = 20$, below which we consider the disk model to be appropriate. Without adequate sampling, a larger fraction of our sample would be incorrectly considered to be rotationally supported.</figcaption></figure>

(Figure 2) as well as the difficulty of fitting a disk model in such circumstances.

## 5. METAL GRADIENTS

Given the capability of the current integral field instruments, the only practical method to measure metallicity in high-redshift galaxies is to use strong line metallicity calibrators. The direct measurement of metallicity ($T_e$ method) is not feasible for individual galaxies at this redshift since it relies on the measurement of emission lines that are too faint to detect. In this work, we use the [N II] $\lambda6584/H\alpha$ ratio (N2) and the ([O III] $\lambda5008/H\beta$)/([N II] $\lambda6584/H\alpha$) ratio (O3N2) with the calibrations from Pettini & Pagel (2004, PP04), Steidel et al. (2014, S14), and Maiolino et al. (2008, M08):

$$12 + \log(O/H) = 8.90 + 0.57 \times N2 \quad (\text{N2, PP04})$$

$$12 + \log(O/H) = 8.62 + 0.36 \times N2 \quad (\text{N2, S14})$$

$$12 + \log(O/H) = 8.73 - 0.32 \times O3N2 \quad (\text{O3N2, PP04})$$

$$12 + \log(O/H) = 8.66 - 0.28 \times O3N2 \quad (\text{O3N2, S14})$$

where $N2 \equiv \log([\text{N II}]\lambda6584/H\alpha)$ and $O3N2 \equiv \log([\text{O III}]\lambda5008/H\beta)/([\text{N II}]\lambda6584/H\alpha))$.

The PP04 relations were primarily calibrated from the $T_e$ method based on nearby extragalactic H II regions. The M08 relations were based on local galaxies that cover a wider range of metallicity compared to the PP04 relations. The S14 relations were based on a subset of the same nearby extragalactic H II regions sample in PP04 relations that are limited to the range of N2 and O3N2 observed in $z \sim 2.3$ sample which have been shown to lower the systematic offset between the two indicators when applied to high-redshift samples. To counter the possible bias of the inferred metallicity due to the
<table><caption>**Table 4**
Metallicity Gradients</caption><thead><tr><th>ID</th><th>[N II]/Hα</th><th>Central 12+log(O/H)<br>(from N2 PP04)</th><th>N2 PP04<br>dex kpc<sup>-1</sup></th><th>N2 S+ 14<br>dex kpc<sup>-1</sup></th><th>N2 M+ 08<br>dex kpc<sup>-1</sup></th><th>O3N2PP04<br>dex kpc<sup>-1</sup></th><th>O3N2 S+ 14<br>dex kpc<sup>-1</sup></th><th>O3N2 M+ 08<br>dex kpc<sup>-1</sup></th></tr></thead><tbody><tr><td colspan="9">Objects Classified as Rotationally Supported Systems (χ<sup>2</sup> &lt; 20)</td></tr><tr><td>cswa11</td><td>0.14 ± 0.01</td><td>8.54 ± 0.06<br>8.60 ± 0.06</td><td>-0.07 ± 0.02<br>-0.11 ± 0.02</td><td>-0.05 ± 0.01<br>-0.07 ± 0.01</td><td>-0.10 ± 0.03<br>-0.16 ± 0.02</td><td>...</td><td>...</td><td>...</td></tr><tr><td>cswa15</td><td>0.03 ± 0.01</td><td>8.16 ± 0.01<br>8.17 ± 0.01</td><td>-0.04 ± 0.01<br>-0.03 ± 0.01</td><td>-0.02 ± 0.01<br>-0.02 ± 0.01</td><td>-0.05 ± 0.02<br>-0.05 ± 0.01</td><td>-0.02 ± 0.01<br>-0.02 ± 0.01</td><td>-0.02 ± 0.01<br>-0.01 ± 0.01</td><td>-0.01 ± 0.01<br>-0.00 ± 0.01</td></tr><tr><td>cswa28</td><td>0.09 ± 0.01</td><td>8.10 ± 0.01<br>8.33 ± 0.12</td><td>0.11 ± 0.03<br>0.04 ± 0.03</td><td>0.07 ± 0.24<br>0.02 ± 0.01</td><td>0.17 ± 0.59<br>0.07 ± 0.04</td><td>...</td><td>...</td><td>...</td></tr><tr><td>cswa159</td><td>0.08 ± 0.01</td><td>8.50 ± 0.01<br>8.49 ± 0.03</td><td>-0.01 ± 0.01<br>-0.01 ± 0.01</td><td>-0.01 ± 0.01<br>0.00 ± 0.01</td><td>-0.03 ± 0.02<br>-0.02 ± 0.01</td><td>0.00 ± 0.01<br>0.00 ± 0.01</td><td>0.00 ± 0.01<br>0.00 ± 0.01</td><td>-0.01 ± 0.01<br>-0.01 ± 0.01</td></tr><tr><td>a773</td><td>0.17 ± 0.07</td><td>8.42 ± 0.10<br>8.43 ± 0.05</td><td>-0.13 ± 0.18<br>-0.05 ± 0.02</td><td>-0.08 ± 0.06<br>-0.03 ± 0.01</td><td>-0.19 ± 0.28<br>-0.07 ± 0.03</td><td>-0.25 ± 0.41<br>-0.03 ± 0.03</td><td>-0.22 ± 0.35<br>-0.02 ± 0.02</td><td>0.01 ± 0.30<br>-0.08 ± 0.02</td></tr><tr><td colspan="9">Objects with χ<sup>2</sup> &gt; 20</td></tr><tr><td>cswa19</td><td>0.09 ± 0.01</td><td>8.30 ± 0.01<br>8.26 ± 0.01</td><td>-0.01 ± 0.01<br>0.02 ± 0.01</td><td>0.01 ± 0.01<br>0.01 ± 0.01</td><td>-0.01 ± 0.07<br>0.02 ± 0.01</td><td>0.12 ± 0.04<br>0.01 ± 0.01</td><td>0.10 ± 0.01<br>0.01 ± 0.01</td><td>0.07 ± 0.05<br>0.03 ± 0.01</td></tr><tr><td>cswa20</td><td>0.05 ± 0.01</td><td>7.89 ± 0.02<br>7.97 ± 0.04</td><td>-0.15 ± 0.05<br>0.05 ± 0.02</td><td>-0.10 ± 0.07<br>0.03 ± 0.01</td><td>-0.29 ± 0.45<br>0.09 ± 0.03</td><td>-0.41 ± 0.19<br>0.03 ± 0.02</td><td>-0.36 ± 0.17<br>0.02 ± 0.02</td><td>-1.00 ± 0.24<br>0.25 ± 0.05</td></tr><tr><td>cswa31</td><td>0.33 ± 0.06</td><td>8.50 ± 0.03<br>8.54 ± 0.01</td><td>0.02 ± 0.01<br>0.00 ± 0.01</td><td>0.01 ± 0.01<br>0.00 ± 0.01</td><td>0.03 ± 0.01<br>0.00 ± 0.01</td><td>...</td><td>...</td><td>...</td></tr><tr><td>cswa128</td><td>0.10 ± 0.02</td><td>8.57 ± 0.03<br>8.51 ± 0.01</td><td>-0.10 ± 0.02<br>-0.04 ± 0.01</td><td>-0.06 ± 0.01<br>-0.02 ± 0.01</td><td>-0.12 ± 0.02<br>-0.05 ± 0.01</td><td>-0.09 ± 0.02<br>-0.04 ± 0.01</td><td>-0.08 ± 0.02<br>-0.03 ± 0.01</td><td>-0.15 ± 0.02<br>-0.06 ± 0.01</td></tr><tr><td>cswa139</td><td>0.17 ± 0.08</td><td>8.25 ± 0.01<br>8.12 ± 0.03</td><td>-0.01 ± 0.01<br>0.04 ± 0.01</td><td>-0.01 ± 0.01<br>0.02 ± 0.01</td><td>-0.02 ± 0.02<br>0.05 ± 0.01</td><td>0.03 ± 0.01<br>0.02 ± 0.01</td><td>0.03 ± 0.01<br>0.01 ± 0.01</td><td>0.02 ± 0.01<br>0.04 ± 0.01</td></tr><tr><td>cswa165</td><td>0.26 ± 0.06</td><td>8.53 ± 0.01<br>8.53 ± 0.01</td><td>0.03 ± 0.01<br>0.01 ± 0.02</td><td>0.02 ± 0.02<br>0.01 ± 0.01</td><td>0.04 ± 0.05<br>0.01 ± 0.04</td><td>0.13 ± 0.01<br>0.01 ± 0.02</td><td>0.11 ± 0.01<br>0.00 ± 0.01</td><td>0.16 ± 0.05<br>0.13 ± 0.03</td></tr></tbody></table>

Note. For each object, the top line contains values derived along each galaxy’s “major axis.” The bottom line contains values derived from radial binning.

higher N/O ratios in high-redshift galaxies than the ratios in local galaxies, we will compare our measured gradients from the N2 calibrators to the gradients from the O3N2 indicators which have been shown to be less dependent on the N/O ratio than the N2 indicator (Steidel et al. 2014).

For each calibrator, we calculate the metallicity gradient using two methods: from the pixels along the major axis and from all pixels binned radially. The galaxy center, position angle, inclination, and pseudo slit along the major axis are the same as those in Section 4. An advantage of the metallicity gradient measured from the pixels along the major axis is that it is independent of the inclination and thus less model-dependent. However, it is based on fewer measurements than the metallicity gradient measured from all pixels and sometimes not feasible when most of the H II regions do not lie along the major axis (e.g., CSWA20). The best-fit gradients and central metallicities are shown in Table 4 and plots are shown in Figure 9.

In detail, we calculate an average N2 ratio (or O3N2 ratio) for each radial bin along the major axis (or along the annuli) for PP04 and S14 calibrators. We assume that metallicity 12 + log(O/H) is a linear function of radius. Hence, the line ratio is an exponential function of radius, metallicity gradient, and central metallicity. We fit the N2 ratios (or O3N2 ratios)

with the weighted least-square regression to find the metallicity gradient and central metallicity for each galaxy.

For the M08 calibrator, we calculate the N2 ratio (and [O III]/Hβ ratio, if possible) for each radial bin and find the metallicity for that radial bin with a brute force maximum likelihood estimation. We then obtain the metallicity gradient and the central metallicity with a linear fit.

We check the consistency between the metallicity gradients derived from pixels along the major axis and radial annuli obtained in Section 4. We found that the results from the two methods are consistent with each other within one standard deviation especially when the rotation is relatively well-described by a simple disk rotation. The deviation is larger when the major axis is short and the measurement along the major axis suffers from a poor sampling e.g., CSWA20 and Abell773.

Metallicity gradients derived from the O3N2 calibrators are consistent with the gradients derived from the N2 calibrators when using the radial binning method and comparing between the same set of calibrators, e.g., PP04 N2 with PP04 O3N2. Although O3N2 calibrators give more accurate metallicity measurements than N2 calibrators, the consistency in the gradients derived from the two calibrators suggests that the metallicity gradients measured with N2 calibrators should not be heavily affected by a possible variation in N/O ratios in
<figure><img src="image_7.png" /><figcaption>**Figure 4.** BPT diagrams for CSWA15 (upper) and CSWA165 (lower). Gray points are from individual pixels in each galaxy. Red points are from pixels whose values lie significantly above the maximum starburst line (green curve, Kewley et al. 2001). The region corresponding to the red points is shown in red in the map at the corner. The blue curve is the best fitted line for the $z \sim 2.3$ samples in Steidel et al. (2014). Median values in uncertainties are shown in the purple error bars. A few pixels have large error bars, e.g., the data point at the top of the lower plot.</figcaption></figure>

high-redshift galaxies. Hence, for the three galaxies that we do not have information on $[O\,III]/H\beta$, the derived metallicity gradients should be relatively reliable. When compared between the S14 and PP04 measurements, the metallicity gradients derived with S14 calibrators are typically flatter than those derived with PP04 calibrators.

In further analysis, we use the metallicity gradients derived from PP04 in radial binning because we can avoid the bias from poor sampling in short major axis cases. The PP04 calibrator is chosen so that we can easily compare with other observed gradients from previous observations.

## 6. ACTIVE GALACTIC NUCLEUS (AGN) CONTAMINATION AND BIAS FROM LOW SIGNAL-TO-NOISE DATA

We now examine whether our sample is affected by line emission from AGNs. Since our metallicity calibrators are only

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Verifying the validity of our derived metal gradients for low signal/noise $[N\,II]$ data. Metallicity gradient residuals are shown for CSWA128 as a function of the amount of synthetic noise added (see the text in Section 6 for technical details). Gradients are derived with PP04 method both along the galaxy “major axis” (purple) and using radial binning (green).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution of metallicity gradient with redshift. Red squares represent lensed galaxies from the present study (L15). Green squares are the four lensed galaxies in Jones et al. (2013, J13). Other data points are measured from a lensed starburst dwarf (Jones et al. 2015, J15), a lensed galaxy (Yuan et al. 2011, Y11), non-lensed galaxies observed with adaptive optics (Swinbank et al. 2012, S12), and Milky Way planetary nebulae (Maciel et al. 2003, M03). Open and filled symbols represent kinematically well-ordered and disturbed systems, respectively. Turquoise lines show predictions from the hydrodynamical simulations of Gibson et al. (2013, G13) emphasizing the strong sensitivity to the incorporated feedback.</figcaption></figure>

applicable to bounded $H\,II$ regions ionized by hot stars, any emission from nuclear activity would give misleading results for central metallicities and gradients. Furthermore the presence or absence of AGNs is interesting in terms of understanding the relationship between growth of galaxies and their central supermassive black holes. Optical spectra of AGNs are characterized by high ratios of collisional to Balmer lines and by broad line widths arising from outflows and/or broad line region kinematics. These signatures have been confirmed in the nuclear regions of many galaxies at $z \simeq 1-2$ with spatially resolved spectroscopy (Wright et al. 2010; Förster Schreiber et al. 2014; Genzel et al. 2014; Newman et al. 2014). Notably, Genzel et al. (2014) show that the fraction of galaxies showing signatures of AGNs is strongly correlated with stellar mass.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Correlation between the observed metallicity gradient for the present sample (L15) (Jones et al. 2013, J13 and Swinbank et al. 2012, S12) and the integrated [N II]/Hα ratio (left panel) and the star formation rate (right panel). Filled symbols represent isolated galaxies; open symbols represent interacting and/or kinematically immature galaxies.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Absence of a correlation between the metallicity gradient and the degree of rotational support. Data points represent the present sample and the gray scale represents the output of the *Illustris* simulation for a sample of $z \simeq 2$ galaxies (Vogelsberger et al. 2014; Nelson et al. 2015).</figcaption></figure>

For each of the eight new galaxies in this paper with [O III]/Hβ information, we construct a “BPT”:[N II]/Hα and [O III]/Hβ diagram from the spatially resolved spectra. Most show a distribution of values entirely consistent with star-forming regions as seen in other $z \sim 2$ galaxies (Steidel et al. 2014; Shapley et al. 2015). We thus only show two examples where the interpretation is less clear (Figure 4). CSWA15 shows a central region whose strong line ratios significantly lie at an average of $0.16 \pm 0.06$ dex above the maximum starburst classification line (Kewley et al. 2001) as shown in the upper panel of Figure 4. However, the central [N II]/Hα = $0.03 \pm 0.01$ and [S II]λλ6716, 6731/Hα = $0.13 \pm 0.03$ are both low and atypical of AGNs. These ratios are at the extreme low end probed by current $z \simeq 2$ samples (e.g., Steidel et al. 2014) and broadly consistent with an extension of the star-forming locus seen in these surveys. Alternatively the line ratios are consistent with theoretical expectations for an AGN with a gas metallicity of $\simeq 0.5$ solar (Groves et al. 2006), but the surrounding (and clearly star-forming) regions of CSWA15 indicate significantly lower metallicity. If the central line ratios are affected by AGNs then the gas-phase metallicity gradient would be significantly steeper than we infer assuming stellar ionization.

As an additional check for nuclear activity, we fit Hα from the central regions of CSWA15 with a double Gaussian profile. The best-fit broad component has a line width $\sigma = 153$ km s$^{-1}$ and comprises 29% of the total flux (c.f. $\sigma = 59$ km s$^{-1}$ for the dominant narrow component). This is similar to stacked spectra of outer disk regions in $z \simeq 2$ galaxies ($\sigma \sim 200$ km s$^{-1}$ and $\sim 40\%$ of flux in the broad component; Newman et al. 2013; Genzel et al. 2014) and is not indicative of AGN-driven outflows ($\sigma \gtrsim 500$ km s$^{-1}$) or broad-line regions ($\sigma \gtrsim 1000$ km s$^{-1}$). The spectrum of CSWA15 is therefore consistent with expectations for emission from star formation and associated outflows, and does not show strong evidence of an AGN.

The incidence of AGNs in our sample is consistent with expectations given their stellar masses. We find no clear signatures of AGNs and only one possible case discussed above out of the 12 galaxies with available BPT diagnostics, including those from Jones et al. (2013). An example of typical BPT diagrams that are compliant with star-forming regions in our sample, CSWA165, is shown in the lower panel of Figure 4. The remaining galaxies in our sample have low central [N II]/Hα and [S II]/Hα ratios indicating no evidence for AGNs. Genzel et al. (2014) find similar results using the same methods: zero secure and only two potential AGNs in a sample of 17 $z \simeq 1-2$ galaxies with stellar mass $\log M_*/M_{\odot} = 9.4-10.3$, comparable to that of our sample. They likewise find an AGN fraction <10% at $\log M_*/M_{\odot} < 10.3$ based on independent X-ray, infrared, and radio indicators. Therefore the absence of AGNs in our sample is not unexpected given the stellar mass range probed.

We then consider the possible bias in measured metallicity gradients due to possibly low signal-to-noise ratio of [N II] emission lines. We select CSWA128 to be a fiducial sample in this analysis due to its moderately steep gradient so that we can examine possible degradation of this gradient as we increase the noise. We adopt the $2\sigma$ detection limit of the [N II] emission line for CSWA128; the limits for other galaxies are similar. We then synthetically add Gaussian noise to the data cube up to twice the detection limit and re-measure the metallicity gradient with the same procedure. We compare the measured gradients when the noise is added to the data cube to the gradients measured from the original data cube in Figure 5. We find that the gradients measured from the noise-added data cubes are
<figure><img src="image_12.png" /></figure>

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

<figure><img src="image_14.png" /><figcaption>**Figure 9.** From left to right: source plane [N II]/Hα map, radial metallicity gradient from N2 calibrators, and radial metallicity gradient from O3N2 calibrators for a subsample of our targets. The gray dots are measurements from each source-plane pixel. Green points and lines are measurements from radial binning. Purple points and lines are measurements along “major axis.”</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 9. (Continued.)</figcaption></figure>

consistent with the original measurement up to doubling the noise or a quarter of the exposure time for CSWA128. Accordingly we can be confident that the measured gradients in other galaxies are not likely to be affected given our [N II] detection limits.

## 7. DISCUSSION

Previous observations (Yuan et al. 2011; Jones et al. 2013, 2015) led to a simple picture whereby isolated galaxies at high redshifts tend to have steep metallicity gradients that gradually flatten with cosmic time, while interacting galaxies have no discernible gradients at all. Conceptually this could be understood if metal gradients only became properly established when galaxies are kinematically well-ordered, with gradients flattening as galaxies grow in size. However, the enlarged data set presented here shows a large diversity of gradients, which is incompatible with this picture. The redshift-dependent behavior of metallicity gradients in our sample is summarized in Figure 6. As in previous work we find (negatively) steeper gradients in rotationally supported systems compared to the merging and dynamically immature systems. However, the gradients in rotationally supported galaxies are much flatter (>-0.1 dex kpc⁻¹) than those previously observed

in Jones et al. (2013) (~-0.3 dex kpc⁻¹), indicating less evolution in the gradient slope at z < 2. We now seek to explain the diversity in observed gradient slopes especially among isolated rotating galaxies.

One possibility for the difference seen between this work and that of Jones et al. (2013) could be procedural. In the present sample, we avoided a possible selection bias toward unusually active or metal-rich sources. In Jones et al. (2013), integrated [N II] and Hα fluxes were first measured for a larger sample with a slit spectrograph and only those four sources with relatively bright emission lines (particularly [N II]) were selected for subsequent study with OSIRIS. Following upgrades to the OSIRIS grating and Keck I adaptive optics system, this pre-screening step was avoided in the current sample (with the exception of CSWA31). To test for this bias, the left panel of Figure 7 correlates the metallicity gradient with the integrated [N II]/Hα line ratio for both isolated or rotationally supported and kinematically disturbed systems. All galaxies with [N II]/Hα < 0.1 have flat gradients (<0.1 dex kpc⁻¹). Moreover, for galaxies with [N II]/Hα > 0.1 there is a divergence between isolated or rotationally supported and merging or dynamically immature systems with the latter showing zero gradients irrespective of the metallicity. The right panel of Figure 7 examines the dependence on the integrated
<figure><img src="image_16.png" /><figcaption>Figure 9. (Continued.)</figcaption></figure>

SFR to check whether the Jones et al. (2013) sample is biased to more intense star-forming systems; however, no strong trend is revealed.

A more likely possibility for the diversity of gradients at $z \simeq 2-2.5$ is variation in feedback strength. Numerical simulations suggest that metallicity gradients are highly sensitive to feedback in the form of outflows and “galactic fountains” (Pilkington et al. 2012; Gibson et al. 2013; Anglés-Alcázar et al. 2014) in the sense that stronger feedback (i.e., higher mass loading factors) leads to more gas mixing and therefore flatter gradients. This is illustrated in Figure 6 which shows the evolution of the metallicity gradient in a simulated galaxy using two different feedback prescriptions. Simple analytical chemical evolution models show that stronger feedback also results in lower gas-phase metallicity (e.g., Jones et al. 2013). The correlation of gradient slope with

integrated metallicity, as shown in the left panel of Figure 7, could therefore also be due to a variation in feedback strength. In fact, the enhanced feedback scheme in Figure 6 was preferred over the normal feedback scheme considering the stellar-mass-halo-mass relationship matching (Stinson et al. 2013).

While Figure 7 provides evidence that gradient slopes are affected to different degrees by feedback, the trend has a large scatter suggesting that other factors are likely important. One obvious possibility is the degree of rotational support compared to random motions which would mix the gas and flatten any gradients. This effect is clearly seen in the sense that merging and interacting systems have flatter gradients compared to isolated galaxies. Among isolated galaxies, we would naively expect those with higher ratios of rotational to random motion (i.e., higher $v/\sigma$) to have stronger gradients in metallicity and perhaps other properties. However, Figure 8 shows very little
dependence on $v/\sigma$. We compare our observational results with $z=2$ galaxies drawn from the *Illustris* simulation (Vogelsberger et al. 2014; Nelson et al. 2015)<fnref n="8" />, a large-volume cosmological hydrodynamical simulation with baryonic feedback. Metallicity gradients in the simulation data were constructed based on a linear fit to all of the star-forming gas within a radius of 10 kpc from each galaxy’s center. The simulated galaxies have gas-phase kinematics and negative metallicity gradients which are similar to our observations and likewise show no clear trend between gradient slope and rotational support. In essence, galaxies with gas disks supported by rotation have radial metallicity gradients which are similar to systems dominated by random motions, in contrast to expectations that bulk velocity dispersion should flatten any gradients. This suggests that gradients are established on timescales of the order of a dynamical time, sufficiently short that prominent random motions are not effective at mixing the gas.

To summarize, we find a large diversity of metallicity gradients among rotating galaxies at $z \simeq 2$. Most rotating galaxies in our sample have negative gradients, while interacting systems have flatter gradients, as is found in local galaxy samples. Surprisingly we find that the degree of ordered rotation versus random bulk motion has no discernible effect on the gradient slope even among isolated galaxies, indicating that gradients are formed and destroyed on timescales comparable to the galaxy dynamical time (~$10^8-10^9$ years). The steepest observed gradients with slopes $\sim -0.3$ dex kpc$^{-1}$ are less common than found in our previous work suggesting that this is a relatively rare or short-lived phase at high redshift. It would therefore be interesting to verify the time-dependence of steep gradients in the numerical simulations. Finally, we find a weak correlation between galaxy-integrated metallicity and gradient slope which we interpret as an effect of feedback: the same feedback that reduces the gas-phase metallicity evidently also flattens the gradients, as found in cosmological zoom-in simulations. The diversity of metallicity gradients at $z \simeq 2$ is therefore likely caused in part by differences in the recent feedback history. As we argued in Jones et al. (2013), the evolution of metallicity gradients is evidently a sensitive probe of feedback in the form of gaseous outflows and galactic fountains.

## 8. SUMMARY

We present spatially resolved kinematic and gas-phase metallicity gradient measurements for a total of 15 gravitationally lensed star-forming galaxies at $z \sim 2$ based on the analysis of strong emission lines of Hα, [N ii], [O iii], and Hβ. Eleven new sources were observed with the Keck laser-assisted adaptive optics systems and the upgraded OSIRIS IFU spectrograph to which we add earlier data on four sources presented in Jones et al. (2013). With the aid of gravitational lensing, the typical source plane spatial resolution for each source is <500 pc which is considerably better than for other systems studied at similar redshifts.

We found the following key results:

1. High spatial resolution is crucial in diagnosing the kinematic properties and dynamical maturity of $z \simeq 2$ galaxies. We compare our observed velocity map of

interacting pairs of merging galaxies with those simulated at the lower resolution typical in non-lensed surveys and demonstrate how easily such sources can be mistaken to represent rotationally supported disks in poorer sampled data. Even for galaxies with no clear morphological sign of interaction, we often find significant large deviations in the velocity field compared to rotating disk models which would not be apparent in seeing-limited data. As a result we observe a significantly lower fraction (36%) of rotationally supported systems in our sample than has been claimed ($\simeq 74\%$) from larger kinematic surveys undertaken with lower spatial resolution.

2. We find a much higher fraction of $z \simeq 2$ galaxies have weak or non-existent metallicity gradients than in previous studies of smaller samples observed in the same redshift range. It seems unlikely that such a change arises as the result of a bias in the earlier sample which pre-selected sources with stronger [N ii] lines and hence increased metallicity. We observe only a weak correlation between the presence of a gradient and the metallicity and none with the degree of rotational support, the latter being consistent with the recent predictions of the *Illustris* hydrodynamical simulations. We argue that variations in gas and metal mixing due to feedback most likely play the dominant role in modifying metal gradients and thus can explain the sizable scatter we see in our enlarged sample. The sensitivity of the observed metal gradient to the various modes of feedback indicates it will remain a promising tool for understanding galaxy assembly.

The authors recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. We acknowledge useful discussions with Chuck Steidel, Phil Hopkins, Xiangcheng Ma, Drew Newman, and Paul Torrey, and thank Paul Torrey for his *Illustris* metal gradient predictions tailored to the context of our observations. Support for A.Z. was provided by NASA through Hubble Fellowship grant \#HST-HF2-51334.001-A awarded by STScI. T.A.J. acknowledges support from NASA through Hubble Fellowship grant HST-HF2-51359.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

## REFERENCES

Anglés-Alcázar, D., Davé, R., Özel, F., & Oppenheimer, B. D. 2014, ApJ, 782, 84

Auger, M. W., Treu, T., Brewer, B. J., & Marshall, P. J. 2011, MNRAS, 411, L6

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Belli, S., Jones, T., Ellis, R. S., & Richard, J. 2013, ApJ, 772, 141

Bresolin, F. 2011, ApJ, 730, 129

Brewer, B. J., Lewis, G. F., Belokurov, V., et al. 2011, MNRAS, 412, 2521

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Cresci, G., Mannucci, F., Maiolino, R., et al. 2010, Natur, 467, 811

Dahle, H., Gladders, M. D., Sharon, K., et al. 2013, ApJ, 773, 146

Elíasdóttir, Á., Limousin, M., Richard, J., et al. 2007, arXiv:0710.5636

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

Förster Schreiber, N. M., Genzel, R., Lehnert, M. D., et al. 2006, ApJ, 645, 1062

Förster Schreiber, N. M., Genzel, R., Newman, S. F., et al. 2014, ApJ, 787, 38

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

<footnote n="8">⁸ *Illustris* data are available through http://www.illustris-project.org.</footnote>
Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

Genzel, R., Förster Schreiber, N. M., Rosario, D., et al. 2014, ApJ, 796, 7

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

Genzel, R., Tacconi, L. J., Eisenhauer, F., et al. 2006, Natur, 442, 786

Gibson, B. K., Pilkington, K., Brook, C. B., Stinson, G. S., & Bailin, J. 2013, A&A, 554, A47

Groves, B. A., Heckman, T. M., & Kauffmann, G. 2006, MNRAS, 371, 1559

Heckman, T. M., Sembach, K. R., Meurer, G. R., et al. 2001, ApJ, 554, 1021

Jones, T., Ellis, R., Jullo, E., & Richard, J. 2010a, ApJL, 725, L176

Jones, T., Ellis, R. S., Richard, J., & Jullo, E. 2013, ApJ, 765, 48

Jones, T., Wang, X., Schmidt, K. B., et al. 2015, AJ, 149, 107

Jones, T. A., Swinbank, A. M., Ellis, R. S., Richard, J., & Stark, D. P. 2010b, MNRAS, 404, 1247

Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, NJPh, 9, 447

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

Larkin, J., Barczys, M., Krabbe, A., et al. 2006, Proc. SPIE, 6269, 62691

Law, D. R., Steidel, C. C., Erb, D. K., et al. 2009, ApJ, 697, 2057

Limousin, M., Ebeling, H., Richard, J., et al. 2012, A&A, 544, A71

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

Livermore, R. C., Jones, T. A., Richard, J., et al. 2015, MNRAS, 450, 1812

Maciel, W. J., Costa, R. D. D., & Uchida, M. M. M. 2003, A&A, 397, 667

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463

Mieda, E., Wright, S. A., Larkin, J. E., et al. 2014, PASP, 126, 250

Nelson, D., Pillepich, A., Genel, S., et al. 2015, A&C, 13, 12

Newman, S. F., Buschkamp, P., Genzel, R., et al. 2014, ApJ, 781, 21

Newman, S. F., Genzel, R., Förster Schreiber, N. M., et al. 2013, ApJ, 767, 104

Noeske, K. G., Faber, S. M., Weiner, B. J., et al. 2007, ApJL, 660, L47

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59

Pilkington, K., Few, C. G., Gibson, B. K., et al. 2012, A&A, 540, A56

Queyrel, J., Contini, T., Kissler-Patig, M., et al. 2012, A&A, 539, A93

Richard, J., Jones, T., Ellis, R., et al. 2011, MNRAS, 413, 643

Richard, J., Smith, G. P., Kneib, J.-P., et al. 2010, MNRAS, 404, 325

Shapley, A. E., Reddy, N. A., Kriek, M., et al. 2015, ApJ, 801, 88

Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65

Sharples, R., Bender, R., Agudo Berbel, A., et al. 2013, Msngr, 151, 21

Stark, D. P., Auger, M., Belokurov, V., et al. 2013, MNRAS, 436, 1040

Stark, D. P., Swinbank, A. M., Ellis, R. S., et al. 2008, Natur, 455, 775

Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165

Stinson, G. S., Brook, C., Macciò, A. V., et al. 2013, MNRAS, 428, 129

Stott, J. P., Sobral, D., Swinbank, A. M., et al. 2014, MNRAS, 443, 2695

Swinbank, A. M., Sobral, D., Smail, I., et al. 2012, MNRAS, 426, 935

Troncoso, P., Maiolino, R., Sommariva, V., et al. 2014, A&A, 563, A58

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Natur, 509, 177

Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209

Wisnioski, E., Glazebrook, K., Blake, C., et al. 2011, MNRAS, 417, 2601

Wizinowich, P. L., Le Mignant, D., Bouchez, A. H., et al. 2006, PASP, 118, 297

Wright, S. A., Larkin, J. E., Graham, J. R., & Ma, C.-P. 2010, ApJ, 711, 1291

Yang, C.-C., & Krumholz, M. 2012, ApJ, 758, 48

Yuan, T.-T., Kewley, L. J., & Rich, J. 2013, ApJ, 767, 106

Yuan, T.-T., Kewley, L. J., Swinbank, A. M., Richard, J., & Livermore, R. C. 2011, ApJL, 732, L14

Zitrin, A., Fabris, A., Merten, J., et al. 2015, ApJ, 801, 44