Chemical evolution of local post-starburst galaxies: Implications for the mass-metallicity relation

We use the stellar fossil record to constrain the stellar metallicity evolution and star-formation histories of the post-starburst (PSB) regions within 45 local post-starburst galaxies from the MaNGA survey. The direct measurement of the regions' stellar metallicity evolution is achieved by a new two-step metallicity model that allows for stellar metallicity to change at the peak of the starburst. We also employ a Gaussian process noise model that accounts for correlated errors introduced by the observational data reduction or inaccuracies in the models. We find that a majority of PSB regions (69% at $>1\sigma$ significance) increased in stellar metallicity during the recent starburst, with an average increase of 0.8 dex and a standard deviation of 0.4 dex. A much smaller fraction of PSBs are found to have remained constant (22%) or declined in metallicity (9%, average decrease 0.4 dex, standard deviation 0.3 dex). The pre-burst metallicities of the PSB galaxies are in good agreement with the mass-metallicity relation of local star-forming galaxies. These results are consistent with hydrodynamic simulations, which suggest that mergers between gas-rich galaxies are the primary formation mechanism of local PSBs, and rapid metal recycling during the starburst outweighs the impact of dilution by any gas inflows. The final mass-weighted metallicities of the PSB galaxies are consistent with the mass-metallicity relation of local passive galaxies. Our results suggest that rapid quenching following a merger-driven starburst is entirely consistent with the observed gap between the stellar mass-metallicity relations of local star-forming and passive galaxies.


INTRODUCTION
Since the advent of the first large-scale galaxy surveys such as the 2dF Galaxy Redshift Survey (Colless et al. 2001) and the Sloan Digital Sky Survey (York et al. 2000), galaxies have been observed to fall into a bimodal distribution in photometric colours in the local Universe (Strateva et al. 2001;Baldry et al. 2004;Bell et al. 2004;Gavazzi et al. 2010).The two sub-populations are found to exhibit different distributions across many other properties, including total stellar mass (Vulcani et al. 2013), star-formation history (SFH) (Kauffmann et al. 2003), kinematics (Graham et al. 2018), stellar metallicity (Gallazzi et al. 2005;Peng et al. 2015), radial concentration (Hogg et al. 2002), and environment (Balogh et al. 2004;Gavazzi et al. 2010).The red sequence consists of quenched, mostly dispersiondominated galaxies, whilst the blue cloud consists of star-forming, mostly rotationally-supported galaxies.The former also have higher stellar metallicity at a given stellar mass than the latter, which can be used to understand the origin of galaxy bimodality by probing the ★ E-mail: hhl1@st-andrews.ac.uk mechanisms of galaxy formation and quenching (Peng et al. 2015;Trussler et al. 2020).
Metallicity is the measurement of the mass of all elements heavier than hydrogen and helium, relative to the total mass of baryons.The vast majority of metals are produced through stellar processes, including a combination of stellar nucleosynthesis, type Ia and core collapse supernovae (for a review, see Nomoto et al. 2013 and more recently Maiolino & Mannucci 2019).These metals are then released into a galaxy's inter-stellar medium (ISM) through mass loss during the red giant phase in lower mass stars (≈ 2−8M ⊙ ) and supernovae in higher mass stars (⪆ 8M ⊙ ).In a closed box system (e.g.Tinsley 1980) the recycling of this gas into new stars leads to the next generation of stars formed having a higher stellar metallicity than the previous.However, the closed box model is an unrealistic approximation of galaxies, as interactions with the medium outside the galaxy through inflows and outflows are omitted.
Inflows from the galaxy's circum-galactic medium (CGM) bring in metal-poor gas, diluting the gas reservoir and lowering both the gas-phase and subsequently stellar metallicity.Outflows remove gas, slowing down star formation to produce fewer metals.Additionally, outflows that originate from stellar feedback might preferentially re-move high metallicity ISM gas from systems, further strengthening the role of outflows in lowering metallicity, particularly in lower mass galaxies (Chisholm et al. 2018).Therefore, the stellar metallicity of a galaxy is a result of the net sum of three processes: enrichment through stellar processes, inflows, and outflows.These processes are key components of the baryonic cycle in galaxies, which is intrinsically linked to mechanisms that cause galaxy properties to vary with time, including the quenching of star formation.
A key piece of the puzzle to understand the baryonic cycle and the evolution of galaxies is provided by higher redshift galaxy surveys such as UltraVISTA (McCracken et al. 2012).The surveys found that red quiescent galaxies grow in both number and total stellar mass since  = 4 (Ilbert et al. 2013;Muzzin et al. 2013), implying star-forming blue cloud galaxies must shut down (quench) their star formation to form quiescent red-sequence galaxies.However, the demographics of red and blue galaxy populations alone are unable to inform on the timescales of these quenching events: the steady growth in quenched galaxies could arise from the average over many individual galaxies with a wide range of different quenching timescales.
As stars form in molecular clouds, the quenching of star formation can be achieved in two ways.The first is the complete consumption of gas following the (likely gradual) termination of the supply of cold gas into the regions of star formation.The second is the sudden heating and/or disruption of the molecular clouds due to disruptive events originating from either within or outwith the galaxy.These two processes are expected to act on different timescales (e.g.Schawinski et al. 2014), which is consistent with observational findings that quenching of star formation occurs over varying timescales, ranging from > 5 Gyr to < 1 Gyr (Heavens et al. 2004;Pacifici et al. 2016;Carnall et al. 2018;Rowlands et al. 2018a).
Mechanisms proposed for the slow termination of star formation include natural depletion of gas reservoirs over time through the gradual locking up of gas into stars, the "maintenance" of hot gas reservoirs by active galactic nucleus (AGN) feedback preventing cooling of the CGM (Croton et al. 2006), morphological quenching due to the stabilising force of a central spheroid (Martig et al. 2009;Ceverino et al. 2010), shock heating of higher mass halo gas preventing cooling of gas onto galaxies (Dekel & Birnboim 2006), the inhibition of radial inflows of cold gas by the increase in angular momentum of accreted gas due to disc growth (Renzini et al. 2018;Peng & Renzini 2020) and the restriction and/or stripping of galaxy gaseous envelopes by tidal forces in clusters (Balogh et al. 2000;Boselli & Gavazzi 2006).Peng et al. (2015) and Trussler et al. (2020) have argued that slow quenching mechanisms are the main driver of intermediate and low stellar mass ( * < 10 11  ⊙ ) galaxy quenching at  < 1 due to the higher metallicity of quenched galaxies compared to star-forming galaxies in the local Universe.In this model, the slow decrease in cold gas supply leads to gradual quenching, which allows for star formation to continue with the remaining gas in the galaxy while a lack of continued inflow of low metallicity CGM gas brings reduced dilution effects.The combined effect enhances the metallicity of quenched galaxies with respect to star-forming galaxies.Trussler et al. (2020) further concluded that, although the decrease in gas supply is the main driver for quenching, a continuous secondary contribution from gas ejection through outflows is required to match the star-formation rates (SFRs) of local passive galaxies particularly at lower stellar masses.
On the other hand, studies that analysed large scale cosmological hydrodynamical simulations have found an important contribution to the build up of the red sequence from rapidly-quenched galaxies (SIMBA, ≈ 50% contribution of total stellar mass at  ∼ 1: Ro-dríguez Montero et al. 2019;Zheng et al. 2022; IllustrisTNG, ≈ 40% of galaxies over all redshifts: Walters et al. 2022).Suggested mechanisms that could lead to this rapid quenching of star formation include feedback in the form of violent ejection of gas from the central regions of a galaxy powered by AGN outflows (Feruglio et al. 2010;Cicone et al. 2014).Stellar sources such as supernovae and stellar winds could similarly provide substantial feedback, particularly in dense starburst regions (Martin 1998(Martin , 2005;;Bolatto et al. 2013;Molero et al. 2023).In clusters, infalling star-forming satellites can experience processes such as ram pressure stripping, thermal evaporation and viscous stripping, which may be powerful enough to remove cold gas directly from star-forming regions (Boselli & Gavazzi 2006).
Several approaches have been used to measure the relative importance of various quenching mechanisms observationally.This includes, but is not limited to, fitting for the SFHs of quiescent galaxies to obtain their quenching timescales (e.g.Pacifici et al. 2016), identifying star-forming galaxies with unusually low molecular gas fractions and short depletion times (e.g.Gómez-Guĳarro et al. 2022), and the aforementioned difference in mass-metallicity (MZ) relations between star-forming and quenched galaxies (Peng et al. 2015;Trussler et al. 2020).Despite the substantial work in recent years, the various approaches lead to conflicting results in the relative importance of fast and slow quenching mechanisms.
One promising avenue towards resolving this confusion in the literature is the study of post-starburst (PSB) galaxies, which have experienced a recent (< 2 Gyr), rapid drop in star formation activity (e.g.Wild et al. 2020).Studying the prevalence and properties of such objects has the potential to constrain both the contribution of rapid quenching to the growth of the red sequence, as well as the physical mechanisms responsible for such rapid quenching events (e.g.Wild et al. 2009;Rowlands et al. 2018b;Davis et al. 2019;Li et al. 2019;Zheng et al. 2020).Historically these were first identified as "E+A" or "K+A" galaxies due to their strong Balmer absorption lines and a lack of nebular emission lines (Dressler & Gunn 1983).As a result of their SFH, PSBs exhibit an abundance of A and F type stars, while the shorter-lived O and B stars are largely absent, allowing the pre-burst stellar population to not be heavily outshone (see French 2021, for a recent review).PSBs typically display compact morphologies, in both the local Universe and at higher redshifts (e.g.Almaini et al. 2017;Chen et al. 2022).Some studies have suggested that high redshift starburst galaxies such as sub-millimetre galaxies are progenitors of high-redshift PSBs (Toft et al. 2014;Wild et al. 2016Wild et al. , 2020;;Wilkinson et al. 2021) and that low redshift starburst or ultraluminous infrared galaxies (ULIRGs) are progenitors of lowredshift PSBs (Hopkins et al. 2008;Cales et al. 2011;French et al. 2015;Pawlik et al. 2018).
The initial quenching that transitions PSBs away from the starburst phase is expected to be mainly driven by stellar feedback (see e.g.Wild et al. 2009), but current-generation simulations require AGN mechanical feedback (outflows) to completely halt star formation and sustain the reduced SFR after the starburst (e.g.Zheng et al. 2020).Although PSBs account for only a minor < 1% of the galaxy population at redshift  ∼ 0 (Pawlik et al. 2016), the short visibility window of the spectral features means that a considerable fraction of all quenched galaxies could have passed through a PSB phase, particularly at higher redshift (Wild et al. 2009(Wild et al. , 2016(Wild et al. , 2020;;Whitaker et al. 2012;Belli et al. 2019;Taylor et al. 2023).Therefore, PSBs provide a key testing ground to study the effects of fast quenching mechanisms.
Measuring the gas-phase metallicity of PSBs is challenging due to the weakness of nebula emission lines and contamination with AGN, shock or diffuse interstellar excitation mechanisms, and can only be achieved in some cases (see Rowlands et al. 2018b;Boardman et al. 2024).However, we might expect substantial chemical evolution to occur during such extreme changes in star formation rate.Given the negative radial metallicity gradients of star forming galaxies (e.g.Matteucci & Francois 1989;Zaritsky et al. 1994), the inflow of gas required to drive the centralised starburst common to many PSBs might be expected to pull in lower metallicity gas from the outskirts of the galaxies, reducing metallicity.On the other hand, the very high star formation rates over a short period of time will lead to repeated recycling of gas released from evolved stars and a rapid build up in metals.Given the higher metallicity of quiescent galaxies than star-forming galaxies at given stellar mass (Gallazzi et al. 2005;Peng et al. 2015), which of these processes dominate in reality has important implications for how significantly PSB galaxies, and rapid quenching more generally, can contribute to the build-up of the quiescent population.
A systematic characterisation of the stellar metallicity evolution of PSBs has not been attempted previously to our knowledge.In this study, we aim to measure this by taking advantage of the fact that both the pre-burst and starburst stellar population are visible in PSBs' integrated light spectrum.To draw a more direct comparison with simulations that focus on the chemical evolution in the cores of starburst galaxies, we focus this study on analysing galaxies with PSB-like centres.In Section 2, we describe our data and sample selection.In Section 3, we present our method of spectral fitting of the optical continuum through stellar population synthesis models.We test the method with both "self-consistent" and simulation-based parameter recovery in Section 4, to verify we can recover the SFH and chemical history of PSBs.We then apply the method to MaNGA galaxies, present the results in Section 5, and discuss them in Section 6.Where necessary, we assume a cosmology with Ω  = 0.3, Ω Λ = 0.7 and ℎ = 0.7.All magnitudes are in the AB system (Oke & Gunn 1983).We assume a Kroupa (2001) stellar initial mass function (IMF), and take solar metallicity  ⊙ = 0.0142 (Asplund et al. 2009).We re-scale all metallicity measurements quoted from the literature to this solar metallicity for direct comparison.Throughout, we denote lookback time as  and ages of the Universe as  ′ , such that  ′ =   − where   is the age of the Universe.

DATA
MaNGA (Bundy et al. 2015) is an integral field spectrograph (IFS) survey of ≈ 10000  * > 10 9  ⊙ galaxies (11273 datacubes) in the local  < 0.2 neighbourhood, a part of the fourth-generation Sloan Digital Sky Survey (SDSS-IV, Blanton et al. 2017) that ran from 2014 to 2020.It used the Sloan Foundation Telescope at Apache Point Observatory (Gunn et al. 2006) to collect spatially-resolved spectra by using hexagonal bundles of 19 to 127 optical fibres, depending on the apparent sise of the target.The output BOSS spectrographs (Smee et al. 2013) provide high quality spectra in the wavelength range 3622 − 10354Å at a spectral resolution of  ∼ 2000 1 .We access MaNGA data through both the web interface and the python package Marvin (Cherinka et al. 2019).
For all MaNGA galaxies in the full data release DR17 (Abdurro'uf et al. 2022), we obtain redshift from the MaNGA data reduction pipeline (Law et al. 2016(Law et al. , 2021) ) and galaxy stellar mass from the NASA-Sloan Atlas (NSA_ELPETRO_MASS, a K-correction fit to elliptical Petrosian fluxes, see Blanton et al. 2011).We obtain spectral 1  = /Δ FWHM indices along with other necessary properties from the MaNGA data analysis pipeline (Westfall et al. 2019;Belfiore et al. 2019).We adjust the stellar masses from NSA for Hubble constant ℎ = 0.7.Other stellar mass estimates from SDSS-MPA/JHU2 and the Wisconsin method (Chen et al. 2012) were also considered, but provided no qualitative changes to the conclusions.We limit the sample to  < 0.06 in favour of local PSBs with good spatial resolution, leaving 7971 galaxies.
Within each MaNGA galaxy's datacube3 , spaxels4 marked with NOCOV, LOWCOV or DONOTUSE flags are removed.To identify PSB spaxels, we broadly follow the methods in Chen et al. (2019), specifically requiring the spaxels' median spectral SNR > 8 per pixel, strength of the H Balmer absorption line after accounting for emission infilling H A > 3Å (Worthey & Ottaviani 1997), equivalent width of the H nebular emission line after accounting for underlying absorption W(H) < 10Å5 , and log W(H) < 0.23 × H A − 0.46.
To select our galaxy sample, we first selected galaxies with a PSB spaxel fraction > 0.05 among all classifiable spaxels (spaxels not marked with the previous flags nor the SNR threshold that we impose).Next, we sliced the galaxies into 3 elliptical annuli with 0 < /  < 0.5, 0.5 < /  < 1 and 1 < /  < 1.5, where   is the -band elliptical-Petrosian effective radius, using the elliptical polar distance of each spaxel from the galaxy centre.Our galaxy sample is selected to have > 50% of the inner annulus spaxels classifiable, and > 50% of these spaxels to be classified as a PSB, yielding 54 candidates.This sample selection is qualitatively similar to the Chen et al. ( 2019) selection of MaNGA galaxies with "central" PSB regions.
We form a stacked PSB-only spectrum for each galaxy, only including spaxel contributions from the PSB-classified spaxels at all radial distances.To ensure spaxel quality, we remove spaxels marked with quality flags DEADFIBER or FORESTAR in MaNGA's H emission line maps.Spectra are summed unweighted, while uncertainties are summed in quadrature.The stacking of many spaxels for each galaxy allows for a very high SNR to be reached across the full MaNGA wavelength range, with the mean SNR per pixel ranging from 95 to > 1200.The SNR of the sample is listed in Table 1.After correcting for Milky Way dust reddening, we further mask major nebular emission lines, residuals of strong skylines (central flux > 5 × 10 −16 erg s −1 cm −2 Å −1 in Hanuschik 2003) and Balmer infilling.Since Pawlik et al. (2018) have previously shown that stellar population synthesis models based on the MILES stellar library (Falcón-Barroso et al. 2011) showed improved recovery of the SFH of PSBs compared to models based on other libraries, we limit the spectra to rest frame  < 7500Å to be fully within the MILES range when fitting. Figure 1 demonstrates the stacking process with two PSBs as examples.
Limiting the spectra to rest frame  < 7500Å potentially loses valuable constraining power on the older stellar population.Spectral information from longer wavelengths can form a longer wavelength baseline to minimise any age-dust-metallicity degeneracy (see Sections 5.1 and 6.1 of Conroy 2013).Hence, the flux in the portions with observed frame wavelength 7500 <  < 9500Å was summed into a single, artificial photometric data point, passed jointly with the trimmed spectra to the fitting framework (Section 3).However, no significant differences in the estimated stellar and dust properties were observed with or without the photometric data point, therefore, we limit our analysis to the trimmed spectra for the rest of this study.

OPTICAL CONTINUUM SPECTRAL FITTING
To fully utilise the fossil record stored in the high-quality MaNGA spectra, we employ the fully Bayesian spectral energy density (SED) fitting code Bagpipes (Carnall et al. 2018(Carnall et al. , 2019a)).In this section, we describe in detail our spectral fitting procedure of the optical continuum, which includes the assumed parametric SFH model for PSBs from Wild et al. (2020) (Section 3.1).Motivated by a suite of gas-rich binary major merger simulations that create PSB signatures (Zheng et al. 2020), we introduce a novel two-step metallicity model which decouples metallicity before and after the starburst, allowing for any change in stellar metallicity during the starburst to be recovered (Section 3.2).Additionally, we employ a Gaussian process (GP) correlated noise model as an additive term to the physical model's predicted spectrum to account for correlated observational uncertainties and imperfect spectral models (Section 3.3).The sampling of the posterior surface is done using the MultiNest nested sampling algorithm (Feroz & Hobson 2008) and its python interface (Buchner et al. 2014).As shown in Section 4 below, our two-step metallicity model also recovers SFH-related parameters more accurately.
Within Bagpipes, we utilise the Bruzual & Charlot (2003) stellar population synthesis models (2016 version), and assume the initial mass function from Kroupa (2001).We apply the two-component dust attenuation law from Wild et al. (2007) and da Cunha et al. (2008), with a fixed power-law exponent  = 0.7 for the interstellar medium (ISM).The dust law asserts that stars younger than 10 Myr have a steeper power-law exponent  = 1.3 and are more attenuated than older stars by a factor  (= 1/ in Wild et al. 2007;da Cunha et al. 2008), as they are assumed to be surrounded by their birth clouds.
Overall, our model has 18 parameters as listed in Table 2, 3 fixed and 15 free to be estimated.As we follow the Bayesian paradigm, prior distributions are placed on the 15 parameters.It is important to also be aware of the imposed prior probability densities on derived physical properties, for example, specific SFR (sSFR) and massweighted formation age ( M ), as they can impact the estimated galaxy properties and their uncertainties (Carnall et al. 2019b).These are shown alongside SFH draws from the SFH prior in Figure 3 of Wild et al. (2020).

The star-formation history model
The SFH traces the rate of star formation in a galaxy and all of its progenitors back in time, typically expressed in lookback time.
To model both the recent starburst and the underlying older stellar population expected in most local PSBs, we adopt the two-component parametric SFHs of Wild et al. (2020), which provides a good fit to combined spectra and photometry of  ∼ 1 PSBs: This is made up of the older, exponential decay component   and the double power-law starburst component  burst , both a function of lookback time .The lookback time when the older population began to form is denoted as  form , while the time since the peak of the starburst is denoted as  burst .The fraction  burst controls the proportion of mass formed during the starburst.The two components have the forms: (3) All times in Equations 2 and 3 are in ages of the Universe, therefore unlike  burst ,  ′ burst in the starburst component's function represents the age of the Universe at the peak of the starburst.  is the older population's exponential decay timescale, while  and  control the declining and increasing timescales of the burst respectively, with larger values corresponding to steeper slopes.The usage of the fraction  burst instead of parameterizing the stellar mass formed in the components individually allows for an easier application of a flat prior over  burst .This allows for not only SFH shapes with a strong starburst, but also rapid quenching events of the existing star formation when  burst ∼ 0.  We investigated allowing the rising starburst slope  to vary freely, with a similar prior to .However, parameter recovery tests performed using SNR = 100, at the lower end of our observations, showed that  is poorly constrained in older starbursts ( burst > 1 Gyr).Therefore, we fix  = 250, consistent with the typical value found from fits to younger starbursts.
A common alternative SED fitting method avoids assuming parametric forms for the star-formation history and instead allowing stars to form in fixed or variable bins in time (e.g.Cid Fernandes et al. 2005;Tojeiro et al. 2007;Iyer & Gawiser 2017;Johnson et al. 2021).In general these models do well with smooth SFHs, but are less well suited to galaxies which undergo a rapid change in SFR, due to the need for adaptive variability of the number of time bins.However, both Pawlik et al. (2018) and Suess et al. (2022) have successfully employed such methods, often referred to as non-parametric, to fit PSBs.Suess et al. (2022) increased the number of time bins around the time of the starburst, successfully recovering the rapid rise and fall in SFR of mock PSBs.While this can provide more flexibility in theory, in practice the need to define time bins and in some cases the inclusion of some form of regularisation to smooth between time bins makes the method more model dependent than it first seems.Additionally, no code currently exists which can implement both non parametric SFHs and a GP model to account for correlated noise, which we found crucial for our fitting (see Section 3.3).Therefore, we opt for a parametric SFH approach, noting that the GP noise component is able to account for any slight imperfections in the assumed SFH.

Two-step metallicity: insight from PSB merger simulations
During integrated light SED fitting, stellar metallicity is often assumed to be constant (e.g.Onodera et al. 2012;Gallazzi et al. 2014;Carnall et al. 2018;French et al. 2018;Wild et al. 2020;Suess et al. 2022).This is done mainly to limit the dimensionality of the problem, by sacrificing the second-order effects of chemical evolution on observations when compared to that from varying SFH, especially for broad-band photometry.This work aims to explore whether this simplification can be removed, and the chemical evolution of PSBs recovered.
To propose a simple yet representative metallicity evolution model for PSBs, we consult the suite of gas-rich binary major merger smoothed-particle hydrodynamics (SPH) simulations that create PSB signatures in Zheng et al. (2020).The simulations were performed using the SPH code SPHGal (Hu et al. 2014;Eisenreich et al. 2017), which is an updated version of the Gadget-3 code (Springel 2005).SPHGal implements sub-resolution astrophysics models from Scannapieco et al. (2005Scannapieco et al. ( , 2006)), updated by Aumer et al. (2013), and includes gas cooling rates following Wiersma et al. (2009).Chemical evolution and stellar feedback from type Ia and type II supernovae, and AGB stars are accounted for (for details, see Section 3.1 of Zheng et al. 2020).The merger progenitor galaxies were set up following Johansson et al. (2009) with modifications in the SFR adapted from Lahén et al. (2018), and initial orbital configurations following Naab & Burkert (2003).The AGN feedback models are from Choi et al. (2012) and Choi et al. (2014).The galaxy models have a baryonic particle mass of 1.4 × 10 5 M ⊙ for both gas and stars, and a gravitational softening length of 28 pc for all baryonic particles.
For our fiducial model we use the retrograde-prograde orbit merger simulation of two identical progenitor galaxies with initial gas mass fractions of  gas = 0.22 (2xSc_07), simulated with mechanical black hole feedback but no radiative feedback, because it results in strong PSB spectral features.Figure 2 plots the stellar metallicity of simulation particles against their lookback times of formation, together with the simulated SFH.When the merger-triggered starburst occurs at ∼ 550 Myr in lookback time, the newly formed stars have significantly higher stellar metallicity than previous star formation due to rapid recycling of gas to form many generations of stars, and the trend settles on more than twice the pre-burst metallicity after the starburst ends.Similar patterns are seen in other gas-rich merger simulations (Perez et al. 2011;Torrey et al. 2012).
We approximate the rapid metallicity increase with a step function and introduce a two-step metallicity model with the time of transition fixed at the peak of the starburst  burst : Both  and  burst are in lookback times.The two metallicity levels  old and  burst are independent and have identical priors, to ensure the model is equally able to fit an increase, decrease or no change in stellar metallicity during the starburst.
We experimented with several more complex metallicity evolution models: a three-step model (pre-burst, during burst, after burst); a gradual increase in metallicity prior to the burst; a two-step metallicity with scatter in the metallicity of coeval stars, following a lognormal or exponential distribution.None provided significantly improved parameter recovery, and given that we do not expect the simulations to be a perfect representation of the real Universe, we felt that any additional model complexity was not justifiable.

Treatment of correlated errors
When fitting photometric data, it is safe to assume the observational uncertainties in the bands are uncorrelated, due to individual photometric bands being observed at different time points, with different instrument set ups.However, when working with spectra consecutive pixels are not independent, due to the many processing steps involved in translating the raw spectroscopic observations into 1D spectral arrays.Following the methods in Carnall et al. (2019a, see Section 4 for a detailed discussion regarding the treatment of spectroscopic uncertainties), we introduce an additive, GP correlated noise component.As well as allowing for correlated uncertainties that stem from the standard data reduction of spectra, this component also serves to account for model-data mismatch that originates from assumptions and approximations involved at all stages of stellar population synthesis: isochrones, stellar spectral templates, SFH, chemical evolution and dust models (see Conroy 2013, for a review).
A GP can be visualised as a series of random variables along one or more continuous axes that represents some physical property.It is a general technique, that has been used to model data in various sub-fields of astronomy, including light curves of X-ray binaries and AGNs (Kelly et al. 2014), asteroseismic data of stars (Brewer & Stello 2009;Foreman-Mackey et al. 2017), exoplanet transits (Barclay et al. 2015;Foreman-Mackey et al. 2017;Chakrabarty & Sengupta 2019) and radial velocity measurements (Czekala et al. 2017), and the cosmic microwave background (Bond et al. 1999).In the case of spectral fitting, the random variables model a series of spectral flux densities along an array of wavelengths, which forms an SED.Each variable is modelled with a Gaussian distribution, such that for a dataset with  values, an N-dimensional Gaussian distribution is constructed.Before the variables are conditioned on the observed data, the prior mean of the Gaussian distributions is typically set as a vector of zeros.This is also adopted in this study.The covariance matrix describes the relationship between each one of the random variables with all other random variables.Each covariance is described by a kernel function that depends on the separation between two observations considering their physical properties.For an in-depth description of GP modelling, see Rasmussen & Williams (2006).
For the fitting of spectra, the GP's covariance matrix allows us to quantify the correlated noise between the measured flux density of any wavelength bin with all other bins.This is useful since it can account for correlated noise on different wavelength scales, where measurements at close-by wavelength bins are expected to correlate more strongly than measurements separated by longer distances.Hence, the close-to-diagonal terms of the covariance matrix will likely have a larger magnitude than off-diagonal terms.
To reduce computational time, we replace the squared exponential kernel used in Carnall et al. (2019a) with a stochasticallydriven damped simple harmonic oscillator (SHOTerm), implemented through the celerite2 python package (Foreman-Mackey et al. with parameters  = (, , ), where  scales the observational uncertainties  , on the SED fluxes,  is the amplitude of the correlated noise and  is the lengthscale of the squared exponential kernel in units of wavelength.  and   are the wavelengths at indices  and , and    is the Kronecker delta function.The first term allows for scaling of the uncorrelated input observational noise while the second term is the GP kernel function for correlated noise.In this study, we replace the second term with the celerite SHOTerm kernel function , which is a sum of exponential terms: where  = (a, c), with a and c vectors with elements   and   respectively.For a single exponential term of this form, the corresponding inverse covariance matrix is tri-diagonal, which can be computed with a small number of evaluations (Rybicki & Press 1992;Kelly et al. 2011), facilitating a reduction in computation time.
To allow for easier usage of the kernel function, we follow Foreman- Mackey et al. (2017) to take the Fourier transform of equation (6), with the power spectral density where  0 is the frequency of the undamped oscillator,  is the quality factor of the oscillator, and  0 is proportional to the power of the oscillator at  =  0 : ( 0 ) = √︁ 2/ 0  2 .To make the function more intuitive, celerite2 allows for the option to swap frequency  0 with period  via the relationship  = 2/ 0 , such that period  is proportional to a typical lengthscale at which fluxes at   and   correlate by a standard degree.celerite2 also allows for swapping  0 with the standard deviation of the GP realisations  via the relationship  = √︁  0  0 , such that this amplitude parameter is independent of the other parameters  0 and .
Aiming to emulate the behaviour of the squared exponential kernel, which works well on spectral fitting problems, we match its autocorrelation curves with those from the SHOTerm kernel.This process is described in Appendix A. This replacement of kernels allowed for a ∼ 100 fold reduction in computational time.

TESTING OF FITTING METHODS
In this section we demonstrate that the combination of the 2component SFH model with the two-step metallicity model can recover relevant galaxy parameters when presented with spectra of PSBs, with low systematic biases.We perform two types of parameter recovery tests: a "self-consistent" test is described in Section 4.1, and a smooth particle hydrodynamic (SPH) simulation test is described in Section 4.2.

Self-consistent parameter recovery
The "self-consistent" test involves generating mock PSB spectra using the functional forms for all parameters, including SFH, metallicity evolution and dust, then fitting the mock spectra with the same functional forms and spectral templates used for mock spectra generation.This setup ensures there is no model-data mismatch nor correlated errors.If the parameter recovery is successful across a large range of input values, it indicates the fitting process can recover the required parameters with minimal degeneracies when the model can perfectly describe the data.
We generate spectra using Bagpipes, with identical wavelength range and spectral resolution as our real MaNGA spectra, perturbing the generated spectra with Gaussian uncorrelated errors assuming SNR = 100 similar to the minimum SNR of our observed spectra.Typical dust and dispersion values were assumed, based on the results from our observed sample (  = 0.6,  = 3,  disp = 80 km/s).Since   (bottom), showing the input truth (black line), posterior median (solid orange line) and its 1 region (shaded), and 15 random draws from the posterior fit (dashed lines).Right: Violin plots showing posterior distributions of total stellar mass formed (log 10  * / ⊙ ), extinction in V band (  ), velocity dispersion ( disp ), burst mass fraction (  burst ), age of the burst ( burst ) and metallicity levels.The height corresponds to the distribution's density, while the central box plot marks its median (white dot), 1 region (thick brown bar) and 2 region (thin brown line).The vertical black lines indicate the input truths.In the lower right panel, the lighter and darker shaded violins correspond to the posterior of the burst and older metallicities, respectively.All parameters are estimated with accuracy within 2, and the metallicity change is recovered.
we do not inject correlated errors, there is no need to include the GP noise component during fitting.
Figure 3 shows the recovery performance of a self-consistent test with mock input parameter values similar to the SPH-simulated PSB in Figure 2, using the two-step metallicity model.The left panels demonstrate that we are able to recover the input SFH to within 1 for nearly all lookback times.In the top left panel, the apparent mismatch between the posterior median SFH (solid orange) and input SFH (solid black) before  = 1 is partly a plotting artefact.Since each posterior sample is an exponential decay function with an abrupt increase in SFR at  form , the median SFR includes a steadily decreasing fraction of models with no star formation.Hence we also plot the SFHs of 15 posterior samples in the same panel, and show the cumulative fraction of the total stellar mass formed against the age of the Universe in the bottom left panel as an alternative visualisation.In the cumulative SFH, it is easier to see that the discrepancy between the fitted median and input curves is < 1.
The right panels of Figure 3 show violin representations of the posterior distributions for seven key parameters, demonstrating they are recovered to within 2 of the input truths (solid bars).Particularly, we are able to recover the difference between the pre-burst and starburst stellar metallicities, with a high degree of confidence.

Sample of self-consistent model tests
To understand whether the offsets observed between true values posterior median estimates in Figure 3 are systematic, we repeated the recovery test 100 times with randomly drawn input values 6 from the priors in Table 2. Variations in dust and velocity dispersion are 6 The total stellar mass formed and redshift are fixed, as these do not alter the shape of the spectrum, and varying them will not provide additional insight.omitted due to computational time limitations (although see below for a comparison when these are included).We only fit mock spectra that are classified as PSBs under our selection criteria using H A and W(H) (Section 2).The mean and standard deviation of the offset (median of posterior − input truth), and fitting uncertainty for all tests are listed in Table 3. Identifying parameters where the mean offset is greater than the mean uncertainty, we find a very slight average overestimation in burst age.However, this is two orders of magnitude smaller than the range of our sample's estimated burst ages (Section 5), thus this does not impact our main results.We verify that there is no bias for finding a change in metallicity when the true metallicity remains constant, and incorrectly assigning an increase in metallicity when the true metallicity falls, and vice versa.
In addition to the test shown in Figure 3, five self-consistent parameter recovery tests with dust and velocity dispersion are performed based on randomly drawn input values, including   ,  and  disp .Comparing the five test results to the 100 above that did not have the dust component and added velocity dispersion, a ∼ 40% increase in estimation uncertainty is seen across individual SFH properties and metallicity.Despite this, with dust and velocity dispersion, the recovered values remain within 2 of the input truths (see Figure 3).
These tests show that, in the absence of model-data mismatch and correlated noise, we can recover the input parameters of the two-step metallicity model, via integrated light MaNGA-like spectra, for a wide range of PSBs with varying stellar and metallicity histories.

SPH Simulation parameter recovery
The second parameter recovery test involved generating mock PSB spectra from the stellar particles of the SPH simulations in Zheng et al. (2020), and fitting them with our assumed models to see whether we can recover the underlying galaxy properties.Unlike in the "selfconsistent" tests above, the star formation and chemical evolution history of the SPH simulations is complex, and cannot be perfectly described by the simple functional forms of our model.Additionally, stars formed coevally in the simulation can have a range of metallicities, which is not possible in our model.Thus, the mock spectra are created from galaxy properties that do not exist within the prior model space, so parameters can only be approximately recovered.Any inaccuracies and biases found during this test allows for conclusions to be drawn concerning the models' performances when tasked with real data, which will exhibit similar characteristics.
While investigating the cause of the small number of selfconsistent parameter recovery tests with large discrepancies between estimated and true values, we discovered that many occur as the mock galaxy's  form approaches the age of the Universe.The rate of change in the flux of a galaxy spectrum with time decreases with increasing stellar age, hence, errors on  form increase for the oldest galaxies.This issue is not due to the changing metallicity, as it is seen in the parameter recovery tests of both constant and two-step metallicity models.Unfortunately, all the SPH simulations in Zheng et al. (2020) were initialised with analytic templates that began their star formation at age of the Universe  = 0.5 Gyr.Therefore, to enable a better insight into the recovery of star formation and metallicity histories in PSBs, we scale the age of all simulated stellar particles down by 15%, preserving the total stellar mass formed and the shape and chronology of the SFH.The shift away from simulated SFHs that began at very low age of the Universe does not impact our results, as the typical estimated age of the Universe when star formation began for our sample ( ′ form ) is > 3 Gyr.Mock spectra are generated exactly as described in Section 4.1, with the same dust properties, velocity dispersion, SNR and uncorrelated Gaussian noise.Due to model-data mismatch caused by the simulations' parameters lying outside of the prior space, the GP noise component is included when fitting.
Figure 4 compares the results of fitting spectra constructed from the binary merger simulation shown in Figure 2 with the constant and two-step metallicity models.Due to the model no-longer existing within the model prior space, we no longer expect perfect parameter recovery.The top left panels show that the two-step metallicity model outperforms the constant model when recovering the SFH of simulated PSBs that underwent changes in stellar metallicity.The bottom left panel shows the fitted two-step metallicity model closely follows the simulation's metallicity evolution 7 . 7The sudden drop in metallicity of the simulation at 10.5 Gyr is a result As there is no direct input "truth" corresponding to many of the fitted parameters, in the right hand violin plots we instead compare the fraction of mass formed within  < 1.5 Gyr (  young ), mass-weighted mean stellar age within lookback time  < 1.5 Gyr ( M,young ), and the mass-weighted mean stellar metallicity throughout the entire simulation, as well as before and after the peak SFR of the starburst.In all cases, the two-step metallicity model substantially outperforms the constant metallicity model in recovering the underlying galaxy properties.
In the bottom right panel, we see that the fitted metallicity of the constant metallicity model is > 5 higher than the true overall mass-weighted metallicity.The over-estimation of the older stellar population's metallicity results in a redder old-star spectrum, leading to a younger fitted  form .The failure to recover the light from old stars (formed before 6 Gyr), leads to an underestimation of total stellar mass formed by > 0.2 dex.On the other hand, the underestimation of the burst population's metallicity results in a bluer young-star spectrum, leading to an overestimation of the burst age to compensate.The flexibility of the two metallicity levels allows for these problems to be mitigated.As a result, the violin plots show a significantly more accurate recovery of all parameters displayed.
To verify that the two-step metallicity model also enables good recovery of input true values when metallicity declines during the starburst, we artificially flip the metallicity of stellar particles in the simulation, to simulate a decrease in metallicity.We found the twostep model again results in a superior recovery of the SFH compared to the constant model.

Sample of simulation recovery tests
To investigate the possible bias on recovered parameters, we expand the simulation parameter recovery test to a suite of 46 tests performed on PSB spectra predicted from the simulations in Zheng et al. (2020).All tests assumed the same dust properties, velocity dispersion, SNR and perturbation of the generated spectrum as previous tests.We only use the simulation runs that resulted in an obvious starburst followed by rapid quenching i.e. prograde-prograde and retrograde-prograde simulations, more gas rich progenitors, and those with mechanical black hole feedback, but without radiative feedback.The inclusion of radiative feedback was found to be too effective at suppressing the increased star formation after the merger, leading to no/very weak PSB signatures in the resulting galaxy (see Zheng et al. 2020).Specifically, these were simulations Sa_Sc_00, Sa_Sd_00, 2xSc_00, Sc_Sd_00, 2xSd_00, Sa_Sc_07, Sa_Sd_07, 2xSc_07, Sc_Sd_07 and 2xSd_07.The initial gas mass fractions of the progenitors are 0.17, 0.22 and 0.31 for Sa, Sc and Sd, respectively.
From each of the 10 SPH simulations, we extract 10 post-burst spectra equally spaced in time from the peak of the starburst to the end of the simulation.We do this by discarding star particles formed after each time point, input the remaining particles' stellar metallicity and ages (shifted to the new time of observation), into Bagpipes which then constructs the integrated spectrum from SSPs in the same way as our models.We then measure the H A and W(H) from the integrated spectra, and check whether they pass our selection criteria (Section 2).This results in 46 simulated PSB spectra at 0.11 − 0.71 Gyr since the peak SFR of the starburst.Similar to the example shown in Figure 2, all chosen simulations exhibit rapid stellar metallicity increase during the starburst, leading to a much of the switch from analytic progenitor galaxies to SPH simulation and is not physical.2).Most panels are as per Figure 3. Panel (C) additionally presents the stellar metallicity evolution; lines and shaded regions have the same meaning as for the SFH panels.The vertical magenta line marks a lookback time of 1.5 Gyr, before and after which we calculate the fraction of mass formed within  < 1.5 Gyr (  young ) and the mass-weighted mean stellar age within  < 1.5 Gyr ( M,young ) shown in panels (G) and (H).In panel (I), the top-most black vertical line indicates the mass-weighted metallicity throughout the simulation, while the following two are the mass-weighted metallicity of stars formed before and after the peak of the starburst.The two-step metallicity model out-performs the constant model in the recovery of all parameters.
Table 4.The offsets (Δ, median of posterior -input truth) from 46 parameter recovery tests fitting SPH-derived PSB spectra with the constant and two-step metallicity models.We list here the mean and standard deviation of offsets and fitting uncertainty averaged across all 46 simulated spectra.Parameters (1), ( 7) and (8) follow the meanings in Section 3, while the other parameters are (2) mass-weighted mean stellar age before lookback time  = 1.5 Gyr, (3) mass-weighted mean stellar age after lookback time  = 1.5 Gyr, (4) fraction of mass formed within  < 1.5 Gyr, (5) mass-weighted mean stellar metallicity before the corresponding simulation's peak SFR of the starburst, and (6) mass-weighted mean stellar metallicity after the corresponding simulation's peak SFR of the starburst.For a wide range of SFHs, metallicity evolution and dust properties, the two-step metallicity model shows smaller offsets than the constant model.
The mean offset and fitted uncertainty for both constant and twostep metallicity models are presented in Table 4.The two-step metallicity model achieves less bias in SFH-related parameters (total mass formed,  M,young and  young ),   and  disp than the constant model, for a wide range of PSBs with varying SFHs and chemical evolution, ages and scatter.The two-step metallicity model returns a small mean offset in both metallicity measurements, which indicates the model is able to accurately recover the metallicity change for a broad range of simulated PSBs.
Tables 3 and 4 demonstrate that the standard deviation of the offset between true and fitted parameters is greater than their respective mean 1 uncertainty.This is expected, due to the model being an imperfect representation of the data, and illustrates why mock recovery tests with semi-realistic data are necessary.We are careful to take account of this in the interpretation of our results below, but also note that the dynamic range of the parameters of interest is sig-nificantly larger than even the standard deviation of the offset in the mock recovery tests.
We note that among the suite of simulation recovery tests there are several outliers with larger offsets in metallicity estimated with the two-step metallicity model than with the constant model.These are limited to models with a starburst peaking recently (0.2 − 0.4 Gyr in lookback time).For these to be selected as PSBs, they have a correspondingly rapid quenching time (e-folding timescale  ∼ 50 Myr).In this case, the two-step metallicity model can suffer from a degeneracy between the true solution and a slightly older starburst with higher burst mass, longer quenching timescales and a declining stellar metallicity.In most cases where the two-step model fails to recover the correct metallicity evolution, the constant model also suffers from a similar older, more massive starburst degeneracy, albeit to a less severe degree.PSBs with such recent starbursts are not found in our observed sample (Table 5).Therefore, we do not expect this to be a significant concern for our results.
To understand the spectral features that drive the better parameter recovery by the two-step metallicity model, Figure 5 compares the fitted spectra from the two metallicity models.The lower three plots in each panel of Figure 5 show the fitting residuals, the residuals smoothed by a mean boxcar of width 11 pixels and the fitted spectra's GP noise contributions.Vertical coloured regions mark the wavelengths of the Lick indices (Worthey et al. 1994;Worthey & Ottaviani 1997) and indices CNB and H+K from Brodie & Hanes (1986).The root mean square residual for both models is noted in the lower panels.
At first sight the fitted spectra from both metallicity models appear to be very well matched with the mock spectrum.However, the smoothed residual reveals differences at all wavelengths.The two-step model's smoothed residual is smaller particularly within the calcium H+K lines and the iron and magnesium metallicity indices, which contributes to the two-step model's slightly lower root mean square residual.This is consistent with the better fit to stellar metallicities obtained when using the two-step model.
The constant metallicity model's GP noise component led to a best fit model with a prominent slope (the red-end is 2.51 +0.38  −0.36 × 10 −16 erg s −1 cm −2 Å −1 higher than the blue-end), while the physical spectrum is significantly bluer than the input mock spectrum.This could arise from a combination of incorrectly estimated dust attenuation curve, incorrect metallicity or SFH properties, leading to the larger offsets for all properties in Table 4.The two-step metallicity model's GP noise component has a much smaller amplitude at all wavelengths (amplitude , small enough to be hidden behind the dashed line), and no overall slope.This indicates that the two-step model's higher degree of flexibility allows for a better approximation of the mock spectrum generated from simulation, and only minor corrections from the GP noise are required.For completeness, in Appendix B we include a SPH simulation recovery test without the use of the GP noise component to investigate the impact of the component's correction particularly for the constant metallicity model's fit.Qualitatively similar conclusions were made.
To summarise, the two-step metallicity model allows for metallicity evolution during the starburst to be traced without significant estimation biases or reduction in fitting precision due to the increased complexity, while allowing for better recovery of the SFH and its parameters of the galaxy.

RESULTS
We fit all 50 PSBs with the two-step metallicity model described in Section 3, including the new GP correlated noise model.5/50 (10%) resulted in fitted GP noise components showing obvious trends across the fitted spectral range, with amplitudes much larger than the observational uncertainty scaled by the posterior median uncorrelated noise amplitude ( in Equation 5).This can potentially indicate additional sources, such as AGN or foreground/background sources in their fitted spectra, or complex stellar/dust properties that cannot be adequately fit with the model.Due to these considerations, their results are excluded from further analysis (Plate-IFU: 7965-1902, 8080-3702, 9494-3701, 9495-3702, 9507-12704).All remaining 45 were found to have clear PSB SFHs, where rapid quenching follows an episode of strong starburst.In general, the PSB regions underwent a starburst ∼ 1 Gyr before observation, with the youngest starburst occurring ≈ 0.45 Gyr ago.The starbursts have a wide range in mass fraction (≈ 0.06−0.94).The fitted SFH properties, metallicity levels,   and reduced chi-squared 8 of the maximum likelihood posterior sample spectrum are reported in Table 5 and discussed in the following sections.All fitted SFHs and metallicity evolution are plotted in Appendix C. Figure 6 shows an example fit.In the top panel, the posterior best fit spectrum (orange) provides a visibly good fit to the observed PSB-region-stacked MaNGA spectrum (black), as is also seen by the small and normally distributed fitting residuals in the middle panel and the near unity reduced chi-squared.The posterior spectrum contains the median physical model spectrum (blue) and the additive GP noise component (bottom panel).The latter does not exhibit an overall slope and has amplitude comparable to the observational uncertainty scaled by the posterior median uncorrelated noise amplitude (dark blue dashed line), which are signs of a well behaved model.In the top panel, the physical spectrum is further separated to show dust-attenuated light contributions from the starburst (lime) and the older components (red).The split is placed at the time of minimum SFR between the old population and the starburst.This PSB has a burst mass fraction of  burst = 0.24 +0.05 −0.03 , and burst age of  burst = 0.61 +0.12 −0.03 Gyr, leading to a light contribution from the starburst that dominates marginally over the more massive older population in the red end of the optical spectrum, but more significantly at bluer wavelengths.
Since estimating properties of the older stellar population is usually more challenging, the pre-burst stellar metallicity tends to have larger ) Figure 5. Comparing the spectral fitting performance of the constant (blue) and the two-step (orange) metallicity models applied to a mock PSB spectrum generated from an SPH merger simulation (Figure 4).The spectrum is split into two rows for clarity.In each row, the main panel shows the input mock spectrum (black) and the best-fit spectra from the two metallicity models (blue and orange).The lower panels show the fitting residuals, residuals smoothed by a mean boxcar of width 11 pixels, the GP noise contribution to each metallicity models' best-fit spectral model, respectively.The orange GP noise curve has a much smaller amplitude compared to the blue curve and largely lies along the black dashed line.All y-axes have the same units.The vertical grey bars indicate masked regions where emission lines would appear in the MaNGA data.Coloured bars mark the wavelengths of common spectral indices.Root mean square residuals of the maximum likelihood posterior sample spectrum (including GP noise) for both fits are also shown.The two-step model achieves a better fit of the SPH PSB mock, particularly in the spectral regions of the calcium H+K lines and iron and magnesium metallicity indices.
uncertainty.This uncertainty further increases where PSBs are found to have high burst light fractions (  burst,L ) due to heavy outshining of the older population's light contribution.We have therefore separated the sample at  burst,L = 0.9 (calculated by integrating over the full fitted wavelength range) in Figure 7.The objects with high burst light fraction are seen to cluster around the line of constant metallicity (thick dashed line) but with large uncertainty, consistent with no change in metallicity being the a priori assumption given by the two-step metallicity model with independent and identical priors on metallicities.Excluding these 13, the proportion of PSBs that experience a net positive change in stellar metallicity at > 1 is 82% (27/33).

A recovered mass-metallicity relation, both before and after starburst
As reviewed in Section 1, the mass-metallicity (MZ) relation has been used in the literature to infer the dominant quenching timescales that build up the red-sequence.It is interesting to compare the stellar mass and metallicity properties of our sample of PSBs to the MZ relation of local star-forming and quiescent galaxies, as the nature of PSBs being found soon after rapid quenching can provide insight into the impact of these quenching processes on chemical evolution.
The top left panel of Figure 8 shows our PSB sample's stellar mass-metallicity relation before the starburst occurred.Also shown are the MZ relations from three studies of local galaxies: star-forming  1).Top: The stacked observed spectrum created by combining only spaxels classified as PSB (black), the posterior best fit spectrum and its distribution (orange line and orange shaded region), which includes contribution both from the physical model and GP noise.The physical model (blue) can be separated to show light contributions from the starburst (lime and lime shaded region) and the older components (red and red shaded region).The reduced chi-squared value of the maximum likelihood posterior sample spectrum (including GP noise),  2  , is shown.Middle: The fitting residual (orange), defined as the observed stacked spectrum (black curve)posterior best-fit spectrum (orange curve).The light blue line and the blue dashed line show the input observational uncertainty before and after scaling by the fitted noise scaling factor , respectively.An increase of around ×3 − 5 is typically required.Bottom: The fitted GP noise component and its distribution in orange, with blue curves as above.The majority of the fitted GP noise flux lies below the scaled observational uncertainty (blue dashed) and there is no obvious global trend, thus, this galaxy is recognised as a good fit.Note that y-axes have the same units, but the three panels vary in scaling.In all panels, vertical grey bands indicate regions masked due to skyline residuals, strong nebular emission lines or Balmer infilling.

Table 5.
Posterior estimated properties of 50 PSBs from the spectral fitting of stacked MaNGA spaxels.5 PSBs marked by dashes were poorly fit and are not considered in further analysis.Columns are (1) MaNGA Plate-IFU, (2) stellar mass within the stacked PSB spaxels, (3) ISM dust attenuation at 5500Å ( band), (4) the log 10 SFR within the stacked PSB spaxels averaged over the last 100 Myr, (5) time since the peak of the starburst, (6) fraction of mass formed during the starburst, (7) SFR halving timescale of the starburst, (8) stellar metallicity before the burst, (9) stellar metallicity during and after the burst, (10) change in metallicity, and (11) reduced chi-squared value of the maximum likelihood posterior sample spectrum.SDSS galaxies, mass-weighted (Panter et al. 2008); SDSS, all types, light-weighted (Gallazzi et al. 2005); star-forming and passive SDSS galaxies, light-weighted (Peng et al. 2015)  9 .Our PSBs broadly follow the known MZ relation where metallicity increases with mass, especially when we consider only the more reliable lower burst light fraction galaxies (magenta dots).This indicates that prior to the starburst, the PSB progenitors are consistent with being drawn from the 9 We note that the Peng et al. ( 2015) metallicity estimates are from Gallazzi et al. ( 2005), but split into star-forming and passive populations.underlying star-forming population, exhibiting no atypical chemical properties.
The top right panel of Figure 8 shows our PSB sample's stellar metallicity during the starburst, which shows no observable correlations with total stellar mass, suggesting starbursts disrupt the MZ relation.In the bottom left panel, we show the overall mass-weighted stellar metallicity of the PSBs, which exhibit a remarkable agreement with the light-weighted mass-metallicity relation of passive galaxies from Peng et al. (2015).It is remarkable since local PSBs, being only recently quenched quiescent galaxies, might not be expected to be representative of the galaxy population at  = 0.The difference between mass-weighted and light-weighted MZ relations should not affect our conclusions, since the difference is minor for quiescent galaxies (Trussler et al. 2020).
In the bottom right panel of Figure 8 we compare the difference between overall mass-weighted metallicity and the pre-burst metallicity, with the difference between passive and star-forming galaxies from Peng et al. (2015).In both cases the differences decrease with increasing  * .The matching trends point towards the PSB phase as a valid process that can create the large gap found between the star-forming and passive MZ relations reported in the literature.The implications are discussed in Section 6.

DISCUSSION
Our results show that most PSBs in our sample underwent an increase in stellar metallicity during the starburst phase, some very substantially.This indicates that the effect of stellar enrichment from the rapid recycling of gas during multiple rapid generations of star formation usually outweighs the combined effects of metal dilution from gas inflow and metal removal via outflow.In this section, we draw together studies from the literature to explain the metallicity changes we observe in PSB galaxies, and discuss the implications of our results for the role of PSB galaxies in galaxy evolution more generally.

On the implications of our results for the origin of PSBs
Galaxy mergers have been a popular suggested trigger for local PSBs, supported by the large fraction of faint tidal features or companion galaxies (Zabludoff et al. 1996;Chang et al. 2001), high shape asymmetry (Pawlik et al. 2016;Wilkinson et al. 2022), neural network determined post-merger classification (Wilkinson et al. 2022) and unusual features in high spatial resolution images (Sazonova et al. 2021).Recent mergers also exhibit a higher fractions of PSBs than non-merging galaxies (Ellison et al. 2022;Li et al. 2023).Observationally, a lowering of gas-phase metallicity is observed in both starbursts and pairs of interacting galaxies (Kewley et al. 2006;Rupke et al. 2008;Ellison et al. 2008), apparently in contradiction to the results in this study.We explore this further below.
Simulations demonstrate how the disruption of the gravitational potential caused by a galaxy merger leads to strong torques that can drive a rapid gas inflow to the central regions, compress it and fuel a strong starburst (Barnes & Hernquist 1991, 1996).Forward modelling of these simulations have consistently shown this to be a reliable way to create galaxies with PSB spectral features and morphologies (Bekki et al. 2005;Wild et al. 2009;Snyder et al. 2011;Davis et al. 2019;Pawlik et al. 2019;Zheng et al. 2020).Comparatively fewer studies have focused on the chemical evolution of galaxies during gas-rich merger-induced starbursts.Since the outer regions of a star-forming galaxy are typically more metal-poor (e.g.Matteucci & Francois 1989;Zaritsky et al. 1994), the inflow of substantial gas driven by the disrupted gravitational potential would lead to a net decrease in central stellar metallicity, so long as the impact of stellar enrichment from the starburst is sufficiently weak.Initial hydrodynamic simulation studies found that gas funnelling events initially decrease the central gas-phase metallicity by diluting the existing relatively high metallicity gas, smoothing the negative radial metallicity gradients common to most star-forming galaxies (Perez et al. 2006;Rupke et al. 2010;Lahén et al. 2018) in agreement with the lowered gas-phase metallicity observed in local interacting and starburst galaxies.Torrey et al. (2012) conducted ensembles of SPH simulations of major merging pairs of star-forming galaxies, finding that the change in stellar metallicity during the resulting starburst depends on the gas fractions of the progenitor galaxies: progenitors with low gas mass fractions tend to decrease stellar metallicity due to strong dilution from inflowing metal-poor gas during the merger, while progenitors with higher gas mass fractions tend to increase in stellar metallicity due to the stronger starburst and greater PSBs with a less dominating burst light fraction (posterior median  burst,L < 0.9) and a higher burst light fraction (posterior median  burst,L ≥ 0.9) are marked with magenta dots and dark green triangles, respectively.All values are plotted against the estimated stellar mass of the whole galaxy from the NSA catalogue.Stellar mass-metallicity relations from the literature are also plotted for comparison, as indicated in the legends.The dashed black lines mark the 16th and 84th percentiles from Gallazzi et al. (2005).The bottom right panel compares the difference between overall mass-weighted and pre-burst metallicity and the difference between passive and star-forming relations from the literature.The PSB pre-burst metallicity are found to agree with the Peng et al. (2015) star-forming relation, while the overall mass-weighted metallicity agrees with the Peng et al. (2015) passive relation, suggesting the PSB phase can explain the wide gap between the two literature relations.
stellar enrichment.Perez et al. (2011) found that gas-rich mergers drive a net increase in gas-phase metallicity of comparable magnitude to Torrey et al. (2012), though mainly caused by rapid increases in SFR due to fragmentation of the discs before merger-driven inflow occurs.
On the other hand, the more modern simulated major mergers from Zheng et al. (2020) used in this work produce metallicity increases with only a weak trend with gas mass fractions.Orbits that produce stronger starbursts induce stronger metallicity enhancements, and minor mergers require higher gas fractions to achieve the same strength of starburst and therefore metallicity enhancement.Results are sensitive to the AGN feedback prescription used, as this impacts the strength of the starburst.Evidently there is scope for further simulation work on the chemical evolution of galaxies during gas-rich merger induced starbursts, to understand the impact of resolution, code type, AGN and chemical evolution modelling on the final prop-erties of the descendants.The further development of semi-analytic models (e.g.Molero et al. 2023), for the specific case of low redshift elliptical formation, may also prove fruitful.
Gas-rich galaxy mergers are not the only plausible cause of PSB galaxies in the local Universe.Ram-pressure stripping in dense galaxy clusters can compress the cold interstellar gas reservoir (Lee et al. 2017), potentially leading to enhanced star formation in affected galaxies (Vulcani et al. 2020;Roberts et al. 2022), followed by rapid quenching (Poggianti et al. 2019;Werle et al. 2022).Although initially identified in clusters (Dressler & Gunn 1983), it is important to note that PSBs are predominantly located in field environments (e.g.Quintero et al. 2004;Blake et al. 2004;Goto 2005;Wild et al. 2009;Pawlik et al. 2018).The precise enhancement in the fraction of PSBs in dense clusters is still debated and may depend critically on redshift, stellar mass, and cluster and PSB selection methods (e.g.Poggianti et al. 2009;Vergani et al. 2010;von der Linden et al. 2010;Socolovsky et al. 2018;Paccagnella et al. 2019;Wilkinson et al. 2021).
Interestingly, lower stellar mass galaxies that are undergoing rampressure stripping are found to have elevated gas-phase metallicities compared to galaxies in both field and cluster environments of the same mass (Franchetto et al. 2020), which might be a result of the increased stellar enrichment from ram-pressure compression without a significant dilution effect from metal-poor gas inflow.This would produce a rise in metallicity after a starburst, at least qualitatively similar to the metallicity increase seen in the majority of our PSB sample.
To investigate further, we cross-matched our final sample of 45 well-fitted PSBs with the GEMA-VAC cluster catalogue (Argudo-Fernández et al. 2015, version DR17), finding 11/45 to be members of rich clusters (member galaxies > 100).This is around twice as high as a control sample (12.9%) of galaxies from MaNGA DR17 matched in total stellar mass and D 4000 stellar index at 1R e (MaNGA-Pipe3D, Lacerda et al. 2022;Sánchez et al. 2022; for D 4000 , see Bruzual 1983).However, we do not find any observable difference in the metallicity change and post-burst metallicity distributions of PSBs that are within rich clusters, compared to those that are not.Additionally, the PSBs and controls showed no significant difference in their distribution of local density as defined by the projected density to the fifth nearest neighbour (2-sample KS test  = 0.77).Therefore, the importance of environmental processes are unclear from our present sample.
While we find that the majority of PSBs in our sample have undergone significant increases in metallicity, a small number have experienced a metallicity drop.The most straightforward cause of a metallicity drop is strong inflow by metal-poor gas, with the inflow triggering a starburst that either produces not enough metals to counteract the effects of dilution, or the metals produced were preferentially expelled by outflows (Chisholm et al. 2018).However, in Figure 8 the galaxies that experienced a metallicity drop are also the most massive.Since outflows are more effective at suppressing chemical enrichment at lower stellar mass (Tremonti et al. 2004;Chisholm et al. 2018;Maiolino & Mannucci 2019), the expelling of metals by outflows is likely not the main explanation for the metallicity decrease.We verified that there was no systematic correlation between the change in metallicity and size of the PSB regions, either in absolute size or relative to   , as might occur if the different patterns in metallicity evolution were caused by notably different merger types, or different processes entirely.The uncertainty in the simulations means that we cannot rule out mergers as a plausible trigger for either sets of PSBs.

On the implications of our results for quenching pathways
The evolution of stellar metallicity of galaxies has the potential to provide insight into the relative importance of different proposed quenching mechanisms.Peng et al. (2015) measured the massmetallicity (MZ) relation of star-forming and passive galaxies from SDSS, finding passive galaxies to be significantly more metal rich than star-forming galaxies with the same total stellar mass, with the gap widening for lower mass galaxies (see Fig. 8, lower right).The authors conclude the large MZ gap rejects quenching mechanisms that act on short timescales as a major contributor in quenching.They argue this is because rapid quenching would prevent significant increase in stellar metallicity as galaxies become passive, predicting little to no MZ gap, which is inconsistent with their observations.Instead, they favour slower mechanisms such as the strangulation of gas inflow, which allows for quenching galaxies to increase their metallicity from the star-forming MZ to the passive MZ through stellar enrichment, given enough time.Trussler et al. (2020) largely agrees with Peng et al. (2015), but additionally proposes ejection of gas through outflows to play a minor role in quenching.
However, Figure 8 shows a good agreement between the metallicity evolution of our sample of PSBs and the star-forming-passive MZ gap.This indicates that the PSB phase, with relatively short starburst and rapid quenching that follows, is sufficient to provide the observed metallicity enhancement as galaxies move from the blue cloud onto the red sequence.Our results suggest long-term processes such as starvation are not the only viable pathways to explain the MZ gap as has previously been suggested.
This result has global implications for galaxy evolution.Studies of the evolving stellar mass function of the red sequence have found it to grow rapidly during  > 1 (e.g.Ilbert et al. 2013;Muzzin et al. 2013), but the growth slowed or stalled by  = 1 (Ilbert et al. 2013;Rowlands et al. 2018a).Therefore, a large fraction of the presentday red sequence likely quenched prior to  = 1, Evidence from both observations (e.g.Wild et al. 2009Wild et al. , 2016Wild et al. , 2020;;Whitaker et al. 2012;Belli et al. 2019;Taylor et al. 2023) and simulations (e.g.Rodríguez Montero et al. 2019;Zheng et al. 2022;Walters et al. 2022) have suggested that PSBs and rapidly-quenched galaxies could contribute significantly to the growth of the red-sequence at  > 0.5.Hence, a significant fraction of the local red sequence may have arrived there following a PSB phase.Making the leap that our local PSBs are likely undergoing similar chemical evolution to that experienced by PSBs at  > 1, we thus conclude that short lived starbursts followed by rapid quenching might be a significant contributor to the observed MZ gap in local galaxies.
Aside from slow-quenching through strangulation, other explanations for the gap in the local MZ relations between star-forming and passive galaxies have been proposed.Employing a "leaky box" analytical model of chemical evolution, Spitoni et al. (2017) modelled the two MZ relations and gas mass fractions of the star-forming galaxies to suggest that galaxies on the passive MZ relation might have assembled earlier, have shorter formation timescales and experienced weaker or less efficient outflows than star-forming galaxies of the same stellar mass.Trussler et al. (2020) noted that a faster gas accretion timescale in passive galaxies could lead to a result similar to the shut-off of gas inflow in their strangulation hypothesis, supporting their slow-quenching picture.Our results do not directly provide evidence for or against the idea that the shorter assembly time of slow-quenching galaxies could lead to enhanced stellar metallicity once they enter quiescence, when compared to star-forming galaxies of the same stellar mass.However, our results do show that the above galaxy evolution pathway is not the only valid pathway resulting in passive galaxies that are more metal rich.
More recently, Zibetti & Gallazzi (2022) drew connections to the local stellar mass surface density ( * ) to explain the observed gap in the MZ relations through investigating spatially resolved properties of galaxies in the CALIFA IFS survey.The study identified a single positive relation between stellar metallicity and  * , and concluded that the higher global MZ relation in passive galaxies resulted simply from a lack of low mass surface density regions in these galaxies.These results are qualitatively similar to those found in an independent MaNGA sample from Neumann et al. (2021), and are in agreement with relationships found between  * and stellar metallicity (Barrera-Ballesteros et al. 2017) and between  * and gas-phase metallicity (González Delgado et al. 2014).
The picture presented in Zibetti & Gallazzi (2022) is consistent with our results and previous work on PSBs.The central starburst that precedes quenching in most PSBs likely enhances the local stellar mass surface density, leading to more compact structure (e.g.Wellons et al. 2015).Indeed, PSBs are generally found to be more compact than star-forming and quiescent galaxies of the same redshift and stellar mass (Yano et al. 2016;Almaini et al. 2017;Maltby et al. 2018;Wu et al. 2018;Chen et al. 2022;Setton et al. 2022).In addition, enhanced stellar surface mass density is expected to enhance local stellar metallicity through a combination of higher gravitational potential (Barone et al. 2018(Barone et al. , 2020) ) and a larger integral of the local SFH (Zibetti & Gallazzi 2022), but the relationship's dependence on the shape of the SFH is unclear.Combining with our results showing that the PSB phase leads to a significant increase in stellar metallicity in most galaxies, a period of central starburst followed by rapid quenching can lead to the different distributions in stellar surface mass density and local stellar metallicity between younger and older regions found in Zibetti & Gallazzi (2022).

On the implications of our results for the chemical evolution of starbursts
Previous detailed theoretical work on the impact of bursty star formation on metallicity, and chemical abundance patterns more generally, has focused on local and relatively weak fluctuations in star formation rate, as might have occurred within regions of the Milky Way.
On small scales, such as within the solar neighbourhood or within dwarf galaxies such as the Magellanic clouds, periods of increased efficiency of star formation (i.e. = SFR/M gas ) will lead to an increase in metallicity due to gas recycling and stellar enrichment (e.g.Weinberg et al. 2017;Johnson & Weinberg 2020).However, on global galaxy-wide scales, evidence for substantial enhancements in star formation efficiency in starbursts is still unclear, with inferred differences potentially driven by uncertainties in which CO-to-H 2 conversion factor to assume in different galactic environments (see Kennicutt &Evans 2012 for a review, andTacconi et al. 2018 for a summary of relevant results).
We might expect the super-solar metallicity starbursts that we infer occurred in the recent past history of our PSB sample to be visible in analyses of the gas-phase mass-metallicity relation.The SFR dependence of the mass-metallicity relation for star-forming galaxies has been much debated in the past decade.Barrera-Ballesteros et al. (2017) use a sample of spatially resolved MaNGA galaxies to argue for no dependence of the MZ relation on star formation rate, and in particular there is no noticeable increase in metallicity at high sSFR in their data.However, our sample will represent < 1% of the star forming population so will not be captured in large numbers in blind surveys.Previous studies have suggested extreme LIRGs or ULIRGs in the local Universe as progenitors to local PSBs, with LIRGs and ULIRGs having similarly low number densities (Hopkins et al. 2008;Cales et al. 2011;French et al. 2015;Pawlik et al. 2018).The metallicity of such extreme starbursts is very difficult to estimate due to dust obscuration.A recent study by Chartab et al. (2022) used mid IR strong line metallicity diagnostics to show that gas in local ULIRGs is not metal deficient as previously reported using standard optical line diagnostics.The difference arises due to dust obscuring the more metal enhanced star forming regions, and places ULIRGs firmly on the local MZ relation.Further work is clearly needed to verify whether super solar gas can be identified robustly in extreme starburst regions of local galaxies.
We searched for correlations between the stellar metallicity evolution and SFH in our PSB sample, which could further elucidate any relations between starburst properties and chemical evolution.However, the potential for sample selection effects to impact observed relations made it difficult to draw firm conclusions, and we therefore leave this to future work.

Caveats and additional checks
There are a number of caveats that are worth keeping in mind with regards to our study.The most important are the fact that we fit only the spatially resolved pre-selected PSB regions of the galaxies, and selection effects imposed by our selection of PSB spaxels.
We chose to fit only the PSB regions in the galaxies in order to simplify the SFH of the integrated spectrum, improving the accuracy of our results for these regions.By selection these regions are centrally located, and therefore represent the majority of light and mass in the galaxy, but some are more inclusive of the entire galaxy than others.Systematic correlations between the spatial extent of the PSB regions and a number of fitting galaxy properties are found in our sample.In particular, galaxies with larger PSB regions tend to have lower burst mass fractions but more rapid quenching, while those with smaller PSB regions have an even spread across the whole prior range.Further work is needed to explore the relation of these regions to the wider galaxy, and whether there are correlations between chemical evolution and the PSB spatial extent, particularly in the context of local stellar mass surface density as discussed in Section 6.2.
To consider selection effects introduced by our PSB classification scheme (Section 2), we calculate the theoretical selection fraction of galaxy models within our prior space as a function of parameters of interest.This is done by first randomly drawing 10 6 mock galaxy models from the assumed SFH and metallicity priors, constructing mock spectra and measuring the spectral features H A and W(H).The fraction of mocks classified as PSB through our classification scheme is taken as the theoretical selection fraction.Although we found slight selection trends in a variety of physical properties (e.g. both older, weaker bursts and younger, slower decay bursts are less likely to be selected as PSBs), they were not consistent with causing any of the metallicity results presented here.
Minor degeneracies between SFH properties and post-burst metallicity are observed but do not qualitatively affect our results.We observed that for around half of the sample, post-burst metallicity is degenerate to a small extent with burst age ( burst ) and burst mass fraction (  burst ), where fitted models with an earlier, stronger starburst but lower post-burst metallicity, or models with a later, weaker starburst but higher post-burst metallicity have high posterior probabilities.However, the degeneracy is weak and are accounted for in the uncertainties of metallicity estimations shown in Figures 7 (contours) and 8 (error bars).Therefore, these degeneracies do not qualitative affect our conclusions.We do not see any degeneracies between any metallicity measurements with dust properties, velocity dispersion and redshift.
While alpha enhanced SSPs are available for older stellar populations (e.g.ALF, Conroy & van Dokkum 2012;Conroy et al. 2018), and have been directly compared to Bagpipes (Carnall et al. 2022), these models are not suitable for young or intermediate age (≲ 1 Gyr) stellar populations as found in our galaxies.The stellar population synthesis models we used do not include alpha-enhanced SSPs, and in this work we rely on the GP noise component to model out systematic uncertainties from spectral features related to alpha elements.To investigate this further we compared Mg and Fe absorption line indices measured from the stacked observed spectrum with those measured from our best-fit physical models (without GP noise), for all 45 well-fitted PSBs.We found that [/Fe] sensitive indices were substantially different in some cases, while [Fe/H] sensitive indices were relatively unchanged.This provides confidence that our metal-licity estimates are not impacted by the lack of alpha enhancement in our models.In addition, we observe that PSBs with the strongest measured alpha enhancements were fitted with stronger GP noise corrections in the Mg absorption features, indicating that the GP noise component is effectively removing the model-data mismatch stemming from alpha element spectral features.We conclude that the lack of alpha enhanced SSPs do not affect our measurements of stellar metallicity.
We note that our stellar metallicity estimates are limited by the metallicities of the stellar population synthesis models, which has an upper limit of  * / ⊙ ∼ 3.52.A number of our PSB's postburst metallicity estimates are close to this upper limit (see Figure 7).Although this is a potential issue for the accuracy of properties estimated for these PSBs, since the upper limit act to restrict the magnitude of metallicity increase during the starburst and subsequent rapid-quenching, our metallicity results cannot be caused by this limit.
As shown in Figure 1, we do not apply any cut on inclination during sample selection and both edge-on and face-on PSBs are included in our sample.To verify the effect of inclination on our results, we extract the 2D Sersic fit axis ratio (b/a) from the NSA catalogue (Blanton et al. 2011), and found insignificant systematic correlations with our fitted galaxy properties in all cases ( > 0.05, Spearman ranked correlation test).

SUMMARY AND CONCLUSIONS
Through selecting and stacking the PSB regions of 50 central PSB galaxies from the MaNGA integral field spectrograph survey, we fit the resulting high SNR > 100 stacked spectra with the Bayesian spectral energy density fitting code Bagpipes.Taking inspiration from a suite of binary gas-rich merger simulations that created mock PSBs, we implemented a two-step metallicity evolution model where stars formed before and during the starburst are allowed independent metallicities.We reduced the computational time to fit the high SNR spectra by a factor of 100, by replacing the original GP kernel used in Bagpipes with a stochastically-driven damped simple harmonic oscillator (SHOTerm), implemented through the celerite2 code.After careful verification of our fitting procedure through ensembles of "self-consistent" and simulation-based parameter recovery tests, we applied our model to the stacked spectra of MaNGA PSB regions to obtain 45 well-fitted results, where for the first time, the metallicity evolution of PSB galaxies with rapid SFH changes can be directly measured.Our results lead to the following main conclusions: (i) A majority (31/45, 69%) of the PSB regions of galaxies formed significantly more metal-rich stars during the starburst than before (average increase = 0.8 dex with standard deviation = 0.4 dex), while a smaller number of PSB regions formed stars of equal or lower metallicity (Figure 7).This suggests mechanisms that substantially raise stellar metallicity play important roles in the origin of PSBs: the effects of metal enrichment through stellar recycling outweigh those from dilution by gas inflow and metal removal by outflows.
(ii) This rise in metallicity during the starburst is consistent with simulations of gas rich mergers, agreeing with previous results that mergers are the leading cause of low redshift PSBs.However, we note that there is some disagreement on the impact of mergers on chemical enrichment in simulations, and more work needs to be done to corroborate the results from the Zheng et al. (2020) simulations used here.
(iii) A good agreement is found between the PSBs' pre-burst metallicity and star-forming mass-metallicity relations from the literature (Figure 8, top left).This is consistent with PSBs being drawn from the underlying population of star-forming disk galaxies as expected.
(iv) The PSBs' final mass-weighted mass-metallicity relation matches the local passive mass-metallicity relation.This suggests that the stellar metallicity evolution caused by rapid quenching following a starburst is entirely consistent with the observed gap in the stellar mass-metallicity relations between local star-forming and passive galaxies.Our results further validate the idea that rapid quenching following a starburst phase may be an important contributing pathway to the formation of the local quiescent galaxy population.
In this study we have focused on galaxies with central PSB features.Further work will be required to understand the importance of these features' spatial extent and how they compare to galaxies with other PSB spatial distributions (e.g.ringed and irregular PSBs, Chen et al. 2019).The measurement of alpha enhancement in PSBs can allow for more precise timing of their starburst and quenching.Although difficult to obtain for recently quenched systems, alpha enhancement might be detectable in PSBs with older starbursts, for instance through the methods of Conroy et al. (2018).Lastly, further simulation work on the chemical evolution of galaxies during starbursts and rapid quenching is required, to understand the effects of AGN, shocks, stellar feedback, mergers/interactions and environments on chemical evolution.
Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Figure 1 .
Figure 1.Two typical PSBs from our sample.The top represents PSBs with the vast majority of classifiable spaxels classified as PSB, while the bottom represents PSBs with only a core PSB region.The left panels show the SDSS 3-colour image with the galaxy's Plate-IFU marked on the top right corner.The MaNGA field of view is marked as the pink hexagon.The middle panels show the spaxel selection (broadly following Chen et al. 2019), displaying regions with no/faulty observations (transparent), with median spectral SNR < 8 too low to be classified (grey), classified as PSB (blue) and classified as non-PSB (red).The right panels show the stacked observed-frame spectrum of the PSB classified spaxels (black), the stacked 1 observational uncertainty (red, multiplied by 10× to make visible) and spectral ranges masked during the fitting process (grey bands), including major nebular emission lines, skyline residuals and Balmer infilling.The resulting stacked spectra have a mean SNR of 274 and 482 respectively.

Figure 2 .
Figure2.Star-formation history (bottom) and stellar metallicities ( * / ⊙ ) of the stars formed (top) in the binary gas-rich major merger simulation 2xSc07 that creates PSB signatures inZheng et al. (2020).The full SFH is shown in the inset panel, where the shaded region indicates the assumed SFH of the progenitor galaxies.Stellar metallicity increases from ∼solar levels to more than twice solar not long after the peak of the starburst at ∼ 550 Myr lookback time.

Figure 4 .
Figure 4. Simulation-based parameter recovery test, comparing the results of fitting the constant (blue) and two-step (orange) metallicity models to star particles generated from an SPH merger simulation (from Zheng et al. 2020, see Figure2).Most panels are as per Figure3.Panel (C) additionally presents the stellar metallicity evolution; lines and shaded regions have the same meaning as for the SFH panels.The vertical magenta line marks a lookback time of 1.5 Gyr, before and after which we calculate the fraction of mass formed within  < 1.5 Gyr (  young ) and the mass-weighted mean stellar age within  < 1.5 Gyr ( M,young ) shown in panels (G) and (H).In panel (I), the top-most black vertical line indicates the mass-weighted metallicity throughout the simulation, while the following two are the mass-weighted metallicity of stars formed before and after the peak of the starburst.The two-step metallicity model out-performs the constant model in the recovery of all parameters.

Figure 7 .
Figure7.Pre-burst and post-burst median posterior stellar metallicities of PSB regions in MaNGA galaxies.The right panel is a zoomed-in view of the region bound by red dashed lines in the left panel.PSBs with a less dominating burst light fraction (posterior median  burst,L < 0.9) are shown in magenta dots, otherwise in dark green triangles.The contours correspond to 1 regions (enclosing the top 39.3% of marginalised posterior probability), highlighting estimation uncertainties and degeneracies.The dashed black diagonal line marks constant stellar metallicity, while the dotted lines mark a 5× and 0.2× change in metallicity.Most PSB regions are found to increase in stellar metallicity during the starbursts.

Figure 8 .
Figure 8. Stellar mass-metallicity relations of PSB regions both before (upper left), during the starburst (upper right) and overall mass-weighted (bottom left).PSBs with a less dominating burst light fraction (posterior median  burst,L < 0.9) and a higher burst light fraction (posterior median  burst,L ≥ 0.9) are marked with magenta dots and dark green triangles, respectively.All values are plotted against the estimated stellar mass of the whole galaxy from the NSA catalogue.Stellar mass-metallicity relations from the literature are also plotted for comparison, as indicated in the legends.The dashed black lines mark the 16th and 84th percentiles fromGallazzi et al. (2005).The bottom right panel compares the difference between overall mass-weighted and pre-burst metallicity and the difference between passive and star-forming relations from the literature.The PSB pre-burst metallicity are found to agree with thePeng et al. (2015) star-forming relation, while the overall mass-weighted metallicity agrees with thePeng et al. (2015) passive relation, suggesting the PSB phase can explain the wide gap between the two literature relations.

)Figure B2 .
Figure B2.Comparing the spectral fitting performance of the constant (blue) and the two-step (orange) metallicity models applied to a mock PSB spectrum generated from an SPH merger simulation (FigureB1), without the GP noise component.All symbols follow the same meaning as in Figure5, except the GP noise panels are removed.The constant metallicity model remains a poorer fit at most metal indices.The residual at calcium H+K lines worsened.

Figure C1 .
Figure C1.Fitted SFHs (top panel of each subplot) and metallicity evolution (bottom panel of each subplot) of 45/50 successfully fitted PSBs.The solid black lines, grey and light grey shaded regions mark the posterior median, the 1 region and 90% confidence interval of the fitted SFHs or metallicity histories.Ten random samples drawn from the posterior distributions are shown using grey dashed lines.The vertical orange lines and orange shaded regions mark the posterior median and 1 region of the time of peak of starburst ( burst ).

Table 2 .
Model priors used for fitting PSB SEDs.The parameter symbols are described in Sections 3 to 3.3, or otherwise have their usual meanings.Some parameters have prior shape log 10 uniform, which indicates a flat prior in uniform space log() ∼  (log(), log( ) ).
Self-consistent parameter recovery test with the two-step metallicity model.Left: SFH (top) and fractional cumulative SFH

Table 3 .
The mean and standard deviation of offsets (Δ, median of posterior − input truth), and mean uncertainties from 100 self-consistent parameter recovery tests using the two-step metallicity model.The input values are randomly drawn from the priors given in Table2, but were then checked to ensure the resulting spectrum satisfied our PSB selection criteria.We list here the mean offset and fitting uncertainty averaged across all 100 tests.All symbols follow the definitions in Section 3.
The full table is available as supplementary online material.
(Waskom 2021)kom 2021)For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.