Discovery of Strongly Inverted Metallicity Gradients in Dwarf Galaxies at z ∼ 2

We use the diffraction-limited spatially resolved slitless spectroscopy, obtained using the HST Wide Field Camera 3 (WFC3) NIR grisms (G102 and G141), acquired by the Grism Lens-Amplified Survey from Space (GLASS, Schmidt et al. 2014; Treu et al. 2015). GLASS observes the distant universe through 10 massive galaxy clusters as natural telescopes, exposing 10 orbits of G102 (0.8–1.15 $\mu {\rm{m}}$, R ∼ 210) and four orbits of G141 (1.1–1.7 $\mu {\rm{m}}$, R ∼ 130) per sightline. This amounts to a total of ∼22 ks of G102 and ∼9 ks of G141, as well as ∼7 ks of F140W+F105W direct imaging for astrometric alignment and wavelength/flux calibrations per field. These exposures are distributed over two separate pointings per cluster with nearly orthogonal orientations, designed to help disentangle spectral contamination from neighboring objects. So two sets of G102+G141 spectra are obtained for each source, covering an uninterrupted wavelength range of 0.8–1.7 $\mu {\rm{m}}$ with almost unchanging sensitivity, and reaching a 1σ surface brightness of $3\times {10}^{-16}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{arcsec}}^{-2}$ across the entire spectral range. The GLASS collaboration has made the catalogs of their redshift identifications in the 10 fields, based on visual inspections of emission line (EL) features, publicly available at https://archive.stsci.edu/prepds/glass/.

To explore the chemical properties of galaxies at the peak epoch of cosmic chemical enrichment, we select from these catalogs a parent sample consisting of ∼300 galaxies with secure redshifts (i.e., redshift quality ≥3 in the publicly available catalogs as described by Treu et al. 2015) in the range 1.2 ≲ z ≲ 2.3. This range is chosen for the detection of multiple nebular ELs¹⁴ —in particular the Balmer lines, [O iii], and [O ii]—enabling the metallicity measurements, as in our earlier work (Jones et al. 2015b; Wang et al. 2017). The GLASS data for these ∼300 galaxies are reduced using the Grism Redshift and Line Analysis software (Grzili;¹⁵ G. Brammer et al. 2019, in preparation). Grzili presents an end-to-end processing of the paired grism and direct exposures. The procedure includes five steps: (1) preprocessing of the raw grism exposures, (2) full field-of-view (FoV) grism model construction, (3) 1D/2D spectrum extraction, (4) solving for best-fit redshift from spectral template fitting (see Appendix A for more details), and (5) refining the full FoV grism model and extractions of source 1D/2D spectrum and EL stamps. In step (1), the preprocessing consists of hot-pixel/persistence masking, cosmic-ray flagging, flat-fielding, astrometric alignment, sky background subtraction, and extraction of visit-level source catalogs and segmentation maps. In step (5), the EL stamps are drizzled onto a grid with a pixel scale of 006, Nyquist sampling the WFC3 point-spread function. We apply an additional step on the Grzili output products to obtain pure 2D maps of [O iii] λ5008 and Hβ, clean from the partial contamination of [O iii] λ4960, due to the limited grism spectral resolution and extended source morphology. Our procedure properly combines EL maps at multiple orientations, preserving angular resolution and accounting for EL blending.

In addition to the deep NIR spectroscopy, there exists a wealth of ancillary imaging data with equally high spatial resolution on the 10 GLASS fields, which encompass all six clusters in the Hubble Frontier Fields (Lotz et al. 2017) and four from the Cluster Lensing and Supernova Survey with Hubble (Postman et al. 2012). This broadband photometry, covering observed wavelengths of ∼0.4–1.7 $\mu {\rm{m}}$, can help constrain the properties of stellar populations (especially ${M}_{* }$) in our selected ∼300 galaxies at sufficient confidence. We use the images sampled with 006 pixel size, and apply kernel convolutions to match the angular resolution of all images to that of the F160W filter. We subtract contamination from intracluster light using established procedures (Morishita et al. 2017). Since our targets have rest-frame optical ELs with high equivalent widths (EWs), we subtract the nebular flux contribution from the broadband photometry to obtain the stellar continuum flux. As given in Appendix A, we model the nebular ELs as Gaussian profiles centered at the corresponding wavelengths. The nebular contribution to the broadband photometry is thereby estimated by convolving the filter throughput with the best-fit Gaussian profiles, and then subtracted off.¹⁶ We then fit the spectral energy distribution (SED) of the stellar continuum with the stellar population synthesis models of Bruzual & Charlot (2003) using the software FAST (Kriek et al. 2009). We assume the initial mass function (IMF) of Chabrier (2003), constant star formation history, stellar dust attenuation in the range ${A}_{{\rm{V}}}^{{\rm{S}}}$ = 0–4 with an extinction curve from Calzetti et al. (2000), and age ranging from 5 Myr to the Hubble time at the redshifts of our targets. Stellar metallicity is fixed to one-fifth solar and we verify that this assumption affects the results by no more than <0.05 dex on ${M}_{* }$.

Out of the parent sample of ∼300 galaxies, we are able to secure accurate (i.e., at sub-kpc resolution) radial metallicity gradients on 81 sources with suitable spatial extent and nebular emission with high signal-to-noise ratio (S/N). These extended sources typically have half-light radii R₅₀ ≳ 015. In the range $7\lesssim \mathrm{log}\left({M}_{* }/{M}_{\odot }\right)\lesssim 10$ given by the analyses in Section 2.3, our sample probes much lower ${M}_{* }$ than other surveys of spatially resolved line emission at similar redshifts (Wuyts et al. 2016; Förster Schreiber et al. 2018), thanks to the enhanced resolution from lensing magnification, and high sensitivity of the HST NIR grisms. We have previously described the properties of 10 galaxies in our sample from the cluster MACS 1149.6+2223 (X. Wang et al. 2017); results for the full sample are in preparation.

In most cases, we find that metallicity gradients are approximately flat (i.e., consistent with zero given the typical σ = 0.03 dex kpc⁻¹) or slightly negative. A minority (10/81) of our sample show positive (i.e., "inverted") gradients, which are of interest because they pose a challenge to standard models of galactic chemical evolution (e.g., Mollá & Díaz 2005). We have selected the two best examples with strongly inverted gradients for further study in this paper. The two sources are ID 03751 ($z=1.96,\,{M}_{\ast }=1.12\times {10}^{9}\,{M}_{\odot }$) in the prime field of A370, and ID 01203 ($z=1.65,\,{M}_{\ast }=2.55\times {10}^{9}\,{M}_{\odot }$) in the prime field of MACS 0744.9+3927.

Table 1 presents their properties. Figure 1 shows the color-composite HST images of these two galaxies and their 2D spatially resolved maps of the nebular ELs. Remarkably, they have ${M}_{* }$ considerably lower—by one order of magnitude—than those of previous positive gradients measured at similar redshifts (see, e.g., Cresci et al. 2010; Queyrel et al. 2012; Stott et al. 2014; Troncoso et al. 2014). To complement the low-dispersion grism spectra, we have obtained adaptive optics (AO)-assisted kinematic data on our sources using a ground-based integral-field unit (IFU) spectrograph when available. The observation of source ID 01203 is presented in Appendix B. The full data analysis is presented in Hirtenstein et al. (2019) in detail.

Zoom In Zoom Out Reset image size

Table 1.  Measured Quantities of the Two Dwarf Galaxies

ID	03751	01203
Cluster	A370	MACS 0744.9+3927
R.A. (deg.)	39.977361	116.197585
Decl. (deg.)	−1.591636	39.456698
zspec	1.96	1.65
μa	${6.35}_{-0.58}^{+0.40}$	${2.25}_{-0.03}^{+0.04}$
Observed emission line fluxes
f[O iii] [10−17 erg s−1 cm−2]	111.41 ± 0.84	117.66 ± 1.17
fHβ [10−17 erg s−1 cm−2]	17.68 ± 0.68	17.46 ± 1.06
f[O ii] [10−17 erg s−1 cm−2]	29.57 ± 0.51	34.00 ± 0.96
fHγ [10−17 erg s−1 cm−2]	7.21 ± 0.67	7.06 ± 1.00
Rest-frame equivalent widths
EW[O iii] [Å]	466.22 ± 3.52	797.14 ± 7.95
${\mathrm{EW}}_{{\rm{H}}\beta }$ [Å]	73.98 ± 2.83	118.29 ± 7.18
EW[O ii] [Å]	79.14 ± 1.37	123.91 ± 3.50
${\mathrm{EW}}_{{\rm{H}}\gamma }$ [Å]	30.18 ± 2.82	25.73 ± 3.68
Estimated physical parameters
${M}_{* }$ [109 M⊙]b	${1.12}_{-0.14}^{+0.14}$	${2.55}_{-0.04}^{+0.04}$
$12+\mathrm{log}({\rm{O}}/{\rm{H}})$ c	${8.08}_{-0.12}^{+0.11}$	${8.10}_{-0.11}^{+0.11}$
${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})/{\rm{\Delta }}r$ [dex kpc−1]	0.122 ± 0.008	0.111 ± 0.017
SFR $[{M}_{\odot }\,{\mathrm{yr}}^{-1}]$ b	25.39 ± 2.19	48.86 ± 3.04
${A}_{{\rm{V}}}$	0.84 ± 0.13	0.90 ± 0.16
${t}_{\mathrm{age}}$ [107 yr]	7.93 ± 0.88	3.98 ± 0.51
${M}_{\mathrm{gas}}$ [109 ${M}_{\odot }$]b	4.07 ± 1.27	23.85 ± 7.33
${f}_{\mathrm{gas}}$ d	0.56 ± 0.24	0.86 ± 0.35
B/Te	0.36 ± 0.14	0.14 ± 0.07
Reff [kpc]b	1.53 ± 0.12	1.66 ± 0.17
Gas kinematics
σ [km s−1]	...f	73 ± 3
V/σ	...f	1.3 ± 0.1
Measurements of the gaseous outflows at the central 1 $\mathrm{kpc}$
λ	49.9 ± 14.7	52.1 ± 20.2
Ψ $[{M}_{\odot }\,{\mathrm{yr}}^{-1}]$	311.3 ± 96.2	1700.9 ± 681.9

Notes.

^aThe magnification estimates are obtained from the Sharon & Johnson version 4 model of A370 (Johnson et al. 2014) and the Zitrin PIEMD+eNFW version 2 model of MACS 0744.9+3927 (Zitrin et al. 2015), for the two sources respectively. ^bValues presented here are corrected for lensing magnification. ^cValues represent global metallicity, inferred from integrated line fluxes. ^dHere the gas fraction is calculated according to Equation (3), using surface densities of the stellar and gas components, the latter of which is given by inverting the extended Schmidt law (Equation (4)) and the former from spatially resolved SED fitting. We caution that stellar mass estimates from resolved photometry can be systematically higher (by factors of up to 5) than spatially unresolved photometry, as elaborated in Sorba & Sawicki (2018). In our cases, the total stellar masses derived from adding up stellar surface densities given by resolved photometry amount to (3.20 ± 1.92) × 10⁹ ${M}_{\odot }$ and (3.88 ± 2.76) × 10⁹ ${M}_{\odot }$, for ID 03751 and ID 01203, respectively. ^eIn the bulge–disk decomposition, we fix the Sérsic index as n = 4 (i.e., de Vaucouleurs) for the bulge component, and n = 1 (i.e., exponential) for the disk component. ^fGround-based Keck OSIRIS follow-up observations targeting Hα or [O iii] gas kinematics for this source are not feasible due to significantly low atmospheric transmission at the corresponding wavelengths.

Download table as:  ASCIITypeset image

Since we infer metallicity from strong-line flux ratio diagnostics, calibrated by either empirical methods or theoretical methods, or a hybrid of both, it is essential to make sure that the line emission is not contaminated by ionization from an active galactic nucleus (AGN) or shock by excitation. As shown in Figure 2, we verify that our targets have a low probability (<10%) of being classified as AGNs according to the mass–excitation diagram (Juneau et al. 2014). Their individual radial annuli also have excitation and ionization states compatible with H ii regions, as revealed in their loci in the "blue" diagnostic diagrams of f_{[O iii]}/f_Hβ versus f_{[O ii]}/f_Hβ (the middle panel of Figure 2) and O₃₂ = (f_{[O iii] λ5008} + f_{[O iii] λ4960})/f_{[O ii]} versus R₂₃ = (f_{[O iii] λ5008} + f_{[O iii] λ4960} + f_{[O ii]})/f_Hβ (the right panel of Figure 2) (Lamareille et al. 2004; Rodrigues et al. 2012; Jones et al. 2015a). Notably, the intrinsic ratios R₂₃ and O₃₂ decrease with galactocentric radius for both of our sources, indicative of positive radial gradients of metallicities, opposite to the normal trend, which is confirmed with direct metallicity measurements in nearby disk galaxies (see, e.g., Bresolin et al. 2009; Berg et al. 2015; Croxall et al. 2015, 2016). In the source ID 01203 covered by our follow-up OSIRIS observations (see Appendix B), its integrated f_{[N ii]}/f_Hα (≲0.1 at 3σ) also shows no sign of AGN or shocked gas emission.

Zoom In Zoom Out Reset image size

Our measurements of radial metallicity gradients largely follow the procedures described in our previous work (Wang et al. 2017). We use a Bayesian approach to jointly infer metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$), nebular dust extinction (${A}_{{\rm{V}}}^{{\rm{N}}}$), and dereddened Hβ flux (f_Hβ). We explore the parameter space using the Markov chain Monte Carlo sampler Emcee (Foreman-Mackey et al. 2013). The likelihood function is given by $L\propto \exp (-{\chi }^{2}/2)$ with

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

where ${\mathrm{EL}}_{i}$ corresponds to each available EL: [O iii], Hβ, [O ii], and Hγ. As shown in Figures 8 and 9 in Appendix A, H$\delta$ is only weakly detected for both of our galaxies, with S/N ≲ 5 from their spatially integrated grism spectra combined from two orientations, and thereby not included in our Bayesian analysis. On the other hand, our source spectra show marked detection of the EL [Ne iii] λ3869. However, due to the relatively low wavelength dispersion of the HST NIR grism channels (i.e., ${R}_{{\rm{G}}141}\sim 130$ and R_G102 ∼ 210), [Ne iii] λ3869 is heavily blended with He i λ3889 + H8 on the red side and H9 on the blue side. As a result, we exclude [Ne iii] from our subsequent analyses as well. ${f}_{{\mathrm{EL}}_{i}}$ and ${\sigma }_{{\mathrm{EL}}_{i}}$ are the flux and uncertainty of ${\mathrm{EL}}_{i}$. R_i is the flux ratio between ${\mathrm{EL}}_{i}$ and Hβ, with ${\sigma }_{{R}_{i}}$ being the intrinsic scatter at fixed physical properties. In the case ${\mathrm{EL}}_{i}={\rm{H}}\gamma$, R_i is given by the Balmer decrement ${f}_{{\rm{H}}\gamma }/{f}_{{\rm{H}}\beta }=0.47$. For ${\mathrm{EL}}_{i}$ ∈ {[O ii], [O iii]}, R_i and ${\sigma }_{{R}_{i}}$ are given by the strong-line metallicity diagnostics (f_{[O iii]}/f_Hβ and f_{[O ii]}/f_Hβ) calibrated by Maiolino et al. (2008). The calibrations of Maiolino et al. (2008) combine the direct electron temperature measurements from the Sloan Digital Sky Survey (SDSS) in the low-metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$ ≲ 8.35) branch (Nagao et al. 2006) and the predictions of the photoionization model in the high-metallicity ($12+\mathrm{log}({\rm{O}}/{\rm{H}})$ ≳ 8.35) branch (Kewley & Dopita 2002), providing a continuous and coherent recipe over a wide metallicity range.

Our metallicity inference presented here primarily originates from the flux ratios between the oxygen collisionally excited lines (CELs) (i.e., [O iii], [O ii]) and Hβ (see the recent review by Maiolino & Mannucci 2019 for pros and cons of all strong-line ratio diagnostics in the literature). Despite their bimodal relationship with metallicity, R₂₃ and f_{[O iii]}/f_Hβ are found to produce the best accuracy among all the strong-line ratios used in z ∼ 2 studies (Patrício et al. 2018). Notably, multiple independent analyses have shown that there is a clean metallicity sequence in the diagnostic diagram spanned by f_{[O iii]}/f_Hβ and f_{[O ii]}/f_Hβ (Andrews & Martini 2013; Gebhardt et al. 2016; Jones et al. 2015a; Curti et al. 2017, also see Figure 2). Since [O iii] and [O ii] are among the brightest and thus most accessible strong ELs in high-z galaxies, the flux ratios involving oxygen CELs will remain the most promising diagnostics of O/H, since they trace directly the emissivities and abundances of the oxygen ions.

In Appendix C, we perform a comparative study of measuring metallicity gradients with the pure empirical calibrations by Jones et al. (2015a) and Curti et al. (2017), based upon the metallicity measurements using electron temperature from the data sets of DEEP2 (at z ∼ 0.8, Newman et al. 2013) and the SDSS (at 0.027 < z < 0.25, Abazajian et al. 2009). We find that the systematic differences between the absolute gradient slopes derived assuming different strong-line calibrations can be as high as 0.06 dex kpc⁻¹. However, a positive radial gradient slope can always be recovered in each of our sources at high statistical significance, regardless of the calibrations used (see Appendix C for more details). We thus verify that there is no significant bias from the strong-line calibrations adopted in our gradient measurements. The same process is applied both to galaxy-integrated fluxes and to fluxes contained in individual spatial pixels (spaxels).

To obtain the correct intrinsic deprojected distance scale for each spaxel, we conducted full morphological reconstruction of our sources in the source plane. We ray-trace the image of each galaxy to its source plane using up-to-date lens models for each cluster: the macroscopic model of Sharon & Johnson version 4 for A370 (Johnson et al. 2014), and the Zitrin version 2 model for MACS 0744.9+3927 (Zitrin et al. 2015). Other lens models are available for these clusters (e.g., Diego et al. 2018; Strait et al. 2018) and we verified that the morphology of each source is robust to the choice of model.

For each source, we fit the SED of individual spaxels using the procedures described in Section 2, obtaining the map of 2D stellar surface density (${{\rm{\Sigma }}}_{* }$) shown in Figure 1. Then we reconstruct the ${{\rm{\Sigma }}}_{* }$ map in the source plane by de-lensing the surface densities according to the deflection field given by the macroscopic lens models. To minimize the stochasticity in stellar population synthesis (Fouesneau & Lançon 2010; Eldridge 2012), we make sure that the source-plane resolution elements during this reconstruction contain enough stellar masses ($\gtrsim {10}^{5}\,{M}_{\odot }$) to be representative of complete stellar populations. The axis ratios, inclinations, and major axis orientations are determined from an elliptical Gaussian fit. This procedure provides the intrinsic lensing-corrected morphology, and in particular, the galactocentric radius at each point of the observed images. The radial scale shown by black contours in all figures is used to establish the absolute metallicity gradient slope (i.e., in units of dex per proper kiloparsec). From the source-reconstructed morphology, we measure the effective radius where the enclosed mass reaches half the total mass of the source. The measurements are represented by R_eff in Table 1.

Figure 3 shows the 2D maps of metallicity of our selected two dwarf galaxies at z ∼ 2. Clearly, the outskirts of our galaxies display highly elevated oxygen abundance ratios. In particular, the outskirts of ID 03751 are more metal-enriched by ∼0.4 dex (i.e., a factor of 2.5) than its center, and more metal-rich by ∼0.2 dex than the value inferred based on the fundamental metallicity relation (FMR) given its integrated ${M}_{* }$ (Mannucci et al. 2010, 2011). Note that our metallicity measurements extend beyond the source effective radius to cover a large enough dynamic range, but not into the region where a plateau/flattening in metallicity (i.e., at $R\gt (2-2.5){R}_{{\rm{eff}}}$, Sanchez et al. 2014; Sánchez-Menguiano et al. 2016) is likely to occur, which might bias the overall gradient determination. For the first time, we are able to detect strongly inverted metallicity gradients in z ∼ 2 dwarf galaxies at unprecedentedly high confidence: 0.122 ± 0.008 dex kpc⁻¹ for ID 03751 (∼15.2σ) and 0.111 ±0.017 dex kpc⁻¹ for ID 01203 (∼6.5σ).

Zoom In Zoom Out Reset image size

The question is thus: what caused these dwarf galaxies to have such strongly inverted gradients? First of all, our sources show no evidence of major mergers, supported by their regular morphology displayed in the 2D maps of ${M}_{* }$ and EL surface brightness in Figure 1. For source ID 01203 with OSIRIS data, this statement is further strengthened by the kinematic evidence of orderly rotation in the disk. Second, the fact that the outskirts of our sources show elevated metallicity as compared to the FMR expectations indicates that there are more metals in the outer regions than could be produced by the stars in those regions. This discourages any explanations involving solely low-metallicity gas inflows, not limited to those induced by mergers. In the subsequent sections, we thus gather all available pieces of observational evidence to further investigate the possible cause.

To understand the cause of the strongly inverted metallicity gradients seen in these dwarf galaxies, we combine their EL maps with HST broadband photometry to derive 2D maps of ${M}_{* }$, SFR, stellar population age, and gas surface density for each galaxy. The SFR is derived from extinction-corrected Balmer emission line flux. Maps of Hβ and Hγ emission are shown in Figure 1. The Hβ/Hγ line ratio provides a measurement of nebular extinction although it is limited by the modest S/N of Hγ. We obtain more precise results from HST photometry, by converting ${B}_{435}$ − ${I}_{814}$ color maps to spatial distributions of stellar reddening E_S(B − V) (Daddi et al. 2004). Nebular reddening E_N(B − V) is then calculated following Valentino et al. (2017). The nebular reddening maps of both our galaxies show lower dust attenuation in centers than that in outskirts, consistent with the inverted metallicity gradients shown in Figure 3.

We calculate extinction in Hβ adopting a dust extinction law of Cardelli et al. (1989) (with ${R}_{{\rm{V}}}$ = 3.1) and assuming Case B recombination with Balmer ratios appropriate for fiducial H ii region properties (i.e., Hα/Hβ = 2.86). Finally, we convert intrinsic Hα luminosity to SFR through the commonly used calibration (Kennicutt 1998a),

$$\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{ 2 }$$

appropriate for the IMF of Chabrier (2003). This provides the instantaneous star formation rate on ∼10 Myr timescales; we note that the ultraviolet continuum probed by HST photometry is sensitive to recent SFR over a longer time span (∼100–300 $\mathrm{Myr}$). The short timescales probed by Balmer emission are most relevant for determining outflow physical properties, which are highly dynamic on small spatial scales, e.g., at sub-kpc level.

Next we derive maps of average stellar age, using the spatial distribution of EL EWs as the primary constraint. We calculate Hβ rest-frame EWs from our maps of the emission line flux and stellar continuum flux density. Stellar continuum maps are corrected for emission line contamination as described in Section 2. We correct for stellar Balmer absorption, which we estimate to be rest-frame EW ∼ 3 Å in Hβ based on the derived galaxy properties (Kashino et al. 2013). Maps of Hβ EW are then converted to average stellar age using a series of Starburst99 stellar population synthesis models (Leitherer et al. 1999; Zanella et al. 2015) and assuming one-fifth solar metallicity and constant star formation history.

We also compare the age estimates given by our SED fitting (Section 2) and Hβ rest-frame EW using the method described above. The median values given by the former practice are systematically larger than those of the latter by ∼0.5 dex, but we note that the uncertainties by the SED fitting are usually much larger due to the absence of prominent continuum spectral age indicators, e.g., ${D}_{n}(4000)$ and ${\rm{H}}{\delta }_{A}$ (Kauffmann et al. 2003). Hence, we adopt the results from Hβ rest-frame EW as the average age for stellar populations throughout our paper, because we consider this a more reliable estimate.

Finally, we calculate the gas fraction defined as

$$\begin{eqnarray}&&{f}_{\mathrm{gas}}={{\rm{\Sigma }}}_{\mathrm{gas}}/\left({{\rm{\Sigma }}}_{\mathrm{gas}}+{{\rm{\Sigma }}}_{* }\right).\end{eqnarray} \tag{ 3 }$$

Since we do not directly observe the bulk of interstellar gas, we instead estimate gas surface density ${{\rm{\Sigma }}}_{\mathrm{gas}}$ by inverting the Kennicutt–Schmidt (KS) law (Schmidt 1959; Kennicutt 1998b), i.e., ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{N}$ together with our measurements of ${{\rm{\Sigma }}}_{\mathrm{SFR}}$ described above. We adopt the more robust extended version of the KS law developed by Shi et al. (2011, 2018), which is especially useful in low density regimes:

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-2}}={10}^{-4.76}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{* }}{{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{0.545}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{gas}}}{{M}_{\odot }{\mathrm{pc}}^{-2}}\right)}^{1.09}.\end{eqnarray} \tag{ 4 }$$

This extended KS law has been tested in numerous ensembles of galaxies as well as regions of low surface brightness in individual galaxies, and is shown to have relatively small scatter (∼0.3 dex) over a large dynamic range of gas and SFR surface densities. We have combined in quadrature this systematic uncertainty of 0.3 dex in our estimates of ${{\rm{\Sigma }}}_{\mathrm{gas}}$.

Figure 4 shows the derived 2D maps of SFR, average stellar age, and gas fraction. In general, we observe centrally concentrated star formation, with the most actively star-forming regions having surface densities $\gtrsim 10\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{kpc}}^{-1}$. On average, the central regions also have older stellar populations and smaller gas fractions than the outskirts, indicating that the outer regions are still in the early stages of converting their gas into stars. These features together indicate that we are witnessing the rapid build-up of galactic disks through in situ star formation and strongly support an inside-out mode of galaxy growth (Jones et al. 2013; Nelson et al. 2014).

Zoom In Zoom Out Reset image size

We compare our radially averaged ${f}_{\mathrm{gas}}$ and metallicity measurements against the predictions from the simple chemical evolution model developed by Erb (2008). To separate the effects of gas inflows and outflows, we compute two extreme sets of models, one being pure gas accretion (i.e., with no outflows, ${f}_{o}={\rm{\Psi }}/\mathrm{SFR}=0$) and the other corresponding to the leaky box model (i.e., with no inflows, ${f}_{i}={\rm{\Phi }}/\mathrm{SFR}=0$). The results are shown in Figure 5. We note that for the pure gas accretion scenario, ${f}_{\mathrm{gas}}$ cannot decrease beyond a certain value, i.e.,

$$\begin{eqnarray}&&{f}_{\mathrm{gas}}^{\min }=1-\displaystyle \frac{1-R}{{f}_{i}-{f}_{o}},\end{eqnarray} \tag{ 5 }$$

where f_o = 0 and R is the instantaneous return fraction. This ${f}_{\mathrm{gas}}^{\min }$, implicitly imposed by Equation (11) of Erb (2008), physically indicates that galaxies cannot exhaust their gas reservoir to below a certain amount without the help of outflows, under the equilibrium condition with steady gas accretion (see Section 3.3 when this equilibrium assumption is relaxed). Therefore, the pure gas accretion scenario cannot explain the observed gas fractions in the central regions of our sources (at ≲2 $\mathrm{kpc}$), where metallicities are also lower. The leaky box model, on the other hand, provides a plausible explanation for our observation such that the outflow rate tends to increase toward the galaxy center. However, we stress that in reality both gas outflows and inflows are acting together to redistribute metallicity. This test using simple chemical evolution models just clearly shows that using gas accretions alone cannot explain our spatially resolved measurements.

Zoom In Zoom Out Reset image size

The application of simple chemical evolution in Section 3.2 is enlightening but depends on strong assumptions, such as that the azimuthal variations are negligible and galaxies exist in equilibrium. In reality, these conditions might not be valid, e.g., due to rapid gas flows. To gain a more precise understanding of the physics of galactic winds and the role of gaseous outflows in shaping the observed spatial distribution of metallicity, independent of those assumptions, we can turn to a more advanced framework for galaxy chemical evolution: the gas regulator model (Lilly et al. 2013; Peng & Maiolino 2014). This model provides an informative and coherent view of the full baryon cycle, involving the accretion of underlying DM halos, as well as the instantaneous regulation of star formation by a time-variable gas reservoir. A key feature of this model is that it does not assume that galaxies exist in an equilibrium state, where the total gas mass remains constant. The non-equilibrium flexibility is especially important for applying this model to spatially resolved regions within a galaxy, where gas may be transported radially from one region to another. Chemical evolution within the gas regulator model is described by the equations

$$\begin{eqnarray}\begin{array}{rcl}{Z}_{\mathrm{gas}} & = & \left[{Z}_{0}+y{\tau }_{\mathrm{eq}}\epsilon \left(1-\exp \left(-\tfrac{t}{{\tau }_{\mathrm{eq}}}\right)\right)\right]\\ & & \times \left[1-\exp \left(\displaystyle \frac{-t/{\tau }_{\mathrm{eq}}}{1-\exp (-t/{\tau }_{\mathrm{eq}})}\right)\right],\\ {\tau }_{\mathrm{eq}} & = & \displaystyle \frac{1}{\epsilon (1-R+\lambda )}.\end{array}\end{eqnarray} \tag{ 6 }$$

Here we adopt the convention of symbols itemized in Table 1 of Peng & Maiolino (2014): ${Z}_{\mathrm{gas}}$ is the mass fraction of metals in the gas reservoir (determined from the observed $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ as in Peeples & Shankar 2011), t is the average stellar population age, ${\tau }_{\mathrm{eq}}$ is the timescale on which the baryon cycle reaches equilibrium, is the efficiency (defined as $\epsilon \equiv \mathrm{SFR}/{M}_{\mathrm{gas}}={{\rm{\Sigma }}}_{\mathrm{SFR}}/{{\rm{\Sigma }}}_{\mathrm{gas}}$), and λ is the mass loading factor (defined in terms of the mass outflow rate Ψ, such that $\lambda ={\rm{\Psi }}/\mathrm{SFR}$ ¹⁷ ). We adopt a stellar nucleosynthesis yield y = 0.003 (Dalcanton 2007) with R = 0.4 estimated from stellar population models of Bruzual & Charlot (2003). Finally, we assume that gas inflows are pristine (Z₀ = 0).

For each spatial region where we have estimated the metallicity, SFR, gas surface density, and age (Figures 3 and 4), we solve the above equations for the mass loading factor λ and subsequently calculate the mass outflow rate Ψ. The 2D distribution of Ψ is displayed in Figure 6. Taking the gradient field of this gaseous outflow map, we obtain the net direction of the outflowing mass flux on sub-galactic scales, projected along the line of sight, denoted by the red arrows in Figure 6. The results demonstrate that strong galactic winds transport mass from the center to the outskirts, with the net radial transport of heavy elements causing the inverted gradients observed in our targets.

Zoom In Zoom Out Reset image size

The distribution of mass loading factors λ within each of our targets is also shown in Figure 7, revealing higher λ (and therefore a higher fraction of metals lost) in the central regions. This preferential removal of metals from the center, and subsequent deposition at larger radii, gives rise to the strong, positively sloped metallicity gradients evident in Figure 3. The high values of λ have important implications for the role of feedback in galaxy formation. Most fundamentally, our results support feedback as a solution to the "overcooling" problem in galaxy formation, by ejecting gas and preventing overly condensed baryonic regions at high redshifts (White & Rees 1978; Dekel & Silk 1986). Such strong outflows are also expected to suppress the formation of stellar bulges from gas with low angular momentum (Governato et al. 2010; Brook et al. 2012). This is consistent with low bulge fraction in these two galaxies measured from high-resolution HST imaging (Table 1).

Zoom In Zoom Out Reset image size

A key feature in the λ distribution is that neither of the wind modes, driven by momentum or energy conservation, can explain the behavior of the mass dependence of λ alone, within individual galaxies. Outflows are typically parameterized by either a momentum-driven (Oppenheimer & Davé 2006, 2008) or an energy-driven (Springel & Hernquist 2003) wind mode, both of which are physically well motivated (Murray et al. 2005). The energy-driven wind scenario assumes that outflows are launched by the thermal pressure of supernova explosions and/or winds from massive stars. A portion of this thermal energy provides the outflow kinetic energy, i.e., ${\rm{\Psi }}\times {v}_{\mathrm{wind}}^{2}\,\sim \mathrm{SFR}$, where the wind speed v_wind can mimic the escape velocity from the DM halo, i.e., ${v}_{\mathrm{esc}}\sim {M}_{{\rm{h}}}^{1/3}$ given by the virial theorem. This results in the scaling relation $\lambda \propto {M}_{* }^{-2/3}$, assuming the linear correlation between the mass constituents of stellar and dark components. The energy-driven wind model is found to be successful in explaining the low abundance of satellite galaxies in the Milky Way (Okamoto et al. 2010). The momentum-driven wind model instead relies on the momentum injection deposited by radiation pressure from supernova explosions and/or massive stars, leading to ${\rm{\Psi }}\times {v}_{\mathrm{wind}}\sim \mathrm{SFR}$ and $\lambda \,\propto {M}_{* }^{-1/3}$. In this scenario, v_wind is proportional to ${M}_{* }$ and SFR, broadly consistent with some observational results (Martin 2005). The transition from energy- to momentum-driven winds is typically thought to be a galaxy-wide phenomenon, resulting in the steepening of the mass–metallicity relation below ${M}_{\ast }\simeq {10}^{9.3}\,{M}_{\odot }$ at $z\simeq 2$ (Henry et al. 2013b). However, our analysis indicates that a single mode is not sufficient to describe spatially resolved data within one galaxy, and it is highly likely that the transition from energy- to momentum-driven winds occurs on sub-galactic scales, governed by local gas and star formation properties in addition to the global gravitational potential.