REQUIEM-2D: A diversity of formation pathways in a sample of spatially-resolved massive quiescent galaxies at z~2

REQUIEM-2D (REsolving QUIEscent Magnified galaxies with 2D grism spectroscopy) is comprised of a sample of 8 massive ($\log M_*/M_\odot>10.6$) strongly lensed quiescent galaxies at $z\sim2$. REQUIEM-2D combines the natural magnification from strong gravitational lensing with the high spatial-resolution grism spectroscopy of \emph{Hubble Space Telescope} through a spectrophotometric fit to study spatially resolved stellar populations. We show that quiescent galaxies in the REQUIEM-2D survey have diverse formation histories manifesting as a gradient in stellar ages, including examples of (1) a younger central region supporting outside-in formation, (2) flat age gradients that show evidence for both spatially-uniform early formation or inside-out quenching, and (3) regions at a fixed radial distance having different ages (such asymmetries cannot be recovered when averaging stellar population measurements azimuthally). The typical dust attenuation curve for the REQUIEM-2D galaxies is constrained to be steeper than Calzetti's law in the UV and generally consistent with $A_V<1$. Combined together and accounting for the different physical radial distances and formation time-scales, we find that the REQUIEM-2D galaxies that formed earlier in the universe exhibit slow and uniform growth in their inner core, whereas the galaxies that formed later have rapid inner growth in their inner core with younger ages relative to the outskirts. These results challenge the currently accepted paradigm of how massive quiescent galaxies form, where the earliest galaxies are thought to form most rapidly. Significantly larger samples close to the epoch of formation with similar data quality and higher spectral resolution are required to validate this finding.

we find that the REQUIEM-2D galaxies that formed earlier in the universe exhibit slow and uniform growth in their inner core, whereas the galaxies that formed later have rapid inner growth in their inner core with younger ages relative to the outskirts. These results challenge the currently accepted paradigm of how massive quiescent galaxies form, where the earliest galaxies are thought to form most rapidly. Significantly larger samples close to the epoch of formation with similar data quality and higher spectral resolution are required to validate this finding.

INTRODUCTION
One of the key empirical findings in observational extragalactic astronomy is the existence of a starformation main sequence (SFMS) manifesting as a population of galaxies with a relatively tight correlation between their star-formation rate (SFR) and stellar mass (e.g. Elbaz et al. 2007;Daddi et al. 2007;Noeske et al. 2007;Whitaker et al. 2012Whitaker et al. , 2014Speagle et al. 2014, see also Figure 1). It is also shown observationally that there is a population of massive galaxies with significantly lower SFR than what is expected from the SFMS, known as quiescent galaxies (e.g. Brammer et al. 2009;Kriek et al. 2006;Muzzin et al. 2013). Understanding the existence of the SFMS (e.g., Kelson 2014) and its outliers is one of the major challenges of extragalactic astronomy.
Spectrophotometric observations, that include UV to far IR data, are extensively used to study quiescent galaxies at z 1 (e.g., Geier et al. 2013;Hill et al. 2016;Belli et al. 2017;Toft et al. 2017;Abramson et al. 2018;Ebeling et al. 2018;Newman et al. 2018;Morishita et al. 2018Morishita et al. , 2019Belli et al. 2019Belli et al. , 2021Estrada-Carpenter et al. 2019Jafariyazani et al. 2020;Setton et al. 2020;Man et al. 2021;Tacchella et al. 2021;Whitaker et al. 2013Whitaker et al. , 2021. These studies imply a diversity among the progenitors, including compact galaxies that form stars centrally (e.g., Williams et al. 2014) and rotationally supported larger star-forming galaxies (e.g., Wisnioski et al. 2015). Quiescent galaxies are thought to form via either "slow" or "fast" channels, leading to different colors, formation time-scales and the ages of the stellar populations (e.g., Belli et al. 2019). The cold gas content of quiescent galaxies has also been studied to some extent, revealing a medley of values for the inferred depletion timescales and gas to stellar-mass ratios (e.g., Belli et al. 2021;Bezanson et al. 2019;Man et al. 2021;Spilker et al. 2018;Suess et al. 2017;Whitaker et al. 2021;Williams et al. 2021). These early studies find cases that support both fast depletion of the cold gas and/or reduced star-formation efficiency associated with the quiescent population.
There is an equally diverse set of potential mechanisms for explaining the existence of quiescent galaxies in cosmological simulations (e.g., Dekel & Birnboim 2006;Tacchella et al. 2015;Wellons et al. 2015;Zolotov et al. 2015;Tacchella et al. 2016), motivating observa-tional studies in order to constrain the imprints of the corresponding physical processes on their stellar populations. Specifically, constraining the spatial variation of stellar age within the galaxies can shed light on their formation pathways (e.g., Wellons et al. 2015;Tacchella et al. 2016). However, as a stellar population gets older, the variations in the age-sensitive spectrophotometric features decrease, making it very difficult to detect age gradients for an older age baseline. This limits the ideal candidates to "recently" 1 quenched quiescent galaxies, which peak in number density at z 1 (e.g., Whitaker et al. 2012;Cassata et al. 2013).
We present the REQUIEM-2D survey in this paper, studying the stellar populations of quiescent galaxies that cover a range of global ages, stellar masses and starformation rates (see Akhshik et al. 2020, and Section 2). This survey is designed to use a specific feature of the HST Wide Field Camera 3 (WFC3) grism, its ability to preserve the spatial information perpendicular to the dispersion angle, to study spatially-resolved stellar populations. When combining the high spatial-resolution of HST with a natural boost from strong gravitational lensing, the REQUIEM-2D survey provides a unique opportunity to constrain the spatially-resolved ages and star-formation histories of the inner cores of 8 massive quiescent lensed galaxies in order to gain further insights into their formation pathways (see also Abramson et al. 2018, for a similar study).
Herein, we present the analyses of the 7 targets from the REQUIEM-2D survey; a detailed study of the REQUIEM-2D's pilot target, MRG-S0851 is presented elsewhere (Akhshik et al. 2020Caliendo et al. 2021). We summarize the data reduction in Section 2, and present the analysis in Sections 3.1 and 3.2, referring the reader to Akhshik et al. (2020) for a more in depth discussion of the adopted methodology. In Sections 4 and 5 we discuss the implications of our results. The details of the previously unpublished lensing model for MRG-P0918 are presented in Appendix A.
In this paper, we assume a standard simplified ΛCDM cosmology with Ω M = 0.3, Ω Λ = 0.7 and H 0 = 70 kms −1 Mpc −1 and the Chabrier (2003) initial mass  Momcheva et al. 2016). The stellar populations' parameters are constrained using Prospector-α Leja et al. 2017Leja et al. , 2020. The red line shows the best fit to the SFMS at z = 1-3 following Leja et al. (2021); Whitaker et al. (2014). The REQUIEM-2D galaxies are shown with the circles. For the REQUIEM-2D galaxies, the stellar masses are the total stellar mass formed in all spatial bins and the SFR is sum of the SFR of individual bins (see Sections 2-3.2). MRG-M2129 is on the SFMS because it has a dust-obscured star-forming component (see Section 4). Both parameters are corrected assuming a linear gravitational magnification.
function. All magnitudes are reported in the AB system.
We use the Mikulski Archive for Space Telescopes (MAST) to access all publicly available HST Wide Field Camera 3 (WFC3) and Advanced Camera for Surveys (ACS) data (see Table 1) along with the REQUIEM-2D Cycle 26 HST program (HST-GO-15663). All HST and Spitzer data is reduced using the "Grism redshift & line analysis software for space-based slitless spectroscopy", or Grizli . Grizli performs astrometric calibrations of the WFC3-IR and WFC3-UVIS images using the Pan-STARSS (Flewelling et al. 2016) and the Gaia-DR2 catalog (Gaia Collaboration et al. 2018;Lindegren et al. 2018). The final drizzled mosaics are constructed using AstroDrizzle(Avila & Hack 2012) with a pixfrac of 0.33 and a pixel size of 0. 1.

Constraining Photometric Fluxes and Morphologies
After reducing HST and Spitzer data, we measure photometric fluxes that are major components of our joint spectrophotometric fitting. We first construct HST point spread functions (PSFs) using Grizli. In this step, Grizli registers predefined empirical PSFs, retrieved from MAST (Anderson 2016), of each filter for the existing mosaics using the spatial resolution and orientation angles of the contributing exposures. We also construct Spitzer PSFs in a similar way to the HST PSFs. To build our Spitzer PSFs, we first utilize the rich data set in the GOODS-S field and construct the 'base' empirical PSFs in each filter in a 0. 1 pixel scale. Then, we mimick the visits of the telescope over the RE-QUIEM fields, building our position-dependent PSFs on a 5 × 5 grid, in a box with (a minimum of) 30 on a side around the target. 1 Figure 2. Reduced HST observations of the REQUIEM-2D sample. The color images (left) are constructed using HF160W as red, either YF110W or JF125W as green, and UF390W, or IF814W or VF555W as blue. The green contours show one of the flexible banana apertures that we use to measure photometric fluxes. The blue lines in each sub-panel on the left define the five central bins used in the spatially-resolved analysis. The two right sub-panels show data for the two grism dispersion angles used to observe each target in order to minimize contamination. The wavelength is increased from left to right. The blue lines in the right sub-panels represent the equivalent extent of the 5 central bins.
Using drizzled mosaics and weight maps, we measure the Spitzer and HST photometric fluxes and uncertainties using a methodology combining earlier works by Akhshik et al. (2020), Wuyts et al. (2010), and Skelton et al. (2014). For all targets, we construct a Galfit (Peng et al. 2011) model in the H F160W band for a 20 × 20 box around each target, modeling all galaxies that can potentially contaminate our target (for MRG-M0138 the box is 37. 5 × 19 to accommodate its extended morphology). In each box, we have several light profiles, modeled using different Galfit components, and they are used to clean the contamination in all other HST bands after matching the PSF, allowing them to have different overall flux normalizations. We finally use these models to remove the light profiles of all surrounding objects from the images.
Conventional circular-aperture-photometric methods fail to provide a reliable estimate for most targets in the REQUIEM-2D survey due to their extended and somewhat irregular morphologies. We therefore use a variation of "object apertures" following Wuyts et al. (2010): the axis of symmetry is determined by fitting a polynomial function to isophotes of each galaxy, re-quiring that the distances of opposing points from the axis of symmetry be minimum. We inspect each fit to make sure that the axis of symmetry is satisfactory. Apertures with increasing radii are centered on this axis defining banana-shaped apertures. These apertures are then used to measure the flux densities (see Figure 2).
It is a standard practice to measure the flux density using a smaller 0. 7 diameter aperture and implement a correction factor to total (e.g., Skelton et al. 2014), as using larger apertures does not optimize the signal to noise ratio. Thus, we pick an aperture with a radius of 0. 35 √ µ, where µ is the average linear gravitational magnification of the different images of the galaxy should it be multiply-imaged 2 . The aperture correction is then defined to be the ratio of the flux measured using the 2 Using smaller apertures are not suitable when the galaxy has color gradients. The radius here is only measuring distances perpendicular to the axis of symmetry and the full size of the banana aperture is indeed larger (see Figure 2). We also use the average gravitational magnification to avoid making the aperture large perpendicular to the axis of symmetry in case of extreme magnifications such as MRG-M1341. same aperture with a radius of 0. 35 √ µ on the noise-free Galfit model to the total flux.
To estimate the uncertainties of the HST aperture photometry, we first define a set of empty apertures with increasing radii for the different filters of each target, measuring their fluxes. We then fit a power-law function to the width of the flux distribution and use this scaling relation to estimate the noise following Whitaker et al. (2011). For the Spitzer IRAC channels, we use the square root of the sum of the inverse weight-map pixels as the uncertainty.

Analysis
We analyze spectrophotometric data using the requiem2d code 3 , referring the reader to Akhshik et al. (2020) for more details. The fully Bayesian model is able to constrain the spatially resolved ages and star-formation histories (SFHs) of massive distant galaxies using spectrophotometric observations that assume a nonparametric SFH. We choose the grism dispersion angle such that it is perpendicular to the axis of symmetry while also minimizing the contamination from nearby objects. We define 7 spatially resolved bins per galaxy using the H F160W image for our analyses following Akhshik et al. (2020) (see Figure 2). Each one of the five central bins are defined to be 3-4 pixels wide in the image plane (0. 3-0. 4) noting that this resolution is greater than the FWHM of the PSF of 0. 18. The two outer bins include all of the remaining pixels within the segmentation map on each side.
One important difference in the analyses here and that of Akhshik et al. (2020)  , whereas we find that this approach in-troduces some artifacts in the recovered SFHs that are produced by different built-in FSPS assumptions. One notable example is the age when dust1 is added to the model: all stars younger than this time scale, with a default value of 10 7 years, are attenuated by both dust1 and the usual dust2 (Conroy & Gunn 2010). It can be shown that at around this time scale, some slight enhancements/depressions manifest themselves in the SFH (e.g., Figure 1 in Akhshik et al. 2021, around 0.01 Gyr). We explicitly check this effect for MRG-P0918 by changing this time scale from 0.01 to 0.015 and 0.02 Gyr and refit the model. The artifacts move in a consistent manner appearing at the assumed dust1 time scale. To avoid adding artifacts into the recovered SFH, we instead define 20 time bins for all samples in the model 4 . We however note that the requiem2d code is fully capable of fitting the finest lookback time resolution (∼ 70-90) produced by FSPS if required.
For those targets with ALMA 1.3 mm detections, MRG-M0138 and MRG-M2129, we carefully examine the distribution of the dust continuum emission. In both sources, the dust emission is unresolved in the ∼1-1. 5 ALMA data . While this resolution is significantly coarser than our spatially resolved bins, the fact that the emission is pointlike implies that there is some information about the spatial distribution of the 1.3 mm flux in the existing data. We performed a number of tests to determine the level of dust emission in each resolved bin that would be consistent with the existing ALMA data, jointly fitting for the 1.3 mm flux density of a series of unresolved point sources with positions fixed to the centers of each of our resolved bins. This analysis suggests that for both MRG-M0138 and MRG-M2129, the 1.3 mm emission is associated with the three central bins, though the distribution of flux among these bins is essentially unconstrained. In other words, the existing ALMA data constrain the total dust emission in these sources and localize it to the central three bins in each source. Therefore, in our fitting procedure we sum the predicted 1.3 mm flux density of the central 3 bins and compare this to the observed 1.3 mm emission, while for the exterior bins we use only upper limits on the millimeter flux. We implement the ALMA semi-resolved flux in the Bayesian model using the integrated likelihood (e.g., Section 4.3 of Stan Development Team 2020), i.e. for those bins with only an upper limit detection, the integral of the flux likelihood up that up-per limit is included in the Bayesian fitting. For other targets that are undetected in ALMA 1.3 mm, we can in principle use the upper limit but as these targets are not particularly dusty (see Figure 5), adding the upper limit would not be additionally constraining.
We use the Noll et al. (2009) parameterization for dust, assuming a Kriek & Conroy (2013) dust model (e.g., Akhshik et al. 2020). The dust model has three main parameters: (1) dust1 that describes the extra attenuation for stars younger than 10 7 years, (2) dust2, describing the attenuation for stars older than 10 7 years, and (3) dust index, that is the index of the power-law correction of the attenuation curve.
For each target we perform a separate stellar population fit using Johnson et al. (Prospector-α;; Leja et al. (Prospector-α; to the semi-resolved photometric data, and subsequently use these posteriors in the joint spectrophotometric fit (Akhshik et al. 2020). We divide the metallicity and dust2 posterior (its 3σ width) into 16 boxes, generating a series of FSPS templates for each box (Akhshik et al. 2020). By refitting the weight of each box along with the SFHs using requiem2d, we practically use these Prospector-α posteriors as weakly informative priors, and we do not expect this choice of priors affect our age estimates (see Figure 15 in Akhshik et al. 2020, for an example of starting from a different Propsector-α posterior). We adopt a slightly different approach to model MRG-M2129 and MRG-M0138. Because of the ALMA detections, we have non-trivial Prospector-α posteriors for the dust emission parameters (duste umin, duste qpah, duste gamma in FSPS) along with dust2 and metallicity. We therefore perform a principal component analysis for the full Prospector-α posteriors and use two first principal components instead of dust2 and metallicity to generate more representative FSPS templates. We also allow the dust attenuation curve to be fit separately, to ensure that the model is flexible enough to deviate from the Prospector-α posterior to fit the spectrophotometric data. This step is implemented by defining half-normally distributed hyper-parameters, and using them as the standard deviations of the log-normally and normally distributed priors for dust1, dust2 and dust index. Each one of these priors are centered at the corresponding mean values in each one of the 16 boxes and they are used to generate the attenuation curve.
The No-U-Turn Sampler (NUTS; Homan & Gelman 2014), that is wrapped in pymc3 (Salvatier et al. 2016), is used to sample the full posterior. We generate two Monte Carlo chains, drawing 2000 samples from the posterior in each one (with the exception of older galaxies, MRG-M0138 and MRG-M1341, where we draw 3000 samples for each chain). We discard the first half of each chain as burn-in phase and compare the two chains to check for divergences usingR (Gelman & Rubin 1992). After ensuring convergence, the two chains are combined to construct a total of 2000 draws from the posteriors (3000 for MRG-M0138 and MRG-M1341).

Fits
The requiem2d code provides the posterior of different parameters as well as the posterior predictive dis-tributions. We show the spectral energy distributions (SEDs), the 1D collapsed grism spectra, and the posterior predictive distributions of each REQUIEM-2D target in Figure 3 here and Figures 8 to 13 in Appendix B.
The SED fit to the REQUIEM-2D sample is shown in the left panels of Figures 3 and 8-13. We include the Spitzer IRAC channels 1 and 2 as constraints by imposing likelihoods for the global fluxes that are normally Figure 4. The SFH and the evolution of stellar mass, median age and U-V colors for different redshift snapshots. In the SFH panels, the central bin is shown by thick black line. The blue and red lines are the outskirt bins on each side, with the dashed line being the bin right next to the center bin. r/rc indicates the bin's distance from the center normalized to the distance that encloses half-light (see the text, Section 3.2 and Table 4).
distributed with the observed global values and their uncertainties as means and standard deviations and they are modeled by the summed fluxes of the individual bins. In other words, these data points are not spatially resolved like the HST photometric bands. We therefore demonstrate their corresponding wavelength ranges with a grey shade in the SED plots.
The right panels of Figures 3 and 8-13 show the 1D collapsed spectra of the five central bins for each target, extracted following Horne (1986). While we fit the two outer bins for all targets, the fit is not as reliable because of the fainter wings of the galaxy and/or contamination by nearby sources (e.g., Akhshik et al. 2020, also see Figure 2). We therefore do not consider the outer bins when analyzing the radially-varying trends.
For each target, we demonstrate the important absorption lines that drive the age and SFH fit using vertical lines. One notable exception is Na D absorption line in MRG-M1341 (Figure 3). We mask this line since it seems that the stellar absorption cannot fully explain its strength, especially at the center, potentially implying further absorption by the ISM. This is consistent with earlier observations of the Na D line in MRG-M0138 (Jafariyazani et al. 2020) 5 . The masked region is shown with a shaded grey box in Figure 3, right panel; the spectral model in this shaded box is not a fit but is rather an extrapolation. The tension between the data and the best fit model clearly demonstrates that the stellar absorption assumed for the rest of the model cannot fully explain the Na D strength.
The distances reported in Figures 3 and 8-13 are measured simply by finding the distance between the centers of corresponding adjacent bin in arcseconds. This distance is then translated into physical units using the standard cosmology, and finally it is divided by √ µ to estimate the source plane distance. The linear magnification is fairly uniform for all of the targets. We also calculate the distance where half of the total flux is enclosed (r c in Figure 4, see Table 4 for the measured values and their uncertainties) when moving along the axis of symmetry (see Section 2.2) while adding all flux in the perpendicular direction. We prefer to use r c instead of the usual half-light radius, r eff because it better represents our image-plane bins, that are in turn defined according to the grism observational constraints (Akhshik et al. 2020). We nevertheless note that the values of r c are consistent within the uncertainties with the published half-light radii (Man et al. 2021;Newman et al. 2018), justifying this decision. We also use the term "inner core" to refer to r < r c .
To better visualize the evolution of the REQUIEM-2D sample as a function of lookback time, we generate a series of animations that shows the evolution of the different parameters at different redshift/lookback time snapshots (Figure 4). In each snapshot, we use the SFH to generate the SEDs and to derive different parameters. Two main assumptions are made here: (1) we assume that there is no significant evolution of dust, and (2) we assume that stars and dust in the galaxies follow non-evolving distributions. To first order, we know these assumptions do not hold, especially in the case of mergers or starbursts, or if the galaxy experiences a major episode of dust formation/destruction (e.g., Li et al. 2019). However, it is safe to assume that these animations are accurate representations of each galaxy in snapshots closer to the observed redshift; we include the earlier epochs for visual purposes. The stellar population maps are superimposed over the color image of each target in the three bottom panels.

DISCUSSION
The goal of the REQUIEM-2D survey is to measure age and SFH gradients for the eight lensed quiescent galaxies through a spectrophotometric fit. We discuss our fitting algorithm in the preceding Sections, and here we present and discuss the implications of the measured gradients in the stellar populations.
We show the gradients of the different stellar population parameters in Figure 4 and 5. While Figure  4 shows the evolution of parameters such as median age and sSFR for different redshift snapshots, Figure  5 demonstrates the values at the redshift of observation. We can identify three different categories of galaxies based on the detected age gradient patterns when noting the difference in the physical distances: (1) younger in the center, (2) flat, and (3) non-axisymmetric gradients. We choose mass-weighted age as the important parameter because it summarizes the average behavior of SFH. Also, it is suggested that age gradients can be connected to the formation pathway mechanisms of the galaxies (e.g., Wellons et al. 2015).
MRG-S1522 is the only galaxy in the REQUIEM-2D sample which is clearly younger in the inner core (r/r c 1). Its SFH in Figure 4 reveals that the center bin remains more star-forming relative to the outer bins between lookback times of ∼300 Myr and ∼1 Gyr, explaining the younger age at the center. Overall, this galaxy is consistent with an outside-in formation.
Both MRG-M0138 and MRG-M1341 have flat age gradients in the inner core. Even though the age gradients rate across the galaxy, the median age does not show a statistically significant gradient 6 . 6 A careful observer, however, notices that the average medianage is lower in the outskirts consistent with the SFH pattern even if it is not statistically significant.  Figure 6. Measured spatially resolved age differences for the REQUIEM-2D galaxies versus the distance and normalized distances to the half-light radii. δt50 is the difference between the age of the corresponding bin and the center for each target in Gyr, i.e. δt50 = t spatial bin − tcenter and the median age, t50, is defined to be 0.5 = where t0 is the observing lookback time. rc denotes the half light radius, and the black dashed line and shaded region are the best fit and the 1σ width of the linear fit to all data points. The red dashed line is the same fit excluding MRG-S1522 data points, which is the only target showing a significant age difference beyond the half-light radius. β denotes the slope of the linear fit and its 1σ posterior uncertainty with a prior of β ∼ N (0, 100). Data points are color coded with the quenching time-scale of the central bin, defined following Carnall et al. (2018). similar distances. This could indicate that a merger or interaction with another galaxy played a role in their formation. For three galaxies, MRG-M0150, MRG-M0454 and MRG-P0918, the SFH on one side dips below the other side for lookback times of ∼100 Myr to ∼1 Gyr, leading to an older age. The asymmetrical features are not uncommon in spatially-resolved studies of quiescent galaxies (e.g., Li et al. 2015;Setton et al. 2020, Figures  8 and 4, respectively). However, this non-axisymmetric trend vanishes if we azimuthally combine the data points given their distances from the center, consistent with Setton et al. (2020).
The non-axisymmetric age gradient in MRG-M2129 is caused by the centrally concentrated 1.3 mm ALMA flux. One bin remains star-forming all along in MRG-M2129 in order to produce the significant measured 1.3 mm flux. It is interesting that this feature does not exist in the one other galaxy with a 1.3 mm ALMA detection in the REQUIEM-2D sample (MRG-M0138). This is likely because the ratio of the ALMA flux relative to other wavelengths such as H F160W is significantly smaller than that of MRG-M2129. The ALMA detection can be roughly attributed to the 3 central bins. The fitting itself prefers the adjacent to central bin as the dominant producer of the ALMA flux in MRG-M2129. We fit another model assuming that 1.3 mm ALMA flux can only be attributed to the central bin, and the resulting SFH changed to significant star formation at the center. As the ALMA beam size is ∼ 1. 5 × 1. 1 , we cannot distinguish between these two cases. We therefore caution that the age gradient pattern in MRG-M2129 is largely affected by the assumption for the spatial location of the 1.3 mm ALMA detection. The ALMA detection also leads to a different pattern of SFR gradients comparing to Toft et al. (2017).
With the exception of MRG-M0138 8 , we measure metallicities consistent with sub-solar to solar for all of our targets, broadly in agreement with previous measurements (Toft et al. 2017;Newman et al. 2018;Jafariyazani et al. 2020;Man et al. 2021). While we do not detect any stellar metallicity gradients, there are slight hints of a relatively metal-rich center for MRG-M0454 and MRG-M1341, but our measurements are not statistically significant.
Taken at face value, the measured flat metallicity gradient in MRG-M0138 is consistent with Jafariyazani et al. (2020). However, for the flat metallicity gradients in other targets, we note an important caveat. WFC3/G141 has a lower spectral resolution relative to the ground-based spectroscopic measurements with instruments such as Keck/MOSFIRE, it may not help with constraining potential metallicity gradients, given the estimated metallicity uncertainties here (∼ 0.3-0.5 dex). These orders of magnitude uncertainties are consistent with similar studies (e.g., Estrada-Carpenter et al. 2019). For all targets, we use the same priors for dust and metallicity for all bins, and for a half of the targets (MRG-P0918, MRG-M0150, MRG-M1341 and MRG-S1522) we get a metallicity posterior consistent with our flat prior. Therefore, our measurements cannot conclusively reject flat metallicity gradients for the REQUIEM-2D targets. Future measurements with the James Webb Space Telescope (for example with the NIR-Spec/IFU) should be suitable to answer this important question.
We find that the median stellar masses of the REQUIEM-2D galaxies (Figure 1) are ∼0.1-0.5 dex higher than previously published by Man et al. (2021) and Newman et al. (2018) for most targets, although for some targets these measurements are consistent within the 1σ uncertainties. This offset can be partially explained by the fact that we use non-parametric SFHs, unlike Man et al. (2021) and Newman et al. (2018), and this leads to higher stellar masses and older ages (e.g., Leja et al. 2019). However, the most extreme offset of ∼0.5 dex in MRG-M1341 is also related to an offset in the HST photometry. We use banana-shaped apertures (Section 2.2) that efficiently enclose fluxes even for unconventional morphologies such as that of MRG-M1341 (Figure 2), whereas Source Extractor (Bertin & Arnouts 1996) and slit-shaped apertures are used in Man et al. (2021) and Newman et al. (2018). Therefore, we measure higher fluxes for some targets in general, leading to higher stellar masses in our SED fits. Specifically, our measured magnitude for MRG-M1341 in HST WFC3/F140W band is ∼0.2 dex lower (brighter).
We find dust index < 0 for the majority of targets and A V < 1, consistent with the earlier studies of lowredshift galaxies (e.g., Salim et al. 2018). A negative dust index corresponds to a steeper attenuation curve relative to the Calzetti et al. (2000) law, leading to more UV attenuation and less IR attenuation compara-bly (Kriek & Conroy 2013). We do not detect gradients in either of these dust parameters for the majority of our sample. The only exceptions are MRG-M0138 and MRG-M2129, both having a centrally concentrated 1.3 mm ALMA detection. A V is only greater than 1 magnitude for the central region of MRG-M2129 that appears to host a dust-obscured star-forming region, producing the significant observed 1.3 mm ALMA flux. It is perhaps surprising that one data point can have such a significant impact on the fit. We explore this in more detail in a future paper (Popescu et. al., in prep).
We are able to constrain the ages and the SFHs of the REQUIEM-2D galaxies at different physical radii thanks to the morphological diversity of the targets. The variation of the magnification-corrected physical distances across the galaxies can be used to combine the age gradients for the full sample. In Figure 6, we combine the measured age gradients relative to the center for the REQUIEM-2D galaxies as a function of radial distances normalized using their half-light radii, r c (see Section 3.2, for the exact definition). We color code the age differences using the quenching time-scale, τ center , defined following Carnall et al. (2018) 9 . The age difference, δt 50 , is defined to be the difference between the age of a given bin relative to the central bin for each target in Gyr. It is evident that most measured age differences for the REQUIEM-2D galaxies are only sampling the inner core, r < r c , and are consistent with 0, i.e. a flat age gradient. However, when we consider measurements at r r c , especially for MRG-S1522, there seems to be an upward trend. We perform two separate linear fits to the REQUIEM-2D age gradients with and without MRG-S1522. Only when including MRG-S1522 do we have a statistically significant positive slope. So while the combined age gradients for the REQUIEM-2D sample seem to be consistent with younger stellar populations in the inner core, this effect is driven by a single galaxy. To robustly confirm such a conclusion, a larger sample with more independent data points at or beyond the half-light radius is required.
We find a correlation between the quenching timescale at the center, τ center , and the age gradient patterns presented in Figure 6. The REQUIEM-2D galaxies that have smaller τ center , that is their cores are forming fast, also have a younger stellar population in their inner core relative to their outskirts. Conversely, those galaxies  Carnall et al. (2018), versus the global formation time-scale defined as the ratio of the age of the universe when 50% of the total stellar mass is formed to the age of the universe at the redshift of observation. The global and central SFHs of all galaxies are plotted in the smaller panels. The portions of SFH that are used to define τcenter are highlighted with blue, black and red. The vertical dashed lines in the right panels show the global measurement of t50 (i.e., the x axis in the left panel). Quiescent galaxies that form their inner core slower (higher on the y axis, left panel) are also formed earlier in the universe (to the right on the x axis, left panel) and vice versa.
that form their inner cores slowly also do so in a uniform manner, as they have flat age gradients in r < r c .
To further investigate this correlation in Figure 6, we plot the quenching time-scale at the center versus the global formation time-scale t(t 50 )/t univ (i.e., the ratio of the age of the universe when 50% of the total stellar mass formed to the age of universe at observation). Small values for this formation timescale ratio represent galaxies that have already formed half of their stellar mass early in the universe, whereas larger values are those that form later. Interestingly, galaxies that form their cores faster (low τ center ), also seem to form later in the universe (high t(t 50 )/t univ ) and vice versa. We therefore conclude that the REQUIEM-2D galaxies that are formed earlier in the universe are also forming their inner cores slowly and uniformly, whereas the galaxies in our sample that are formed later in the universe are forming their inner cores faster and have a younger center in general. The existence of fast/slow quenching channels is generally consistent with the unresolved observational studies of quiescent galaxies (e.g., Belli et al. 2019), but our spatially-resolved study confirms it for the inner core and it supports that the fast quenching channel should lead to age gradients. The existence of age gradients in quiescent galaxies is also broadly consistent with cosmological simulations (e.g., Wellons et al. 2015;Tacchella et al. 2016).
We show that the quiescent galaxies that form later in the universe build up their cores faster and vice versa. At face value, this result is in tension with previous chemical abundance analyses that imply the opposite (e.g., Thomas et al. 2005). Our observations challenge the current paradigm where the earliest massive galaxies form most rapidly. However, the chemical abundance analyses of local early-type galaxies in Thomas et al. (2005) that support the current framework are unresolved and naturally probe long timescales owing to the epoch of observation. Over ten billion years, dynamical processes could in principle homogenize the central formation time-scale between z ∼ 2 and z ∼ 0 and effectively wash out the signatures that we capture at high redshift. Our observations could be consistent with a physical picture in which z ∼ 2 quiescent galaxies are formed via two different channels: central starbursts often triggered by gas-rich mergers at z ∼ 2 − 4 or gradual, accretion-throttled formation at early cosmic times when the Universe was denser Wellons et al. (2015). Future spatially-resolved chemical abundance studies of galaxies closer to the epoch of formation at z 2 with higher spectral resolution, potentially with the James Webb Space Telescope, are necessary to understand the origin of this tension. Our understanding of how massive galaxies formed the bulk of their stars thus may be incomplete and/or incorrect.
We finally note that the correlation between τ center and t(t global 50 ) does not trivially hold as the measured time-scales in principle represent different properties of the SFHs. τ center measures the time difference between forming 50% of the stellar mass and when the SFR drops to 10% of its average, t quench , normalized by t(t quench ) for the center, whereas the t(t global 50 ) measures the age of the universe where 50% of the global stellar mass is formed. For example, we have two galaxies in the sample, MRG-M1341 and MRG-S1522, with almost the same t(t global 50 )/t univ but very different quenching timescales at their centers. As it is clear from Figure 7, they also have distinct SFHs at their centers while having pretty similar global SFHs. However, we acknowledge that a sample of eight is not adequate to sample the full parameter space possible in this diagnostic plot. More data of similar quality would help build a more representative picture of central quenching relative to global formation timescales for the overall quiescent population at cosmic noon.

SUMMARY
In this paper, we present the novel analysis of deep 5-15 orbit grism spectroscopy from the REQUIEM-2D survey of eight strong-lensed quiescent galaxies ranging from z=1.59 to z=2.92. By combining the grism spectroscopy with HST and Spitzer imaging, we perform a spatially resolved spectrophotometric analysis, developing and publicly releasing the requiem2d software package. With this methodology, we can robustly constrain the star formation histories in the inner cores of distant galaxies in order to constrain signatures of formation and quenching. Our conclusions from this analysis are detailed as follows.
• The REQUIEM-2D galaxies are not dusty (A V < 1) and have attenuation curves that are steeper than the Calzetti et al. (2000) law in the UV and flatter in the IR. We cannot reject flat metallicity gradients in the REQUIEM-2D survey.
• Despite having only eight targets in the REQUIEM-2D survey, we already observe a diversity of formation pathways.
• Our measurements only include the 5 central bins (Section 3.1 and Figure 2), sensitive to the inner ∼1-2 kpc of the galaxy cores. For these bins, corresponding to different radial distances given the morphology and gravitational magnification (see the x-axis in Figures 4-6), we have MRG-M0138 and MRG-S0851 with flat-age gradients in their inner core, consistent with an early and almost spatially-uniform formation scenario (Akhshik et al. 2020). The flat age gradient in these galaxies are consistent with the results of Jafariyazani et al. (2020) and Setton et al. (2020). We speculate that if there was a gradient at some point during their evolution in the inner half-light radius, it has since diminished, for example as a result of mergers. MRG-M1341 has a flat age gradient, too, but its formation pathway seems different, as the SFH implies an inside-out quenching (e.g., Nelson et al. 2021;Tacchella et al. 2016;Wellons et al. 2015).
• We detect a statistically significant age gradient in MRG-S1522, with the galaxy looking younger in the inner half-light radius. This result is consistent with Gobat et al. (2017), and it may imply non-uniform quenching (Wellons et al. 2015).
The recovered SFH of MRG-S1522 clearly indicates that this galaxy quenches in the outer bins first.
• For the remaining four REQUIEM-2D galaxies, we find non-axisymmetric age gradients. Similar asymmetries are reported for both local and high redshift galaxies (e.g., Li et al. 2015;Setton et al. 2020). While this asymmetry may have significant implications for quenching, we note that by combining data points azimuthally, it washes out any statistically significant gradients. We also note that some massive rotationallysupported star-forming galaxies at z ∼ 2 show slight disturbances in their morphology. (Girard et al. 2018).
• While the combined age gradients for the REQUIEM-2D sample are consistent with younger stellar populations in the inner half-light radius, only one galaxy (MRG-S1522) samples the age beyond the half-light radius. A larger sample with more independent data points at or beyond the half-light radius are required to reach a robust conclusion.
• By investigating the age gradient patterns, the global formation time-scales, and quenching time-scales at the center, we conclude that the REQUIEM-2D galaxies that are forming their inner cores slowly and uniformly are consistent with an early formation scenario, whereas the galaxies that are forming their inner half-light radius faster are younger at center and are generally formed later in the universe.
• Our main result is in tension with previous unresolved chemical abundance analyses of local earlytype galaxies (e.g., Thomas et al. 2005), supporting the current paradigm where the earliest massive galaxies form most rapidly. While dynamical processes can alleviate this tension, future spatially-resolved chemical abundance studies closer to the epoch of formation with higher spectral resolution are necessary to understand its origin.
Future high spatial and spectral resolution observations of lensed quiescent galaxies, such as those possible with the James Webb Space Telescope, can provide further insights to the metallicity, age and chemical abundance gradients in the quiescent population to better constrain their formation pathways. For example, using NIRSpec/IFU, one should be able to constrain metallicity gradients due to NIRSpec's superior spectral resolution relative to the G141 grism spectroscopy used herein. Also, by using the IFU, one can avoid the complications of the contamination by nearby objects in grism spectroscopy. This analysis serves as a jumping off point, proposing a new diagnostic of formation pathways (i.e., Figure 7) that can be tested with larger samples, requiring high spatial resolution and enough signal-to-noise in the outskirt of galaxies.

APPENDIX
A. GRAVITATIONAL LENSING MODEL MRG-P0918 is a highly magnified, singly-imaged galaxy lensed by the massive galaxy cluster PSZ1 G295.24-21.55 at z = 0.61 (Planck Collaboration et al. 2014). We make use of the HST images to locate three multiple lensed systems (numbered 1 to 3), composed of 3 images each. We present these multiple images in Table 2, adopting the usual notation X.Y where X identifies a system of multiple images and Y is a running number identifying all images for a given system. System 1 has a known spectroscopic redshift with z = 2.146. This redshift is measured using Grizli and its automated redshift fitting pipeline leveraging grism spectroscopy (observed as a part of HST-GO-15663). In particular, the [O III] doublet at rest-frame wavelengths of 4960/5008Å is used to constrain the redshift with a reduced χ 2 of 1.04. However, systems 2 and 3 have an unknown redshift.
We model the mass distribution of PSZ1 G295.24-21.55 using the latest version (v7.1) of the public Lenstool software 10 (Jullo et al. 2007). Following previous work on cluster lens models (e.g., Richard et al. 2014), we assume the total mass distribution to be a combination of double Pseudo Isothermal Elliptical mass profiles parametrised with a velocity dispersion σ * , a core radius r core and a cut radius r cut (Suyu & Halkola 2010). The HST images show that the cluster is dominated by a central Brightest Cluster Galaxy (BCG), and we therefore model the smooth cluster-scale mass distribution using a single dPIE mass distribution. Galaxy-scale mass components are added as individual dPIE profiles, where we select galaxies based on the cluster red sequence from HST color F555W-F814W. To reduce the number of parameters following previous works (e.g., Richard et al. 2010), we assume that the mass component parameters associated with the cluster members follow a scaling relation with respect to σ * and r * cut for a L * galaxy. We fix r * cut = 45 kpc to remove degeneracy between these parameters. This does not affect the predictions from the lens model for MRG-P0918.
Lenstool uses a Monte Carlo Markov Chain (MCMC) to sample the posterior probability distribution of the model, expressed as a function of the likelihood defined in Jullo et al. (2007). In practice, we minimise the distances in the image plane: with θ (i,j) obs and θ (i,j) pred representing the observed and predicted vector positions of the multiple image j in system i, respectively. Furthermore, σ pos is a global error on the position of all multiple images, which we fix at 0. 5.  A. and ∆Dec.), dPIE shape (ellipticity and orientation), core and cut radii and velocity dispersion. The final row is the generic galaxy mass at the characteristic luminosity L * , which is scaled to match each of the cluster member galaxies. Parameters in square brackets are fixed a priori in the model.  Using all three systems from Table 2 as constraints, we obtain a best fit rms of 1.4 for this mass model. The best fit parameters for the mass distribution and the optimised redshifts of systems 2 and 3 are presented in Tables 3 and  2 respectively. Based on the best fit mass model of the cluster, we predict a magnification of µ = 7.69 ± 0.86 for MRG-P0918. We also summarize all the measured magnifications and half-light radii in Table 4.

B. JOINT SPECTROPHOTOMETRIC FIT
We discussed the joint spectrophotmetric fit in Section 3.2 and present the results for one of the targets, MRG-M1341, there ( Figure 3). We show the similar plots for the rest of the targets here in Figures 8-13.
Despite a careful selection of the grism dispersion angles, it is impossible to fully remove the contamination from the outer fainter wings of these galaxies because the targets are dispersed perpendicular to their longer axis of symmetry and are also located in the crowded fields of views (lensed by clusters). Therefore, the contamination of nearby sources in WFC3/G141 grism data remains significant for three targets: MRG-P0918, MRG-S1522 and MRG-M0454. We show the collapsed 1D spectrum of the contamination in green. The contamination is modeled and removed before fitting following Akhshik et al. (2020), but this procedure is not perfect and therefore the shape of the 1D collapsed spectra is affected to some extent for the mentioned targets (Figures 10, 11 and 13, right panels). We finally note that the full requiem2d model is constructed and fit in the native 2D space of the WFC3/G141 spectra whereas the 1D collapsed spectra are constructed only for visual purposes.  The resolved extracted 1D grism spectra. The green line in the spectra denotes the collapsed 1D contamination, affecting the shape of the spectra of the outer bins significantly. Figure 11. Left: The SED of MRG-S1522. The blue line shows 200 draws from the FSPS model. The shaded regions in the SED coincides with the unresolved Spitzer IRAC channels 1 and 2 band passes. Right: The resolved extracted 1D grism spectra The green line in the spectra denotes the collapsed 1D contamination, affecting the shape of the spectra of the outer bins significantly.