Introduction

As stars form and evolve in galaxies, they pollute their surroundings with metals. One might expect that the radial distribution of gas-phase metals (or, the abundance ratio of metals to nonmetals, i.e., the metallicity) will follow the radial distribution of past generations of star formation. At 1  z  3, galaxies on average have negative radial gradients in star formation rate surface density (e.g., Nelson et al. 2016a). In a simple closed-box with no radial transport of metals and no external sources acting to enrich or dilute the interstellar gas, local enrichment of the interstellar medium (ISM) through star formation should lead these galaxies toward negative radial gradients in gas-phase metallicity too-i.e., higher metallicities in their centers than in their outskirts.

This description illustrates an important point, but is incomplete. Galaxies form and evolve in a complex ecosystem, with gas cycling in, out, and around galaxies on rapid timescales (Somerville & Davé 2015). This cycle includes (1) stellar feedback driving metal-enriched winds out into the circumgalactic medium with the potential for re-accretion at later times (Tumlinson et al. 2017), (2) metal-poor gas accretion onto galaxies from intergalactic filaments (Kereš et al. 2005), and (3) the rapid shuffling of metals through galaxy mergers (Rupke et al. 2010). These, and other actors, can effectively redistribute metals on galaxy and halo scales. Indeed, by the present day, galaxies are estimated to have retained only ∼20% of the metals they produced over their lifetime (Peeples et al. 2014)-the rest presumably lost to the circum-and intergalactic medium.

These processes are more prevalent at earlier times. At the peak of cosmic star formation at z ∼ 2 (Madau & Dickinson 2014), the rates of star formation of galaxies were 10 times higher (Whitaker et al. 2012), the rates of accretion onto galaxies were 10-30 higher (Scoville et al. 2017), and the rates of galaxy-galaxy mergers were 10 times higher (Rodriguez-Gomez et al. 2017) than they are today.

Numerical simulations reveal that these processes can flatten gas-phase galaxy metallicity gradients-through, e.g., star formation and strong stellar feedback (Gibson et al. 2013;Ma et al. 2017), fountain flows (Grand et al. 2019), and mergers (Rupke et al. 2010;Perez et al. 2011;Torrey et al. 2012). As discussed at the outset, we expect that the simple continuous shedding of metals into the local ISM through stellar evolution will continuously lead galaxies toward negative metallicity gradients. Observations of galaxies with flat or positive gradients (i.e., excursions from this simple expectation) lend insight into the prevalence and timescales of the processes that redistribute metalrich and metal-poor gas listed above. The demographics of gasphase metallicity gradients over time (e.g., the population mean and scatter) provide an important benchmark for theoretical models of galaxy formation.

With the introduction of sensitive near-infrared spectrographs on the Hubble Space Telescope and 10 m class groundbased facilities, the past decade has seen a rapid increase in observations of gas-phase metallicity gradients in high-redshift galaxies (Maiolino & Mannucci 2019). Such measurements are now available for hundreds of galaxies at intermediate redshifts (0 < z < 1; e.g., Frye et al. 2012;Stott et al. 2014;Gupta et al. 2016;Carton et al. 2018;Patrício et al. 2019) and high redshifts (1 < z < 4; e.g., Cresci et al. 2010;Yuan et al. 2011;Queyrel et al. 2012;Swinbank et al. 2012;Jones et al. 2013Jones et al. , 2015;;Leethochawalit et al. 2016;Wuyts et al. 2016;Molina et al. 2017;Wang et al. 2017;Förster Schreiber et al. 2018;Curti et al. 2019;Wang et al. 2019Wang et al. , 2020;;Gillman et al. 2021)-with samples that are large enough to form meaningful conclusions about galaxy populations.

Up to z ∼ 2, high-mass galaxies (M * ∼ 10 10 -10 11 M e ) tend to have slightly negative or flat metallicity gradients (e.g., Wuyts et al. 2016;Carton et al. 2018;Förster Schreiber et al. 2018). Until recently, however, little information on the lowmass galaxy population (M * < 10 10 M e ) above z > 1 was available. The main limitation at these masses and redshifts is the small angular sizes-poor resolution leads to artificially flat metallicity gradients (Yuan et al. 2013;Carton et al. 2017;Acharyya et al. 2020). This poses a significant challenge for ground-based seeing-limited spectrographs, where the typical resolution is ∼ 0 5-1 0.

In the past few years, however, the low-mass population at these redshifts has started to be explored with large numbers (N ∼ 100) through deep galaxy surveys taking advantage of the magnification afforded through gravitational lensing. These include recent efforts with the Hubble Space Telescope/Wide Field Camera 3 slitless spectrograph (Wang et al. 2017(Wang et al. , 2019(Wang et al. , 2020) ) and the Very Large Telescope/KMOS integral-field spectrograph (Curti et al. 2019). These surveys reveal a zoo of flat, positive, and negative gradients -with a majority of galaxies having flat gradients. These highredshift galaxies are in marked contrast with today's star-forming galaxies, of which the majority show declining gas-phase metallicity gradients (e.g., Belfiore et al. 2017).

In this paper, we use deep Hubble slitless grism spectroscopy from the Wide Field Camera 3 to study the metallicity gradients of 238 galaxies at 0.6 < z < 2.6 over a stellar mass range of 8.5 < log M * /M e < 10.5-effectively doubling the low-mass sample size at these redshifts. Our goal is to provide a statistical understanding of the gas-phase metallicity gradients at high redshift, and to examine the population scatter as a function of stellar mass.

Our outline for this paper is as follows. In Section 2, we detail the observations, data reduction, and sample used in this paper. In Section 3, we describe the measurements of the physical properties of our sample, including the metallicity gradients. Next, in Section 4, we examine the gas-phase metallicity gradients as a function of galaxy stellar mass and assess the intrinsic scatter of the population. We also explore correlations between the metallicity gradient and stellar population properties. We use these results to place constraints on predictions of correlation strengths from theoretical models. In Section 5, we discuss our findings and, finally, in Section 6 we summarize the results of the paper.

Throughout this paper, we adopt a ΛCDM cosmology with the Planck 2015 cosmological parameters (h, Ω m , Ω λ ) = (0.67, 0.31, 0.69). For relevant derived quantities (e.g., stellar masses, star formation rates), we assume a Chabrier (2003) initial mass function (IMF). The FWHM spatial resolution of the Hubble Space Telescope observations used in this paper span 0 12-0 15. The resolution increases toward the redder end of the near-infrared spectral range of the Wide Field Camera 3 (0.8-1.7 μm). This corresponds to a physical resolution of ∼1-1.3 kpc at the median redshift of our sample (z = 1.6).

Observations, Data Reduction, and Sample Selection

In this section, we describe the observations (Section 2.1), data reduction/spectral fitting (Section 2.2), and sample selection (Section 2.3) used in this paper.

Hubble Slitless Spectroscopy

The Hubble Wide Field Camera 3 (HST/WFC3) slitless spectra used in this paper are taken from a collection of programs. The root program is the Cycle 23 CANDELS Lyα Emission at Reionization survey (CLEAR; GO-14227, PI Papovich). The CLEAR footprint covers 12 pointings in the GOODS-S deep and GOODS-N deep CANDELS fields (Grogin et al. 2011;Koekemoer et al. 2011). Each pointing includes 10-12 orbit depth (10 for GOODS-N and 12 for GOODS-S) HST/WFC3 G102 grism spectroscopy and companion HST/WFC3 F105W direct imaging. Each pointing is observed at 3 orients, separated by >10°-helping relieve source confusion in the grism spectra. A primary motivation for the CLEAR survey was to constrain the evolution of the distribution function of Lyα emission at z > 6.5. Such constraints from CLEAR data will be presented in a forthcoming paper (I. Jung et al. 2021, in preparation). Previous work has used the CLEAR spectra to study the metallicities, ages, and formation histories of massive high-redshift galaxies (Estrada-Carpenter et al. 2019, 2020), and to appraise Paschen-β as a star formation rate indicator in low-redshift galaxies (Cleri et al. 2020).

Given their location in the well-studied GOODS-S and GOODS-N fields, the CLEAR observations are supported with extensive ancillary photometry spanning the ultraviolet to the nearinfrared. We use an augmented version of the v4.1 3D-HST GOODS-S and GOODS-N photometric catalogs (Skelton et al. 2014)-adding Y-band photometry measured from the CLEAR F105W imaging and archival HST/WFC3 F098M and/or F105W imaging. Y-band imaging is available over the majority of the GOODS-S and GOODS-N fields, and the entire CLEAR footprint. The new Y-band photometry is measured in a manner that is consistent with the existing 3D-HST photometric catalog (see details in Estrada-Carpenter et al. 2019 andSkelton et al. 2014). We use the Python version of Eazyfoot_0 (Brammer et al. 2008) to re-derive the photometric zero-points of the augmented catalog.

To extend the spectral coverage and maximize the depth of our grism observations, we query the Mikulski Archive for Space Telescopes (MAST) for all publicly available HST/WFC3 G102 and G141 grism observations (and their associated direct imaging) that overlap with the root CLEAR footprint. In doing so, we retrieve a total of 52 additional orbits with the WFC3/G102 grism and 76 orbits with the redder WFC3/G141 grism. Of note, we retrieve an ultra-deep 40 hr G102 pointing in the GOODS-S/UDF from the Faint Infrared Grism Survey (FIGS; Pirzkal et al. 2017). The FIGS/UDF field comprises the deepest grism data in this collection. Hereafter, we refer to the full collection of grism observations as "CLEAR+."

The archival G141 observations were acquired by the following programs: GO-11600 ("AGHAST"; PI Weiner), GO-12461 ("SN COLFAX", PI Reiss), GO-13871 (PI Oesch), GO/DD-11359 ("ERS", PI: O'Connell), GO-12099 ("GEORGE, PRIMO", PI Reiss), and GO-12177 ("3D-HST", PI van Dokkum). The archival G102 observations are from: GO-13420 (PI Barro), GO/DD-11359 ("ERS", PI O'Connell), and GO-13779 ("FIGS", PI Malhotra). The full collection of G102+G141 observations used in this paper is packaged in https://doi.org/10.17909/t9ctff-wx60.

In Figure 1, we show the redshift windows over which several strong rest-frame optical lines are accessible with the G102 (0. ) can be accessed over just 1.2 < z < 1.3 with the G102 grism and 2.0 < z < 2.4 with the G141 grism. With joint WFC3/G102+G141 spectral coverage, the same redshift windows are significantly widerthe full R 23 complex can be accessed over 1.2 < z < 2.4 when both grisms are employed. Limiting to any two lines of the R 23 complex, the redshift window increases to 0.6 < z < 2.6 with both grisms employed. Furthermore, the spectral overlap of the G102 and G141 grisms leads to continuous coverage over the full spectral range.

Grism Data Reduction and Spectral Extractions

To process the set of G102+G141 grism observations described in the previous subsection, we use the grism redshift and line analysis software Grizlifoot_2 (Brammer 2019). Grizli performs full end-to-end processing of HST imaging and slitless spectroscopic data sets, including: retrieving and pre-processing the raw observations, astrometric alignment, modeling contamination from overlapping spectra, extracting 1D and 2D spectra, fitting full continuum+emission-line models, and generating emission line maps.

First, we use Grizli to retrieve and pre-process the observations described in Section 2.1 from the MAST archive. The raw WFC3 data are reprocessed with the calwf3 pipeline with corrections for variable sky backgrounds as described by Brammer (2016). Cosmic rays and hot pixels not flagged by the pipeline are identified with the AstroDrizzle software (Gonzaga et al. 2012). As in Momcheva et al. (2016), the grism exposures are flat-fielded using the F105W and F140W calibration images for the G102 and G141 grisms, respectively, and grism sky subtraction is performed using the "Master Sky" provided in Brammer et al. (2015). Relative astrometric corrections are applied to the processed data through an alignment to the deeper F140W HST mosaic galaxy catalog from the 3D-HST survey (Skelton et al. 2014).

A contamination model of each CLEAR+ pointing is created using a forward-model of the HST Y-band full-field mosaic. For each spectrum, this model is used to subtract the contamination from adjacent spectra. This is an iterative process. The first pass contamination model is created for all objects in the Y-band segmentation map brighter than m F105W = 25, assuming spectra that are flat in units of f λ flux density and normalized at F105W given the image segment. Next, a refined continuum model is created for objects brighter than m F105W = 24 by fitting thirdorder polynomials to the spectrum of each source after subtracting the models of contaminating sources. These steps result in a full contamination model of the detector for each visit of each observational program.

We use Grizli to extract the 2D grism spectra (i.e., "beams") of objects in the field of view to a limiting magnitude of m F105W < 25. This is performed for each grism visit of each object. The extractions carry the full description of the WFC3 detector and the contamination model. The grism exposure times of the CLEAR+ objects range from 2-30 hr in G102 and 2-10 hr in G141. In total, 6048 objects are extracted from the CLEAR+ data. Of these, 533 have sole coverage with the G102 grism, 808 have sole coverage with the G141 grism, and 4707 have joint coverage with both grisms.

Redshift fits are carried out using the grism spectra and available multiwavelength photometry. The spectra are scaled to the photometry using a first-order polynomial correction. A basis set of template Flexible Stellar Population Synthesis models (FSPS; Conroy et al. 2009;Conroy & Gunn 2010), including emission line complexes with fixed ratios, are used in the fit. The FSPS templates are constructed to span a broad diversity of galaxy types following the methodology described by Blanton & Roweis (2007) and Brammer et al. (2008). At a given trial redshift, the redshifted templates are both convolved with the photometric filter bandpasses and projected to the space of each extracted 2D spectral "beam" using the direct Yband image to define the spatial morphology. This approach takes into account the unique morphological broadening of each galaxy due to the limited spectral resolution of the grism. A final model is determined from a nonnegative linear combination of the model template spectra and the goodness of fit is determined from the combined χ 2 of all photometry and 2D spectral pixels using the uncertainties from the photometric catalogs and exposure-level noise model, respectively. The "best" redshift is taken to be where this χ 2 is minimized across a grid of trial redshifts from 0 < z < 12.

Emission line fluxes are fit at the best-fit redshift using the FSPS basis templates described above, but now including separate components for each line species without fixing their line ratios. 14 These fits are carried out using the forwardmodeling technique described above.

The Grizli-derived redshifts and emission line fluxes are consistent with previous measurements from the 3D-HST survey (Momcheva et al. 2016)-their differences divided by the sum of their uncertainties have a near standard normal distribution in both cases. Comparing to ground-based spectroscopic measurements of sources in the CLEAR fields, we find an overall redshift precision of σ NMAD = 0.0024 in Δz/(1 + z).

Emission line maps are derived for several lines (notably the strong lines listed in Figure 1) by drizzling the contaminationand continuum-subtracted 2D spectral beams using the astrometry of the slitless exposures projected along the spectral trace to the wavelength of the redshifted line center. The line maps are created with a pixel scale of 0 1 and are convolved with a 2 pixel boxcar kernel to increase their signal-to-noise ratio. The effective resolution of the line maps is ∼ 0 16. The uncertainties on the line maps are calculated using the drizzle weights from the constituent beam pixels. Example line maps are shown in the top panels of Figure 4.

Sample Selection

To choose galaxies we can reasonably derive gas-phase metallicity maps for, we select objects with a 5σ integrated detection in at least two of the strong lines in the R 23 complex:

This selection leads to a total of 486 galaxies. Of these, 19% (N = 86) have integrated detections of the full R 23 complex, 39% (N = 190) have detections of only the R 3 complex, less than 2% have detections of only the R 2 complex (N = 9), and 41% have detections of only the O 32 complex (N = 201).

We match the CLEAR sample to the Chandra Deep Field X-ray point-source catalogs (Xue et al. 2016;Luo et al. 2017) and remove galaxies harboring X-ray bright active galactic nuclei (AGNs). AGN ionization is not accounted for in the photoionization models. It will act to increase the O 32 ratio and artificially decrease the inferred metallicities in the centers of galaxies. This selection removes 28 galaxies from the sample.

As discussed in the following section, the sample is further culled through a selection on the radial extent of the detected metallicity signal. This leaves 238 galaxies in our final science sample. The redshift distribution of this sample is shown in the top panel of Figure 1. The selection criterion constrains the redshifts to 0.6 < z < 2.6. The majority of the sample are at the high-redshift end of this range-11%, 37%, and 52% of the sample are at 0.6 < z < 1.0, 1.0 < z < 1.5, and 1.5 < z < 2.6, respectively.

Measurements of Physical Properties

In this section, we describe the measurements of the physical properties of the CLEAR+ sample. In Section 3.1, we discuss derivations of their stellar masses, star formation rates, and dust extinctions, as well as measurements of their sizes. In Section 3.2, we derive HST-resolution metallicity maps of the CLEAR+ galaxies using pixel-by-pixel fits to the emission line maps, and fit the radial gradients of the gas-phase metallicity.

Stellar Mass, Star Formation Rates, Dust Extinction, and Sizes

Stellar masses, star formation rates, and dust extinctions are measured from the ancillary multiwavelength photometry of the full CLEAR+ sample using Eazy-py. The photometry spans 0.3-8 μm (as described in Table 3 of Skelton et al. 2014), including the new HST/WFC3 F105W photometry described above. The fits are carried out with the "fsps_QSF_12_v3" spectral energy distribution (SED) template set available in the Eazy library. These templates are created using FSPS, described in Section 2.2, and assume a Chabrier (2003) IMF.

The circularized effective radii (R eff ) of the CLEAR+ galaxies were measured in van der Wel et al. (2012). We adopt R eff as measured from the CANDELS HST/WFC3 F125W+F160W imaging, with the redshift-dependent correction outlined in van der Wel et al. (2012).

UV+IR-derived star formation rates are taken from Whitaker et al. (2014). These are measured using a conversion from total UV+IR luminosity to total star formation rate (Bell et al. 2005), assuming a Chabrier (2003) IMF. A template-based (Dale & Helou 2002) conversion is used to translate the observed Spitzer/MIPS 24 μm flux density of each galaxy to its total IR luminosity. The total UV luminosity is linearly scaled from the rest-frame 2800 Å luminosity derived from the best-fit SED template. UV+IR star formation rates are available for ∼83% of our sample. For galaxies where it is not available, we use the SED-derived star formation rates described above. We note that the SED-derived star formation rates are, on average, 0.09 dex lower than the UV+IR-derived star formation rates. This systematic difference is less than our adopted uncertainties on the UV+IR star formation rates (0.2 dex), and does not affect the conclusions of this paper.

In Figure 2, we show the distribution of stellar mass, star formation rate, circularized effective radius, and V-band dust extinction for the sample of 238 galaxies. The 16th-50th-84th percentiles for stellar mass are *  M M log : 9.3-9.7-10.2. The same percentiles for the V-band dust extinctions, star formation rates, and effective radii are: 0.09-0.29-0.64, 3-9-34 M e yr -1 , and 0 13-0 24-0 42. In Figure 3, we plot the star formation as a function of stellar mass for this sample. We also include the star formation mass sequence (SFMS) fits from Whitaker et al. (2014). The galaxies in our sample lie along the SFMS at their respective epochs.

Metallicity Maps

We carry out pixel-by-pixel fits of metallicity and ionization parameter to the grism-derived emission line maps using the Bayesian photoionization fitting code IZIfoot_4 (Blanc et al. 2015). IZI uses a grid of outputs from a user-specified photoionization model to fit observed dust-corrected emission line fluxes. We dust-correct the observed emission line maps using the Eazy-derived A V value of each galaxy, assuming a Calzetti et al. (2000) extinction law. The distribution of A V for our galaxy sample is shown in Figure 2-it is highly skewed to low A V .

We do not account for differences between stellar and nebular extinction at these redshifts (Price et al. 2014) and assume no radial dependence on A V . The dust gradients in lowmass galaxies (<10 10 M e ) at the redshifts of our sample (which comprise the majority of our sample) have been measured to be relatively flat on average (Nelson et al. 2016b). We test potential biases introduced by this latter assumption using a suite of simulated profiles. These simulated profiles span a range of intrinsic gas-phase metallicity gradients, with a fixed central A V = 1.5 (the most extreme case for the galaxies in our sample), and an intrinsic gradient in A V of -0.08 dex kpc -1 . The latter is chosen to match the broadband-derived measurements of high-mass galaxies (>10 10 M e ) at z = 2 (Tacchella et al. 2018). Using the same measurement techniques we use for the real galaxies in our sample (outlined below), we find that excluding a dust gradient does not significantly bias the recovered metallicity gradients-there is a ∼5%-30% bias depending on the intrinsic metallicity gradient and the metallicity diagnostics available. The bias is such that the recovered gradient is more positive than the intrinsic gradient for both intrinsically positive and intrinsically negative input gradients. The bias is highest (∼30%) for galaxies with steep intrinsic metallicity gradients, and for the R 23 diagnostic. Taken as a whole, and without uniform empirical constraints for our sample, we consider the choice of a flat dust gradient justified.

We use the MAPPINGS-IV photoionization models (Dopita et al. 2013), with the parameter describing the distribution of electron energies κ set to 20. The MAPPINGS-IV models adopt a functional dependence of the nitrogen-to-oxygen ratio and the oxygen abundance, based on observations of local galaxies (van Zee et al. 1998). Because the gradient is a relative measure, the results presented herein are less sensitive to the normalization differences between photoionization libraries. A comparison between the predictions of the photoionization models and observed line ratios from the CLEAR+ sample appears elsewhere (C. Papovich et al. 2021, in preparation).

For each object, we create a segmentation map where the [O II] line map (or Hβ, if [O II] is not observed) exceeds an S/N of 1. While this S/N threshold is relatively low, the S/N of the Hβ and [O II] line fluxes are generally the lowest of the line fluxes that are used in the photoionization fit. Furthermore, the metallicity gradients (described later) are fit using a large number of pixels and are not generally susceptible to low S/N outliers. Pixels inside this segmentation region are fit with IZI, and pixels outside are masked. For each 0 1 × 0 1 pixel in the unmasked CLEAR line maps, a posterior inference of the metallicity and ionization parameter is recovered. IZI is flexible-it allows for upper limits as well as line sums (from, e.g., unresolved line complexes in the grism). We provide IZI the strong emission lines and line sums listed in Figure 1 The R 2 , R 3 , and R 23 ratios have degenerate solutions with metallicity-a low-metallicity branch and a high-metallicity branch. An ancillary measurement of other line ratios, e.g., O 32 , can help break this degeneracy (see C. Papovich et al. 2021, in preparation, for the application to the CLEAR+ data), but high S/N detections of the [O II] line are not always available. Without such a constraint, the branch favored by the posterior is generally a strong function of the metallicity prior adopted. To avoid internal "branch-switching" in our fits-i.e., a subset of galaxy pixels marginally favoring the lower branch and the remaining pixels marginally favoring the upper branch-we set a prior that favors the high-metallicity branch of the R 2 , R 3 , and R 23 MAPPINGS-IV models. Specifically, we adopt a top-hat prior over 12+log O/H = [8.5, 9.5]. Our results are insensitive to the specific choice of the lower bound of this prior down to 12+log O/H ∼ 8.2-below which the MAPPINGS R 23 and R 3 models turn over to the low-metallicity branch. The choice to favor the high-metallicity branch is supported by the locus of the global gas-phase mass-metallicity relation (MZR) at the extreme ends of our sample-the highest redshifts and lowest masses-where very few objects have integrated gas-phase oxygen abundances that fall on the lower branch derived from the integrated R 23 and O 32 emission lines (C. Papovich et al. 2021, in preparation). 16 This is generally consistent with existing measurements of the MZR over 0.6 < z < 2.6. At the modal stellar mass of our sample (M * ∼ 10 9.5 M e ), the MZR exceeds the lower-branch turnover of the MAPPINGS-IV models (12+log O/H ∼8.2) across our full redshift range: at z ∼ 2.3 (e.g., Sanders et al. 2018Sanders et al. , 2021)), z ∼ 1.5 (e.g., Henry et al. 2013;Wang et al. 2017;Gillman et al. 2021;Topping et al. 2021), and z ∼ 0.6 (e.g., Gillman et al. 2021).

Example emission line maps, and derived metallicity maps, are shown in Figure 4. We adopt the luminosity-weighted center of the direct F105W image as the center of each galaxy, and measure the radial position and gas-phase metallicity of each pixel. We do not de-project the galactic coordinates, as the emission line kinematics of galaxies at this redshift indicate that the emission line structure is rarely disky (Kassin et al. 2012;Simons et al. 2016Simons et al. , 2017)). As such, a de-projection using, e.g., the continuum axis ratio, has an uncertain meaning. Finally, we remove galaxies from the sample if their metallicity is not recovered beyond 0 3. This ensures that the radial profiles in our final sample extend to at least ∼2 HST/WFC3 point-spread function FWHM resolution elements (0 12-0 15).

For each galaxy, we carry out a regression to the radial metallicity profile using a simple line of the form 12 log

where R is the projected radial distance from the galaxy center, m is the metallicity gradient, and b is the central metallicity.

We expect a level of unknown correlated noise in the radial profiles at the level of the spatial resolution of the emission line maps. Ideally, this should be accounted for in the fit, and ultimately reflected in the uncertainties of the fitted model. With that motivation, we use the Gaussian Processes fitting routine Georgefoot_6 (Ambikasaran et al. 2015). George allows for a level of unknown covariance in the radial profile and models the structure of that covariance matrix as a Gaussian Process. We define the shape of the covariance matrix using a stationary exponential-squared kernel, such that the covariance of any two points in the radial profile depends on their relative distance from each other. The length scale of the kernel is set to 0 2, which is the size of our boxcar smoothing kernel and the scale at which we expect nearby pixels to be strongly correlated. We fit the full suite of galaxies in our sample with simple least-square linear regression as well. We find that there is not a systematic difference between the Gaussian Processes and least-squares regression fits, but the Gaussian Processes fits have on average 45% higher uncertainties.

Example fits are shown in Figure 5. A truncated table of the best-fit parameters is shown in Table 1. A figure set (the Appendix; Figure 12) is provided with the online version of this article that includes for each galaxy: (1) the emission line maps and their associated S/N maps, (2) the integrated galaxy spectrum with the emission lines marked, (3) the derived metallicity map and its upper and lower uncertainty maps, and (4) the radial profile of the metallicity with its best-fit linear model.

Metallicity Gradients

In this section, we present the main results of the paper. In Section 4.1 and Figures 6 and7, we show the gas-phase metallicity gradients of the galaxies in our sample as a function of stellar mass. We compare that with the z ∼ 0 galaxy population and highlight the inferred population evolution. In Section 4.2 and Figure 8, we assess the intrinsic scatter (and its mass dependence) of the metallicity gradients. Finally, in Section 4.3 and Figures 9 and 10, we assess correlations between the metallicity gradient and various galaxy properties.

As a Function of Stellar Mass

In Figure 6, we show the gas-phase metallicity gradients of the galaxies in our sample as a function of their stellar mass.

In the left panel, we show individual galaxies from the CLEAR+ sample. Across all masses, the CLEAR+ galaxies are generally consistent with a flat or slightly positive gradient. Roughly 82% of the sample (197/238 galaxies) are 3σ consistent with a flat gradient (we define a "flat" gradient as objects where R log O H 0 ( ) D D = is within 3σ) and 16% of the sample (38/238 galaxies) are 3σ consistent with a positive gradient. Taken together, nearly 99% of the sample (235/238) are consistent with a flat or positive gradient and only 1% of the sample (3/238 galaxies) are consistent with a negative gradient. With 1σ and 2σ confidence, respectively, 87% and 95% of the sample are consistent with a flat or positive gradient.

We also include a collection of recent measurements from the literature from the following surveys at 0.8  z  2. ) slightly favor a positive gradient. As discussed later, the population scatter marginally (but with statistical significance) increases toward lower stellar masses.

In Figure 7, we compare our median population trends of the metallicity gradient (normalized by their circularized effective The examples include a galaxy that is more metalrich in its center than in its outskirts (i.e., a galaxy with a negative metallicity gradient; left) and a galaxy that is more metal-poor in its center than in its outskirts (i.e., a galaxy with a positive metallicity gradient; right). A gray x marks the center of each galaxy and a 0 3 scale is included for reference. The white bar extends 0 3 in all panels. The resolution of the emission line maps is indicated with the white circle. The panels span 2″ × 2″. A figure set of the emission line maps and derived metallicity maps for the full sample are provided in the online version of this article (the Appendix; Figure 12). As an important note, the z ∼ 1.5 CLEAR+ galaxies are generally star-forming and thus still growing in stellar mass. For instance, the galaxy population with log M * /M e ∼ 9 at z = Example fits to the radial gas-phase metallicity profile for three galaxies in the CLEAR sample are shown, including: a galaxy with a flat metallicity gradient (left column), a galaxy with a positive gradient (middle column), and a galaxy with a negative gradient (right column). An HST/ACS WFC3 color map is shown in the top row. The gas-phase metallicity map is shown in the middle row, with the strong line diagnostics (including detections and upper limits) that are used in the photoionization fit indicated. The radial metallicity profile is shown in the bottom row. In the top and middle rows, a gray square marks the center of the galaxy, as defined by the center-of-light of the HST/WFC3 F105W and F125W image, respectively. A white bar extending 0 3 is shown for scale. Random draws from the posteriors of the linear fits to the metallicity profile are shown in black in the bottom row. A figure set of the metallicity maps and radial metallicity profiles for the full sample are provided in the online version of this article (the Appendix; Figure 12). galaxy with log M * /M e = 10 at z = 1.5 can double its effective radius from 2.1 kpc to ;4 kpc by z = 0 (e.g., Papovich et al. 2015). As a consequence, the true evolutionary tracks of the z ∼ 1.5 galaxy populations include an increase in stellar mass coincident with a decline in metallicity gradient (i.e., toward the lower right in Figure 7)-the inferred evolution is stronger than is indicated by the population differences at fixed mass.

Intrinsic Population Scatter

In Figure 8, we assess the intrinsic scatter of our sample and the collated literature sample as a function of stellar mass. To do so, we first measure the observed scatter using a running standard deviation of the sample. The running variable is stellar mass and we use a running width of 0.5 dex. We report the result every 0.25 dex. The conclusions below are relatively insensitive to the specific intervals chosen. To calculate the excess scatter, we measure the running median of the observational uncertainties of the sample and subtract it from the running standard deviation. We use bootstrap resampling to measure the standard error on the difference. At face value, this technique recovers scatter that cannot be attributed to the observational uncertainties.

At all masses, we measure nonzero excess scatter, which we attribute to the intrinsic scatter of the population. The intrinsic scatter ranges from ∼0.4-0.8 dex kpc -1 and it continuously increases with decreasing stellar mass. This is consistent with Wang et al. (2020), which also reports an inverse correlation between stellar mass and intrinsic scatter at these redshifts.

To assess the slope of the relation between the intrinsic population scatter and stellar mass, we carry out a least-squares linear fit. We do this using: (1) the CLEAR+ sample, (2) the collated literature sample collected in the previous subsection, and (3) a combined CLEAR+ and literature sample. We report 3.9σ evidence for a slope in the CLEAR+ sample Note. Linear fits to the radial profiles of the gas-phase metallicity for a subset of the galaxies used in this paper. The fits take the simple form log(O/H) (R) = m × R + b, where R is the projected radial distance from the galaxy center. The IDs are matched to the 3D-HST photometric catalog (Skelton et al. 2014).

(This table is available in its entirety in machine-readable form.) (m = -0.0120 ± 0.0031 kpc -1 ) and 2.2σ evidence in the literature sample (-0.0255 ± 0.0114 kpc -1 ). In the combined literature and CLEAR+ sample, the evidence for a slope increases to 5.2σ (-0.0252 ± 0.0049 kpc -1 ). In summary, we report excess (intrinsic) scatter in the gas-phase metallicity gradients at all masses (8.5 < log M * /M e < 10.5) and find strong evidence for a dependence of the intrinsic scatter on stellar mass-wherein higher-mass galaxy populations have lower intrinsic scatter.

Correlations with Galaxy Properties

To better understand the physical mechanisms responsible for the results in the previous two subsections, we assess correlations of the observed metallicity gradients with various galaxy properties.

In Figure 9, we use our sample to explore the relationships of metallicity gradient and several stellar population properties. Specifically, we compare the metallicity gradients of our sample to their circularized effective sizes, star formation rates, and average star formation surface densities (which we define as SFR/π R eff 2 ). We also compare the metallicity gradients to the star formation rate per galaxy gravitational potential energy, which is defined as SFR/( * M 2 R eff ). This latter quantity is a proxy for the ratio of available kinetic energy (from feedback associated with the SFR) to the gravitational binding energy. It allows us to test if this ratio impacts the ability of a galaxy to redistribute metals (produced from star formation) to the rest of the galaxy, eject them into the circumgalactic medium, or unbind the gas entirely from the galaxy. The quantities are relevant because there are observed correlations between the velocities of star formation driven galaxy winds (usually metal enriched) and the local star formation surface density (Kornei et al. 2012;Heckman et al. 2015). The presence (or absence) of a correlation with the last two parameters lends insight into the impact of star formation winds on metallicity gradients, and/or the timescales over which they leave an imprint in the observations.

To remove the mass and redshift dependence of these properties, we adopt their value relative to the population average at their mass and redshift. To do so, we use the masscircularized size relation of late-type galaxies from van der Wel et al. ( 2014) and the SFMS from Whitaker et al. (2014).

We measure the Pearson coefficient for each correlation and report the value (and an interpretation) in the top right of each panel. We report no evidence for a correlation between metallicity gradient and any of these mass-normalized physical properties.

As an aside, we do note that the galaxies with the highest star formation rate surface densities (those that are a dex above the population average) tend to have more positive metallicity gradients than the rest of the sample (bottom left panel, Figure 9). A plausible physical explanation that can account for this result and the measured null correlation is that the star formation rate surface densities need to reach a certain threshold before they are able to effectively launch metal-rich winds (see, e.g., Heckman et al. 2015). If such winds are preferentially launched from the galaxy center, then they could presumably drop the central metallicity and push the metallicity gradient in the positive direction. We report significant evidence for an increase in intrinsic scatter with a decrease in stellar mass. The distribution of stellar mass for both the CLEAR+ sample and the collection of galaxies adopted from the literature (Swinbank et al. 2012;Jones et al. 2013;Wang et al. 2017Wang et al. , 2019Wang et al. , 2020;;Curti et al. 2020) is shown in the top panel.

Figure 9 shows the results for the galaxies in our sample. Again, we find no statistically significant evidence for a correlation between the metallicity gradients and these galaxy properties. However, even these null/weak correlations place important constraints on predictions from theoretical models. We discuss this in more detail below.

In Figure 10, we estimate the intrinsic correlations that are ruled out by the CLEAR observations-specifically, those that are inconsistent with the null observed correlations between the metallicity gradient and each of the galaxy properties above.

To do so, we simulate "CLEAR-like" realizations of a galaxy population assuming there is some intrinsic correlation between the metallicity gradient and each galaxy parameter. The correlations are defined by two parameters: the intrinsic slope of the correlation, and the intrinsic vertical scatter of the correlated variables. For each simulated intrinsic correlation (i.e., each pixel in the Figure 10 heat maps), we generate 100 sample realizations and then add noise to match the quality of the CLEAR+ measurements. Specifically, we draw a mock sample using the intrinsic correlation parameters with the same dynamic range of metallicity gradients as the CLEAR+ sample, and the same number of galaxies. We add simple uncorrelated Gaussian noise to mimic the observational uncertainties. Finally, we measure the Pearson correlation coefficient for each realization.

In Figure 10, we show the fraction of simulated realizations that are as (or more) uncorrelated than the observed CLEAR sample. To interpret these diagrams, the parameter space where the fractions are high are consistent with the CLEAR observationse.g., a fraction of 80% means that 80% of the realizations are as uncorrelated as the observed correlations. The parameter space where the fractions are low are strongly disfavored by the CLEAR results. The region of the diagram outside the 5% contour shown in the diagrams is ruled out at 2σ confidence or more. The results shown in Figure 10 offer a direct constraint on theoretical predictions as our data rule out these regions of the parameter space at a 95% confidence. . Correlations between gas-phase metallicity gradient and various galaxy properties: the galaxy star formation rate (top left), the circularized effective radius (top right), the galaxy star formation rate surface density (bottom left), and a proxy for the star formation rate per unit gravitational potential energy (bottom right). All of the stellar population measurements have been normalized by the population average at the mass and redshift of the galaxy. We report null correlations with all of these galaxy properties. The typical uncertainty is denoted by the error bar shown in each panel.

Discussion

Metals are a dye for galaxies. They are formed in stars, deposited around young star-forming regions, and swept up in supernovae (SNe) or massive star driven winds. As gas moves around galaxies, so too will its entrained metals (or, equivalently, its lack of metals). If the spatial distribution of the gas-phase metallicity strays from the distribution of a galaxy's stars (an integrated record of past star formation and metal enrichment), it indicates one (or both) of the following:

(1) the metals that formed in the galaxy were re-distributed from their birthplace or (2) the ISM of the galaxy was unevenly diluted through metal-poor gas accretion. In either case, the evolution of the metallicity gradient is intimately linked with galaxy-and halo-scale gas flows.

In this section, we develop and test a toy model to highlight the rapid evolutionary timescales implied by the observed metallicity gradients (Section 5.1 and Figure 11) and discuss various processes known to flatten metallicity gradients (Section 5.2).

A Toy Model: Timing the Development of Metallicity

Gradients at z ∼ 1

The fact that a large fraction of galaxies have flat and positive metallicity gradients at 0.6 < z < 2.6 (Section 4.1)

Figure 10. Correlation parameters favored-and regions of the parameter space that are "ruled out"-by the CLEAR data. In each region (pixel) of the plot, we simulate data sets of a fixed intrinsic correlation slope and vertical scatter. We then add noise to match the quality of the CLEAR data (as described in the text) and measure the fraction of those realizations that are less correlated than the real data. The purple region outside of the contour is strongly disfavored by the CLEAR data (at 95% confidence: less than 5% of the realizations fall in these regions). Specifically, this region is inconsistent with the null correlations measured between the metallicity gradients and these galaxy properties in Figure 9. strongly favors a scenario in which metallicity gradients are destroyed in galaxies on short timescales-quicker than the short time it takes for star formation to develop an (observable) declining gradient. In Figure 11, we show the results of an empirical toy model that illustrates this point.

Constructing the Model. We adopt a set of empirical prescriptions to evolve the radial metallicity profiles of mock galaxy populations. For simplicity, these models assume metal production from star formation, no radial mixing of metals, and an instantaneous replenishment of the gas used to form stars. The empirical prescriptions include (1) the average radial profiles of star formation rate surface density at z ∼ 1 measured from the 3D-HST survey (Nelson et al. 2016a), (2) the global MZR at z ∼ 2.3 (Sanders et al. 2018), and (3) the Kennicutt-Schmidt (KS) law (Kennicutt 1998) relating the star formation and gas mass surface densities. The first dictates the rate at which stars form per galactic annulus and, thus, the rate at which metals are deposited into the local ISM. The second dictates the global normalization of the metallicity at the start of the simulation. The third dictates the differential increase in metallicity per unit of star formation. If the KS relation has a power-law slope that exceeds 1, it means that stars are forming more efficiently at higher gas densities. As a result, star formation would more efficiently increase the local ISM metallicity at higher gas densities. This dependence of the rate of change of metallicity (a subtle difference from the absolute metal production rate) on gas density is key to the model.

The KS law ( n SFR gas S µ S ) at high surface densities at these redshifts has a slope between unity and 2 (Narayanan et al. 2011;Tacconi et al. 2013). We adopt a probabilistic broken power law for the KS relation in our model, reflecting a regular and starburst sequence. The break occurs at 100 M e pc -2 , above which the relation assumes a flatter slope. The slope above and below the power-law break are drawn from a Gaussian distribution with a mean of n = 1.1 and n = 1.4, respectively, and a 1σ width of 0.1. We adopt a probabilistic metal yield from type II SNe, using a Gaussian with a mean of y Z,ii = 0.03 and a 1σ width of 0.005. We simulate four galaxy populations of log M * /M e = [9. 0-9.5, 9.5-10.0, 10.0-10.5, 10.5-11.0]. Lastly, we assume that the ratio of oxygen to total metals is constant, such that the change in the oxygen abundance gradient is equivalent to the change in the metallicity gradient. For each population, we simulate 100 realizations-each using a random draw from the distributions of KS slope and metal yield. We evolve each simulation for 1 Gyr. The initial metallicity gradient is set as flat with a global metallicity matching the z ∼ 2.3 MZR (Sanders et al. 2021).

This is a simple empirical model, and should not be confused with more sophisticated chemical evolution models (e.g., Spitoni & Matteucci 2011;Mollá et al. 2019;Henriques et al. 2020). With that said, it illustrates an important point.

Toy Model Predictions. Given the observed radial profile of star formation and the KS relation at z ∼ 1, Figure 11 indicates that the rate of change of metallicity due to star formation should be higher in the centers of galaxies than in their outskirts.

Furthermore, if metals were stationary (in a radial sense), we anticipate galaxies should develop detectable (i.e.,

log O H R 0.05 ( ) D D < -) metallicity gradients rapidly-on the order of 10-100 Myr in galaxy subpopulations below log M * /M e < 10.5. To reiterate, the fact that we observe a large fraction of flat and positive gradients indicates that declining metallicity gradients must be destroyed in galaxies on timescales shorter than the short time it takes for them to establish.

In Figure 11, the higher-mass galaxy populations exhibit a shallower evolution than the lower-mass galaxy populations. This is a consequence of the broken power-law form of the KS relation. The lower-mass galaxies tend to draw more from the low Σ gas branch of the KS relation (i.e., the branch with the steeper n ∼ 1.4 slope), while the higher-mass galaxies tend to draw more from the high Σ gas branch (i.e., the branch with the shallower n ∼ 1.1 slope). The shallower the slope in the KS relation, the smaller the gradient that develops. At n = 1, the rate of change of the metallicity is the same everywhere-no gradient would develop.

Flattening and Inverting Metallicity Gradients

The majority of galaxies over the redshift (0.6 < z < 2.6) and mass (8.5 < log M * /M e < 10.5) ranges probed by the CLEAR + sample have flat or positive metallicity gradients (Figure 6). This is the key result of the paper, and is in qualitative agreement with previous observations of galaxies spanning similar masses and redshifts (Stott et al. 2014;Wuyts et al. 2016;Wang et al. 2017Wang et al. , 2019;;Curti et al. 2019;Wang et al. 2020;Gillman et al. 2021). At z = 0.6-2.7, Wuyts et al. (2016) report for the KMOS 3D survey (using the near-infrared VLT/KMOS IFU) that 90% of the galaxies in their sample (log M * /M e > 10) are within 3σ of a flat gradient. The galaxies in the KMOS 3D sample rest at the highmass end of the CLEAR+ sample. If we restrict CLEAR+ to a similar mass range (log M * /M e = 10-11), we find that 89% of the ) in star-forming galaxies at z ∼ 1. The toy model is governed by the observed radial profiles of the star formation rate density in z ∼ 1 galaxies (Nelson et al. 2016a) and a broken power-law Kennicutt -Schmidt relation (KS; Kennicutt 1998). The toy model starts the simulation at t = 0 with a flat gradient and a metallicity matching the z ∼ 2.3 massmetallicity relation (Sanders et al. 2021). The shaded regions show the results from this toy model. Without additional metal redistribution, we should expect galaxies to rapidly develop declining metallicity gradients through star formation and stellar evolution-and to do so faster in lower-mass galaxies. The ubiquity of flat gradients in the observations implies that gradients are being disrupted on short timescales (∼10-100 Myr). The shading of each line represents the scatter of 100 realizations with varying KS slopes and metal yields (see the text for details).

sample are within 3σ of a flat gradient-in strong quantitative agreement with Wuyts et al. (2016). Using a sample of gravitationally lensed galaxies spanning similar masses and redshifts to the full CLEAR+ sample, Curti et al. (2019;KLEVER survey; using VLT/KMOS) and Wang et al. (2020; GLASS survey; using the near-infrared HST/WFC3 grism spectrographs) report that 89% (within 3σ) and 67% (within 2σ) are consistent with a flat gradient, respectively. As noted in Section 4.1, these are generally consistent with the CLEAR+ population measurements.

On average, we find that galaxies have more positive metallicity gradients at lower stellar masses-albeit with significant population scatter. Moreover, the medians of the lowest mass galaxy populations (log M * /M e < 10) are generally consistent with a positive slope. This is qualitatively consistent with several studies reporting a negative correlation between metallicity gradient and stellar mass at 0 z 1 (e.g., Carton et al. 2018) and 1 z 3 (Wuyts et al. 2016;Curti et al. 2020;Wang et al. 2020;Gillman et al. 2021). Wang et al. (2019) was the first to discover strong positive (also called "inverted") gradients in lowmass (log M * /M e ∼ 9) galaxies at this redshift. Positive gradients are physically incompatible with the closed-shell models constructed in Section 5.1, and may suggest that galactic-scale gas flows are re-distributing metal-rich gas from the centers of galaxies to their outskirts. This scenario is more explicitly explored in Wang et al. (2019). Wang et al. (2019) use maps of star formation survey density, stellar age, and gas fraction derived for the two galaxies in their sample to construct maps of inferred mass outflow rates. They find that the implied outflow rates are highest in the centers of their galaxies. They argue that the net impact of these flows is the transportation of metal-rich gas radially outward in galaxies. As such, galaxy winds will act to increase (i.e., lead to a flat or positive gradient) the gas-phase metallicity gradient of these galaxies. The fact that the median gradient of low-mass galaxies at this redshift is positive (Figure 6) might indicate that such flows are prevalent in this galaxy population.

In Figure 7, we compare the medians of the galaxy populations probed by our sample at 0.6 < z < 2.6 with the galaxy populations today measured from the MaNGA survey (Belfiore et al. 2017), as a function of mass. We infer a strong massdependent evolution in the metallicity gradient, with the steepest evolution needed for the higher-mass galaxy population.

This has an enticing parallel with the observed kinematic evolution of star-forming galaxies from z ∼ 2 to now. At z ∼ 2, the velocity dispersions of the ionized gas in galaxies are several factors higher than they are today (Kassin et al. 2007;Wisnioski et al. 2015;Simons et al. 2016;Johnson et al. 2018;Gillman et al. 2019;Molina et al. 2019;Übler et al. 2019;Wisnioski et al. 2019;Förster Schreiber & Wuyts 2020;Price et al. 2020). High gas velocity dispersions are likely accompanied by efficient radial mixing-acting to redistribute metal-enriched gas from the centers of galaxies to their outskirts. With time, star-forming galaxies gradually increase in rotational support and decline in dispersion support (Kassin et al. 2012;Simons et al. 2017). This evolution is a strong function of mass-with less massive galaxy populations on average having higher contributions from dispersion support at all times.

These kinematic parallels extend down to z ∼ 0. In the local universe, the majority of massive star-forming galaxies ( *  M M log 9.5 >

) have disk-like gaseous kinematics that are dominated by rotation (e.g., Barat et al. 2019), whereas only a fraction of low-mass star-forming galaxies ( *  M M log 9.5 < ) have formed a disk (Simons et al. 2015). The lingering dispersion support in the ionized gas in today's low-mass galaxies may be indicative of large-scale turbulent motions. Such turbulent motions could act to redistribute metals and sustain the observed average flatness of the gas-phase metallicity profiles at these masses (Belfiore et al. 2017).

Using the FIRE galaxy formation simulations, Ma et al. (2017) conclude that strong (declining) metallicity gradients only appear in galaxies with a well-formed disk (i.e., a galaxy with a rotation velocity higher than its local velocity dispersion)-but not all well-formed disks have a strong metallicity gradient. Furthermore, they find that highly perturbed non-rotating (mostly post-merger) galaxies tend to have flat gradients.

In Section 4.2, we report evidence for an increase in the intrinsic scatter of the metallicity gradients toward lower masses. This is consistent with recent results at this redshift (Wang et al. 2020), and with observations of the local universe (Ho et al. 2015;Bresolin 2019). The scatter in metallicity gradients provides a key benchmark for galaxy evolution models.

Specifically, the mean and scatter of the metallicity gradients have been shown to be a sensitive probe of stellar feedback. Using closed-box chemical evolution models, Ho et al. (2015) demonstrate that the population mean and scatter of metallicity gradients is sensitive to mass accretion rates and mass-loading factors-with high rates of both producing a narrow distribution centered around flat gradients (see also Belfiore et al. 2019). Gibson et al. (2013) compared the metallicity gradients of simulated galaxies from two simulation suites-one run with a "conservative" feedback model and one with an "enhanced" feedback model that can more efficiently drive hot SNe winds. Galaxies forming in the simulation with enhanced feedback always have flat metallicity gradients, while galaxies forming in the simulation with conservative feedback always rapidly build up metallicity gradients at high redshift and flatten at late times as the galaxy grows. By contrast, in the FIRE simulations, Ma et al. (2017) recover a wide scatter in metallicity gradients all with the same feedback model. They attribute this, in part, to the burstiness of star formation in the FIRE galaxies, and the ability of the subkiloparsec feedback model to switch between favorable and unfavorable conditions for driving metal-enriched outflows. Hemler et al. (2021) study the gas-phase metallicity gradients in the Illustris-TNG simulations, and find that the Illustris-TNG galaxies generally have steeper gradients than those found in the FIRE simulations. These differences are attributed, at least in part, to the differences in the feedback models-the Illustris-TNG model leads to less bursty and disruptive behavior than the FIRE feedback model.

It is clear that statistical measures of metallicity gradients at high redshift serve as an important benchmark for galaxy formation feedback models. Gradient demographics provide a unique way to confront the physical models underpinning these numerical simulations.

Another important actor for setting and disrupting metallicity gradients is gas accretion-either through metal-poor accretion from intergalactic filaments, (re-)accreted material from the circumgalactic medium, or galaxy mergers. Each channel acts in a unique manner. The impact of metal-poor accretion from cold filaments strongly depends on where that material is deposited. If deposited directly into the centers of galaxies (or radially mixed on quick timescales), it should act to dilute the central metallicity and flatten metallicity gradients. If the metal-poor gas is deposited on the outskirts of galaxies, it should push the gradient negative. For (re-)accreted material, a general flattening of metallicity gradients will occur as metals ejected into the CGM rain back onto the outskirts of galaxies (Fu et al. 2013;Grand et al. 2019). For mergers, numerical simulations indicate that they can rapidly redistribute metals around galaxies and flatten metallicity gradients (Rupke et al. 2010;Perez et al. 2011;Torrey et al. 2012). In general, accretion and merger rates are expected to be factors of 10-30 higher at z ∼ 2 than they are today (Rodriguez-Gomez et al. 2015;Scoville et al. 2017).

As the prevalence of all of these "flattening mechanisms" (star formation, accretion, mergers) declines with time, it is reasonable to assume that there will be a corresponding decline in the rate in which metals are (re-)distributed around galaxies and the ISM is diluted by metal-poor accretion. At late times, in a more hospitable universe, galaxies will have the chance to develop long-lived declining metallicity gradients through star formation.

Conclusions

We report on the gas-phase metallicity gradients in 238 galaxies over 0.6 < z < 2.6. The observations include deep near-infrared HST/WFC3 G102 grism spectra taken through the CANDELS Lyα Emission at Reionization (CLEAR) survey, as well as publicly available WFC3 G102+G141 grism spectra overlapping the CLEAR footprint. The combined G102+G141 spectral coverage allows for simultaneous coverage of multiple strong line metallicity indicators (notably [O II], [O III], and Hβ); and the high spatial resolution of the grism enables for resolved maps at the resolution of HST.

We summarize our conclusions as follows:

1. The majority of the galaxies (96%) in our sample have flat or positive radial gas-phase metallicity gradients. This result is generally consistent with previous results in the literature Curti et al. 2020;Wang et al. 2020), and is in marked contrast with the z ∼ 0 galaxy population (Belfiore et al. 2017). We produce an empirical toy model to demonstrate that flat and positive metallicity gradients require gas-phase metals to be (re-)distributed, or central metallicities to be diluted from metal-poor gas accretion, on rapid ∼10-100 Myr timescales. 2. We detect an intrinsic population scatter-i.e., scatter in excess of that expected from the observational uncertainties-across our full mass range (8.5 < log M * /M e < 10.5). We find strong evidence for a mass dependence to the scatter (such that there is excess scatter at lower masses) in the CLEAR sample (3.9σ) and even stronger evidence when we combine with existing measurements in the literature (5.2σ).

3. We explore correlations between the gas-phase metallicity gradient and stellar population properties at fixed stellar mass, including: star formation, sizes, star formation surface density, and star formation per potential energy. We find no statistically significant correlations between the metallicity gradient and any of these properties. We place constraints on the parameter space of intrinsic correlations ruled out by the CLEAR data. These constraints can be readily confronted with predictions from theoretical models.

These results strongly favor a scenario in which metals are re-distributed around galaxies on faster timescales than the short time it should take for star formation and stellar evolution to lead to a declining metallicity gradient. Despite rapid progress on both the observational and theoretical front, a full description of the processes responsible for this (re-)distribution-specifically, their timescales and relevance as a function of galaxy mass and redshift-remains outstanding.