NGDEEP Epoch 1: Spatially Resolved Hα Observations of Disk and Bulge Growth in Star-forming Galaxies at z ∼ 0.6–2.2 from JWST NIRISS Slitless Spectroscopy

We study the Hα equivalent width (EW(Hα)) maps of 19 galaxies at 0.6 < z < 2.2 in the Hubble Ultra Deep Field using NIRISS slitless spectroscopy as part of the Next Generation Deep Extragalactic Exploratory Public Survey. Our galaxies mostly lie on the star formation main sequence with stellar masses between 109 and 1011 M ⊙, characterized as “typical” star-forming galaxies at these redshifts. Leveraging deep Hubble Space Telescope and JWST images, spanning 0.4–4.8 μm, we perform spatially resolved fitting of the spectral energy distributions for these galaxies and construct specific star formation rate (sSFR) and stellar-mass-weighted age maps with a spatial resolution of ∼1 kpc. The pixel-to-pixel EW(Hα) increases with increasing sSFR and with decreasing age. The average trends are slightly different from the relations derived from integrated fluxes of galaxies from the literature, suggesting complex evolutionary trends within galaxies. We quantify the radial profiles of EW(Hα), sSFR, and age. The majority (84%) of galaxies show positive EW(Hα) gradients, in line with the inside-out quenching scenario. A few galaxies (16%) show inverse (and flat) EW(Hα) gradients, possibly due to merging or starbursts. We compare the distributions of EW(Hα) and sSFR to star formation history (SFH) models as a function of galactocentric radius. We argue that the central regions of galaxies have experienced at least one rapid star formation episode, which leads to the formation of the bulge, while their outer regions (e.g., disks) grow via more smoothly varying SFHs. These results demonstrate the ability to study resolved star formation in distant galaxies with JWST NIRISS.


INTRODUCTION
Investigating spatially resolved properties of galaxies offers important constraints on the mechanisms by which galaxies form stars and grow in stellar mass (e.g., Brinchmann et al. 2004;Papovich et al. 2005).Studies on the structure of galaxies up to z ∼ 2.5 have converged to a coherent picture that star-forming galaxies (SFGs) are on average larger and disk-dominated systems, while quiescent galaxies are more compact and bulge-dominated (e.g., Shen et al. 2003;Wuyts et al. 2011;Weinzirl et al. 2011;van der Wel et al. 2014;van Dokkum et al. 2015).These studies seem to imply that as SFGs evolve and quench their star formation, they may experience bulge assembly and/or compaction (Dekel & Burkert 2014;Papovich et al. 2015;Zolotov et al. 2015;Tacchella et al. 2016b), which results in an inside-out growth/quenching pattern that the central regions are formed at earlier times and therefore show older age and lower specific star-formation rates (sSFR ≡ SFR/M * ) than the outer disc.This inside-out quenching has been supported by extensive observational evidence from spatially-resolved studies at low redshift that the sSFR profiles appear to be centrally suppressed in massive galaxies log(M * /M ⊙ ) > 10.6 (e.g.Spindler et al. 2018), in green valley galaxies (e.g.Belfiore et al. 2018), and in galaxies below the star formation main sequence (SFMS; e.g., Ellison et al. 2018).
At higher redshift (z ≳ 1), studies have used Hα to measure spatially-resolved star formation in galaxies using the integral field units (IFU; Tacchella et al. 2015Tacchella et al. , 2018;;Wilman et al. 2020) and the HST Grisms (Nelson et al. 2012(Nelson et al. , 2016a;;Matharu et al. 2022).These studies have found evidence that ongoing star formation occurs in disks that are more extended than those occupied by existing stars in SFGs at z ∼ 0.5 − 2.7, and the interpretation is that the disks form "inside-out".Tacchella et al. (2018) focused on ten massive galaxies at z ∼ 2 and found flat sSFR radial profiles for galaxies with M * = 10 10.5−11 M ⊙ , and found suppressed central sSFR radial profiles for galaxies with M * > 10 11 M ⊙ .They interpreted this as galaxies being quenched from the inside out.This scenario is also supported by the study of the gas reservoirs in a massive star-formation galaxy at z ∼ 2.2 that shows that both molecular gas fraction and molecular gas depletion time are lower in the central region compared to the outskirts (Spilker et al. 2019).However, these studies are limited to massive and brighter galaxies or are limited to lower redshifts because HST can only study Hα to z < 1.6, or rely on "stacking" to get a sufficient signal-to-noise ratio.
Recently, there has been an increasing number of studies focusing on the 2D and radial properties of high-redshift galaxies using JWST data (Pérez-González et al. 2023;Giménez-Arteaga et al. 2023;Wang et al. 2022;Abdurro'uf et al. 2023).In this paper, we build on these using measurements of the spatially resolved Hα derived from slitless grism spectroscopy from JWST (Gardner et al. 2006(Gardner et al. , 2023) NIRISS (René et al. 2023) for 19 galaxies at 0.6 < z < 2.2, using the first epoch of the deep JWST NIRISS slitless spectroscopy from the Next Generation Deep Extragalactic Exploratory Public (NGDEEP) Survey (Bagley et al. 2022).We focus on the equivalent width (EW) of Hα which provides an independent method for determining the sSFR, as the EW(Hα) is a ratio of Hα flux, an SFR indicator, and the continuum flux at λ rest = 6563 Å which traces the stellar mass.Additionally, EW(Hα) is less susceptible to the presence of dust compared to the rest-frame UV.As a comparison, we explore the sSFR and stellar mass-weighted age maps of these galaxies derived from a spatially resolved SED fitting method on similar spatial resolution using 9 bands from Hubble Space Telescope (HST) ACS and WFC3 imaging from CANDELS, and 11 bands from JWST imaging from NGDEEP, the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2023;Rieke & the JADES Collaboration 2023) and the JWST Extragalactic Medium-band Survey (JEMS, Williams et al. 2023).We compare these spatially resolved properties and their radial profiles to reveal the spatial distribution of stellar mass, star formation, and age of these galaxies.Thanks to the high spatial resolution, deep imaging, and slitless spectroscopy, we can perform these spatially-resolved analyses on individual galaxies, which provide details on the star formation, such as those contributed from off-center clumps that are mostly erased when stacking.

The NGDEEP survey
We use the JWST NIRISS Wide Field Slitless Spectroscopy (WFSS) observations of galaxies as part of NGDEEP (Proposal ID #2079; Bagley et al. 2023).NGDEEP targets the Hubble Ultra Deep Field (HUDF) with NIRISS slitless spectroscopy.The science aims of NGDEEP include measuring metallicities and SFRs for lowmass galaxies in the redshift range 0.5 < z < 4. In parallel, NGDEEP targets the HUDF-Par2 parallel field with NIRCam with the science aim of discovering galaxies up to z > 12 and constraining the slope of the faint-end of the rest-UV luminosity function (Leung et al. 2023).
The first epoch of NGDEEP observations was taken over 2023 Jan 31 -Feb 2. These include observations with NIRISS with the GR150R and GR150C grisms and the F115W, F150W, and F200W filters.The total exposure times on grisms are 95 ks in F115W, 43 ks in F150W, and 32 ks in F200W.In addition, the observations include NIRISS direct imaging in the same filters with total integration times of 5.4, 1.7, and 1.7 ks, respectively.
We reduced the direct images using v1.9.4 of the JWST Pipeline 1 (Bushouse et al. 2023).From the direct images with NIRISS, we measured the 3σ depth of F115W, F150W, and F200W imaging with 0. ′′ 2 radius aperture to be 29.1, 28.6, and 28.6 AB mag, respectively.For the grism spectroscopic imaging, we reduced the NIRISS data using two methods, including the Grism Redshift & Line software (GRIZLI; Brammer et al. 2022) and EMission-line two-Dimensional (EW2D, Pirzkal et al. 2018).The latter method is described in Pirzkal et al. (2023).For the analysis of emission-line maps here, we used the images from GRIZLI (see details in Section 2.3 below).From the one-dimensional extracted spectra, we found that the grism data are sensitive to line fluxes to f lim ≈ 2 × 10 −18 erg s −1 cm −2 (3σ, Pirzkal et al. 2023).

Available Optical/NIR Imaging and Photometry
We utilized a vast array of deep imaging taken with HST and JWST available in the HUDF field.This includes HST ACS F435W, F606W, F775W, F814W, F850LP, and WFC3/IR F105W, F125W, F140W and F160W, and we use the latest reductions processed as part of the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS, Koekemoer et al. 2011;Grogin et al. 2011), see Finkelstein et al. (2022) for more details.Besides the three JWST NIRISS direct-imaging bands from NGDEEP, we include the nine JWST NIRCam bands: F090W, F335M, F356W, and F444W from the DR1 release from JADES (Eisenstein et al. 2023;Rieke & the JADES Collaboration 2023) and F182M, F210M, F430M, F460M and F480M from the DR1 release from JEMS (Williams et al. 2023).
We adopt the HST photometry catalog from Finkelstein et al. (2022), and use the same method to extract the photometry in JWST imaging.To measure consistent photometry across all bands, we align each image astrometrically and on the same pixel scale (0. ′′ 06).We then match each image to the PSF of the F160W image.The PSFs of HST images are constructed by stacking selected stars (see Finkelstein et al. 2022 for details).However, there are no isolated, unsaturated stars available in the NGDEEP NIRISS field with which to construct an empirical PSF.Thus, we constructed a PSF for each JWST bandpass using WebbPSF (Perrin et al. 2012(Perrin et al. , 2014)).For the PSF in each bandpass, we generated 100 different PSFs at different position angles, which we then averaged.
We constructed kernels to match each PSF to that in the HST WFC3/IR F160W band using the pypher 2 Python routine.Each image was then convolved with its respective kernel.To examine the accuracy of this PSF matching process, we measured the growth curve of four, compact elliptical galaxies in the PSF-matched images.The average ratio of the growth curve of these elliptical galaxies in the PSF-matched JWST images and the HST F160W image spans a range of 0.99 to 1.10 at a radius of 1 ′′ .This was one of the motivating factors to include a 10% additional error (in quadrature) to the uncertainties of the fluxes when we performed the integrated and spatial resolved SED fitting.
We use the F160W image as the detection image and cycle through every aligned and PSF-matched JWST image as the measurement.We adopt the small Kron elliptical apertures and correct to total using an aperture correction.We estimate photometry uncertainties empirically to account for the Poisson correlation following Finkelstein et al. (2022).

Grism Data Reduction and Spectral Extractions
For the analysis of emission-line maps here, we processed the data using the latest version of the GRIZLI.GRIZLI performs full end-to-end processing of NIRISS imaging and slitless spectroscopic data sets, including retrieving, preprocessing the raw observations for cosmic rays, flat-fielding, sky subtraction, astrometric corrections, alignment, modeling contamination from overlapping spectra, extracting 1D and 2D spectra, and fitting full continuum+emission-line models.
From the data products, GRIZLI derives the emissionline maps by drizzling the contamination-and continuumsubtracted 2D spectral beams back to the imaging plane.It uses the JWST/NIRISS imaging to scale the spectroscopic data to match the total fluxes measured in the direct images.We created emission-line maps with a pixel scale of 0. ′′ 1, similar to the FWHM of PSF in our NIRISS filters.The uncertainties on the line maps are calculated using the drizzle weights from the constituent beam pixels.For additional details on GRIZLI and its data products, we refer the reader to Estrada-Carpenter et al. ( 2019

Sample Selection
Here, our goal is to study the emission-line maps for galaxies with high signal-to-noise ratios (SNR).To construct this sample, we selected galaxies that have Hα line maps where the total Hα SNR is > 20σ, and the galaxies have data in at least 50 pixels each with Hα SNR > 3σ.The latter criterion is to ensure that we have enough pixels to derive the EW(Hα) radial gradient.For example, a circular area with ≥50 pixels would have a radius of ≥4 pixels, providing at least two annular bins each with a width of at least two pixels to facilitate a measurement of the gradient.We visually examined these Hα line maps and excluded those with any contamination across the source possible due to the failure of the continuum subtraction.
Our final sample includes a total of 19 galaxies.Fig. 1 displays the false color images and Hα line maps for this sample.We note that the current NIRISS grism data has 2 position angles (GR150R and GR150C grisms), providing good spatial coverage of the emission-line features.The full NGDEEP NIRISS data will double the data and have four position angles, which should largely reduce the fraction of contamination and increase the number of galaxies with secure line maps.
The GRIZLI-derived redshifts are consistent with previous ground-based and space-based spectroscopic measurements of our sources, where we find a median offset of only ∆z/(1+ z) = 0.002.Therefore we are confident in the identification of the Hα emission line for these galaxies.
We identify four galaxies in our sample as AGN by crossmatching to the latest AGN catalog from Lyu et al. (2022).Two of them (NGDEEP_01837 and NGDEEP_02071) are classified as AGN based on their near-to mid-IR excess emission from SED fitting and the other two (NGDEEP_01520 and NGDEEP_00923) are classified as AGN based on their X-ray to radio luminosity ratio higher than expected for stellar processes in a galaxy.For these galaxies, to remove any AGN contamination, we mask their central 3 × 3 pixels (0. ′′ 3 × 0. ′′ .3,larger than the PSF FWHM of F160W image), where the center is derived from the peak of the F160W image.
Figure 2 shows the stellar mass versus SFR of our sample using quantities obtained from SED fitting to the integrated emission (see Section 3.1).We compared our galaxies to the "star formation main sequence" (SFMS) derived from larger samples from Whitaker et al. (2014) in the redshift range of our sample.Our sample generally follows the SFMS with a median of log(SFR/SFR MS ) = 0.05 with a scatter from -0.05 to 0.31 dex from the 16th/84th percentile.Therefore, our sample represents galaxies along the main sequence of star formation at these redshifts.Although we preferentially select galaxies with extended Hα, we do not expect our sample to be biased toward galaxies with more star formation, because no significant correlation is found between SFR and effective radius at fixed stellar mass for galaxies at 0.5 < z < 2.5 and stellar mass within 9.0 ≤ log(M * /M ⊙ ) < 11.0 (Lin et al. 2020).

SED fitting
We employed the SED fitting Code Investigating GALaxy Emission (CIGALE) (Boquien et al. 2019;Yang et al. 2020) to derive constraints on the integrated properties of the stellar populations of our galaxies using the 21 bands photometry catalog spanning the wavelength range of 0.43 -4.8 µm (see details in Section 2.2) and the spatially-resolved properties using the 21-band fluxes in each Voronoi bin (see Section 3.2).The parameters for these fits are listed in Table 1.In detail, we adopted a delayed exponential star formation history (SFH) allowing the τ and stellar age to vary from 0.05-20 Gyr and 0.02-10 Gyr, respectively.We assumed a Chabrier (2003) IMF and the stellar population synthesis models presented by Bruzual & Charlot (2003) with solar (Z ⊙ ) metallicity.The dust attenuation follows the extinction law of Calzetti et al. (2000) allowing the dust attenuation in emission lines from nebular regions E(B −V ) l to vary from 0 to 1.1, and a lower dust attenuation in stellar continuum with a fixed dust attenuation ratio (E(B − V ) factor = 0.44) between emission lines and stellar continuum.This corresponds to allowing dust attenuation in the stellar continuum to vary from 0 to 0.5.The amplitude of the absorption UV bump feature produced by dust at 2175 Å is allowed to range from 0 to 3, and the slope of the dust power law is allowed to vary −0.5 to 0.
For the integrated properties, e.g., the mass, SFR, and age, we also test the robustness of these properties by including mid-and far-IR data from IRAC/Spitzer, MIPS/Spitzer from the GOODS-Herschel catalog (Elbaz et al. 2011; see discussion in Section 5.2).Therefore, to fit the dust emission, we also adopt dust emission templates from Dale et al. (2014), which model the star-forming component as dM d (U) ∝ U −α dU, where M d is the dust mass, U is the ra- False color images and Hα maps for the galaxies in our galaxy sample.For each galaxy, we show the RGB images using the HST/ACS F814W, JWST/NIRISS F150W and JWST/NIRCam F444W (left), and the Hα emission line map constructed from the NIRISS slitless spectroscopy (right).The sample is ordered by increasing redshift.Each image is 2 ′′ × 2 ′′ , corresponding to 13.4×13.4kpc 2 and 16.6×16.6kpc 2 at z = 0.6 and z = 2.2, respectively.The text inset lists the galaxy ID, redshift, and whether galaxies are identified with AGN.
Figure 2. The stellar mass versus the SED-derived SFR of our final sample.Galaxies are colored by their redshift.The two galaxies with X-ray-detected AGN are marked by red diamonds, and the two with AGN implied by their near-/mid-IR excess are marked with orange squares.For comparison, the star-formation-main-sequence (SFMS) for galaxies at 1.0 < z < 1.5 and 1.5 < z < 2.0 from Whitaker et al. (2014) are shown in blue and orange, respectively, as labeled.
diation field intensity.We allow the α to vary from 0.25 -4.
We adopt quantities and associated errors for each parameter from the posteriors from the CIGALE output catalog.In particular, we use the stellar mass (M * ), SFRs averaged over 100 Myr (SFR 100Myr ), and the mass-weighted-stellar age, and their associated errors.The sSFR is calculated as sSFR ≡ SFR 100Myr /M * .

Spatially-resolved photometry and SED fitting
To measure the spatially-resolved photometry, we use the PSF-matched imaging and reproject them to a pixel scale of 0. ′′ 1, to match the Hα emission line map.To place stronger constraints on the resolved properties, we applied an adaptive binning algorithm 3 based on the Weighted Voronoi Tessellation (MVT) method described in (Cappellari & Copin 2003;Diehl & Statler 2006).We chose 2 pixels to be the minimum size of the Voronoi bin, which is approximately the size of the FWHM of the PSF of F160W (FWHM = 0. ′′ 18).We run CIGALE using the 21-band fluxes in each Voronoi bin, and the same modules and parameters as described in Section 3.1.

Hα Equivalent Width Maps
We use the Hα maps from GRIZLI reconstructed from the NIRISS WFSS data for each galaxy.From these maps, we calculated the rest-frame EW(Hα) maps as follows, where we take the continuum flux f continuum to be, In these equations, f Hα is the flux from the Hα emission-line map, f direct−image is flux from the NIRISS direct imaging in the bandpass that contains the Hα line, ∆λ is the FWHM of the direct-image bandpass, and (1 + z) is a factor to correct to the rest-frame.The uncertainty on EW(Hα) is obtained following the propagation of uncertainty using the associated error maps.To obtain secure EW(Hα) measurements, we include only those pixels where the Hα flux has significance > 3σ, and the underlying continuum flux from the direct image has significance > 5σ.
Due to the low spectral resolution of the NIRISS grism, the Hα and [N II] lines are blended.Therefore, in this work, EW(Hα) includes the sum of the Hα and [N II] lines (For simplicity, we use EW(Hα) in the text, but use "EW(Hα+[N II])" in the figure for clarity).For the stellar mass and redshift of the galaxies in our sample, we expect [N II]/(Hα+[N II]) to be ∼ 0.2 (Faisst et al. 2018).Because we are looking at radial gradients in the EW(Hα), this has no impact on our results so long as there are no strong changes in Hα+[N II] throughout the galaxies (see further discussion in Section 5.2.2).Another effect on the EW(Hα) measurements is impacted by the presence of AGN, which could have higher Hα+[N II] values (Baldwin et al. 1981).This was the motivation to remove the inner 3 × 3 pixels of galaxies that show AGN activity (see above).
We also compare the EW(Hα) from the spatially resolved measurements against the EW(Hα) measured from the bestfitted 1D integrated spectrum.The 1D spectrum is extracted from the stacked 2D spectra using an optimally-weighted method following Horne (1986).The integrated EW(Hα) is calculated using where F Hα is the best-fitted Hα emission line profile, and F continuum is obtained by fitting a linear function to the bestfitted spectra in two wavelength ranges 6510Å < λ rest < 6550Å and 6620Å < λ rest < 6660Å.These two wavelength windows are chosen to be close to the Hα+[N II] regions but are not affected by the Balmer absorption.The uncertainty is obtained from a Monte-Carlo where we perturb each bin using the full covariance matrix of the best-fitted spectra.
The average EW(Hα) from the emission line flux maps is calculated using where f i,Hα are the Hα fluxes in each pixel and f i,continuum are the continuum fluxes in each pixel from the direct imaging.The sum is over all i pixels that have EW(Hα) measurements.The uncertainty on ⟨EW (Hα)⟩ is obtained using the associated error maps.We note that the uncertainties on the integrated EW(Hα) and for ⟨EW (Hα)⟩ from the maps do not include uncertainties for the background subtraction and contamination uncertainties.The former should be negligible assuming the relative uncertainties are small on the size of our galaxies.Our inspection of the images suggests the latter is also negligible, but we can test this with the second epoch of NGDEEP NIRISS data.
The results are shown in the right panel of Figure 3.We see the majority of data follows the 1-to-1 line, with a few outliers.The median difference is 0.02 dex with the 16th/84th percentiles of the distribution of differences are −0.05/0.04dex.
We consider here some possible reasons that could cause the differences between the integrated EW(Hα) and that calculated from the resolved maps.The integrated EW(Hα) is calculated from the 1D spectra, and the extraction from the stacked 2D spectra uses the optimally weighted method.We calculate the ⟨EW (Hα)⟩ from the resolved maps using an unweighted average over the pixels that have EW(Hα) detections.Secondly, the 1D spectra are stacked from spectra taken in two diagonal PAs, (i.e., GR150R and GR150C), there can be issues where resolved structures are poorly modeled by Grizli, which models all features as a single-Gaussian model.This could lead to an underestimate of the fluxes in some galaxies.Thirdly, for the ⟨EW (Hα)⟩ from the resolved maps, it is also possible that for more dusty and redder galaxies where the continuum is rising toward redder wavelengths, our assumption of a flat continuum assumption could potentially bias the EW(Hα) and depend on the observed-frame wavelength of Hα relative to the bandpass.In addition, there are differences in how the stellar absorption is accounted for in the EW measurements.GRIZLI includes an estimate of the stellar continuum when it generates the emission-line flux images Here, we then use the direct image to model the underlying continuum (Equation 2).Because the direct imaging includes the effects of the stellar continuum, we are doublesubtracting the continuum from the emission line fluxes.This could underestimate the underlying continuum, and overestimate the resolved the ⟨EW (Hα)⟩.We tested the effects of stellar absorption in the continuum when measuring the EW by comparing the EW(Hα) measured using the continuum correction from GRIZLI compared to using a linear fit to the continuum (see above).We measure the differences in the range of 0.06 -0.08 dex in the integrated EW(Hα) between using a linear function as an underlying continuum and using a continuum with absorption.This is an upper limit because the continuum from the direct image is averaged over a longer wavelength, instead of averaging over the Hα region as in the integrated EW(Hα).Therefore, we take this as an upper limit on the systematic uncertainty in the EW arising from how we model the continuum.In addition, the effect of other emission lines falling in the same filter is negligible.After removing other emission lines, the ⟨EW (Hα)⟩ is lower by a range of 0.01 -0.03 dex.We note that these differences could affect our results when we compare them to integrated studies (Section 4.1), but this does not change our results on the EW(Hα) gradients.The analysis of the integrated emission of galaxies has shown strong correlations between the equivalent widths of various nebular emission lines, the specific SFR (sSFR), and the age of the stellar population from z ∼ 0 to 6 (Schaerer et al. 2013;Reddy et al. 2018).Here, we use the spatially resolved information in the JWST and HST data to study these correlations within individual galaxies.
Figure 4 shows the 2D histograms of all pixels from the spatially resolved maps that include sSFR and age versus EW(Hα).We also show the average (median) values in bins along the abscissa.The median and associated median error are calculated with a bootstrap method following Shen et al. (2023).Each median and its error are calculated using approximately a similar number of data points (∼200 individual pixels) per median.
Focusing on the spatially resolved results, the EW(Hα) increases with increasing sSFR and decreases with increasing age.We adopt the non-parametric Spearman's ρ test to assess these correlations4 .The returned p-values are much smaller than 0.01, therefore, we reject the null hypothesis that the EW(Hα) and SSFR and age are uncorrelated.Furthermore, we note that the EW(Hα) and sSFR (and EW(Hα) and age) relations are independent measurements because the EW(Hα) is obtained from the NIRISS grism data independent of the SED fitting results.In contrast, the sSFR and age are derived from SED fitting to the broad-band and mediumband photometry.
We compare our trends derived from the spatially resolved maps to those derived from integrated measurements.For the integrated properties of galaxies, we use the sSFR (and age) versus EW(Hα+[N II]) relations for star-forming galaxies at z = 1.4 − 3.8 from the MOSFIRE Deep Evolution Field (MOSDEF) survey (Reddy et al. 2018).Our median sSFR versus EW(Hα) and mass-weighted age versus EW(Hα) trends are offset from those of MOSDEF by a median of 0.11 dex and 0.18 dex, respectively.This is comparable to the 1σ scatter reported in MOSDEF by Reddy et al. (2018).This offset could be due to the different star formation histories adopted in the SED fitting that MOSDEF assumed a constant star formation history.
Furthermore, we see a clear curvature in our resolved median sSFR-EW(Hα) and age-EW(Hα) values.Specifically, the slope of the relations is "flatter" for high sSFR and young ages.This is also apparent in the integrated relations.In addition, the sSFR versus EW(Hα) and age versus EW(Hα) relations show a large scatter.Quantitatively, the scatter is 0.30 dex and 0.29 dex in the sSFR versus EW(Hα) relation and age versus EW(Hα) relation, respectively.This scatter is larger than that reported from the MOSDEF derived from the integrated properties of galaxies, where Reddy et al. (2018) found 0.12 dex and 0.23 dex, respectively.We do not consider these to be significant given the small number of objects in our sample, and the different methods used by MOSDEF and here.Instead, we suspect that the scatter and the shape of our resolved sSFR and age versus EW(Hα) reflect complex evolutionary trends inside galaxies and galaxy-to-galaxy variation, while the trends of integrated measures only reflect the latter.We will discuss this further in Section 5.
Regardless of the systematic offsets and scatter, these comparisons indicate a consistent relation between the sSFR and age versus EW(Hα) relations for both the integrated emission from galaxies and the spatially-resolved properties.This means that these trends are driven by star formation and stellar populations within galaxies on spatial scales on the order of kiloparsecs.

Radial Gradients in EW(Hα), sSFR and Age
The radial profiles of EW(Hα), sSFR, and age can test the stellar build-up and quenching processes.For example, different radial trends are expected if these processes proceed in an "inside-out" or "outside-in" fashion (e.g., Ellison et al. 2018).To quantify the radial gradient, we fit for each galaxy a linear relation to the logarithm of EW(Hα), sSFR, and age as a function of the proper distance to the luminosity-weighted center (the "galactocentric radius").We use a Gaussian Mixture Model method (LINMIX, Kelly 2007) to fit the relation, which allows for uncertainties on the dependent variables (EW(Hα), sSFR, and age).The best-fit parameters are the median of the fitted parameters (the slope and intercept) from 400 random draws from the posterior, and the associated errors are the 16th and 84th percentiles of each parameter.We reiterate that for galaxies with AGNs, their central pixels (0. ′′ 3 × 0. ′′ 3) are not included in the fitting.However, if these pixels are included, we do not find a significant change in the measured gradients.In Figure 5, we show two examples of the EW(Hα), sSFR, and age maps, and as a function of galactocentric radii with their best-fitted lines.The radial gradients of EW(Hα), sSFR, and age for our galaxies are listed in Table 2.In the Appendix we show the EW(Hα), sSFR, and age maps, and their radial profiles for the rest of the galaxies in our sample.
For the sSFR and age gradients, there are 10 and 16 galaxies, respectively, that have positive sSFR gradients (sSFR increases with galactocentric radius) and negative age gradients (age decreases with radius).The median sSFR and age gradients are 0.003 +0.049 −0.038 dex/kpc and −0.04 +0.04 −0.03 dex/kpc.There is a correlation between the strength of the EW(Hα) gradient and that of the sSFR gradient in the galaxies in our sample.The Spearman's ρ test applied to the data in Figure 7 returns a p-value of 0.003 which we interpret as evidence (> 2.5σ) that a correlation exists.While there is an  anti-correlation between the strength of the EW(Hα) gradient and that of the age gradient (right panel of Figure 7, the Spearman's ρ test (p-value=0.07)result is less conclusive.This could be a result of different factors, which we discuss below in Section 5, but at least some of the weakness in the correlations is because of one galaxy that "bucks the trend", which we discuss in the next paragraph.
There is one galaxy, NGDEEP_00922, that shows a large offset from the relationship in both EW(Hα) versus sSFR and EW(Hα) versus age gradients in Figure 7.This galaxy has a high EW(Hα) gradient of 0.22 dex kpc −1 but very negative sSFR gradient of −0.14 dex kpc −1 .The reason for this is that this galaxy's sSFR and age maps show a region of high sSFR and low age region near the center (the top row of Figure 8), while this feature does not appear with an excess in the EW(Hα) map.This galaxy may have either very recently quenched the star formation in the core (leading to a drop in EW(Hα)) that has yet to appear in the broad-band SED.It is also possible that the degeneracy between the age and dust attenuation in the broad-band SED fits have misrepresented the physical conditions in this bin (leading to an erro-neously high sSFR).This could be tested by using spatially resolved mid-IR imaging (e.g., JWST/MIRI), which will be available in the future.Regardless, if we remove this one object from our fits, the correlation between the EW(Hα) gradient and sSFR gradient becomes even stronger, where the Spearman ρ test returns a p-value of 10 −7 .Similarly, the anti-correlation between the EW(Hα) gradient and age gradient becomes stronger, where this test returns a p-value of 0.004.
Previous studies have argued that the shape of the sSFR radial profile correlates with a galaxy's distance from the center of the SFMS.This has been reported at low-redshift (Ellison et al. 2018), at z ∼ 1 (Nelson et al. 2021) and at z ∼ 2 (Tacchella et al. 2018).We therefore consider if our galaxies follow this trend.We color-coded galaxies based on the distance from the SFMS, using the ratio between the measured SFR and SFR on the SFMS relation at their stellar mass and redshift (log(SFR/SFR MS )).This is shown in Figure 7 where positive (negative) values indicate galaxies with SFRs above (below) the SFMS.Qualitatively, the gradients of EW(Hα) and sSFR tend to be more positive for galaxies below the Here the gradients are the slope measured for the radial gradients for each galaxy (see Figure 5).In each plot, the histograms show the distribution of EW(Hα), sSFR, and age gradients for the individual galaxies.The markers in the scatter plots are color-coded by the distance of the galaxies from the SFMS (positive/negative values are above/below the SFMS).Galaxies identified with AGN are marked by open red diamonds and orange squares.The best-fit linear relations are shown as the purple dash lines along with 400 random draws from the posterior (faint grey lines).main sequence, and more negative for galaxies above the main sequence.We do not observe such a trend in the age gradients.However, we cannot make any quantitative conclusions based on the small sample of 19 objects that span a wide range of stellar mass, SFR, and redshift.We plan to return to this question once larger samples are available.

DISCUSSION
Previous studies have seen that the integrated (total) specific sSFRs and gas fractions are correlated with galaxy structure (e.g., Papovich et al. 2015;Whitaker et al. 2015;Freundlich et al. 2019).This has been observed also in studies of spatially resolved galaxies (e.g, Wuyts et al. 2011Wuyts et al. , 2013)).The implication is that galaxy quenching and star-formation activity are linked.Here we see that this process is tied to the sSFR and EW(Hα) values within a galaxy.
The majority of our sample (16 out of 19 galaxies) show positive EW(Hα) gradients where "positive" means the EW increases with radius.The median EW(Hα) gradient is 0.09 dex kpc −1 .This implies there is active build-up of the stellar content in the outer regions of galaxies, compared to the inner-most regions.This result supports an "inside-out" scenario of galaxy disk formation, where the central component is built first (and quenches first) and the disk structure grows radially, forming stars at preferentially larger radii (e.g., Cole et al. 2000;van den Bosch 2002;Papovich et al. 2005;Aumer et al. 2013;Somerville & Davé 2015).Cen-Figure 8.The EW(Hα), sSFR, and age spatial maps for galaxies with atypical gradients.The top row shows NGDEEP_00922, which shows a large offset from the relationship in both EW(Hα) -sSFR and -age gradients (Section 4.2).The middle and bottom rows show the two galaxies with negative EW(Hα) gradients (NGDEEP_01831 and NGDEEP_01948; Section 5.3) tral suppressed star formation has been seen in the interpretation of EW(Hα) gradients in massive high-redshift galax-ies with M * > 10 10 M ⊙ at 0.7 < z < 1.5 via stacking (Nelson et al. 2016a) and in an M31-sized progenitor at z ∼ 1.25 with M * ∼ 10 11 M ⊙ (Nelson et al. 2019).These results are also consistent with studies of the Hα morphologies, which generally find that the effective radii of galaxies measured in Hα are larger than that of stellar continuum for star-forming galaxies at z ∼ 0.5 -3 (Wilman et al. 2020;Matharu et al. 2022).In addition, this is in line with studies where they found a reduced specific SFR in the centers of massive galaxies at z ∼ 1 (Nelson et al. 2016a), z ∼ 2 (Tacchella et al. 2018) and at z ∼ 4 (Jung et al. 2017;Ji et al. 2023).
In section 4.2, we found that of the 16 galaxies with positive EW(Hα) gradients, 10 also show positive sSFR gradients.The median EW(Hα) gradient is more positive (0.09 dex kpc −1 ) than both the median sSFR gradient (0.003 dex kpc −1 ) and the median age gradient (−0.04 dex kpc −1 , see Figure 7).It is interesting then, that it is not universally true that galaxies with positive EW(Hα) gradients also have positive gradients in sSFR nor negative gradients in age (see Section 4.1).
There are multiple reasons we may expect the gradient in EW(Hα) and the sSFR to behave differently.Firstly, this can be a result of the complex star-formation history within galaxies.Because the Hα emission responds to changes in the SFR on shorter physical timescales compared to the SFR derived from the SED fitting to the broadband data.We discuss the implication of this on the relation between the EW(Hα) and sSFR gradients in the context of understanding the spatially-resolved star formation in Section 5.1.While the choice of the parameterization of the star-formation history is important, as we argue below this is unlikely to change our overall interpretation (Section 5.1).
Secondly, other effects could contribute to the gradient in EW(Hα) and the sSFR.There could be spatial variations in the dust attenuation and/or metallicity that impact our measurements.We consider this in Section 5.2 and argue that these factors are less likely to be driving the differences between EW(Hα) and sSFR that we observe.
Finally, there are exceptions.Two interesting galaxies in our sample show negative EW(Hα) and sSFR gradients, and are outliers compared to the rest of our sample.Understanding how these connect to the rest of the sample is important, and we discuss them in Section 5.3.

Spatially Resolved Star Formation Histories in
Galaxies: Evidence of Bulge Formation Comparing the spatially resolved EW(Hα) and sSFR gradients in galaxies can provide insights into different timescales in the recent star-formation history (SFH).For example, the Hα-to-UV ratio has been used to quantitatively evaluate SFH burstiness (e.g., Guo et al. 2016;Faisst et al. 2019;Emami et al. 2019;Asada et al. 2023).Here we can use these results to study the spatially resolved SFHs in galaxies.
The Hα emission is sensitive to ionization from O-type stars with lifetimes ∼5-10 Myr.The SFR derived from SED fitting primarily responds to changes in the rest-UV continuum, which is sensitive to changes on ∼100 Myr timescales (the main-sequence lifetime of B-type stars).The massweighted ages from the SED fitting probe the ages on longer timescales and are therefore sensitive to stars with lifetimes of 100-1000 Myr.Thus, as shown in Figure 7, the relation between the EW(Hα) and sSFR gradients shows a tighter correlation compared to that of the EW(Hα) and age gradients, where the latter also shows larger scatter.
Considering that EW(Hα) and the SED-derived sSFR trace the star formation on different timescales, we expect these differences should manifest in the position of galaxies in the sSFR -EW(Hα) diagram.To study this, we used the BAGPIPES code (Carnall et al. 2018) to generate SEDs from stellar population models with simple, delayed, exponentially declining star-formation histories, where SFR ∼ t × exp(−t/τ ), where we allow both the age (t) and the e-folding timescale (τ ) to vary between 0.1-2 Gyr.We also assumed solar metallicity, a fixed nebular emission with log(U) = −2 and a fixed Calzetti et al. (2000) dust attenuation law with A(V ) = 0.4 mag.We then computed the EW(Hα) and sSFR from these models.The left panel of Figure 9 overplots these models on the 2D histogram of sSFR versus EW(Hα) for all pixels from all the galaxies in our sample.Inspection of the figure shows that these models with ages of 0.5 to 2 Gyr and τ values of 0.1 to 2 Gyr broadly span the parameter space seen in the distribution of EW(Hα) and sSFR within galaxies in our sample.Therefore, variations in the star-formation histories can account for the distribution of EW(Hα) and sSFR we observe.
We further explored the effects of the star-formation histories on the EW(Hα) and sSFR values as a function of galactocentric radii.We separated the galaxy pixels based on their normalized galactocentric radii.Here, we normalized the galactocentric radius of each galaxy by the galaxy's effective radius (r eff ) measured in the bandpass closest to rest-frame 1.2 µm. 5The effective radii are measured using GALFIT (Peng et al. 2002(Peng et al. , 2010) ) with a single Sérsic profile, applied to the image for each galaxy.We then divided the pixels from each galaxy into four bins of normalized radial distance, r/r eff < 0.5, 0.5 < r/r eff < 1, 1 < r/r eff < 1.5, and r/r eff > 1.5, and plot their distribution in sSFR versus EW(Hα).The right four panels of Figure 9 show these results.In these panels of Figure 9, we see that as the galactocentric distance increases, the 2D distributions of EW(Hα)-sSFR shift vertically upwards toward higher values of EW(Hα), albeit with large scatter.The median EW(Hα) increases by 0.24 dex from the innermost regions (r/r eff ≤ 1) to the outermost regions (r/r eff > 1).
Comparing to the BAGPIPES models, the innermost regions (r/r eff ≤ 0.5) prefer models with relatively young ages (∼500 Myr) but short e-folding times, (τ ∼100 Myr).In contrast, the outermost regions (1 < r/r eff ≤ 1.5 and r/r eff > 1.5) favor models with large τ (∼2 Gyr) and a wide range of age (∼0.5-2Gyr).Therefore, if the differences in EW(Hα) and sSFR trace differences in the star-formation history, then these results imply that the innermost regions of the galaxies have experienced a faster build-up of stellar mass (with a star-formation history represented by a short τ , i.e., a starburst), followed by the rapid decay of quenching.This is evident by the ratio of the age/τ ∼ 5, which implies the SFRs in the central regions of these galaxies have declined by a factor of exp(−age/τ ) ∼ exp(−5) ≈ 1/150.Meanwhile, the outer regions of galaxies follow a more steady, slowly varying starformation history represented by a longer τ and prolonged ages, where the SFR has declined from its peak by at most exp(−age/τ ) ∼ exp(−2/2) ∼ 1/3.
We note that the values of EW(Hα) and sSFR from the BAGPIPES models do not change with dust attenuation (i.e., changing A(V )), nor does changing the metallicity have a strong impact on our results.We explore the effects of dust and metallicity in the next subsection.
We also considered how changes to the parameterization of the SFH could impact our interpretation.Our results depend on the strength of the Hα emission and the SEDderived sSFRs, where (as argued above) these are sensitive to changes to the SFH when averaged over ∼10 Myr and ∼100 Myr timescales, respectively.The delayed-τ models, while simplistic, allow for such changes in the relative SFR (i.e., models with short τ values act as bursts or sudden quenching).The impact of more complex SFHs would be minimal when averaged on these timescales.Indeed, Cutler et al. (2023) reconstructed the spatially resolved starformation histories of the center and outskirts of 66 galaxies at z ∼ 2.3 using the "flexible" star-formation histories with PROSPECTOR (Johnson et al. 2021), combining information from HST Spitzer/IRAC, ground-based photometry, and constraints on the gas-phase metallicity from MOS-DEF spectroscopy.They found that for the majority of their galaxies, the outskirts experienced smoother star-formation histories, while the central regions experienced more recent (<100 Myr) bursts.This is similar to our findings and implies our interpretation is qualitatively robust against variations in the assumption of the parameterization of the SFH.
This scenario where the galaxy centroids form in shorterduration star-formation episodes can further account for observations of galaxy sizes measured from spatially resolved imaging in the mid-IR.The mid-IR probes emission from polycyclic aromatic hydrocarbons (PAHs), a tracer of star formation from heating in photo-dissociation regions, (e.g., Ronayne et al. 2023).Early studies with JWST/MIRI (Magnelli et al. 2023;Shen et al. 2023) show that the effective radii of galaxies at 0.2 < z < 2.5 over a range of stellar mass measured in the mid-IR emission are more compact (with smaller effective radii) than either the light in the rest-frame optical/near-IR (which trace stellar light), the near-UV light (which traces the direct continuum from young stars with ∼100 Myr lifetimes), or the Hα emission.One explanation for this is that the mid-IR emission tracks star formation over longer timescales than either the UV or Hα emission, as longer-lived stars (with lifetimes of 500 Myr -1 Gyr) can heat PAHs efficiently (Salim et al. 2009;Kennicutt & Evans 2012;Salim & Narayanan 2020).In this scenario, the mid-IR emission in the central regions of galaxies may be a remnant of the previous starburst episodes, whereas the Hα emission, which traces only the most massive, shortest-lived stars, declines faster.This implies that the galaxy centroids were actively forming in the past ≲ 1 Gyr, and that the SFR has declined since then.If this is correct, then the galaxies in our sample should show similar mid-IR emission concentrated in their centers.This can be tested with studies with JWST/MIRI in the future.
In the outer regions of galaxies, star formation appears to be more steady-state and slowly varying to account for the distribution of EW(Hα)-sSFR seen in our galaxies.Again, this is in line with the "inside-out" formation of galaxy disks, and has been inferred from past studies of spatially resolved Hα sizes versus the sizes of galaxies in their stellar light at the redshifts of our sample (e.g., Tacchella et al. 2015;Nelson et al. 2016a;Wilman et al. 2020;Matharu et al. 2022).Our results support this scenario, with the added detail that this star formation acts over longer timescales than that responsible for the formation of galaxy centroids (see also Cutler et al. 2023).
Moreover, our finding is reminiscent of galaxy bulge formation (e.g, Zoccali et al. 2003;Ballero et al. 2007).The majority of galaxies in our sample have positive EW(Hα) gradients and positive/flat sSFR gradients (Figure 7), implying that we may be witnessing the end of the formation of galaxy bulges in the galaxies in our sample.Theoretical work sug-gests that the formation of dense stellar cores (e.g., bulges) can stabilize gaseous disks, leading to a local suppression of star-formation (e.g., Bournaud et al. 2007;Martig et al. 2009;Genel et al. 2012;Sales et al. 2012;Dekel & Burkert 2014;Tacchella et al. 2019;Semenov et al. 2023).The center regions of galaxies may experience a higher number of discreet starbursts as a result of multiple minor mergers, counter-rotating streams, recycled gas, or tidal compression, which preferentially cause gas to flow toward the galaxy center, triggering star-formation (e.g., Zolotov et al. 2015;Tacchella et al. 2016a,b).The build-up of these galactic cores may ultimately lead to the formation of the bulge and also to a decrease or cessation of star formation.Indeed, the distribution of metallicity in the Milky Way stars shows evidence for such a formation scenario that bulge/thick disk formed through multiple discreet bursts then quench rapidly, where there is more gradual/smooth star-formation in the thin disk (e.g., Zoccali et al. 2003;Ballero et al. 2007;Queiroz et al. 2021Queiroz et al. , 2023)).Therefore, discrete bursts of star formation would account for the EW(Hα)-sSFR and EW(Hα)-age distributions in the innermost galaxy regions in our sample, and possibly indicate bulge formation.

Systematic Effects of Dust and Metallicity
In this section, we discuss other possible effects from dust attenuation, metallicity, and their radial distributions on the measurement of EW(Hα), sSFR, and age and their gradients.

Dust Attenuation
There are two ways that systematics associated with dust attenuation could impact our interpretation.The first is that there is a well-known degeneracy between the stellar population age and dust attenuation in SED fitting.Second, as mentioned in the introduction, while the EW(Hα) is generally not affected by dust, there may exist higher attenuation from the birth cloud on the nebular emission, compared to the stellar continuum.We consider both of these effects here.
The degeneracy between age and dust attenuation in the SED fitting is most problematic when only rest-frame UVoptical data are available (e.g., Papovich et al. 2001;Conroy 2013).That is, the same photometry data can be fitted with either a younger model and more extinction or an older model with less extinction.However, we mitigate against this degeneracy as we use imaging from HST and JWST that covers 0.4 to ∼5 µm.For the galaxies in our sample, this corresponds to coverage through the rest-frame K-band, where the effects of dust and stellar population ages are separable (e.g., see next paragraph).Furthermore, we found no strong correlation between dust attenuation and mass-weighted age from the results of spatially resolved SED fitting, indicating the SED fitting results are largely not affected by such degeneracy.
As an additional check, the age -dust degeneracy can be broken when including mid-/far-IR data.To test this, we repeated the SED fitting using the integrated photometry of the galaxies in our sample including longer-wavelength data in the mid-IR and far-IR from IRAC/Spitzer at 5.8µm and 8 µm, MIPS/Spitzer at 24 µm, and Herschel/PACS at 70, 100, 160 µm.We cross-matched our catalog with the GOODS-Herschel catalog (Elbaz et al. 2011).We identify 9 galaxies in our sample with > 3σ MIPS 24µm detections.For these 9 galaxies, we reran CIGALE SED fitting using the same parameters as in section 3.1, but now including the additional data.We find that when including these mid-IR and far-IR data, the sSFR values decrease slightly, with a median change of 0.07 dex.The mass-weighted ages increase slightly, with a median of 0.03 dex.These changes are within the error budget derived on the same quantities without the mid-IR and far-IR data, and some of these changes can be attributed to the associated uncertainties on the mid-IR and far-IR data.
Because we use these same data to derive the spatially resolved SED fits, we expect that our values for the resolved sSFR and ages are not seriously impacted by uncertainties in dust attenuation.In addition, the spatially resolved dust attenuation profiles are smooth for the majority of our galaxies.To further test the effect of dust attenuation, we refit the SEDs for each Voronoi bin for each galaxy while fixing the dust attenuation law to be the same as that from the fit to the galaxy's integrated SED.On average, this increases the sSFR gradients by a median value of 0.08 dex kpc −1 and decreases the age gradients by a median value of 0.06 dex kpc −1 .These are within our error budget and the correlation between the EW(Hα) and sSFR gradients remains significant with a pvalue of 0.03.
Therefore, within the limitations of our dataset, our results are relatively robust to dust attenuation.We expect future analyses including spatially resolved observations of the Balmer decrement to correct for EW(Hα) and/or modeling that includes JWST/MIRI can probe fully the impact of dust attenuation on our results.
One problematic scenario would be if the Hα emitting regions experience significantly higher dust attenuation than the continuum.This would decrease the Hα emission more, lowering the EW(Hα).
There are also some theoretical predictions that dust could directly affect Lyman-continuum photons, which would also possibly lead to a reduction of Hα (Tacchella et al. 2022).Previous studies find an increase in dust attenuation toward the center of galaxies at 1 ≲ z ≲ 3 (e.g., Nelson et al. 2016b;Wang et al. 2017;Tacchella et al. 2018;Matharu et al. 2023), but this does not necessarily affect EW(Hα) unless the dust law itself is also changing, or if there exists a deeply obscured, optically thick component.In fact, Hemmati et al. (2015) studied the attenuation ratio between nebular and stellar in eight galaxies at z ∼ 0.4 and found a flat radial profile of the attenuation ratio in the ma-jority of their galaxies.Although this result is based on a small sample, future studies of resolved hydrogen recombination lines (e.g., Hα/Hβ) in a large sample would allow constraining the radial profile of attenuation ratio between nebular and stellar and the dust-corrected EW(Hα) radial profile.Furthermore, for the cases which our galaxies have mid-IR and far-IR detections, we observe only minor differences in the resulting dust attenuation and sSFR values from model fits to the integrated SEDs.Therefore, it is unlikely that these galaxies have dominant, obscured components missed by the SED modeling that would impact our interpretation.Nevertheless, future studies of hydrogen recombination lines at longer wavelengths (e.g., Paschen lines) would allow for studies of dust attenuation affecting the nebular regions directly.

Metallicity
Another possible concern is if there are strong metallicity gradients in our galaxies that alter the interpretation of the EW(Hα) and sSFR gradients.However, we expect this to be a minor effect.The gas-phase metallicity gradients of galaxies at 0.6 < z < 2.6 are found to be flat for galaxies with M * ≳ 10 10 M ⊙ and slightly positive metallicity gradients (∼0.05 dex/kpc) for galaxies with lower masses (e.g., Simons et al. 2021) The gradients of gas-phase metallicity could affect the [N II]/Hα ratio across galaxies and affect the EW(Hα+[N II]) gradient.As mentioned in Section 4.1, we measure the EW(Hα+[N II]) as EW(Hα) and assume a constant [N II]/Hα/ across one galaxy.In the case of a positive metallicity gradient of 0.05 dex kpc −1 , the log([N II]/Hα) would increase by 0.09 dex kpc −1 following Pettini & Pagel (2004), which would lead to the EW(Hα+[N II]) gradient increase by 0.02 dex kpc −1 .This change is within the error budget we measure on the EW(Hα) gradients.
In addition, we tested the effects of changing the metallicity in the models, but this had only a minor effect on the results.In Figure 9 we show the shift in the expected EW(Hα) and sSFRs when using models with a gas-phase metallicity of Z/Z ⊙ = 0.6 and a stellar metallicity of Z/Z ⊙ = 0.16 as arrows for models with age = 0.5 Gyr (compared to the canonical Solar metallicity values).These alternative values for the metallicities are adopted based on measurements from galaxies with similar stellar mass and redshift (Sanders et al. 2021 andChartab et al. 2023).The modeled EW(Hα) values generally increase by a median of ∼ 0.06 dex and increase more with higher EW(Hα) in the range of −0.05 to 0.12.These are smaller by more than a factor of two than the median difference (0.24 dex) between the innermost region (r/r eff ≤ 1) and the outermost region (r/r eff > 1).Thus, although the metallicity cannot fully explain the shift in the EW(Hα) from the galaxy's center to the outskirt region, we cannot rule out the possible effect of the metallicities causing some differences in the sSFR-EW(Hα) histogram, at least partially.We expect future analysis on details of metallicity profiles of these galaxies and using more extended grids of stellar evolutionary models such as Mesa Isochrones and Stellar Tracks (MIST, Dotter 2016;Choi et al. 2016) could help in understanding the effect of metallicity.5.3.Galaxies with negative EW(Hα) gradients: the duty cycle of Bulge Formation Interestingly, two galaxies in our sample show a negative EW(Hα) gradient, a negative sSFR gradient, with a flat age gradient.These two galaxies (NGDEEP_01948 and NGDEEP_01831) show irregular morphological signatures, possibly indicative of a merger or that they have starbursting clumps near their centers (see Fig. 1).
To test for mergers, we measured the Gini-M 20 values for our galaxies using STATMORPH (Rodriguez-Gomez et al. 2019) from the JWST imaging corresponding to rest-frame 1.2 µm (as above).Compared to the separation from Lotz et al. (2008), most of our galaxies (16 out of 19) fall well within the region associated with galaxy disks, including one of the irregular galaxies here (NGDEEP_01948).Three of the galaxies fall in the region of the Gini-M20 plot associated with mergers, including the other irregular galaxy here (NGDEEP_01831), but they are all very close to the mergerdisk separation line.Therefore, it is inconclusive if these two irregular galaxies are mergers.
Inspecting the spatially resolved sSFR and age maps of these two galaxies, NGDEEP_01948 and NGDEEP_01831 (see Figure 8), they both have centers with high sSFR and low ages, suggesting rapid growth.Therefore we favor the scenario where these galaxies have recently undergone a nuclear starburst, which could result from a recent gas flow, or as a massive star-forming clump(s), which may form and migrate from the disk (Dekel et al. 2009a,b;Zolotov et al. 2015;Tacchella et al. 2016b).Other galaxies at these redshifts indicate possible central gas accretion events based on "inverted" metallicity gradients, which could be taken as evidence for metal-poor cold-gas inflows from the ISM or circumgalactic medium depositing directly onto the centers of galaxies (e.g., Wang et al. 2019;Simons et al. 2021).The latter are rare (occurring in less than 10% of the samples), which is consistent with the two galaxies we find with central regions with elevated sSFRs and high EW(Hα).
Therefore, these two galaxies may be the "exceptions that prove the rule" in our interpretation of the EW(Hα), sSFR, and age gradients.If we return to the argument from Section 5.1 that bulge formation is episodic, with a short duty cycle of star formation, then we may have caught these two galaxies in a growth phase of their central regions.Taking the ratio of the number of galaxies with forming centroids (two) to the total (19) within the period of cosmic time (4.9 Gyr) spanned by the redshift range of our sample (0.6 < z < 2.2), this duty cycle is approximately ≈ (2/19) ≈ 10%, and the timescale for this period is 10% ×4.9 Gyr = 500 Myr.This is remarkably similar to the ages of the galaxy cores in our sample, supporting our argument that we are witnessing bulge formation in this sample.
Taken together, our data imply that most (16 out of 19) bright, star-forming galaxies at 0.6 < z < 2.2 form stars in disks that grow in an "inside-out" fashion.At the same time, they are forming their bulges through multiple, shortduration bursts.In some cases (two out of 19 in our sample) there are indications that we have caught galaxies during a formation stage, likely as a result of a recent increase in the nuclear gas-accretion rate or via clump migration.

SUMMARY
We study the Hα equivalent width maps of 19 galaxies at 0.6 < z < 2.2 from NIRISS slitless spectroscopy as part of the NGDEEP survey.Our sample is dominated by bright galaxies that lie along the star-formation main sequence with a stellar mass range of 10 9 − 10 11 M ⊙ .We combine these data with ancillary deep HST and JWST imaging from CAN-DELS, JADES, and JEMS, from which we perform spatiallyresolved SED fitting with 21-band fluxes for these galaxies and construct their specific SFR and mass-weighted age maps.The resolved SED fits and Hα maps resolve angular structures in galaxies on scales of 0. ′′ 18, allowing quantitative measurements of nebular emission lines and galaxy properties at a spatial resolution of ∼1 (proper) kpc.
The EW(Hα) and sSFR are independent measurements, derived from the NIRISS grism data and multi-band SED fitting, respectively.They both trace the current star formation activity, and the mass-weighted ages from the SED fitting probe the time since the last major star formation episode.We compare the spatially-resolved EW(Hα), sSFR, and mass-weighted age and their radial profiles.Our results are summarized as follows: • We find strong correlations between EW(Hα) and sSFR, and an anticorrelation between EW(Hα) and the mass-weighted age.This is true for both the integrated emission from galaxies and on a pixel-by-pixel scale.Therefore, spatially within galaxies, the EW(Hα) increases with increasing sSFR and the EW(Hα) increases with decreasing age.
• Quantifying the radial profiles of EW(Hα), sSFR and mass-weighted age, most galaxies show positive EW(Hα), positive sSFR gradients, and negative age gradients.The median values of the gradients in EW(Hα), sSFR, and age gradients are 0.09 +0.03 −0.06 dex kpc −1 , 0.003 +0.049 −0.038 dex kpc −1 and −0.04 +0.04 −0.03 dex kpc −1 .This indicates that, for most galaxies, their outer regions exhibit higher levels of star formation and appear to be younger when compared to their central regions.This is consistent with the inside-out formation of galaxies.
• We compare our measured EW(Hα) and sSFR to stellar population obtained from a set of star formation histories following a delayed-τ model using BAGPIPES.We find that these models roughly span the observed distribution of EW(Hα) and sSFR and EW(Hα) and age.We then compare the distributions of EW(Hα) and sSFR to the models as a function of galactocentric radius (normalized by the galaxy effective radii).The central regions of galaxies (r/r eff ≤ 0.5) favor models with a small age (∼500 Myr) and rapidly declining star-formation histories with τ ∼100 Myr.In contrast, the outer regions (r > r eff ) favor models with slowly evolving star-formation histories with larger τ ∼2 Gyr, and a wide range of age (∼0.5-2Gyr).We argue that this indicates that while most galaxies grow their disks in an "inside-out" fashion, their cores (i.e., bulges) grow through multiple, short-duration bursts.However, there are important systematics that could impact this result, including assumptions about the dust law, metallicity gradients, and parameterization of the star-formation history.We discuss these uncertainties and offer predictions that our results make that can be tested in future studies (in Section 5).
• Two galaxies show negative EW(Hα) and sSFR gradients with a flat (i.e., constant) age gradient.These two galaxies show irregular morphologies, but we do not find strong evidence of a merger from their morphologies (using Gini-M20).We, therefore, conclude these galaxies are experiencing central star-forming episodes, either as a result of recent gas accretion in their cores or through the migration of clumps from their disks.Comparing the ratio of the number of galaxies in our sample with active central star formation (two) to the total (19), gives a short duty cycle of 10%, albeit with a large statistical uncertainty.This suggests that bulge formation may occur at these redshifts in multiple star-formation episodes of < 500 Myr.Based on the age of Galactic bulge stars and the redshift of our sample, it is possible that we are witnessing the end of the period of star formation in the galaxies in our sample.
This paper demonstrates the capability of JWST/NIRISS slitless spectroscopy data to construct the nebular emission line maps of high-redshift galaxies, which can be combined with broad-band imaging to gain useful constraints on the spatially resolved formation of galaxies during the peak of the cosmic SFR density.In the future, we will expand on this work by adding the second half of the NGDEEP NIRISS observations, which will provide a large sample of high-redshift galaxies with secure nebular emission lines (Hα, Hβ, [O II], [O III]) maps for star formation, dust attenuation, and metallicity studies, which can inform a more complete picture of galaxies' growth and quenching.We also expect to expand this work using larger samples of galaxies with broad-band imaging, extending to the mid-IR (with JWST MIRI) and WFSS from NIRISS.
We acknowledge the hard work of our colleagues in the NGDEEP collaboration and everyone involved in the JWST mission.This work benefited from support from the George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University.CP thanks Marsha and Ralph Schilling for generous support of this research.This work is based on observations made with the NASA/ESA/CSA JWST.The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.These observations are associated with program #2079.Some/all the JWST data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.The specific observations analyzed can be accessed via Leung, Gene et al. (2023), Williams, Christina et al. (2023) and Rieke, Marcia et al. (2023).JM is grateful to the Cosmic Dawn Center for the DAWN Fellowship.PGP-G acknowledges support from grant PGC2018-093499-B-I00 and PID2022-139567NB-I00 funded by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033,FEDER, UE.The work of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
Figure1.False color images and Hα maps for the galaxies in our galaxy sample.For each galaxy, we show the RGB images using the HST/ACS F814W, JWST/NIRISS F150W and JWST/NIRCam F444W (left), and the Hα emission line map constructed from the NIRISS slitless spectroscopy (right).The sample is ordered by increasing redshift.Each image is 2 ′′ × 2 ′′ , corresponding to 13.4×13.4kpc 2 and 16.6×16.6kpc 2 at z = 0.6 and z = 2.2, respectively.The text inset lists the galaxy ID, redshift, and whether galaxies are identified with AGN.

Figure 3 .
Figure3.Left and Middle: Comparison between the stellar mass and star-formation rate derived by summing the values from the bins in the spatially-resolved "map(s)" with the same aperture as used for integrated photometry and the "integrated" values derived from integrated photometry.Right: Comparison between the Hα equivalent width average from the EW(Hα) "map(s)" and the "integrated" EW(Hα) measured by GRIZLI from the total 1D spectra.Note that statistical errors on the EW(Hα) values from the maps and from GRIZLI are smaller than the data points.

Figure 4 .
Figure 4. Relation between spatially resolved sSFR and EW(Hα) (left), and spatially resolved age and EW(Hα) (right).The 2D histograms include all pixels in the spatially resolved images for galaxies in our sample.In each panel, the orange dots show the medians of the values and errors derived from the individual pixels.The red crosses show values derived from the integrated photometry.The median error of individual, spatially resolved pixels is indicated by the blue error bar in each panel.For comparison, the dashed lines show the relations derived from MOSDEF galaxies at z = 1.4 − 3.8 (Reddy et al. 2018) (median and 1σ scatter).The median values of our spatially resolved measurements broadly agree with the results from the integrated emission of our galaxies.However, these are offset slightly to lower EW(Hα) compared to the relation from MOSDEF.

Figure 5 .Figure 6 .
Figure5.Example of the EW(Hα), sSFR, and age spatial maps and their radial gradients for galaxies in our sample.The sSFR and age maps use the Voronoi tessellation (see Sect. 3.2).In the left panels, the red crosses mark the luminosity-weighted galaxy centroid derived from the F160W image.Orange contours are 5, 10, 25, 50, 100 σ in the F160W image.The example here shows a galaxy with no AGN indication (NGDEEP_ 01614).The middle panels show the radial gradient of the EW(Hα), sSFR, and age.The purple solid lines show best-fit linear relations, along with 400 random draws from the posterior.Grey dots in the sSFR and age panels show Voronoi bins that lack EW(Hα) measurements, which are excluded in the linear fitting.The right panels show SED fits at three different radial bins ("inner", "middle", and "outer" regions, top-to-bottom), indicated by the red boxes in the middle panels.Each panel shows the observed photometric fluxes with errors (purple), the CIGALE-derived best-fitted model (black curve), and model photometry (red dots).The best-fitted CIGALE model is the sum of the contributions from a dust-attenuated stellar emission and nebular emission (the unattenuated stellar model is shown in blue).A galaxy with an AGN is shown in Figure6.The maps and radial gradients of the remainder of our galaxies are shown in Appendix A.

Figure 7 .
Figure 7.The radial gradients of EW(Hα) versus sSFR (left), and the radial gradients of EW(Hα) versus stellar mass-weighted age (right).Here the gradients are the slope measured for the radial gradients for each galaxy (see Figure5).In each plot, the histograms show the distribution of EW(Hα), sSFR, and age gradients for the individual galaxies.The markers in the scatter plots are color-coded by the distance of the galaxies from the SFMS (positive/negative values are above/below the SFMS).Galaxies identified with AGN are marked by open red diamonds and orange squares.The best-fit linear relations are shown as the purple dash lines along with 400 random draws from the posterior (faint grey lines).

Figure 9 .
Figure 9. Left panel: The grids of sSFR versus EW(Hα) from a set of delayed SFH models overlaid on the 2D histograms of sSFR versus EW(Hα) for all pixels.The diamonds show expectations from SED models with delayed-τ star-formation histories at different ages and with different τ values.The colors of these crosses correspond to models of different ages and τ values.The models with the same age are linked with colored, solid line.The τ values are labeled by each model τ =0.1, 0.2, 5, 1, 2 Gyr.The 2D histogram is the same as Figure 4 but in grey.All models assume Solar metallicity.The arrows on models at an age = 0.5 Gyr show the shift in the models if we adopt a gas-phase metallicity Z/Z⊙ = 0.6 and a stellar metallicity Z/Z⊙ = 0.16 (see Section 5.2.2).The black error bar shows the median error of individual, spatially resolved pixels.Right panels: the 2D histograms for pixels in four bins of galactocentric radius, normalized by the galaxies' effective radii (reff) measured at rest-frame 1.2 µm (see text).The radii of each bin are listed in each panel.The colored lines show the same models as in the left panel.

Table 1 .
Parameters used in the SED fitting with CIGALE.