The Grism Lens-Amplified Survey from Space (GLASS). X. Sub-kiloparsec Resolution Gas-phase Metallicity Maps at Cosmic Noon behind the Hubble Frontier Fields Cluster MACS1149.6+2223

2.1.1. The Grism Lens-Amplified Survey from Space

The Grism Lens-Amplified Survey from Space²¹ (GLASS; Proposal ID 13459; P.I. Treu, Schmidt et al. 2014; Treu et al. 2015) is an HST cycle 21 large general observing (GO) program. GLASS observed 10 massive galaxy clusters with the Wide Field Camera 3 Infrared (WFC3/IR) grisms (G102 and G141; 10 and 4 orbits per cluster, respectively) targeted at their centers and the Advanced Camera for Survey (ACS) Optical grism (G800L) at their infall regions. The exposure on each cluster core is split into two nearly orthogonal position angles (PAs), in order to disentangle contamination from neighboring objects. Data acquisition was completed in 2015 January. Here we focus on the grism data targeted on the center of MACS1149.6+2223 taken in 2014 February (PA = 032^◦) and 2014 November (PA = 125^◦), marked by blue squares in Figure 1 (also see Table 1). GLASS provides ∼10% of the G141 exposures used in this work, and 100% of the G102 exposures. Details on GLASS data reduction can be found in Schmidt et al. (2014) and Treu et al. (2015).

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Grism Spectroscopic Data Used in This Work

PAa	Grism	Exposure Time	Program	Time of Completion
(deg.)		(s)
032	G102	8723	GLASS	2014 Feb
	G141	4412	GLASS
111	G141	36088	SN Refsdal follow-up	2014 Dec
119	G141	36088	SN Refsdal follow-up	2015 Jan
125	G102	8623	GLASS	2014 Nov
	G141	4412	GLASS

Note. Here we only include the grism observations targetted on the prime field of MACS1149.6+2223.

^aThe position angle shown here corresponds to the "PA_V3" value reported in the WFC3/IR image headers. The position angle of the dispersion axis of the grism spectra is given by ${\mathrm{PA}}_{\mathrm{disp}}\approx \mathrm{PA}\_{\rm{V}}3-45$.

Download table as:  ASCIITypeset image

Shallow images through filters F105W or F140W were taken to aid the alignment and extraction of grism spectra. The imaging exposures are combined with the exposures obtained by the Hubble Frontier Fields (HFF) initiative and other programs to produce the deep stacks, released as part of the HFF program (see Section 3 for more details).

2.1.2. The Supernova Refsdal Follow-up Program with HST G141

2.1.3. Grism Data Reduction

The selection of our metallicity gradient sample is based upon the master redshift catalog for²³ MACS1149.6+2223, published by the GLASS collaboration. As described by Treu et al. (2016), redshifts were determined by combining spectroscopic information from GLASS and the SN Refsdal follow-up HST grism programs, ground-based MUSE observations and Keck DEIMOS data.

For the purpose of this study, we compiled an exhaustive list of spatially extended sources with secure spectroscopic redshift measurements in the range of $z\in [1.2,2.3]$, showing strong ELs (primarily $[{\rm{O}}\,{\rm{III}}]$ and ${\rm{H}}\beta$) in the 2D grism spectra. The reason this redshift range is selected is that we require high signal-to-noise ratio (S/N) detections of $[{\rm{O}}\,{\rm{III}}]$, ${\rm{H}}\beta$, and $[{\rm{O}}\,{\rm{II}}]$ ELs in the spectral region where grism sensitivity and throughput do not drop significantly, in order to deliver reliable metallicity estimates from well-calibrated EL diagnostics.²⁴ The selection results in a sample of 14 galaxies. The positions of these 14 galaxies relative to the cluster are denoted by the magenta circles in Figure 1, while the postage stamp images of these galaxy are shown in Figure 2. The full sample is also listed in Table 2. In particular, one source in our sample, i.e., galaxy ID 04054, is one of the multiple images of SN Refsdal's host galaxy (see Figure 3 for its appearances from various perspectives). A steep metallicity gradient (−0.16 ± 0.02 dex kpc⁻¹) has been previously measured on another multiple image of this spiral galaxy (Yuan et al. 2011). Our work is the first metallicity gradient measurement on this particular image, which is least distorted/contaminated by foreground cluster members. We also take advantage of a more precise lens model of both the entire cluster and the SN Refsdal host, from the extensive lens modeling effort summarized by Treu et al. (2016).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Global Demographic Properties of the Metallicity Gradient Sample Analyzed in This Work

ID	R.A. (deg.)	Decl. (deg.)	zspec	μa	${H}_{160}$ Magnitude	SED Fitting		EL Diagnostics
					(ABmag)	log(${M}_{* }$ b/${M}_{\odot }$)	${A}_{{\rm{V}}}^{{\rm{S}}}$ c	$12+\mathrm{log}({\rm{O}}/{\rm{H}})$	SFRb [${M}_{\odot }$ yr−1]	${A}_{{\rm{V}}}^{{\rm{N}}}$ c
01422	177.398643	22.387499	2.28	2.35 [2.09, 2.26]	24.46	${9.29}_{-0.01}^{+0.07}$	<0.01	${8.26}_{-0.13}^{+0.11}$	${15.10}_{-7.11}^{+16.76}$	$\lt 1.07$
02389	177.406546	22.392860	1.89	43.25 [37.47, 50.28]	23.29	${8.43}_{-0.07}^{+0.06}$	<0.01	${7.88}_{-0.15}^{+0.16}$	$\lt 6.37$	${1.36}_{-0.82}^{+0.69}$
02607d	177.413454	22.393843	1.86	3.01 [2.79, 2.98]	22.63	${10.25}_{-0.10}^{+0.06}$	${1.70}_{-0.18}^{+0.42}$	${7.70}_{-0.11}^{+0.13}$	$\gt 36.51$	${2.26}_{-0.58}^{+0.31}$
02806	177.392541	22.394921	1.50	4.82 [3.11, 3.51]	23.35	${8.72}_{-0.01}^{+0.19}$	${0.50}_{-0.01}^{+0.01}$	${7.96}_{-0.16}^{+0.15}$	${1.96}_{-0.13}^{+1.37}$	$\lt 0.18$
02918	177.412652	22.395723	1.78	2.95 [2.72, 2.91]	25.22	${8.37}_{-0.08}^{+0.15}$	${0.10}_{-0.10}^{+0.34}$	${8.42}_{-0.14}^{+0.11}$	$\lt 70.32$	${1.25}_{-0.80}^{+0.72}$
03746	177.391848	22.400105	1.25	3.78 [3.56, 3.83]	22.48	${8.81}_{-0.01}^{+0.03}$	<0.01	${8.11}_{-0.11}^{+0.10}$	${10.95}_{-0.58}^{+1.61}$	$\lt 0.17$
04054	177.403393	22.402456	1.49	3.39 [3.04, 3.28]	21.27	${9.64}_{-0.01}^{+0.05}$	${1.10}_{-0.01}^{+0.01}$	${8.70}_{-0.11}^{+0.09}$	${16.99}_{-2.80}^{+5.71}$	${0.72}_{-0.23}^{+0.23}$
05422	177.382077	22.407510	1.97	3.63 [3.11, 4.41]	25.05	${8.45}_{-0.46}^{+0.15}$	${0.00}_{-0.01}^{+0.56}$	${7.56}_{-0.12}^{+0.16}$	$\lt 38.18$	$\lt 1.85$
05709	177.397234	22.406181	1.68	7.10 [6.74, 7.53]	25.44	${7.90}_{-0.03}^{+0.02}$	<0.01	${8.21}_{-0.21}^{+0.13}$	$\lt 17.97$	$\lt 1.18$
05732	177.415126	22.406195	1.68	1.56 [1.51, 1.61]	23.45	${9.09}_{-0.01}^{+0.01}$	<0.01	${8.41}_{-0.12}^{+0.10}$	$\lt 84.70$	$\lt 1.31$
05811	177.389220	22.407583	2.31	7.50 [6.88, 9.34]	23.79	${8.85}_{-0.10}^{+0.11}$	<0.01	${8.16}_{-0.21}^{+0.14}$	${4.94}_{-2.47}^{+3.13}$	$\lt 0.74$
05968	177.406922	22.407499	1.48	1.84 [1.78, 1.96]	22.24	${9.34}_{-0.03}^{+0.02}$	${0.60}_{-0.01}^{+0.20}$	${8.47}_{-0.08}^{+0.07}$	${27.95}_{-4.41}^{+4.70}$	${0.23}_{-0.13}^{+0.15}$
07058	177.405976	22.412977	1.79	2.13 [1.97, 3.01]	24.28	${8.74}_{-0.15}^{+0.17}$	${0.00}_{-0.01}^{+0.12}$	${8.38}_{-0.19}^{+0.13}$	$\lt 44.99$	$\lt 1.85$
07255	177.385990	22.414074	1.27	2.34 [2.31, 2.99]	22.47	${9.36}_{-0.21}^{+0.01}$	${0.30}_{-0.01}^{+0.21}$	${8.38}_{-0.08}^{+0.08}$	${8.73}_{-3.10}^{+3.53}$	$\lt 0.79$

Notes. The error bars and upper/lower limits shown in the columns of SED fitting and EL diagnostics correspond to 1σ confidence ranges. Note that the errors on the SED fitting results do not include any systematic uncertainties associated with the Bruzual & Charlot (2003) stellar population models.

^aBest-fit magnification values and 1σ confidence intervals. Except for galaxy ID 02389, the magnification results are from the Glafic version 3 model, calculated by the HFF interactive online magnification calculator available at https://archive.stsci.edu/prepds/frontier/lensmodels/webtool/magnif.html. For galaxy ID 02389, we use the Sharon & Johnson version 3 model instead to compute the magnification results and correct for lensing magnification. ^bValues presented here are corrected for lensing magnification. ^cThe superscripts of "S" and "N" refer to the stellar and nebular V-band dust extinction in units of magnitude, respectively. ^dThe EL diagnostic result on this source is not trustworthy, since it is classified as an AGN candidate (see Section 5).

Download table as:  ASCIITypeset image

The full-depth seven-filter HFF photometry was fitted with the template spectra from Bruzual & Charlot (2003) using the stellar population synthesis code FAST (Kriek et al. 2009), in order to derive global estimates of stellar mass (${M}_{* }$) and the dust extinction of stellar continuum (${A}_{{\rm{V}}}^{{\rm{S}}}$) for each source in our sample. We take a grid of stellar population parameters that include exponentially declining star-formation histories with e-folding times ranging from 10 Myr to 10 Gyr, stellar ages ranging from 5 Myr to the age of the universe at the observed redshift, and ${A}_{{\rm{V}}}^{{\rm{S}}}=0\mbox{--}4$ mag for a Calzetti et al. (2000) dust extinction law. We assume the Chabrier (2003) initial mass function (IMF) and fix the stellar metallicity to solar. Table 2 lists the best-fit ${M}_{* }$ values, corrected for lensing magnification. Although the adopted solar metallicity is higher than the typical gas-phase abundances we infer for the sample, this has little effect on the derived stellar mass. Fixing the stellar metallicity to 1/5 solar (i.e., Z = 0.004 or $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.0$, corresponding to the sample median) reduces the best-fit ${M}_{* }$ by only ∼0.03 dex.

Measurements of metallicity gradients require knowledge of the galaxy center, major axis orientation, and inclination or axis ratio. We derive these quantities from spatially resolved maps of the stellar mass surface density following the methodology described in Jones et al. (2015b). Briefly, we smooth the HFF photometric images to a common PSF of 02 FWHM and fit the spectral energy distribution (SED) of each spaxel using the same procedure described in Section 4.1. The resulting ${M}_{* }$ maps in the image plane are shown in Figures 3 and 4. We reconstruct the ${M}_{* }$ maps in the source plane using the adopted lens models and determine the centroid, axis ratio, inclination, and major axis orientation from a 2D elliptical Gaussian fit to the ${M}_{* }$ distribution. This allows us to measure the galactocentric radius at each point, assuming that contours of ${M}_{* }$ surface density correspond to constant de-projected radius (following Jones et al. 2015b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the following analysis is mostly done with the robust image plane data, we reconstruct the derived radius maps in the image plane. The contours of constant source plane galactocentric radius for each galaxy in our sample are shown in Figures 3 and 4. Note that these contours are distorted from ellipses ascribed to the shear of gravitational lensing. For nearly half of our sample where we have seeing-limited IFU data described in Appendix A, we verify that for galaxies with significant rotation support (e.g., IDs 05709 and 05968), the morphological major axis agrees at 2σ with the pseudo-slit that maximizes velocity shear according to our kinematic analysis.

The broad grism wavelength coverage provides spatially resolved maps of multiple ELs, such as $[{\rm{O}}\,{\rm{II}}]$, ${\rm{H}}\gamma$, ${\rm{H}}\beta$, $[{\rm{O}}\,{\rm{III}}]$, ${\rm{H}}\alpha$+$[{\rm{N}}\,{\rm{II}}]$, and $[{\rm{S}}\,{\rm{II}}]$. The integrated fluxes of these ELs for our sample of galaxies are presented in Appendix B. To obtain pure EL maps, we first use the aXe continuum models described in Section 2 to subtract contaminating flux from neighboring sources. The continuum model for the target object is scaled to match the locally observed levels, and then subtracted. We check that the flux residual in regions near the ELs of interest follow a normal distribution with zero mean.

Due to the limited spectral resolution of HST WFC3/IR grisms, ${\rm{H}}\beta$ and the $[{\rm{O}}\,{\rm{III}}]$ λλ4960,5008 doublets are partially blended in the spatially extended sources. We adopt a direct de-blending technique to separate these emission lines following the procedure used in Jones et al. (2015b). First, an isophotal contour is measured from the co-added HST ${H}_{160}$-band postage stamp, typically at the level of ∼10% of the peak flux received from the source. We then apply this contour to the 2D grism spectra at the observed wavelengths corresponding to $[{\rm{O}}\,{\rm{III}}]$ λλ4960,5008 and ${\rm{H}}\beta$, maximizing the flux given grism redshift uncertainty (∼0.01). We use the theoretical $[{\rm{O}}\,{\rm{III}}]$ λλ4960,5008 doublet flux ratio (i.e., $[{\rm{O}}\,{\rm{III}}]$ λ5008/$[{\rm{O}}\,{\rm{III}}]$ λ4960 = 2.98, calculated by Storey & Zeippen 2000) to subtract the flux contribution of $[{\rm{O}}\,{\rm{III}}]$ λ4960 corresponding to the region where $[{\rm{O}}\,{\rm{III}}]$ λ5008 is unblended. We iterate this procedure to remove the $[{\rm{O}}\,{\rm{III}}]$ λ4960 emission completely, resulting in pure maps of ${\rm{H}}\beta$ and $[{\rm{O}}\,{\rm{III}}]$ λ5008.

Note that this direct de-blending will be compromised if the source is so extended that its $[{\rm{O}}\,{\rm{III}}]$ λ5008 and ${\rm{H}}\beta$ are blended. In practice, we fine-tune the isophotal contour level to avoid this contamination. As a result, no cases in our sample have problems with our direct de-blending method, and the resulting pure $[{\rm{O}}\,{\rm{III}}]$ λ5008 and ${\rm{H}}\beta$ maps are uncontaminated.

In the case of ${\rm{H}}\alpha$ and $[{\rm{N}}\,{\rm{II}}]$, the emission lines are separated by less than the grism spectral resolution. Therefore, we cannot use the direct technique to de-blend ${\rm{H}}\alpha$ and $[{\rm{N}}\,{\rm{II}}]$ emission, and instead treat the $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ ratio as an unknown parameter in Section 4.4.

After obtaining pure EL maps, we rotate the maps from each PA to the same orientation and then combine the data, in order to utilize the full depth of the grism exposures. We note that EL maps must first be extracted from grism spectra at each PA before being combined, in order to account for the different rotation and spectral resolution. There are, in total, six different PA-grism combinations as given in Table 1. At certain redshifts, ELs are visible in the overlapping region between the two grisms at λ ∈ [1.08, 1.15] μm resulting in up to six individual maps of the same EL. We use the data taken at PA = 119^◦ as a reference, as this minimizes errors arising from rotating the various PAs into alignment. We first apply a small vertical offset to each EL map to account for the slight tilt of the spectral trace (i.e., known offset between the dispersion direction and the x-axis of extracted spectra; Brammer et al. 2012). Next, we apply a small horizontal shift to correct for uncertainties in the wavelength calibration and redshift. Finally, we rotate each EL map to align it with the reference orientation at PA = 119^◦. The best values of these three alignment parameters are found by maximizing the normalized cross-correlation coefficient (NCCF) in order to achieve optimal alignment of multiple PAs:

$$\begin{eqnarray}&&\mathrm{NCCF}=\displaystyle \frac{1}{n}{{\rm{\Sigma }}}_{x,y}\displaystyle \frac{({S}_{\mathrm{ref}}(x,y)-\langle {S}_{\mathrm{ref}}\rangle )({S}_{\mathrm{PA}}(x,y)-\langle {S}_{\mathrm{PA}}\rangle )}{{\sigma }_{{S}_{\mathrm{ref}}}{\sigma }_{{S}_{\mathrm{PA}}}}.\end{eqnarray} \tag{ 1 }$$

Here S is the surface brightness profile for a specific EL, ${S}_{\mathrm{ref}}$ corresponds to the reference alignment stamp adopted as the PA = 119^◦ data, and n is the number of spaxels in the surface brightness profile. $\langle S\rangle$ and ${\sigma }_{S}$ represent the average and standard deviation of the surface brightness, respectively.

Once the data from each PA are aligned, we vet the EL maps from each PA-grism combination and reject those that show significant contamination-subtraction residuals or grism reduction defects. We combine the remaining maps of a given EL using an inverse variance weighted average in order to properly account for the different exposure times and noise levels. The final combined maps of each EL for each galaxy are shown in Figures 3 and 4.

Generally speaking, two methods exist in deriving gas-phase oxygen abundance in star-forming galaxies at high redshifts. "Direct" determinations based on electron temperature measurements, which require high S/N detections of auroral lines (e.g., $[{\rm{O}}\,{\rm{III}}]$ λ4364), are still extremely challenging beyond the local universe (see Sanders et al. 2016, for a rare example). We therefore rely on the calibrations of strong collisionally excited EL flux ratios to estimate metallicities.

One of the most popular strong EL diagnostics is the flux ratio of $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ calibrated by Pettini & Pagel (2004). However, a large offset (0.2–0.4 dex) between the loci of high-z and present-day star-forming galaxies in the Baldwin–Phillips–Terlevich (BPT, Baldwin et al. 1981) diagram has recently been discovered, indicating that extending the locally calibrated $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ to high-z can be problematic (Sanders et al. 2015; Shapley et al. 2015). The interpretation of this offset is still a topic of much debate. It has been interpreted as the existence of a fundamental relation between the nitrogen–oxygen abundance ratio and ${M}_{* }$ (Masters et al. 2016), a systematic dependence of strong EL properties with respect to Balmer line luminosity (Cowie et al. 2016), a combination of harder ionizing radiation and higher ionization states of the ISM gas at high redshifts. (Steidel et al. 2014, 2016), or an enhancement of nitrogen abundance in hot H ii regions (Pilyugin et al. 2010). Fortunately, the alpha-element BPT diagrams show no offset with z (Sanders et al. 2015; Shapley et al. 2015), and therefore the diagnostics based upon those ELs are more reliable.

In this work, we adopt the calibrations of Maiolino et al. (2008, M08), including the flux ratios of $[{\rm{O}}\,{\rm{III}}]$/${\rm{H}}\beta$, $[{\rm{O}}\,{\rm{II}}]$/${\rm{H}}\beta$, and $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$. Given the potential systematics related to nitrogen ELs, we do not use the $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ calibration of M08. M08 fit the relation between these EL ratios and gas-phase oxygen abundance based upon "direct" measurements (from $[{\rm{O}}\,{\rm{III}}]$ λ4364) for $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ < 8.3, and photoionization model results for higher metallicities. This provides a continuous framework valid over the wide range of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ ∈ (7.1, 9.2). We note that although there are alternative calibrations of metallicity available in the literature and the applicability of these recipes is currently a hotly debated topic (see, e.g., Dopita et al. 2013; Blanc et al. 2015). Here we are primarily interested in relative metallicity measurements. Even though the absolute measurements of metallicity may change if we used a different calibration, the gradients and morphological features in the maps should be more robust to changes in the calibration. We will restrict our comparisons with other samples to only include studies assuming the M08 calibrations.

Unlike previous work in which calibrated relations are applied to various EL flux ratios (see, e.g., Pérez-Montero 2014; Bianco et al. 2015), we design a Bayesian statistical approach that directly uses the individual EL fluxes as input, such that the information from one EL is only used once. Our approach presents several advantages over those based on line flux ratios. First, it properly accounts for flux uncertainties for lines that are weak or undetected. Second, multiple lines can be taken into consideration without the risk of double counting any information.

Our approach simultaneously infers gas-phase metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$), nebular dust extinction (${A}_{{\rm{V}}}^{{\rm{N}}}$), and expected de-reddened ${\rm{H}}\beta$ flux in units of 10⁻¹⁷ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$. Details of the prior assumptions applied to these parameters are explained in Table 3. For values of these three parameters, we can compute the expected line fluxes, and compare with measured values to compute the likelihood and then the posterior probability. The ${\chi }^{2}$ statistic in our inference procedure is calculated as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\,\displaystyle \frac{{({f}_{{\mathrm{EL}}_{i}}-{R}_{i}\cdot {f}_{{\rm{H}}\beta })}^{2}}{{({\sigma }_{{\mathrm{EL}}_{i}})}^{2}+{f}_{{\rm{H}}\beta }^{2}\cdot {({\sigma }_{{R}_{i}})}^{2}}.\end{eqnarray} \tag{ 2 }$$

Here ${\mathrm{EL}}_{i}$ denotes the set of emission lines $[{\rm{O}}\,{\rm{II}}]$, ${\rm{H}}\gamma$, ${\rm{H}}\beta$, $[{\rm{O}}\,{\rm{III}}]$, ${\rm{H}}\alpha$, and $[{\rm{S}}\,{\rm{II}}]$. ${f}_{{\mathrm{EL}}_{i}}$ and ${\sigma }_{{\mathrm{EL}}_{i}}$ represent the observed de-reddened ${\mathrm{EL}}_{i}$ flux and its uncertainty given a value of ${A}_{{\rm{V}}}^{{\rm{N}}}$. We adopt the Cardelli et al. (1989) galactic extinction curve with ${R}_{{\rm{V}}}$ = 3.1. R_i is the dust-corrected flux ratio between ${\mathrm{EL}}_{i}$ and ${\rm{H}}\beta$, with ${\sigma }_{{R}_{i}}$ being the associated intrinsic scatter. Note that the content of R_i varies depending on the corresponding ${\mathrm{EL}}_{i}$ in the following ways.

Table 3.  Sampling Parameters and Their Prior Information

Set	Parameter	Symbol (Unit)	Prior
Vanilla	gas-phase metallicity	$12+\mathrm{log}({\rm{O}}/{\rm{H}})$	flat: $[7.0,9.3]$
	nebular dust extinction	${A}_{{\rm{V}}}^{{\rm{N}}}$	flat: $[0,4]$
	de-reddened ${\rm{H}}\beta$ flux	fHβ (10−17 $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$)	Jeffrey's: $[0,100]$
Extended	$[{\rm{N}}\,{\rm{II}}]$ to ${\rm{H}}\alpha$ flux ratio	$[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$	Jeffrey's: [0, 0.5]
Derived	star-formation rate	SFR (${M}_{\odot }$ yr−1)	⋯

Note. Most constraints presented in this work are obtained under the "Vanilla" parameter set with the extended parameter fixed as $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ = 0.05.

Download table as:  ASCIITypeset image

• 1.  
   For ${\mathrm{EL}}_{i}\in \{{\rm{H}}\alpha ,\,{\rm{H}}\gamma \}$, R_i corresponds to the Balmer decrement. We use the intrinsic Balmer ratios of ${\rm{H}}\alpha /{\rm{H}}\beta =2.86$ and ${\rm{H}}\gamma /{\rm{H}}\beta =0.47$, appropriate for case B recombination and fiducial H ii region conditions (see, e.g., Hummer & Storey 1987). We set ${\sigma }_{{R}_{i}}=0$, assuming that the intrinsic Balmer ratios are fixed with no scatter. These ratios are proxies of the nebular dust extinction ${A}_{{\rm{V}}}^{{\rm{N}}}$.
• 2.  
   For ${\mathrm{EL}}_{i}\in \{[{\rm{O}}\,{\rm{II}}],[{\rm{O}}\,{\rm{III}}]\}$, R_i and ${\sigma }_{{R}_{i}}$ are taken from the M08 calibrations. These ratios are diagnostics of the oxygen abundance $12+\mathrm{log}({\rm{O}}/{\rm{H}})$.
• 3.  
   For ${\mathrm{EL}}_{i}=[{\rm{S}}\,{\rm{II}}]$, we compute the values of R_i and ${\sigma }_{{R}_{i}}$ using the Balmer ratio of ${\rm{H}}\alpha$/${\rm{H}}\beta$ and our new calibration of the ratio $[{\rm{S}}\,{\rm{II}}]$/${\rm{H}}\alpha$, derived following Jones et al. (2015a). Jones et al. (2015a) used the same data as the low-metallicity calibrations of M08. Although M08 do not provide any calibrations for $[{\rm{S}}\,{\rm{II}}]$, we expect our calibration to be self-consistent and reliable for $12+\mathrm{log}({\rm{O}}/{\rm{H}})\lesssim 8.4$. $[{\rm{S}}\,{\rm{II}}]$/${\rm{H}}\alpha$ is primarily useful as a diagnostic to differentiate between the upper and lower branches of the $[{\rm{O}}\,{\rm{III}}]$/${\rm{H}}\beta$ calibration. This is valuable because $[{\rm{N}}\,{\rm{II}}]$ is not directly measured and $[{\rm{O}}\,{\rm{II}}]$ is typically of low signal-to-noise in the lower redshift galaxies for which $[{\rm{S}}\,{\rm{II}}]$ is available, due to decreasing grism throughput at $\lambda \lt 0.9$ μm.
• 4.  
   For ${\mathrm{EL}}_{i}={\rm{H}}\beta$, R_i = 1 with ${\sigma }_{{R}_{i}}=0$. This is just comparing the observed and expected ${\rm{H}}\beta$ flux, corrected for dust extinction.

These EL flux ratios can be computed given a value of $12+\mathrm{log}({\rm{O}}/{\rm{H}})$, using a universal polynomial functional form, i.e., $\mathrm{log}R={\sum }_{i=0}^{5}{c}_{i}\cdot {x}^{i}$, where $x=12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.69$. We summarize the coefficients (c_i) for each ratio diagnostic used in our statistical inference in Table 4.

Table 4.  Coefficients for the EL Flux Ratio Diagnostics Used in Our Bayesian Inference

R	c0	c1	c2	c3	c4	c5
${\rm{H}}\alpha$/${\rm{H}}\beta$	0.4564	⋯	⋯	⋯	⋯	⋯
${\rm{H}}\gamma$/${\rm{H}}\beta$	−0.3279	⋯	⋯	⋯	⋯	⋯
$[{\rm{O}}\,{\rm{III}}]$/${\rm{H}}\beta$	0.1549	−1.5031	−0.9790	−0.0297	⋯	⋯
$[{\rm{O}}\,{\rm{II}}]$/${\rm{H}}\beta$	0.5603	0.0450	−1.8017	−1.8434	−0.6549	⋯
$[{\rm{S}}\,{\rm{II}}]$/${\rm{H}}\alpha$	−0.5457	0.4573	−0.8227	−0.0284	0.5940	0.3426

Note. The value of the EL flux ratio is calculated by the polynomial functional form, i.e., $\mathrm{log}R={\sum }_{i=0}^{5}{c}_{i}\cdot {x}^{i}$, where $x=12+\mathrm{log}({\rm{O}}/{\rm{H}})-8.69$. The ratio of $[{\rm{S}}\,{\rm{II}}]$/${\rm{H}}\alpha$ is calibrated by the work of Jones et al. (2015a), and the Balmer line ratios correspond to the Balmer decrement, whereas the ratios between oxygen lines and ${\rm{H}}\beta$ are from M08.

Download table as:  ASCIITypeset image

In addition, the $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ ratio is included as an additional parameter for galaxies at $z\lesssim 1.5$ where these lines are observed. However, we do not attempt to determine $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ using the M08 calibrations, as locally calibrated diagnostics tend to underestimate the $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ ratio in high-z galaxies whereas oxygen and other α-element line ratios remain constant with redshift (Steidel et al. 2014; Jones et al. 2015a; Shapley et al. 2015). Instead, we either leave this parameter free ("extended" priors), or fixed to $[{\rm{N}}\,{\rm{II}}]/{\rm{H}}\alpha =0.05$ ("vanilla" priors). The vanilla value is typical of galaxies with similar z and oxygen EL ratios as our sample. Fixing the value of $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ provides faster convergence and does not significantly affect the inferred metallicity. We therefore report constraints for the "vanilla" parameter set. For galaxy ID 04054, the SN Refsdal host galaxy, we fix its $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ to be 0.11, as measured by Yuan et al. (2011).

We use the Markov Chain Monte Carlo sampler Emcee (Foreman-Mackey et al. 2013) to explore the parameter space, setting the number of walkers to 100 with 5000 iterations for each walker, and a burn-in of 2000. The autocorrelation times are computed for each parameter and are used to make sure chains have converged. An example constraint result for the SN Refsdal host galaxy's global EL fluxes is shown in Figure 5. While this paper is primarily concerned with gas-phase metallicities, our method simultaneously provides constraints on ${A}_{{\rm{V}}}^{{\rm{N}}}$ and SFR, which can be compared to the results of SED fitting. We then calculate the lensing-corrected de-reddened ${\rm{H}}\alpha$ luminosity from f_Hβ (assuming the fiducial Balmer ratio) and the magnification estimate of the source. Then the instantaneous SFR can be derived, via the commonly used calibration of Kennicutt (1998) and Moustakas et al. (2010), i.e.,

$$\begin{eqnarray}&&\mathrm{SFR}=4.6\times {10}^{-42}\,\displaystyle \frac{L({\rm{H}}\alpha )}{\mathrm{erg}\,{{\rm{s}}}^{-1}}\quad [{M}_{\odot }\,{\mathrm{yr}}^{-1}],\end{eqnarray} \tag{ 3 }$$

which is valid for a Chabrier (2003) IMF. This SFR estimate does not rely on any assumptions about star-formation history, and thus is a more direct measure of ongoing star formation than values given by SED fitting.

Zoom In Zoom Out Reset image size

We apply this inference procedure to both galaxy-integrated line fluxes (to be discussed in Section 5) as well as individual spaxels in the EL maps (which will lead to metallicity gradient measurements described in Section 6).

To estimate the spatial distribution of the constrained parameters, we apply the analysis described in Section 5 to each individual spaxel in the EL maps. In Figure 10, we show the maps of best-fit metallicity and the corresponding conservative uncertainty (the larger side of asymmetric error bars given by the Bayesian inference) measured for galaxies in our metallicity gradient sample. We bin spaxels 2 by 2 to regain the native WFC3/IR pixel scale. In the maps, we only include spaxels with at least one EL (primarily ${\rm{H}}\alpha$ or $[{\rm{O}}\,{\rm{III}}]$) detected with $\geqslant 5\sigma$ significance. The metallicity gradients measured from these individual spaxels are shown in the right column of Figure 10. We also measure metallicity gradients using another independent method, i.e., via radial binning. The range of each radial annulus is determined by requiring its S/N ≳ 15. In order to avoid biasing toward low metallicities, we put the S/N threshold on ${\rm{H}}\alpha$ whenever available (i.e., $z\lesssim 1.5$) instead of on $[{\rm{O}}\,{\rm{III}}]$, since the latter correlates tightly with $12+\mathrm{log}({\rm{O}}/{\rm{H}})$, in the lower branch of $[{\rm{O}}\,{\rm{III}}]$/${\rm{H}}\beta$, where the majority of our sample resides. Based upon our S/N threshold criterion, three sources in our sample, IDs 02918, 05422, and 07058, do not have enough spaxels to constitute a trustworthy metallicity map, and therefore were removed from the spatially resolved analysis. So in total, we present 10 star-forming galaxies with sub-kiloparsec resolution metallicity maps. In order to measure radial gradients, we perform Bayesian linear regressions in the same manner as described in Section 5. We take into account the conservative uncertainties for individual metallicity constraints in each spaxel or annulus. The resulted metallicity gradients are given in Table 5. We note that the metallicity maps are not always symmetric around the center of stellar mass. Therefore, while gradients are a useful statistic to describe the overall behavior and to compare with numerical simulations, they do not contain all the available information about the apparent diversity of morphologies. Thus, it is helpful to keep the 2D maps in mind while interpreting the gradients.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result, seven galaxies have gradients consistent with being flat at 2σ, among which three show almost uniformly distributed metals (IDs 02806, 03746, and 05968), two have marginally positive gradients (IDs 01422 and 05732), and the other two display mildly negative gradients (IDs 02389 and 04054). Apart from these, the other three galaxies in our sample have very steep negative gradients: IDs 07255, 05709, and 05811. In particular, after the lensing magnification correction, the stellar mass of galaxy ID 05709 is estimated to be $\sim {10}^{7.9}$ ${M}_{\odot }$. This is the first time a sub-kiloparsec scale spatially resolved analysis has been done on such low-mass systems at high redshifts.

In order to verify that our results are not contaminated by ionizing radiation from AGNs, we examine spatially resolved BPT diagrams. This approach is only possible for sources at $z\lesssim 1.5$, where ${\rm{H}}\alpha$ is detectable. Furthermore, due to the low spectral resolution of HST grisms, we slightly modify the BPT diagram to plot a flux ratio of $[{\rm{O}}\,{\rm{III}}]$/${\rm{H}}\beta$ as a function of $[{\rm{S}}\,{\rm{II}}]$/${\rm{H}}\alpha$ + $[{\rm{N}}\,{\rm{II}}]$, using the "vanilla" value for $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$. As Figure 11 shows, we do not identify any strong correlation between source plane galactocentric radius and deviation from the H ii region loci, suggesting that there are no AGNs hidden in the center of these galaxies. Moreover, the integrated EL fluxes of all three galaxies lie in close proximity to the extreme starburst model prescribed by Kewley et al. (2006), which confirms our working assumption that AGN contributions are negligible for our sample, except in one case (i.e., ID 02607, see Figure 6).

Zoom In Zoom Out Reset image size

We also test our metallicity estimates using the pure empirical EL calibrations based on a large sample of "direct" metallicity determinations in SDSS local H ii regions recently published by Curti et al. (2016), replacing the M08 hybrid EL calibrations. It is found out that our results on integrated metallicities and metallicity gradients of the galaxies in our sample remain unaffected.

The radial range where metallicity gradients are measured is critical since H ii regions at far outer/inner regions of galaxies are found to show significantly elevated/truncated oxygen abundances (see, e.g., Sánchez et al. 2014). A reasonable radial range for metallicity gradient measurements is believed to be between the disk effective radius (Freeman 1970; Sánchez et al. 2014) and the optical radius (Mollá & Díaz 2005; Mollá et al. 2016c), i.e., roughly twice the size of the disk effective radius. We verified that all of our metallicity gradients are derived in this radial range, which enhances the significance of the comparative studies (of our measurements at least) presented in Sections 6.1.1 and 6.1.2.

6.1.1. Cosmic Evolution of Metallicity Gradients

We plot the evolution of the metallicity gradient across cosmic time in Figure 12. Here we include all existing high-z gradient measurements obtained with sufficiently high spatial sampling, i.e., finer than $\mathrm{kpc}$ resolution in the source plane. Together with the 10 new metallicity gradient measurements from the individual spaxel method presented here, in Figure 12, we show one highly magnified galaxy of an interacting triple at z = 1.85 analyzed in our previous work (Jones et al. 2015b). We also include 11 measurements by Leethochawalit et al. (2015; 5 isolated and 6 merging), 4 by Jones et al. (2013; 3 isolated and 1 merging), 1 sub-${L}_{* }$ post-merger at $z\sim 1$ by Frye et al. (2012), and 7 data points from Swinbank et al. (2012; 5 isolated and 2 interacting). The scatter of these high-resolution gradient measurements is ∼0.12 dex and ∼0.22 dex at $z\sim 1.5$ and $z\sim 2.3$ respectively. The spread of recent gradient measurements from the seeing-limited ${\mathrm{KMOS}}^{3{\rm{D}}}$ survey is also overlaid in Figure 12. Note that the PSF FWHM given by the median seeing of their observations is 06, which results in a 5 $\mathrm{kpc}$ resolution at $z\sim 2$, twice the size of the half-light radius of an ${L}_{* }$ galaxy at this redshift. We thus only focus on interpreting the sub-kiloparsec resolution gradient measurements, in order to avoid possible biases toward null values associated with poor spatial sampling (see, e.g., Yuan et al. 2013). Although corrected for beam smearing, the ${\mathrm{KMOS}}^{3{\rm{D}}}$ gradient measurements still display small scatter than that of the high-resolution results. When possible, we also divide each data set in terms of sources being isolated or kinematically disturbed (i.e., undergoing merging). In part of our sample, signatures of post-merger tidal remnants can be clearly identified, which strongly indicates that gravitational interaction plays a key role in shaping the spatial distribution of metals in these galaxies (IDs 02806, 03746, and 05968).

Zoom In Zoom Out Reset image size

In order to compare theoretical predictions with observations, we incorporate into Figure 12 the predictions for the metallicity gradient cosmic evolution from canonical chemical evolution models (CEMs) based upon the "inside-out" disk growth paradigm given by Chiappini et al. (2001, C01) and Mollá et al. (2016a, 2016b, M16). C01 constructed a two-phase accretion model for galaxy formation. The first infall period corresponds to the formation of the dark matter halo and the thick disk, with an exponentially declining infall rate at a fixed e-folding timescale (∼1 Gyr). In this phase, the infall of pristine gas is rapid and isotropic, giving birth to a relatively extended profile of star formation. The thin disk is formed in the second phase, where gas is preferentially accreted to the periphery of the disk, with the e-folding timescale proportional to galactocentric radius (vis-á-vis the "inside-out" growth). The foundation of M16 is Mollá & Díaz (2005), who proposed their prescription by treating the ISM as a multiphase mixture of hot diffuse gas and cold condensing molecular clouds, and calibrating against different star-formation efficiencies. Popular nucleosynthesis prescriptions for SNe Ia/II and AGB stars are employed by both to construct the yields of dominant isotopes (e.g., H, He, C, N, O, Fe). Improving over what is done in the Mollá & Díaz (2005) multiphase CEMs, M16 adopted newly calibrated gas infall rates from Mollá et al. (2016c) and a variety of different prescriptions for H i-to-${{\rm{H}}}_{2}$ conversion, with the best being the ASC technique (Y. Ascasíbar et al. 2017, in preparation). We thus show this model solution in the black dashed–dotted curve. As shown in Figure 12, the majority of the high-z measurements are generally compatible with both model predictions, with the M16 model showing better consistency, largely ascribed to its more accurate assumptions of ISM structure and more realistic treatment of star formation. Note that these CEMs do not take galactic feedback or radial gas flows into account. So without the consideration of strong physical mechanisms that stir and mix up the enriched gas with relatively unenriched ambient ISM, conventional analytical models are still able to match the metallicity gradients observed at high redshifts.

In addition to the predictions from these analytical models, we also compare our measurements with results from numerical simulations. First, we show the fiducial representative (i.e., the realization g15784) of the McMaster Unbiased Galaxy Simulations (MUGS, Stinson et al. 2010) using the gravitational N-body and SPH code Gasoline. This simulation employed the "conventional" feedback scheme, i.e., depositing ∼10%–40% kinetic energy from SN explosions into heating ambient ISM. Furthermore, it adopted a high star-formation threshold (gas particle density >1 cm⁻³), which makes star-forming activities highly concentrated in the galaxy center at early stages of disk formation. As a result, a steep gradient (i.e., $\gtrsim -0.2$ dex kpc⁻¹) is persisting at $z\gt 1$ and then significantly flattens with time. Although greatly offset from most of the measurements (and other theoretical predictions) as seen in Figure 12, this particular evolutionary track is actually quite consistent with some of the observed steep gradients (galaxies 07255, 05709, and 05811 in our analysis). In comparison, the same galaxy (g15784) is re-simulated with an enhanced feedback prescription of sub-grid physics as part of the Making Galaxies In a Cosmological Context (MaGICC, Stinson et al. 2013) program. A factor of 2.5 difference in the energy output rate from SN feedback diverges the evolutionary track of metallicity gradients given by the MUGS and MaGICC simulations, making the evolution of the MaGICC gradient less dramatic and more similar to the M16 model solution. This demonstrates that amplifying the feedback strength can significantly flatten the metallicity gradients in early star-forming galaxies (Gibson et al. 2013). This is actually adopted by some recent studies (Leethochawalit et al. 2015; Wuyts et al. 2016) to explain their metallicity gradient measurements, which are also shown in Figure 12. However, as we have already mentioned, this is not the only explanation possible.

Lastly, in Figure 12, we also show another set of numerical results, from an ensemble of 19 galaxies in different environments (9 in the isolated field and 10 in loose groups), simulated by the Ramses Disk Environment Study (RaDES, Few et al. 2012) using the adaptive mesh refinement code Ramses. Unlike what is shown for the MUGS and MaGICC simulations (i.e., just one single realization), the full range of the whole suite of RaDES simulations is marked in Figure 12, which provides a statistical comparison sample, though the scatter of the RaDES simulations cannot capture the scatter seen in actual measurements. Generally speaking, the evolutionary trends of the abundance gradients of RaDES galaxies lie somewhat in-between the two tracks of g15784 in MUGS and MaGICC simulations, but are still quite consistent with the majority of the high-z gradient measurements. Because no mass loading is assumed in the RaDES simulations, outflows do not play any role in this observed consistency. Moreover, since a lower star-formation threshold is used (>0.1 cm⁻³), which is similar to what the analytical CEMs assumed, star formation in RaDES galaxies occurs more uniformly in the early disk-formation phase. This strongly suggests that having more extended star formation can also result in galaxies possessing shallower gradients, confirming the tight link between star-formation profile and metallicity gradient (Pilkington et al. 2012).

6.1.2. Mass Dependence of Metallicity Gradients

In Figure 13, we plot a metallicity gradient as a function of stellar mass for a subset of all gradient measurements shown in Figure 12, where reliable ${M}_{* }$ estimates are available. Notably, in the low-mass regime of ${M}_{* }\lesssim {10}^{9}\,{M}_{\odot }$, all the current sub-kiloparsec resolution metallicity gradient measurements at $z\gtrsim 1$ are obtained from our work, and our measurements even extend to below 10⁸ ${M}_{\odot }$. This illustrates the power of combining the space-based grism spectroscopy with gravitational lensing as a means to recover reliable metal abundance distributions in early star-forming disks in the low-mass regime. From these measurements, we can tentatively observe an intriguing correlation between ${M}_{* }$ and metallicity gradient: more massive galaxies tend to have shallower (less negative) gradients. This is consistent with the galaxy formation "downsizing" picture (see, e.g., Brinchmann et al. 2004; Fontanot et al. 2009) such that more massive galaxies are more evolved and their metal distributions are flatter since they are in a later phase of disk mass assembly where star formation occurs in a more coherent mode throughout the disk. This correlation between mass and gradient slope has also been seen in a variety of theoretical studies including the RaDES simulations (Few et al. 2012) and the multiphase CEMs by (Mollá & Díaz 2005).

Zoom In Zoom Out Reset image size

Several groups of authors have reported a characteristic oxygen abundance gradient for nearby disk galaxies, such that the normalized gradient in units of dex per scale radius is approximately invariant with galaxy global properties (${M}_{* }$, morphology, etc.; e.g., Sánchez et al. 2014; Bresolin & Kennicutt 2015; Ho et al. 2015). Consequently, larger and more massive galaxies have shallower gradients in physical units (dex $\mathrm{kpc}$⁻¹; e.g., Ho et al. 2015), as we tentatively see in Figure 13. We verify that 8/10 of our gradient measurements, if expressed in terms of disk effective radii, are consistent with the common gradient slope reported by Sánchez et al. (2014). This common normalized abundance gradient may reflect a more fundamental relation between local surface mass density and metallicity (e.g., Bresolin & Kennicutt 2015). Both of these quantities decrease more rapidly with radius in smaller galaxies, in line with our "downsizing" picture of metallicity gradients. Downsizing is further reinforced by the evolution in the size–mass relation recovered by van der Wel et al. (2014).

As a comparison, we include the linear fit to the low spatial resolution metallicity gradient measurements from the ${\mathrm{KMOS}}^{3{\rm{D}}}$ survey in the similar z range (see Figure 6 in Wuyts et al. 2016). Since the mass completeness limit for the ${\mathrm{KMOS}}^{3{\rm{D}}}$ data set is $\sim {10}^{9.7}\,{M}_{\odot }$, our analysis is highly complementary to the ${\mathrm{KMOS}}^{3{\rm{D}}}$ result in terms of ${M}_{* }$, albeit our metallicity gradients are measured at much higher spatial resolution. However, the ${\mathrm{KMOS}}^{3{\rm{D}}}$ fitting relation deviates systematically from the "downsizing" mass dependency of metallicity gradients and is also in conflict with some sub-kiloparsec resolution gradient measurements obtained by Jones et al. (2013) at an over 3σ confidence level. This can be attributed to possible biases associated with the low spatial resolution of the ${\mathrm{KMOS}}^{3{\rm{D}}}$ observations and the small mass range of their sample.

We also show the 2σ spread of the galaxy metallicity gradients from the FIRE cosmological zoom-in simulations (Ma et al. 2016b). This suite of simulations has been used to study the MZR as we discussed in Section 5 (Ma et al. 2016a). The simulations show a diversity of gradient behaviors at four redshift epochs (z = 0, 0.8, 1.4, 2.0). At z = 2, the scatter of their simulation results can reach 0.2 dex when radial gradients are measured in the central 2 kpc region, quite consistent with the scatter (∼0.22 dex) seen in high-resolution gradient measurements at $z\sim 2.3$. However, their simulations still cannot reproduce some of the steep gradient measurements, especially some of our results in the low-mass regime. According to Ma et al. (2016b), galactic feedback from intense starburst episodes on a 0.1–1 Gyr timescale can effectively flatten metallicity gradients. Hence for the galaxies measured with steep metallicity gradients, galactic feedback must be of limited influence, at least in affecting metal distribution in H ii regions on such a short timescale.

In short, we observe a tentative mass dependence of metallicity gradients that low-mass galaxies have steeper gradients compared with high-mass galaxies at high redshifts. A larger set of accurately measured and uniformly analyzed metallicity gradients well sampled in ${M}_{* }$ (in the regime of ${M}_{* }\lesssim {10}^{9}\,{M}_{\odot }$ in particular) is needed in order to fully explore the validity of the metallicity gradient "downsizing" picture and whether galactic feedback plays a crucial role in altering this picture.

In this section, we highlight some noteworthy sources in our metallicity gradient sample. For reference, all the measurements of global properties, metallicity gradients, and kinematic decompositions in these sources are presented in Tables 2, 5, and 6, respectively.

02389: This is the most highly magnified source in our sample, with $\mu ={43.25}_{-5.78}^{+7.03}$, given by the Sharon & Johnson version 3 model. It is one of two connected double images that straddle the lensing critical curve at a source redshift of z = 1.89, forming a highly sheared arc more than 5 arcsec. This is confirmed by the velocity shear of ∼50 $\mathrm{km}\,{{\rm{s}}}^{-1}$ detected in both components of the fold arc. After lensing correction, the stellar mass of this source is $\mathrm{log}({M}_{* }/{M}_{\odot })={8.43}_{-0.07}^{+0.06}$, consistent with the very low dynamical mass indicated by both the velocity shear and integrated velocity dispersion. The relatively high $V/\sigma \gtrsim 2$ of this source shows that it is a disk galaxy. The metallicity gradient is mildly negative, i.e., −0.07 ± 0.04, consistent with the prediction by analytical CEMs for disk galaxies at similar redshifts.

02806: This is a star-forming galaxy at z = 1.50 with a very flat metallicity gradient of −0.01 ± 0.02. We removed the contamination from a marginally overlapping cluster member (with greenish color in the color composite stamp shown in Figure 4). Note that such a process would not be possible with seeing-limited data quality. The MUSE observation does not spatially resolve this source well, and its $[{\rm{O}}\,{\rm{II}}]$ emission line is detected with low S/N near the red edge of MUSE spectral coverage. A tentative velocity shear is seen at low significance. The faint fuzzy elongation to the southeast of this galaxy suggests a likely recent merging event, which accounts for its gradient being flat. (See sources 03746 and 05968 for more outstanding cases where post-merger features can be identified.) The global metallicity of this source ($12+\mathrm{log}({\rm{O}}/{\rm{H}})={7.91}_{-0.16}^{+0.16}$) is slightly lower than the average given by the MZR at similar ${M}_{* }$, in agreement with the study of Michel-Dansac et al. (2008), suggesting that the infall of low-metallicity gas during mergers reduces the average metallicity.

03746: This is an intense starburst galaxy with $\mathrm{SFR}\,={10.95}_{-0.58}^{+1.61}$ ${M}_{\odot }$ yr⁻¹ at z = 1.25, as revealed by the BPT diagram in Figure 11. It also displays a flat metallicity gradient of −0.03 ± 0.03. Moreover, the EL maps of this source show highly clumpy star-forming regions, which do not have direct correspondence in the broadband photometry. The MUSE observation gives a velocity shear of ∼20 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and an integrated dispersion of ∼60 $\mathrm{km}\,{{\rm{s}}}^{-1}$, which suggests that this system is dispersion dominated ($V/\sigma \lt 1$). There is an extended feature southward to the source center, which is probably a tidal remnant, consistent with this source being in a recent post-merger phase. Besides, given the irregularities seen in the central region of the ${M}_{* }$ map, this system is still long before dynamical relaxation, and we likely capture this merger system almost right after the coalescence. Therefore, we conclude that the flat metallicity gradient recovered in this system is from tidal interactions. See galaxy 05968 below for an even more outstanding post-merger case.

04054: This is one of the multiple images of the SN Refsdal host galaxy, a highly magnified almost face-on spiral at z = 1.49. It is labeled as image 1.3 by Treu et al. (2016).²⁸ Our global demographic analysis yields the following measurements: $\mathrm{log}({M}_{* }/{M}_{\odot })={9.64}_{-0.01}^{+0.05},\mathrm{SFR}={16.99}_{-2.80}^{+5.71}$ ${M}_{\odot }$ yr⁻¹, and $12+\mathrm{log}({\rm{O}}/{\rm{H}})={8.70}_{-0.11}^{+0.09}$. Compared with the previous analysis of this galaxy (Livermore et al. 2015; Rawle et al. 2016, not necessarily the same image though), our work benefits from a more complete wavelength coverage of imaging data (from ${B}_{435}$ to ${H}_{160}$ bands) at a greater depth. In particular, we used the most updated lens model, optimized specifically for this spiral galaxy (Treu et al. 2016), giving a reliable estimate of the magnification. Rawle et al. (2016) measured SFR ∼ 5 ${M}_{\odot }$ yr⁻¹ on image 1.1 of the SN Refsdal host from Herschel IR photometry. However, their value is derived assuming a magnification factor of $\mu =23.0$, given by an older lens model for MACS1149.6+2223 (Smith et al. 2009). If the magnification value of $\mu =8.14$ given by the latest Glafic model is adopted instead, their measurement yields SFR ∼ 14, comparable to our SFR measurement. As a result, our measurements are fairly consistent with the SFMS, the MZR and the FMR at similar redshifts. Given the high ${M}_{* }$ and $12+\mathrm{log}({\rm{O}}/{\rm{H}})$, this galaxy is significantly evolved by $z\sim 1.5$ and with the still relatively high SFR, it will likely turn into a galaxy more massive than our Milky Way at z = 0. According to the kinematic map derived from the MUSE spectroscopy of the $[{\rm{O}}\,{\rm{II}}]$ EL doublet, we obtain the velocity shear to be ∼110 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (see also Figure 8 in Grillo et al. 2016) and the integrated dispersion ∼60 $\mathrm{km}\,{{\rm{s}}}^{-1}$, in good agreement with Livermore et al. (2015) who found $2{v}_{2.2}=118\pm 6$ $\mathrm{km}\,{{\rm{s}}}^{-1}$, and local (spatially resolved) velocity dispersion of 50 ± 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$, based upon AO data.

The measured metallicity gradient for the SN Refsdal host based upon our analysis of image 1.3 is −0.04 ± 0.02 from the individual spaxel method, and −0.07 ± 0.02 from the radial annulus method. This measurement is very close to the usual abundance gradient of local spirals, including our Milky Way, measured from either H ii regions or planetary nebulae (Smartt & Rolleston 1997; van Zee et al. 1998; Henry et al. 2004; Esteban et al. 2015). It is worth noting that our measurement is not as steep as the previous result reported by Yuan et al. (2011), i.e., −0.16 ± 0.02, measured from the AO assisted OSIRIS IFU observation of image 1.1 of the SN Refsdal host, assuming the Pettini & Pagel (2004) $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ empirical calibration. After the discovery of SN Refsdal, Kelly et al. (2015c) obtained a 1 hr MOSFIRE spectra on the nuclear region and the SN explosion site in image 1.1, and measured $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.6\pm 0.1$ and $12+\mathrm{log}({\rm{O}}/{\rm{H}})=8.3\pm 0.1$, respectively, using the same $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ calibration. The spatial offset between these two measurements in terms of galactocentric radius is 8.2 ± 0.5 $\mathrm{kpc}$, given by the Glafic model. As a result, the recent MOSFIRE observations indicate a gradient being −0.04 ± 0.02, highly consistent with our measurements.

Moreover, we obtain a robust constraint on the nebular dust extinction, i.e., ${A}_{{\rm{V}}}^{{\rm{N}}}={0.72}_{-0.23}^{+0.23}$. Compared with the stellar dust extinction retrieved from SED fitting, i.e., ${A}_{{\rm{V}}}^{{\rm{S}}}={1.10}_{-0.01}^{+0.01}$, we see that in this galaxy, the dust reddening of the ionized gas in H ii regions is less severe than that of the stellar continuum (see another source, i.e., ID 05968, from which the same conclusion is drawn). This is different from what has been assumed previously: the color excess of nebular ELs is larger than that of the stellar continuum, i.e., ${E}_{{\rm{S}}}(B-V)=(0.44\pm 0.03){E}_{{\rm{N}}}(B-V)$ (Calzetti et al. 2000). Recently, Reddy et al. (2015) conducted a systematic study of dust reddening using the early observations from the MOSDEF survey and found that only 7% of their sample is consistent with ${E}_{{\rm{N}}}(B-V)\gtrsim {E}_{{\rm{S}}}(B-V)/0.44$. They also discovered that as sSFR (derived from Balmer EL luminosities) increases, ${E}_{{\rm{N}}}(B-V)$ diverges more significantly from ${E}_{{\rm{S}}}(B-V)$, with a scatter of ∼0.3–0.4 dex at sSFR ∼ 3 × 10⁻⁹ yr⁻¹, which can account for the difference we see in the total dust attenuation (${A}_{{\rm{V}}}$) of the line and continuum emissions of this galaxy.

05709: This galaxy has the lowest stellar mass in our metallicity gradient sample ($\mathrm{log}({M}_{* }/{M}_{\odot })={7.90}_{-0.03}^{+0.02}$) after lensing correction.²⁹ This is for the first time any diffraction-limited metallicity gradient analysis and seeing-limited gas kinematics measurement have been achieved in star-forming galaxies with such small ${M}_{* }$ at the cosmic noon. This galaxy shows a steep metallicity gradient, i.e., −0.22 ± 0.05. The velocity shear is measured to be ∼25 $\mathrm{km}\,{{\rm{s}}}^{-1}$, with integrated dispersion being ∼20 $\mathrm{km}\,{{\rm{s}}}^{-1}$. The inferred low dynamical mass is consistent with the result from SED fitting mentioned above. In light of the kinematic result $V/\sigma \gtrsim 1$, this galaxy is probably a clumpy thick disk. More interestingly, observed from a multi-perspective point of view (see Figure 4), this galaxy shows a clumpy star-forming region at the periphery of its disk (having a distance of 3–4 kpc to the mass center), which has a tremendous amount of line emission, yet very little stellar mass, and has low gas-phase metallicity. Also see galaxy ID 07255 for a more dramatic case.

05968: With a tidal tail more appreciably shown to the northeast, and almost perpendicular to the galactic plane, this galaxy (at z = 1.48) is safely classified as a post-merger. However, note that this identification would not be possible and this source would have surely been classified as a thick disk due to its symmetric morphology and kinematic properties, without the ultra-deep HFF imaging data. Compared to the ${M}_{* }$ and EL maps of ID 03746, the ${M}_{* }$ map of this galaxy is more smooth and EL maps are less clumpy, which indicates that merging took place longer before than that in source ID 03746. Still, the star formation is high, i.e., $\mathrm{SFR}={27.95}_{-4.41}^{+4.70}$ ${M}_{\odot }$ yr⁻¹. The global metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})={8.47}_{-0.08}^{+0.07}$) and stellar mass ($\mathrm{log}({M}_{* }/{M}_{\odot })={9.34}_{-0.03}^{+0.02}$) are quite comparable to the MZR at similar redshifts. The metallicity gradient, as expected for systems that have experienced recent mergers (see, e.g., Kewley et al. 2010), is constrained to be robustly flat (−0.01 ± 0.02). The EL kinematic analysis yields a velocity shear of ∼120 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and an integrated dispersion of ∼75 $\mathrm{km}\,{{\rm{s}}}^{-1}$. The velocity field appears typical of intermediate-mass disk galaxies at this redshift, which infers that this post-merger system has been significantly relaxed and regained at least some of the rotation support. This is also one of the few cases where robust constraints on nebular dust extinction can be derived, i.e., ${A}_{{\rm{V}}}^{{\rm{N}}}={0.23}_{-0.13}^{+0.15}$, which is less than the dust extinction of stellar continuum obtained from SED fitting ${A}_{{\rm{V}}}^{{\rm{S}}}={0.60}_{-0.01}^{+0.20}$, as already seen in the case of ID 04054.

07255: This galaxy is highly consistent with the SFMS at z = 1.27 formulated by Speagle et al. (2014), with $\mathrm{log}({M}_{* }/{M}_{\odot })={9.36}_{-0.21}^{+0.01}$ and $\mathrm{SFR}={8.73}_{-3.10}^{+3.53}$ ${M}_{\odot }$ yr⁻¹. By z = 0, it will turn into a sub-Milky-Way-sized galaxy according to Behroozi et al. (2013). Its slightly lower than average metallicity of $12+\mathrm{log}({\rm{O}}/{\rm{H}})={8.38}_{-0.08}^{+0.08}$ can be a result of low-metallicity gas inflow, following the popular interpretation of the FMR at high z. This hypothesis would also explain its slightly enhanced SFR. In the metallicity map, the inflow can be identified with a bright star-forming clump in the lower right corner, which consists of very little ${M}_{* }$, but is dominating in the EL maps, similar to the case of galaxy 05709 discussed above. This indicates that in this system, the metal enrichment timescale is much larger than the star-forming timescale, which is in turn much larger than the gas infall timescale. Unfortunately, at the moment, we do not have EL kinematic observation on this system.

The metallicity gradient for this galaxy is constrained to be very steep, i.e., −0.16 ± 0.03 using the individual spaxel method, and −0.21 ± 0.03 from the radial annulus method. We also dissected this galaxy to separate the star-forming clump from the bulk part and re-do the analysis. It is found that the H ii regions closely related to the clump (with at least 2.5 kpc distance to the galaxy mass center) have significantly lower metallicity and provide more than one-third of the entire star formation, i.e., $12+\mathrm{log}({\rm{O}}/{\rm{H}})={8.04}_{-0.18}^{+0.14}$, ${A}_{{\rm{V}}}^{{\rm{N}}}\lt 0.78$, $\mathrm{SFR}={3.05}_{-1.10}^{+1.21}$, compared to the quantities measured in the bulk of the galaxy (within 2 kpc to the mass center), i.e., $12+\mathrm{log}({\rm{O}}/{\rm{H}})={8.51}_{-0.10}^{+0.08}$, ${A}_{{\rm{V}}}^{{\rm{N}}}={0.63}_{-0.41}^{+0.59}$, and $\mathrm{SFR}\,={5.78}_{-2.47}^{+3.41}$. In both regions, EL ${\rm{H}}\alpha$ is detected at a significance of $\gtrsim 30\sigma$. This again demonstrates that measurements of global quantities can be to some extent misleading and also demonstrates the importance of detailed spatially resolved information. Cases like this galaxy should be exceedingly interesting, since their steep metallicity gradients at such high z are not predicted by current mainstream numerical simulations (equipped with strong galactic feedback) or conventional analytical CEMs. Through studying systems like this, we can also learn how gas is accreted into the inner disk, and more importantly whether the accretion happens alongside star formation or long before the gas clouds can collapse.

In this paper, we took advantage of the deep HST grism spectroscopic data, acquired by the GLASS and SN Refsdal follow-up programs, and obtained gas-phase metallicity maps in high-z star-forming galaxies. The HST grisms provide an uninterrupted window for observing the cosmic evolution of metallicity gradients through z from 2.3 to 1.2 (corresponding to the age of the universe from roughly 2.5 to 5.5 Gyr), including the redshift range of $1.7\lt z\lt 2.0$, which is unattainable from ground-based observations due to low atmospheric transmission. Another advantage of the HST grisms is that multiple ELs are available within a single setup and all the potential issues related to combining data with different setups and atmospheric conditions are eliminated. Our main conclusions can be summarized as follows.

• 1.  
  We presented 10 maps of gas-phase metallicity (as shown in Figure 10), measured at sub-kiloparsec resolution in star-forming galaxies in the redshift range of $1.2\lesssim z\lesssim 2.3$, in the prime field of MACS1149.6+2223. This field is so far one of the deepest fields exposed by HST spectroscopy, at a depth of 34(10) orbits of G141(G102), reaching a 1σ flux limit of 1.2(3.5) × 10⁻¹⁸ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$, without lensing correction.
• 2.  
  We developed a novel Bayesian statistical approach to jointly constrain gas-phase metallicity, nebular dust extinction, and ${\rm{H}}\alpha$-based dust-corrected SFR. In determining these quantities, our method does not rely on any assumptions about star-formation history, nor the conversion of dust reddening from the stellar phase to the nebular phase. Unlike the majority of previous work on deriving metallicity from strong EL calibrations, our method does not compare directly the observed EL flux ratios with the calibrated ones. Instead, we work directly with EL fluxes, which effectively avoids using redundant information and circumvents possible biases in the low S/N regime.
• 3.  
  The metallicity maps we obtained show a large diversity of morphologies. Some maps display disk-like shapes and have metallicity gradients consistent with those predicted by analytical CEMs (i.e., IDs 04054, 01422, 02389, 05732), whereas other disk-shaped galaxies have exceedingly steep gradients (i.e., IDs 07255, 05709, 05811) and bright star-forming clumps at the periphery of their disks. These large-offset star-forming clumps, containing very little stellar mass, are estimated to have lower metallicity than the corresponding global values, which can be interpreted as evidence for the infall of low-metallicity gas enhancing star formation. We also recovered three systems with nearly uniform spatial distribution of metals (i.e., IDs 02806, 03746, 05968). In these systems, we can identify possible signatures of post-merger tidal remnants, which suggests that the observed flat gradients are likely caused by gravitational interactions.
• 4.  
  Collecting all existing sub-kiloparsec resolution metallicity gradient measurements at high redshifts, we study the cosmic evolution and mass dependence of metallicity gradients.

  • (i)  
    We found that predictions given by analytical CEMs assuming a relatively extended star-formation profile in the early phase of disk growth can reproduce the majority of observed metallicity gradients at high redshifts, without involving strong galactic feedback or radial outflows. This confirms the tight link between star formation and metal enrichment.
  • (ii)  
    We tentatively observed a correlation between stellar mass and metallicity gradient, which is in accord with the "downsizing" galaxy formation scenario that more massive galaxies are more evolved into a later phase of disk growth, where they experience more coherent mass assembly and metal enrichment than lower mass galaxies.

  A larger sample of uniformly analyzed metallicity maps well sampled in the full ${M}_{* }$ range (especially in the low-mass regime, i.e., ${M}_{* }\lesssim {10}^{9}\,{M}_{\odot }$) at high redshifts is needed to further investigate the cosmic evolution and mass dependence of metallicity gradients.
• 5.  
  We also compiled a sample of 38 global metallicity measurements all derived from the Maiolino et al. (2008) EL calibrations in the current literature (including 13 galaxies presented in this work). This sample, at $z\sim 1.8$, spans three orders of magnitude in ${M}_{* }$. We measured the MZR and tested the FMR using this sample.

  • (i)  
    The slope of the MZR constrained by this sample rules out the momentum-driven wind model by Davé et al. (2012) at a 3σ confidence level. Our MZR is consistent with the theoretical prediction given by recent cosmological zoom-in simulations, suggesting that high spatial resolution simulations are favorable in reproducing the metal enrichment history in star-forming galaxies.
  • (ii)  
    Given the intrinsic scatter and measurement uncertainties, we do not see any significant offset from the FMR except the subsample of Wuyts et al. (2012). This discrepancy is likely due to the fact that the Wuyts et al. (2012) data set is the only subsample that relies exclusively on $[{\rm{N}}\,{\rm{II}}]$/${\rm{H}}\alpha$ to estimate metallicity. We therefore advocate to avoid nitrogen lines when estimating gas-phase metallicity in high-z H ii regions.

We now turn to the interpretation of our observed metallicity gradients and maps. Apart from the flattening of metallicity gradients through natural evolution of galaxies that brings the abundances of any elements in the ISM asymptotically to the average stellar yields at the end of the evolution everywhere in the galaxy, there are four physical mechanisms on relatively short timescales that could significantly flatten a star-forming galaxy's pre-existing negative metallicity gradient (not exclusively in high-z scenarios).

• 1.  
  Efficient radial mixing by tidal effects from gravitational interactions, as seen in galaxy IDs 02806, 03746, and 05968 in our sample and the triply imaged arc 4 in our previous work (Jones et al. 2015b).
• 2.  
  Rapid recycling of metal-enriched material by enhanced galactic feedback, as seen in the results of the MaGICC numerical simulations, for example.
• 3.  
  Turbulence mixing driven by thermal instability as elucidated by Yang & Krumholz (2012).
• 4.  
  Funnelling of cold streams of low-metallicity gas directly into the galaxy center, as suggested by Dekel et al. (2009).³⁰

Besides pre-existing negative gradients being flattened by these processes, galaxies can also inherit a broadly flat metallicity gradient from a relatively extended star-formation profile in the early disk-formation phase and coherent mass assembly for the disk growth. Because there is strong evidence that the mass assembly for Milky Way progenitors takes place uniformly at all radii, maintaining almost the same ${M}_{* }$ profile from $z\sim 2.5$ to $z\sim 1$ (van Dokkum et al. 2013; Morishita et al. 2015). However, the uncertainties associated with current measurements are still insufficient to distinguish numerical simulations with enhanced feedback prescriptions from conventional analytical CEMs assuming more uniform star formation (see Figure 12).

One perhaps obvious but nevertheless important caveat common to most galaxy evolution work is that we only have access to cross-sectional data, i.e., observations of different objects at one specific epoch. Therefore, any interpretations of plots like Figure 12, which try to tie these data with longitudinal theories (describing the evolution of one specific object), should be made with caution.

Perhaps an even more intriguing interpretation of Figure 12, which is robust against the aforementioned caveat, is that there exist some really steep metallicity gradients at high z, such as our galaxy IDs 07255, 05709, and 05811. Based upon current measurements, the occurrence rate of high-z metallicity gradients more negative than −0.1 dex kpc⁻¹ is 5/34. We can derive the following constraints on the physical processes that can contribute to these steep metallicity gradients observed beyond the local universe:

• 1.  
  galactic feedback (including supernova explosions and stellar winds) must be of limited influence (via being confined for instance),
• 2.  
  star-formation efficiency is low, or at least lower than their local analogs (in simulations, this can be achieved by having a higher star-formation threshold or less concentrated molecular gas, as shown by the MUGS runs),
• 3.  
  star formation at early times is highly centrally concentrated, and the mass assembly of the inner regions is much faster than that of the outer regions, and
• 4.  
  alternatively, steep gradients are the result of low-metallicity gas falling onto the outskirts of their disks, implying that the dynamic timescale is much shorter than the star-forming timescale, which is much shorter than the metal enrichment timescale, as seen in our galaxy 07255.

We need larger samples with robust measurements in order to more accurately quantify the occurrence rate of these steep metallicity gradients in high-z star-forming disks. The existence of a distinct population of galaxies with steep metallicity gradients, if confirmed by future observations constraining their abundance, is a powerful test of galaxy formation models, because this is not predicted by most current theoretical models. Hence, these galaxies with steep metallicity gradients should be the prime targets for new numerical predictions: future JWST observations will constrain their properties in great detail, allowing discriminating tests of current theories.

Furthermore, most current numerical simulations are aimed at reproducing present-day Milky Way analogs, surrounded by a dark matter halo with virial mass $\sim {10}^{12}\,{M}_{\odot }$. It would be useful to have large samples of simulated sub-${L}_{* }$ galaxies (${M}_{* }$ below $5\times {10}^{9}\,{M}_{\odot }$ at $z\sim 2$) in order to investigate how metals are recycled in these less massive systems and provide a more suitable comparison to our metallicity gradient measurements at high z. With improvements on both the observational and theoretical sides, we will be able to answer why in these high-z star-forming disk galaxies, the metallicity gradients can be established so early and why they are so steep compared to those found in local analogs.