Molecular Gas Reservoirs in Massive Quiescent Galaxies at $\mathrm{z\sim0.7}$ Linked to Late Time Star Formation

We explore how the presence of detectable molecular gas depends on the inferred star formation histories (SFHs) in 8 massive, quiescent galaxies at $\mathrm{z\sim0.7}$. Half of the sample have clear detections of molecular gas, traced by CO(2-1). We find that the molecular gas content is unrelated to the rate of star formation decline prior to the most recent 1 Gyr, suggesting that the gas reservoirs are not leftover from their primary star formation epoch. However, the recent SFHs of CO-detected galaxies demonstrate evidence for secondary bursts of star formation in their last Gyr. The fraction of stellar mass formed in these secondary bursts ranges from $\mathrm{f_{burst}\approx0.3-6\%}$, and ended between $\mathrm{t_{end\mbox{-}burst}\approx0-330~Myr}$ ago. The CO-detected galaxies form a higher fraction of mass in the last Gyr ($\mathrm{f_{M_{1Gyr}}=2.6\pm1.8\%}$) compared to the CO-undetected galaxies ($\mathrm{f_{M_{1Gyr}}=0.2\pm0.1\%}$). The galaxies with gas reservoirs have enhanced late-time star formation, highlighting this as a contributing factor to the observed heterogeneity in the gas reservoirs in high-redshift quiescent galaxies. We find that the amount of gas and star formation driven by these secondary bursts are inconsistent with that expected from dry minor mergers, and instead are likely driven by recently-accreted gas i.e., gas-rich minor mergers. This conclusion would not have been made based on $\mathrm{SFR_{UV+IR}}$ measurements alone, highlighting the power of detailed SFH modeling in the interpretation of gas reservoirs. Larger samples are needed to understand the frequency of low-level rejuvenation among quiescent galaxies at intermediate redshifts, and to what extent this drives the diversity of molecular gas reservoirs.


INTRODUCTION
The most massive galaxies in the nearby universe (log M * > 11 M ) formed the majority of their stars rapidly, during the first few billion years of the universe's history (e.g. Thomas et al. 2010). As cosmic time passes, massive galaxies rapidly truncate or "quench" their star formation, and by z ∼ 1 the majority of massive galaxies have transitioned to a quiescent, passive phase (e.g., Bell et al. 2004;Faber et al. 2007;Muzzin et al. 2013a;Tomczak et al. 2014;Davidzon et al. 2017). Decades of study have revealed independent classes of mechanisms by which galaxies might quench, with internal physics linked to mass and external effects due to environment apparently both contributing independently to the growth of the red sequence (e.g., Peng et al. 2010). Despite quenching being a well-established phenomena, one of the most challenging problems in galaxy evolution is understanding the physical processes driving it.
To provide insights into the physical mechanisms responsible for quenching galaxies, we can study the different timescales over which they form and quench (Wu et al. 2018;Belli et al. 2019;Tacchella et al. 2022a;Akhshik et al. 2022).
For example, it has been proposed that a slow path (≈ 1Gyr) to quiescence is from the cessation of cosmic gas supply and subsequent slow consumption of the remaining gas (Feldmann & Mayer 2015). On the other hand, a rapid quenching (≈ 10 8 yr) could be due to gas-rich processes such as major mergers (e.g., Hopkins et al. 2006) or disk instabilities (Yang et al. 2008;Dekel et al. 2009;Dekel & Burkert 2014;Zolotov et al. 2015). To further complicate things, multiple mechanisms could be at play simultaneously. For example, the growth of a stellar bulge may stabilize cold gas against collapse, decreasing star formation efficiency (Martig et al. 2009), but to quench it may require additional processes such as the inability to accrete gas because of shock-heating at the virial radius of massive dark matter halos (Birnboim & Dekel 2003;Kereš et al. 2005;Dekel & Birnboim 2006;Faucher-Giguère et al. 2011).
Different evolutionary pathways may be accompanied by changes in galactic structure. Studies have found that compact quiescent galaxies are older than more extended, normal-sized quiescent galaxies Williams et al. 2017;Lee et al. 2018;Wu et al. 2018;Estrada-Carpenter et al. 2020) and that quiescent galaxies at high redshift are on average more compact than quiescent galaxies at lower redshifts (Daddi et al. 2005;Trujillo et al. 2006;Buitrago et al. 2008;Cimatti et al. 2008;van Dokkum et al. 2008;Cassata et al. 2013;Mosleh et al. 2017;Suess et al. 2019). One of the possible drivers of size evolution is that the stellar density of galaxies scales with the density of the universe at its formation epoch, such that galaxies that form later are larger (progenitor bias; Poggianti et al. 2013;Carollo et al. 2013;Lilly & Carollo 2016), another possibility is that size increases due to accreted stellar envelopes from dry minor merging (Bezanson et al. 2009;Hopkins et al. 2009;Naab et al. 2009;Trujillo et al. 2011;Hilz et al. 2012;Oser et al. 2012;Bluck et al. 2012;Poggianti et al. 2013;Ferreras et al. 2014;Damjanov et al. 2022). While it is possible local samples are biased against identifying compact galaxies (Taylor et al. 2010;Calvi et al. 2014;Szomoru et al. 2013), empirical studies appear to favor dry minor merging, which we discuss in further detail in Section 5.3.
A new avenue for providing constraints on the quenching process involves studying the molecular gas content of quiescent galaxies, because ongoing star formation requires cold gas to collapse into stars. Indeed, this is reflected in global scaling relations such as the correlation between the surface density of star formation and molecular hydrogen (Kennicutt 1998). For galaxies to quench, this cold gas must be unavailable or unusable to form new stars, which can happen because of many different mechanisms such as stabilization against collapse, feedback, or consumption. Several studies have characterized molecular gas in statistical samples of massive quiescent galaxies at z ≈ 0 (Young et al. 2011;Saintonge et al. 2011Saintonge et al. , 2012Saintonge et al. , 2017Davis et al. 2016) through ambitious large-scale surveys, investing significant telescope time to detect the faint emission from relatively small molecular gas reservoirs. They found that, after a sustained period of quiescence, the majority of galaxies have low gas fractions, < 1%, suggesting that additional gas accretion is prevented (although exceptions with large gas reservoirs exist, up to 24%, see Sansom et al. 2019). Other studies have found that recently quenched galaxies (e.g., post-starburst galaxies) in the local universe have a puzzling diversity of molecular gas, with gas fractions up to ∼ 50% (Rowlands et al. 2015;French et al. 2015French et al. , 2018bAlatalo et al. 2016;Smercina et al. 2018;Li et al. 2019;Yesuf & Ho 2020;Smercina et al. 2022;French et al. 2022), suggesting that at least some fraction of galaxies that quenched recently have done so by a different mechanism that does not destroy or deplete gas. To better understand the quenching process, it is important to study the gas content of galaxies closer to the peak epoch of quenching at z ≈ 1.
Beyond the local universe, relatively few studies have made molecular gas mass measurements. Similarly to poststarburst galaxies in the nearby universe, studies at higher redshifts reveal a heterogeneous range in molecular gas properties. At z ∼ 0.7, there is heterogeneity in the gas content of individual galaxies (Suess et al. 2017;Spilker et al. 2018;Bezanson et al. 2022), while most of the studies at z > 1 of individual galaxies find very low gas content, of a few percent or less (Sargent et al. 2015;Bezanson et al. 2019;Caliendo et al. 2021;Williams et al. 2021;Whitaker et al. 2021a). This finding could indicate that the quenching process in high-redshift galaxies is both efficient and complete (Bezanson et al. 2019;Williams et al. 2021), plausibly due to galaxies running out of gas (Whitaker et al. 2021a). However, average properties measured by stacking dust emission in quiescent galaxies at z > 1 have, in contrast, shown evidence for large gas reservoirs, with gas fractions an order of magnitude higher (Gobat et al. 2018;Magdis et al. 2021). These seemingly conflicting results may result from different methodologies, or they may more broadly reflect gas content in quiescent galaxies that are higher than implied by the small sample of single targets, indicating that the population demonstrates a large diversity of gas content. This diversity could be linked to different quenching processes in galaxies, or it may reflect external factors that could impact their later star formation histories (SFHs), such as recently accreted gas.
Therefore, understanding the link between cold molecular gas reservoirs and the detailed SFHs of quiescent galaxies would provide important insights into the quenching process. Inferring the SFHs of galaxies requires high-quality spectra and photometry to probe stellar populations for the subtle diagnostic features that imprint the detailed timescales of star formation, so most studies have been limited to the local universe. Recently, large spectroscopic surveys have produced high S/N spectroscopy of galaxies outside of the local universe, enabling the properties of galaxies out to z ∼ 1 to be characterized to similar quality of SDSS at z ∼ 0. The Large Early Galaxy Astrophysics Census (LEGA-C) Survey (van der Wel et al. 2016;Straatman et al. 2018;van der Wel et al. 2021) has produced 20-hour depth spectroscopy of galaxies from 0.6 < z < 1 to observed λ 8800Å, providing exquisite constraints on rest-frame optical diagnostics that trace ages, SFHs, and velocity dispersions. In addition, a small sample of LEGA-C quiescent galaxies has been observed with the Atacama Large Millimeter/submillimeter Array (ALMA), providing gas mass measurements based on observations of CO(2-1) emission (Spilker et al. 2018), which comprise the sample for this paper.
In this paper, we infer the detailed SFHs of the Spilker et al. (2018) sample while assuming a highly flexible functional form for the SFHs. We are thus able to reconstruct complex and more realistic evolutionary pathways, such as bursts of star formation that cause sharp transitions in the star formation rate, SFR(t), that cannot be captured by traditional parametric models of star formation history. Having both molecular gas measurements and the ability to constrain SFHs, both of which are challenging to measure, opens the door to a better understanding of the quenching process. The only other ALMA survey at comparable redshifts and with similarly excellent quality spectra is SQUIG-GLE Bezanson et al. 2022), where galaxies were selected to be in a post-starburst phase and therefore are a younger population than our sample. Our work thus provides constraints on the mechanisms that cause more early-forming quiescent galaxies to quench their star formation, and it gives new insights into the heterogeneity of their molecular gas mass measurements, complementing studies that target more-recently formed post-starburst galaxies.
In Section 2, we present our sample. In Section 3, we present our spectral energy distribution (SED) fitting methods and models. In Section 4, we discuss our results. In Section 5, we discuss possible causes for our results. We assume a flat Λ cold dark matter (ΛCDM) cosmology with WMAP-9 parameters, Ω m = 0.287, H 0 = 69.3 km s −1 Mpc −1 (Hinshaw et al. 2013), and a Chabrier initial mass function (IMF) (Chabrier 2003 Straatman et al. 2018;van der Wel et al. 2021)) is a 128-night ESO public spectroscopic survey with VIMOS on the Very Large Telescope. The survey includes ∼ 3500 K s -band selected galaxies at redshifts 0.6 < z < 1.0, stellar masses M * > 10 10 M , and observed wavelength range 6300Å ≤ λ ≤ 8800Å in the widearea COSMOS field ). The 20 hour long integrations produce continuum spectra with S/N ∼ 20Å −1 and high-resolution (R≈ 3500; Straatman et al. 2018). We use the spectra from the third data release (DR3: van der Wel et al. 2021).
In addition, the LEGA-C galaxies have excellent ancillary data from the UltraVISTA catalog (Muzzin et al. 2013b) that we utilize in our fitting. The LEGA-C galaxies were selected from this catalog, which is a collection of photometric data across 30 passbands from 0.2 to 24µm. It includes UV imaging from the GALEX satellite (Martin et al. 2005), optical imaging from the Canada-France-Hawaii Telescope (CFHT) and Subaru telescope (Taniguchi et al. 2007;Capak et al. 2007), near-infrared data from the VISTA telescope (McCracken et al. 2012), and mid-infrared data from Spitzer (Sanders et al. 2007;Frayer et al. 2009). The IRAC and MIPS photometry from Spitzer was deblended using a welltested code for many K s selected catalogs (e.g. Labbé et al. 2005;Wuyts et al. 2007;Marchesini et al. 2009;Williams et al. 2009;Whitaker et al. 2011;Marchesini et al. 2012;Labbé et al. 2010Labbé et al. , 2013. See more details from Section 3.5 of Muzzin et al. (2013b).

ALMA
We use ALMA measurements from Spilker et al. (2018) who observed the CO(2-1) transition, which is a tracer of the molecular ISM. The 8 massive, quiescent galaxies were selected to be log M * > 10.8 M and with SFR UV+IR 3 − 10× below the star-forming main sequence as defined by Whitaker et al. (2014). Four of the eight targets have clear detections of CO(2-1) emission. Spilker et al. (2018) use standard assumptions to convert the CO(2-1) luminosities to molecular gas masses. They assume the luminosity ratio between CO(2-1) and CO(1-0) transitions, r 21 = 0.8, which is within the plausible expected range, r 21 = 0.7 − 1.0, based on observations of the Milky Way, nearby quiescent and star-forming galaxies, and highredshift star-forming galaxies (Fixsen et al. 1999;Combes et al. 2007;Dannerbauer et al. 2009;Young et al. 2011;Spilker et al. 2014;Saintonge et al. 2017). In addition, they assume a "Milky Way-like" conversion factor α CO = 4.4M (K km s −1 pc 2 ) −1 (Solomon et al. 1987;Bolatto et al. 2013;Sandstrom et al. 2013). For the 4 CO-undetected galaxies, conservative 3σ upper limits are estimated. The noise in 800 km s −1 wide channels is used for this estimation, which would fully include all of the line emission from detected sources. The integrated line flux upper limits are proportional to √ ∆v, where ∆v is the velocity interval over which the spectrum is integrated. The range of NOTE-zspec and σ are the spectroscopic redshift and velocity dispersion from LEGA-C DR3 molecular gas masses measured for the CO-detected galaxies is 9.82 < log (M H2 /M ) < 10.16 and the limits for the CO-undetected galaxies are 9.72 < log (M H2 /M ) < 9.91. The redshift range for these galaxies is 0.60 < z < 0.73. We note that the range of half-light sizes for our sample of galaxies is 0.4-1.6 , therefore the full extent of the stellar distribution is larger than this. The VIMOS slits widths are 1 and at ∼ 2mm wavelength of the CO(2-1), ALMA has a primary beam Half Power Beam Width of ∼38.8 . Therefore, ALMA is able to probe the entire galaxy and surrounding regions while the optical spectra probe only the inner regions (<1 ) of the galaxies.

SED FITTING
We use the SED fitting code Prospector (v1.1.0), which uses a Bayesian inference framework (Leja et al. 2017;Johnson et al. 2021b) to infer SFHs. The posterior distributions are sampled using the dynamic nested sampling code dynesty (v1.2.3) (Speagle 2020). Dynamic nested sampling is preferred over the more popular Markov Chain Monte Carlo (MCMC) methods because it has precise and flexible stopping criteria and improved evaluation of multimodal solutions .
The spectra were flux-calibrated by matching the shape of the photometric SED (Straatman et al. 2018). This is beneficial for the simultaneous fitting of both spectra and photometry, as it helps correct any wavelength-dependent normalization offset between the spectra and photometry, from imperfect flux calibration of the spectra. We fit 30 passbands from the photometry, including GALEX (NUV and FUV), CFHT (MegaPrime u*), Subaru Suprime-Cam (IA427, IA464, IA484, IA505, IA527, IA574, IA624, IA679,  IA709, IA738, IA767, IA827, B, V, g+, r+, i+, z+), UltraV-ISTA (Y, J, H, Ks), Spitzer IRAC (ch1, ch2, ch3, and ch4) and MIPS24, see Section 2.1 for more information. We enforce a 5% error floor for all bands of photometry, to allow for systematic errors in the physical models for stellar, gas, and dust emission. Some galaxies in our sample do not have one or both of the UV data from GALEX, in which case only the remaining 28-29 passbands were fit.

Stellar Populations
Prospector uses the Flexible Stellar Population Synthesis (FSPS, v3.2) stellar population code (Conroy et al. 2009) via python-FSPS (v0.4.1) (Foreman-Mackey et al. 2014). We employ the MIST stellar evolutionary tracks and isochrones (Choi et al. 2016;Dotter 2016) which utilizes the MESA stellar evolution package (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018. We use MILES for the stellar spectral library (Vazdekis et al. 2015;Falcón-Barroso et al. 2011) Figure 1. Observed data in purple and best fit Prospector model in black for the CO-detected galaxies. A mask was applied to the spectra before fitting to only cover the approximate spectral response limits and to exclude wavelengths below the MILES spectral library limit. The left column shows the spectra while the right column shows the photometry. The smaller panels at the bottom show χ, defined as (data-model)/σ. The middle columns are HST/ACS F814W images of each galaxy, shown as 4 × 4 cutouts. Galaxies are assigned numbers 1-8 and will be labeled as such in all future plots for ease of comparison. The data are well-fit by the model.

Dust Attenuation
For dust attenuation, we tested different models. We fit with a Kriek & Conroy (2013) flexible attenuation curve with the UV bump tied to the slope of the curve. We use the Charlot & Fall (2000) two-component dust model, which assumes a diffuse dust screen across all stellar light (τ dust,2 ) and additional attenuation of stellar light for stars that have formed in the past 10 Myr (τ dust,1 ). The free parameters are the diffuse dust attenuation index (n),τ dust,1 , andτ dust,2 . We use a joint prior on the ratio of the two dust components 0 <τ dust,1 /τ dust,2 < 2.0 as in . We compared the results obtained using this method to those determined with a Calzetti et al. (2000) attenuation curve and found that all of their inferred masses were consistent within uncertainties. In addition, for 6 of the 8 galaxies, their SFHs were consistent within uncertainties. However, for galaxy 1, five of the agebins were not consistent within uncertainties and for galaxy 2, the three most recent agebins were not consistent within uncertainties. The ratio of the Bayesian evidence can be used to determine which model is better-fit by the photometry and spectra; the evidence increases with more accurate models but is penalized for model complexity. For most galaxies, the Bayesian evidence does not strongly fa-  vor either model. However, for galaxies 1 and 2 the evidence strongly supports the more complex dust model. Because the inferred SFHs were not significantly different from each other with the different dust models for all galaxies except galaxy 1 and 2, we use the evidence-supported Kriek & Conroy (2013) flexible attenuation curve with the Charlot & Fall (2000) two-component dust model as our fiducial model for dust attenuation. This is in line with other studies supporting steeper than Calzetti attenuation curves in quiescent galaxies (e.g. Salim et al. 2018;Nagaraj et al. 2022).

Dust Emission
Our dust emission model assumes energy balance, where all starlight is attenuated by dust and re-emitted in the IR (da Cunha et al. 2008). We use the Draine & Li (2007) dust emission templates, which control the shape of the IR SED with three free parameters. U min and γ e control the shape and lo-cation of the thermal dust emission bump and q P AH is the fraction of total dust mass that is in PAHs. These templates are based on the silicate-graphite-PAH model of interstellar dust (Mathis et al. 1977;Draine & Lee 1984). See Appendix C of Leja et al. (2017) for more information.

AGN
Our AGN model adopts AGN templates from Nenkova et al. (2008a) and Nenkova et al. (2008b), which are included in FSPS. This model has two free paremeters: f AGN , which is the AGN luminosity as a fraction of the galaxy bolometric luminosity, and τ AGN , which is the optical depth of the AGN dust torus. A log-uniform prior describes the observed power-law distribution of black hole accretion rates Georgakakis et al. 2017;Caplar et al. 2018), therefore we adopt a log-uniform prior for f AGN . A loguniform prior is adopted for τ AGN , as the SED response to  Figure 2. A comparison between the parametric and fidudical nonparametric models for individual galaxies. Note that the y-axes are not on the same scale between plots. The parametric SFHs for each galaxy are shown in black. The nonparametric SFHs are shown in purple for the CO-detected galaxies and in red for the CO-undetected galaxies, which are the fiducial models for this study. The shaded regions show the 16th and 84th percentile range of the posterior probability. The time of the burst for the parametric model is shown as a solid black vertical line, with the 16th and 84th percentiles shown as dotted black vertical lines. The fraction of mass formed during the burst for both the parametric and nonparametric models are listed in Table 3. The CO-detected galaxies show evidence for enhanced late time star-formation (0-1 Gyr in lookback time) and secondary bursts of star formation compared to the CO-undetected galaxies.  . The nonparametric SFHs shown in comparison, with the bursty prior in black and the continuity prior in purple for CO-detected galaxies and in red for CO-undetected galaxies. Note that the y-axes are not on the same scale between plots.  logarithmic changes in τ AGN is approximately linear. For more details, see . We use nebular marginalization for the AGN emission lines because they are powered by sources not included in our nebular emission model; see next section for more details.

Nebular Emission
The nebular emission is based on CLOUDY model grids , and includes a nebular continuum component. The ionizing continuum from the stellar populations is modeled to be fully absorbed by the gas and then emitted as both line and continuum emission. Both the gasphase metallicity and gas ionization are set as free parameters, and the gas-phase metallicity is decoupled from the stellar metallicity. The line emission from quiescent galaxies may be caused by evolved stars, gas shocks, or AGN (Kewley et al. 2006), while the assumption in FSPS is that emission lines are caused by HII regions. Evolved stars, more specifically extreme horizontal branch or post-asymptotic giant branch stars, have been shown to produce ionizing photons in theoretical models (Greggio & Renzini 1990;Binette et al. 1994;Cid Fernandes et al. 2011;Belfiore et al. 2016;Byler et al. 2019). Maseda et al. (2021) found that 59% of UVJ-quiescent LEGA-C galaxies have [OII] emission (EW> 1.5Å) and that star formation is unlikely to provide the bulk of the ionizing photons. For this reason, we use the nebular emission-line marginalization option in Prospector. This allows for some flexibility in the emission-line modeling; each emission line is modeled with a Gaussian fit with variable width and amplitude instead of coupling the ob-served emission line luminosities to stellar physics, see Johnson et al. (2021b). In other words, this decouples the emission lines from the SFH.

Noise and Calibration Models
Simultaneously fitting high-resolution spectra and photometry can be challenging owing to differing instruments, data processing, and calibration procedures. Therefore, we include both noise and calibration models in our simultaneous SED fitting. Three parameters are used for noise and calibration including f outlier,spec , spec norm, and spec jitter. f outlier,spec is an outlier model that marginalizes over significant deviations from the model, see Hogg et al. (2010) for more details. For our sample, the median fraction of pixels which are determined to be outliers, or "bad" spectral data points, from the fits is 0.001 (0.1% per galaxy).
To look for spectroscopic uncertainties that are not welldetermined, the parameter spec jitter is also used. This term essentially inflates the errors such that the reduced χ 2 ≈ 1, where when spec jitter=1, there is no inflation. The median spec jitter value for our sample is 1.5, showing some minor inflation.
Because we are fitting both spectra and photometry, the overall continuum shape is set using photometry with the continuum of the spectra adjusted to match the shape from the photometry. This is done with the parameter spec norm, where the model spectrum is multiplied by a 10th-order maximum-likelihood polynomial using the ratio of the observed and model spectrum before computing the likelihood. The value for spec norm is then the ratio between the ob-served spectrum and model spectrum. The median value for spec norm for our sample is 1.03. This value is very close to unity, due to the spectra being flux-calibrated to the photometry before fitting. See Table 1 for more details, including priors used, and Johnson et al. (2021b) for a more in-depth discussion on the noise and calibration models.

Star Formation Histories
We infer the SFHs with three different models included in Prospector to test whether we find the same results while making different assumptions. We include a parametric model and two nonparametric models with differing priors.
For the parametric model, we experiment with the commonly used delayed-τ model, with functional form SF R(t) = A t exp(−t/τ ), and we include an additional instantaneous random burst (delta-function) of star formation. This SFH model has the following free parameters: age, τ (timescale of exponentially declining star formation), t burst (time of instantaneous burst), and f burst (fraction of mass formed in instantaneous burst).
The ratio of the Bayesian evidence (often referred to as Bayes factor) quantifies the relative evidence between two models. Jeffreys (1998) developed a classification scheme based on Bayes factor to interpret the decisiveness of the evidence with the following categories: decisive, very strong, strong, substantial, and barely worth mentioning (weak). Using this classification scheme, we find that the nonparametric model is preferred over the parametric model for six of the galaxies in our sample. The nonparametric model shows decisive evidence for galaxy 1, very strong evidence for galaxy 4, strong evidence for galaxy 6, substantial evidence for galaxies 3 and 5, and weak evidence for galaxy 8. The parametric model shows strong evidence for galaxy 2 and very strong evidence for galaxy 7. The nonparametric model is preferred for most galaxies, and where evidence exists for the parametric model, our results would not change, because the parametric model also recovers the late time burst for galaxy 2 and the shapes of both models agree for galaxy 7.
Nonparametric models assume a highly flexible functional form for SFR(t). Leja et al. (2019a) found that the priors chosen for the nonparametric models are the primary determinants for the size of the posteriors (and therefore the uncertainties). We therefore fit with two different priors for the nonparametric SFHs. We assume the continuity prior using the student's t-distribution with σ = 0.3 and ν = 2. This prior fits for smooth transitions in SFR(t) between adjacent time bins and is a flexible model, able to fit both starforming and maximally old stellar populations, see Leja et al. (2019a). We also assume a bursty continuity prior (Tacchella et al. 2022b), using the student's t-distribution with σ = 1.0 and ν = 2, which leads to a more variable SFH by allowing sharper transitions in SFR(t) between adjacent time bins. Both nonparametric SFH models have 10 time bins specified in lookback time: the first two bins are spaced at 0-30 Myr and 30-100 Myr, while the last bin is fixed at 0.15 · t obs before t obs , where t obs is the age of the universe at the galaxy's observed spectroscopic redshift. The remaining seven bins are spaced evenly in logarithmic time between 100 Myr and 0.85t obs .
The nonparametric models allow us to fit for complex SFH trends, such as bursts and other sharp transitions in SFR(t), that would not be allowed by commonly used parametric forms such as the delayed tau-model. Parametric functions, such as the exponentially declining τ -models, intrinsically couple the older SFH shape to the recent SFH (Simha et al. 2014), potentially causing an unphysical link between unrelated events in the SFH. Johnson et al. (2021b) used mock data based on several galaxies from the large Illustris cosmological hydrodynamic simulation of galaxy formation (Vogelsberger et al. 2014;Diemer et al. 2017). They found that the flexibility of a nonparametric SFH model more accurately recovered both the older and more recent SFH shapes of simulated SFHs from Illustris, in comparison to parametric models, which only recovered the recent SFH shape. Similarly, Lower et al. (2020) use model galaxies from the SIMBA cosmological galaxy formation simulations (Davé et al. 2019) run with GIZMO's meshless finite mass hydrodynamics (Hopkins 2015(Hopkins , 2017. They also find that the nonparametric SFHs were able to match the simulated SFHs from SIMBA across time significantly better than the parametric models. For more discussion comparing parametric models to nonparametric models see

Inferred SFHs
In this section, we present the inferred SFHs for our sample of eight massive, quiescent galaxies. In Figure 1 we compare the spectrophotometric data to the fitted Prospector model predictions, finding excellent agreement. We inspect the residuals for both the spectra and photometry and find that χ, defined as (data-model)/σ, is clustered around 0 for both.
In figure 2 we show the SFHs of individual galaxies inferred with both the parametric and fiducial nonparametric model. The top four panels show the SFHs of the COdetected galaxies, while the bottom four panels show the COundetected galaxies. There is generally fair agreement in the overall shape for galaxies 5, 7, and 8. However, for the remaining galaxies, the parametric SFH does not agree well with the fiducial model. The fiducial models in these galaxies have sharp transitions in SFR(t) which tend to be poorly reproduced using simple parametric models. Figure 3 shows the two nonparametric models with differing priors compared with one another. The shaded regions show the 16th and 84th percentile range of the posterior probability. The results using the continuity prior are not significantly different from those using the bursty continuity prior. In addition to being consistent with one another, Figure 3 shows that the overall qualitative shapes of the SFHs are the same for both priors. We adopt the nonparametric star formation history with the continuity prior as our fiducial model.

Early time SFH
To explore differences between the formation pathways of our sample of galaxies, we derive metrics from their fiducial nonparametric SFH models. In particular, we calculate the mass-weighted ages, formation redshifts, and star-formation timescales. The stellar mass-weighted ages (t MW ), in units of lookback time, are defined as: where t obs is the age of the Universe at the observed redshift. The formation redshifts, z form , were calculated by converting the mass-weighted ages from units of lookback time into units of cosmic time and finding the corresponding redshift. The star-formation timescale, τ SF , is measured as the time between when the galaxy formed 20% and 80% of its total formed mass, in Gyrs (Tacchella et al. 2022a). Figure 4 shows the measurements for the formation redshifts and star-formation timescales versus molecular gas mass. There is no significant difference in the formation redshifts or the star-formation timescales between the COundetected galaxies and the CO-detected galaxies, suggest-ing that the presence of molecular gas does not depend on either metric.
Recent studies have found evidence for both fast and slow paths to quiescence (Carnall et al. 2018;Wu et al. 2018;Belli et al. 2019;Tacchella et al. 2022a;Akhshik et al. 2022), suggesting multiple quenching pathways with different associated timescales. This leads to the idea that higher levels of molecular gas in quiescent galaxies could be leftover after a slower quenching timescale. To test if this was the case for our sample, we measured the timescale of decline for the early SFH, which corresponds to lookback times 1 Gyr. We fit an exponentially declining τ -model, with functional form SFR(t)=A exp(−t/τ ), to the center of the first 4 agebins ( 1 Gyr in lookback time) in the early nonparametric SFH. A high (low) τ value would indicate a slow (fast) quenching timescale. Six of the eight galaxies have values in the range τ = 1.0 − 3.0 Gyr, while galaxies 4 and 7 have higher values of 5.5 Gyr and 7.8 Gyr respectively, with typical average uncertainties of 0.7 Gyr. There is no correlation between gas mass and τ , as shown in Figure 5 and this metric shows no significant difference between the CO-detected and CO-undetected galaxies. Therefore, our data do not show evidence that the gas content is leftover after differing quenching timescales in the early SFHs of the galaxies.

Comparison of SFR indicators
Because the molecular gas content is unrelated to the rate of star formation decline in the early SFHs for our sample of galaxies, we look for differences in their late-time star formation activity. Establishing the level of ongoing star formation in quiescent galaxies is difficult because SFR measurements in quiescent galaxies are known to be uncertain (e.g. Leja et al. 2019b;Belli et al. 2021). This is particularly true of infrared based SFR, since mid-and far-infrared flux does not always trace ongoing star formation in older galaxies, but in-stead can come from sources such as AGN, asymptotic giant branch (AGB) stars, and dust heated by older stars (Salim et al. 2009;Fumagalli et al. 2014;Utomo et al. 2014;Hayward et al. 2014). This makes the most common gold standard SFR indicator at higher redshifts, based on UV+IR measurements, unreliable for quiescent galaxies. For this reason, SFR UV+IR values in particular should be considered an upper limit for quiescent galaxies.
To confirm that our late-time SFH (0-1 Gyr in lookback time) measurements from the full SED-fitting are robust, we present independent SFR indicators based on data that are not used in our SED-fitting. This includes SFRs based on flux from UV+IR, [OII], and radio measurements. SFR UV+IR is based on a weighted sum of UV and IR (24µm) fluxes, published in Spilker et al. (2018). To measure SFR [OII] , we first correct the [OII] flux for dust attenuation using the results from our Prospector fits. We use the dust-corrected fluxes to derive SFR [OII] with the conversion factor in Kewley et al. (2004), assuming that the attenuation of gas is the same as the stars. Table 2 shows the estimated SFRs from the dust-corrected [OII] flux, assuming that the emission is entirely from star formation (Kewley et al. 2004). This [OII] emission is weak (for galaxies at these redshifts), corresponding to average < 5M yr −1 .
To calculate SFR radio , we use the 3 GHz flux measurements from Smolčić et al. (2017) and K-correct to 1.4 GHz assuming the convention S ν ∝ ν α with α = −0.7 (Condon 1992). We then use the luminosity to derive SFR radio with the conversion factor from Murphy et al. (2011).
We compute SFR radio for all the galaxies in the sample. We find that galaxies 2 and 6 have high SFR radio values, in agreement with the results of Barišić et al. (2017) indicating that the radio flux is likely dominated by emission from a radio-loud AGN. This classification was determined using the radio luminosity limit from Best et al. (2005). We find that they do not show any strong emission lines that would be characteristic of AGN in our spectra. We check if the AGN are obscured by dust in the center, but the IRAC colors for these two galaxies are not near the region populated by AGNs (Kirkpatrick et al. 2013), and are instead near the region for quiescent galaxies, meaning they are not highly dust-obscured AGN. The results from our Prospector fits infer low AGN luminosity fractions with values of f AGN = 0.02 +0.12 −0.02 % and f AGN = 0.03 +0.29 −0.03 % for galaxies 2 and 6, respectively. We confirm visually that the optical detections are not offset from their radio detections, meaning they are not contaminated by a nearby companion.
We calculate the radio spectral index for these two galaxies identified as radio loud AGN by (Barišić et al. 2017) using their 1.4 GHz flux measurements from Schinnerer et al. (2010) and 3 GHz flux measurements from Smolčić et al. (2017). The spectral index for galaxy 2 is −0.1 ± 0.1 and for galaxy 6 is −0.4 ± 0.2, both of which are consistent with a flat-spectrum AGN (de Gasperin et al. 2018). The high SFR radio values for these galaxies are likely due to added flux from an AGN and the radio flux is not only tracing star formation. We conclude that our measured current SFR (averaged over the past 100 Myr, also presented in Table 2) we derive from the SED fitting is consistent with the radio measurements, as the high SFR radio is likely contaminated by flux from the AGN.
Five of the eight galaxies do not have significant 1.4 or 3 GHz emission. Therefore we present the 3σ upper limit to SFR radio , 14 M yr −1 , based on the sensitivity of radio imaging of the COSMOS field from Smolčić et al. (2017) (assuming the canonical spectral index for radio emission from star formation). The SFR measurements from the SED fitting are well below the upper limits from the radio imaging, confirming their low SFRs measured by our SED fitting. The remaining galaxy 7 has an SFR radio value that is higher than the SFR determined from the SED fitting. Table 2 shows a comparison of the SFR measurements using all four different techniques for the galaxies in our sample, and we find that there are systematic offsets between them (see also Belli et al. 2021). As expected for the reasons discussed above, the SFR UV+IR values are overestimated, by about a factor of 4, compared to our SED-fitting results (SFR 100Myr ). For galaxies 4 and 5, the SFR UV+IR values are 6 times larger than the SFR [OII] values. Galaxy 1 is an exception, as it has a larger SFR measured from the SED fitting and the [OII] than from the UV+IR. The SFRs determined from the SED-fitting and [OII] are consistent for galaxy 5. However, for galaxy 1 and 4 the SFR from dustcorrected [OII] is smaller than the measurement from the SED-fitting, consistent with similar measurements in Belli et al. (2021). They discuss possible causes for this disagreement and conclude that it could be due to an inadequate dust correction because of the inability to measure nebular dust attenuation by fitting stellar emission. We therefore determine that the most robust SFR measurements for our sample of quiescent galaxies are based on full SED-fitting, and we adopt these SFRs for the remainder of this work. SED-based SFRs are less likely to be contaminated by other sources in quiescent galaxies and are in line with accumulating evidence that SED-fitting produces the least biased of SFR indicators for quiescent galaxies Belli et al. 2017), see also Section 6.2.1 of Leja et al. (2021).

Late time SFH
We look for differences between the late-time SFHs of CO-detected and CO-undetected galaxies. To determine if a galaxy's SFH shows evidence of a secondary burst, we calculate the difference between the peak in the late-time SFH (most recent 1 Gyr) and the previous minimum in SFR. We  Figure 4. The formation redshifts and star-forming timescales as a function of gas mass. There is not a significant difference in the formation redshifts or the star-forming timescales between the CO-undetected galaxies and the CO-detected galaxies.
consider the galaxy to have a secondary burst if the increase in the SFR at its peak is larger than the previous minimum by a factor of 2.0 or more, and if the resulting difference in SFR between peak and previous minimum is significant at 1σ given the uncertainties in the SFRs. Using this criteria, the CO-detected galaxies have possible evidence for secondary bursts of star formation in their late SFH (0-1 Gyr in lookback time), while the CO-undetected galaxies do not. While CO-detected galaxy 3 demonstrates a secondary burst using both priors, the uncertainties are large, and therefore a flat SFH would also be consistent with our measurements. Therefore, for galaxy 3, higher S/N is needed to determine if there is a significant secondary burst of star formation. Galaxies 1, 2 and 4 have statistically significant evidence for secondary bursts of star formation.
We therefore explore whether the evidence for the secondary burst in galaxies 1, 2 and 4 exists under different assumptions, namely our assumed parametric form for the SFH of delayed-tau+burst. The bursts for the parametric and nonparametric models agree quite well for galaxies 1 and 2. The parametric delayed-tau+burst model shows an instantaneous burst at a lookback time of 1 Gyr that is consistent with the burst found in the nonparametric model for these two galaxies. The fraction of mass formed in these bursts is 2-3% for galaxy 1 and 6-7% for galaxy 2, consistent between all 3 models. However, for galaxy 4, the parametric model did not identify a secondary burst with consistent time and mass fraction as the nonparametric model. The fraction of mass formed in this burst is relatively small, only 0.3%, which could explain why it was not identified by the parametric model. We note that the instantaneous burst time measured by the parametric model for this galaxy is early in its history, aligning with the peak epoch of the SFH (see black lines in Figure 2), and forms a relatively minor amount of mass.
Comparisons of the fraction of mass formed in the secondary bursts, f burst , are shown in Table 3 for all 3 models. To determine the fraction of mass formed in the bursts for There is no significant difference between the CO-detected and CO-undetected galaxies. Therefore, the detectable gas reservoirs in the CO-detected galaxies are not leftover from their primary epoch of star formation.
the nonparametric models, f burst , we define a burst as an increase in SF. The bin prior to the increase in SF is considered the baseline. The mass formed is then calculated using all the bins following the baseline that have higher SFRs than the baseline.
Given that we see evidence in a few, and some tentative in others, for a late time burst, we explore whether this late time SFH is a plausible explanation for higher gas content. Figure 6 shows the total mass formed in the last 1 Gyr (and last 500 Myr) for each galaxy. We find that CO-detected galaxies formed a higher fraction of their mass in the last Gyr than CO-undetected galaxies. The average fraction of mass formed in the last Gyr and the standard deviation for CO-detected galaxies is 2.6 ± 1.8% and for CO-undetected galaxies is 0.2 ± 0.1%.
To summarize, our data are consistent with the idea that detectable gas reservoirs in quiescent galaxies are linked to late time star formation activity, and we do not see evidence that it is leftover from their primary star formation epoch. We discuss possible causes for this finding in the next section.

DISCUSSION
The main finding in our study is that quiescent galaxies with detectable gas reservoirs have enhanced late time star formation compared to those without significant gas reservoirs. All four of the CO-detected galaxies have plausible evidence for a secondary burst, although in one case the uncertainties are high enough that we cannot reject a flat level of SF in the last 1 Gyr. In contrast to this, the CO-undetected galaxies have declining SFRs in their most recent 1 Gyr. In addition, the CO-detected galaxies formed a higher fraction of mass in their last Gyr (f M 1Gyr = 2.6 ± 1.8%) than the CO-undetected galaxies (f M 1Gyr = 0.2 ± .1%). Here, we explore what physical mechanisms could be influencing the late time star formation in galaxies with higher molecular gas content and discuss how this may be a contributing factor to the heterogeneity of molecular gas content observed in quiescent galaxies.

Relative timescales of gas depletion and star formation decline
The relative timescale over which cold molecular gas is depleted, and over which the star formation declines, is a powerful diagnostic of the quenching mechanisms at play. In principal, the CO-detected galaxies in our sample are quenching over a longer timescale than the CO-undetected galaxies (emulating a slow and fast timescale for quenching). However, we find that this is due to a late-time boost in star formation, possibly influenced by external sources such as recently accreted gas, and not necessarily because they are undergoing a different, slower quenching mechanism (e.g. Belli et al. 2019). With a limited sample such as ours, it is difficult to draw broad conclusions, so in this section we compare to other available timescale constraints in the literature for coeval (at similar redshifts) non-starforming galaxies.
We compare our results to the most comprehensive coeval sample, which are post-starburst (PSB) galaxies from the SQUIGGLE survey. We note that none of the galaxies in our sample fall in the region of the Dn4000 vs. Hδ plane occupied by the SQUIGGLE sample, or by PSBs in general. The Hδ measurements alone exclude our sample from being classified as PSB galaxies by traditional selections, which require Hδ > 3 − 5Å, while the galaxies in our sample have values of Hδ < 3Å. Our sample also have higher Dn4000 than traditional PSBs, with Dn4000 in the range 1.5-1.9Å, while the PSB galaxies in SQUIGGLE have typical . M1Gyr (M500Myr) is the total formed mass in the last Gyr (500 Myr) for each galaxy. The CO-detected galaxies are shown in purple while the CO-undetected galaxies are shown in red. The bootstrapped uncertainties are shown as errorbars. More mass was formed in the last Gyr for galaxies with higher gas masses. This is further evidence that the detectable gas reservoirs are due to an event in their last Gyr, rather than their early SFH. ranges from 1.2-1.5Å. SQUIGGLE galaxies also have ceased star formation, while exhibiting heterogeneous molecular gas reservoirs (Suess et al. 2017;Hunt et al. 2018;Suess et al. 2022;Bezanson et al. 2022). The SQUIGGLE PSB selection targets strong Balmer breaks and blue slopes redward of the break using rest-frame filters, following Kriek et al. (2010). The SQUIGGLE sample is selected at z ≈ 0.6, only slightly lower redshift than this sample. However, their primary burst of SF occurred shortly before observation time, meaning their sample is made up of more-recently quenched galaxies than our sample. Thus, SQUIGGLE PSBs are qualitatively different in that they are ending their primary burst of star formation, while our quiescent galaxies are ending a small secondary burst, having completed their primary growth epoch earlier in cosmic time. Nonetheless, both samples have ceased star formation while retaining signficant gas reservoirs, and together they increase the baseline in galaxy properties for comparing gas and star formation timescales in quenched populations exhibiting heterogeneous gas content. Bezanson et al. (2022) find that 6 of the 13 SQUIGGLE galaxies were CO-detected. Before comparing to our results, we first must convert their molecular gas masses using the same conversion factors assumed in Spilker et al. (2018) with the standard rescaling M H2 × (0.8/r 21 )(α CO /4.4). Bezanson et al. (2022) assume r 21 = 1.0 and α CO = 4.0. We next check for differences in the survey sensitivities. The range of the upper limits on gas masses for the COundetected galaxies in Bezanson et al. (2022) using our conversion factors is log M H2 = 9.1 − 9.4M , somewhat deeper than the limits of the Spilker et al. (2018) data with log M H2 = 9.7 − 9.9M . For both studies, the upper limits were determined using the noise in the 800 km/s wide channel at the 3σ level. Therefore, we can safely make this comparison since our entire sample of CO-detected sources from (Spilker et al. 2018) would be included as detections by the deeper Bezanson et al. (2022) data.
The main finding in Bezanson et al. (2022) is that the galaxies that quenched ≤150 Myr prior to observation have significant gas reservoirs, with an average gas fraction of f H2 ∼ 7% (or f H2 ∼ 5% using our conversion factors), while older post-starburst galaxies do not. Assuming that these individual galaxies are in different phases of an evolutionary trend, this implies that quenching precedes the disappearance of cold gas in their sample, but that molecular gas must be consumed or destroyed shortly thereafter.
To compare to our sample, we calculate the time between the end of the secondary burst and the time of observation (t end-burst ). We find for galaxies 1, 3, and 4 that t end-burst ≈ 0 − 100Myr, which is consistent with the benchmark suggested by Bezanson et al. (2022) for hosting residual gas. Galaxy 2 (which formed the highest fraction of mass in its secondary burst of star formation among our sample) has a secondary burst t end-burst ≈330 Myr prior to observation, despite exceeding the timescale identified by Bezanson et al. (2022) for gas to disappear in their sample of CO-detected galaxies. Based on the overall total fraction of mass formed during the course of the recent bursts (75% of SQUIGGLE galaxies have burst mass fractions >25% compared to only 0.3-6% for our sample), the gas content that drove the recent low-level star formation in our sample is likely much smaller than the total amount of gas that drove the large, primary bursts of star formation in the PSBs. It is thus plausible that the gas consumption may not proceed on similar timescales given the dramatically different scales of the growth, perhaps resulting in slower gas consumption than the clearly rapid decline in PSBs.
There are several physical reasons that would explain a slower consumption at low gas content. Firstly, it is possible that small amounts of gas can be dynamically stabilized to prevent star formation without any gas destruction or consumption (Gensior et al. 2020). This simulation predicts that dynamical stabilization of molecular gas can suppress star formation by as much as a factor of 5 in bulge-dominated systems with relatively low gas fraction ( 5%) without destroying or consuming gas. In fact, galaxy 2 that exhibits a relatively long timescale since burst is bulge-dominated (see Figure 1) and has a low enough gas fraction to be in line with the conditions of dynamical stabilization (f gas = 4 ± 0.2%). While this may also explain the other 3 CO-detected galaxies that had more recent burst times, this physical reason does not appear to be a universal threshold among quiescent galaxies, given the relatively low f gas limits among our non-detected sample, and among others across redshifts (e.g. Davis et al. 2016;Saintonge et al. 2017;Williams et al. 2021;Whitaker et al. 2021a). An alternative explanation for unconsumed gas in quiescent galaxies is a decrease in star formation efficiency due to a lower fraction of dense molecular gas compared to the total molecular gas (e.g. Li et al. 2019;French et al. 2022). Unfortunately, our data do not provide insight into the density of molecular gas.
While gas depletion times can in principal measure any difference in gas decline directly, they are difficult to interpret among galaxies below the main sequence because the value depends strongly on which indicator is used to measure the SFR, which can vary by as much as a factor of 6 (see Section 3; Belli et al. 2021). We note that both the Bezanson et al. (2022) and Spilker et al. (2018) samples exhibit very long depletion times for gas consumption by the residual star formation (> 2 Gyr). 1 In the case of the PSB sample, such long depletion times are actually inconsistent with the inferred rapid < 150Myr timescale for gas disappearance. However, the measure of depletion time implicitly assumes that gas is only destroyed by simple consumption by star formation and neglects processes such as heating or removing gas via AGN feedback. It is possible that in such cases that gas is instead actively destroyed or heated by AGN feedback rather than simply consumed by star formation in PSBs (e.g. the PSB sample studied in French et al. 2018a). Bezanson et al. 2022 also note the possibility that some COdetected PSBs could appear gas rich because they are caught in a phase immediately prior to rejuvenating star formation (e.g. in a secondary burst such as what we observe). We discuss the possible contribution of rejuvenation to the heterogeneity observed among non-starforming galaxies in the next section.

Rejuvenation
50% of our sample of 8 quiescent galaxies show plausible evidence for secondary bursts of star formation along with significant molecular gas content. This is an interesting result, in light of the prevailing theoretical picture that quiescent galaxies form early in cosmic time, and lack appreciable new star formation for 10 Gyr since quenching (e.g. Renzini 2006;Citro et al. 2016). While this implies rejuvenation is not uncommon around z ∼ 0.8, the observed frequency of rejuvenation in quiescent galaxies remains extremely uncertain, owing to the difficulty of characterizing low-level star formation in quiescent galaxies across cosmic time (see Section 3), the short timescale over which such "frosting" might be observable, and the difficulty of measuring detailed SFH outside of the nearby Universe.
The majority of studies characterizing the prevalence of rejuvenated star formation in quiescent galaxies are from the local universe. By z < 0.1 the rejuvenation fraction appears relatively common, with higher rates among quiescent galaxies at lower masses (10< log M * /M <11) (10- 30% Treu et al. 2005;Schawinski et al. 2007;Thomas et al. 2010;Donas et al. 2007;Salim & Rich 2010;Fang et al. 2012;Salvador-Rusiñol et al. 2020;Werle et al. 2020;Jeong et al. 2022), generally in agreement with theoretical predictions (Trayford et al. 2016;Pandya et al. 2017;Behroozi et al. 2019). Recent advances in higher redshift datasets suggest that the rate that high mass quiescent galaxies such as our sample (log M * > 11M ) undergo rejuvenation events at earlier times are by comparison quite low (<1-16%, depending on the threshold of star formation used to define a rejuvenation; Belli et al. 2017;Chauke et al. 2019;Tacchella et al. 2022b). The lower rejuvenation fraction perhaps reflects the shorter available cosmic time, and therefore lower probability, to both quench and experience processes driving rejuvenation such as minor merging. Nonetheless, examples at even higher redshifts z > 1 are now routinely identified thanks to advances in analyses such as the ones we use here (Jørgensen et al. 2014;Williams et al. 2021;Akhshik et al. 2021;Caliendo et al. 2021;Belli et al. 2021).
Generally the criteria used to define rejuvenation identify relatively large secondary bursts: they typically must increase the star formation in quiescent galaxies to within range of the main sequence scatter, with a large range of typical burst masses that can form between ∼ 2 − 25% of the galaxy's mass (Chauke et al. 2019;Tacchella et al. 2022b). We note here that our secondary bursts are small compared to published metrics of rejuvenation in coeval galaxies, with typical mass fractions formed f burst ∼ 2% with a range 0.3 − 6% (see Table 3). Based on our inferred SFHs, one of our CO-detected sources (galaxy 1) meet the stringent criteria for identifying rejuvenations at z ∼ 0.8 (e.g. Tacchella et al. 2022b, although we note that these authors may find a low rejuvenating fraction because they specifically select for rejuvenation that is ongoing at the time of observation). If we exclude the requirement in Tacchella et al. (2022b) that the rejuvenation be ongoing at time of observation, galaxy 4 would be identified using the Tacchella et al. (2022b) criteria. Our secondary bursts are also more minor than the fiducial rejuvenation measured by Chauke et al. 2019, half of which form stellar mass fractions f burst 8% in rejuvenation. Chauke et al. 2019 identify rejuvenations as quiescent galaxies with a prior quiescent phase, with a return to the star forming sequence in between. The star forming and quiescent criteria are based on a sSFR threshold compared to the main sequence (>0.3 dex below), that evolves with redshift as defined by Speagle et al. (2014). Two galaxies (galaxies 1 and 2) meet their criteria using our new SFH measurements (roughly 25% of our sample, slightly higher than their estimate of 16% for the fraction of quiescent galaxies with a rejuvenation). These two galaxies also meet their more-relaxed criteria for a rejuvenation that star formation increase above a fixed limit of log sSFR< −1 Gyr −1 after quiescence, re-sulting in a consistent fraction of rejuvenation (Chauke et al. 2019 find 24±2%).
We thus find that these lower-level rejuvenation events (that we see in our SFH analysis but are missed by previous rejuvenation criteria) may be more common than stringent criteria imply. This may be in agreement with findings of Carnall et al. 2018 that up to 60% of quiescent galaxies may have experienced low level rejuvenation or a merger since z < 1.5 (see also Mancini et al. 2019;Carleton et al. 2020). Some evidence from theory and cosmological simulations imply that a similar fraction of massive galaxies undergo weak rejuvenations (10-40% of massive halos with log M halo > 13), supporting the idea that weak rejuvenations are relatively common (Alarcon et al. 2022). While the rejuvenations we measured are relatively small, they nonetheless contribute to the surprisingly large number of detectable reservoirs of molecular gas. We explore whether the low level secondary bursts we measure are consistent with common quiescent galaxy evolutionary processes such as minor merging in the next section.

Minor Mergers and size evolution
While the prevalence of star formation rejuvenation in high-redshift quiescent galaxies still remains uncertain, it is well known that minor merging must be a significant process in their late time evolution after quenching their primary growth phase. At fixed stellar mass, high redshift quiescent galaxies are 3-5 times smaller than in the local universe (Daddi et al. 2005;Trujillo et al. 2006;Buitrago et al. 2008;Cimatti et al. 2008;van Dokkum et al. 2008), implying that they have increased significantly in size but not mass over the past 10 Gyr. While gaseous mergers can dissipate energy to collapse to form stars centrally in galaxies (serving to decrease size and increase mass) dry minor merging can match observations by increasing size without adding appreciable mass. Empirical studies find that minor merging is likely the dominant mechanism for size growth over 0 < z < 1.5 (van Dokkum et al. 2010;Belli et al. 2014;Hamadouche et al. 2022) and can account for most or all of the size growth at z 1 (Newman et al. 2012). Therefore, due to its ability to explain these empirical trends, and the direct evidence that the process occurs (e.g. the presence of tidal debris composed of old stars, and prevalence of quiescent satellites around quiescent massive galaxies van Dokkum et al. 2010;Ji et al. 2018), dry minor merging has become the prevailing theoretical picture influencing late time quiescent galaxy evolution (Bezanson et al. 2009;Hopkins et al. 2009;Naab et al. 2009;van de Sande et al. 2013;van Dokkum et al. 2010;Oser et al. 2012).
The CO-detected galaxies in our sample tend to have lower stellar densities than the CO-undetected galaxies, which is consistent with having experienced minor merging (Spilker et al. 2018). Given that our data also indicate plausible evidence for accreted gas, perhaps delivered via minor mergers, we explore whether the level of gas measured in the CO-detected sample is consistent with the expectations of dry minor merging. Minor mergers are typically defined as mergers between two galaxies with a mass ratio in the range of 10 − 4 : 1 (Kaviraj et al. 2009(Kaviraj et al. , 2012, and are considered dry if the average gas fraction of both galaxies is gas-poor. To check this possibility, we estimate the limiting quantity of molecular gas allowed in a dry minor merger, assuming our CO-detected galaxies are the larger of the merging pair. For a dry merger, we define a gas-poor galaxy using a very conservative assumption of f gas < 20% (relatively high gas content) following Man et al. (2012, see also Kaviraj et al. 2012. For the minor galaxy, assuming a 10:1 mass ratio, this would correspond to the minor galaxy bringing in gas masses of log M gas ≈ 9.4 − 9.7M (assuming the range of stellar masses of the CO-detected galaxies as the larger of the merger pair). Assuming a 4:1 ratio predicts a range of log M gas ≈ 9.8 − 10.1M . Comparatively, the range in measured gas masses in our CO-detected sample is log M gas = 9.8 − 10.2M .
Thus, we find that a dry minor merger with a 10:1 merger ratio would bring in a lower mass of molecular gas than what is measured for all 4 of our CO-detected galaxies, indicating that our galaxies are inconsistent with the expectations of a dry minor merger with this mass ratio. Assuming a less restrictive 4:1 minor merger ratio, we calculate that a dry minor merger would still be inconsistent with our large measured gas reservoirs in galaxies 3 and 4. However, for galaxies 1 and 2, the limiting gas content calculated for a less restrictive 4:1 dry minor merger ratio are consistent with the gas content that we observe. (We note that since the secondary bursts are nearing their completion at the time of observation, the gas content brought in by the companion galaxy would have to be larger than these estimates). Therefore, we conclude that the detectable gas reservoirs in our CO-detected galaxies are likely too large to result from typical dry minor mergers, favoring instead gas-rich minor mergers as the source of the gas and secondary burst.
While the accreted gas is plausibly from gas-rich minor mergers, Spilker et al. (2018) argue against the idea of an outside gas source such as gas-rich minor mergers because the rotational axes of the stellar and molecular components of 3 of the 4 CO-detected galaxies are consistent with one another. Several studies have suggested that galaxy mergers would cause star-gas misalignments (Bureau & Chung 2006;Crocker et al. 2009;Martel & Richard 2020;Khim et al. 2021). Khim et al. (2021) study galaxies from the Horizon-AGN simulation, which is a large-volume cosmological hydrodynamical simulation (Dubois et al. 2014) run with the adaptive mesh refinement code RAMSES (Teyssier 2002). They find that major mergers are 2.5× more effective in causing star-gas misalignments than minor mergers (with the criterion for misalignments being >30°). While the simulations show that the average lifetime of gas misalignment in quiescent galaxies is long (∼2 Gyr), the scatter in lifetimes is quite large, indicating that misalignments indicative of past mergers can decay even on a few hundred Myr timescale (Khim et al. 2021). This is similar to the decline timescale we observe in our secondary bursts. We conclude that a lack of star-gas misalignment does not necessarily rule out the possibility that enhanced late time star formation in our sample could be due to minor mergers that occurred more than a few hundred Myr prior.
Even though we find that any minor merging in our sample is unlikely dry, it is possible they can still contribute to size evolution because the gas mass accreted in the merger is still relatively small. Sonnenfeld et al. (2014) show that purely dry mergers cannot be the only mechanism driving the size and density evolution of massive early type galaxies. They find that dry minor mergers predict a strong decrease in the mass density profile slope, inconsistent with the nearly constant slope inferred from observations over 0 < z < 1.
To match the observed evolution in both the size and the mass density profile slope, some dissipation is required and their models prefer a gas fraction of f gas ≈ 8% for the minor merging companion, higher than what is measured for 3 of the 4 CO-detected galaxies. While we do not have the data to constrain whether the gas at the time of accretion met this criteria, we can use it as an upper limit on the gas fraction that would be observed after the minor merger, since some amount of accreted gas is consumed as the merger completes. Sonnenfeld et al. (2014) indicate that a gas fraction of 8% in the companion on average results in a 4% increase in mass of new stars formed. This is within the range of f burst values that we measure among our CO-detected sample, f burst = 0.3 − 6%. As such, to the best of the ability of the data, we find that our data are plausibly in line with the theory. Both Mancini et al. (2019) and Ji & Giavalisco (2022) find evidence that higher-mass (log M * > 11 M ) quiescent galaxies quench and then grow due to merging with satellite star-forming (gas-rich) galaxies based on more extended light profiles in the I814 (bluer) band than in the J125 and H160 (redder) bands. In addition, Belli et al. (2017) find lowlevel star formation activity based on weak Hα emission in nine quiescent galaxies that is likely to be fueled by inflowing gas or gas-rich minor mergers. Rutkowski et al. (2014) find recent bursts of star formation that they conclude could be due to gas-rich minor merger activity, as well. Further, latetime star formation from gas-rich minor mergers increases SFR during these secondary bursts, which could slightly decrease the average light-weighted age while also increasing the size. This would be consistent with the size-age trend found in Williams et al. (2017) where younger galaxies are larger, see also Fagioli et al. (2016). This is often attributed to progenitor bias. With our small sample, it is impossible to draw conclusions about the broader massive galaxy population, but we speculate that part of the size-age trend attributed to progenitor bias could be due to gas-rich minor mergers.
Finally, widespread evidence from quiescent galaxies in the local Universe indicates that minor merging is a significant resource for accreted cold gas, fueling very low levels of star-formation in otherwise inactive systems (Kaviraj et al. 2009;Davis et al. 2011;Kaviraj et al. 2011Kaviraj et al. , 2012Kaviraj et al. , 2013van de Voort et al. 2018;Davis & Young 2019). Our data provide evidence that this process also impacts massive galaxies at earlier times, in line with previous findings that gas rich minor merging probably occurred in this population consistently over the past 8 Gyr, and may even be the dominant process for new star formation at z < 1 (e.g. Kaviraj et al. 2008;Vulcani et al. 2016). The extent to which gas rich minor mergers impact the quiescent galaxy population at z > 1 remains unknown, although possible examples do exist (Williams et al. 2021;Akhshik et al. 2021;Caliendo et al. 2021;Belli et al. 2021) and warrants further study. We posit that gas-rich minor mergers could thus be a contributing factor to the heterogeneity of molecular gas reservoirs observed in quiescent galaxies, a process that increases with frequency with cosmic time since z < 1.

CONCLUSIONS
We investigate the nonparametric SFHs of 8 quiescent galaxies with existing cold molecular gas mass measurements from Spilker et al. (2018). Four of the eight galaxies contain signficant gas reservoirs. Our main results are as follows: • We find that the molecular gas content is unrelated to the rate of star formation decline in the early SFH. This suggests that the gas reservoirs are not leftover from their primary SF epoch.
• Quiescent galaxies with detectable gas reservoirs have enhanced late time star formation. We suggest the gas may have accreted via gas-rich minor mergers within the last 1 Gyr based on evidence for a recent low-level rejuvenation (f burst ≈ 0.3 − 6%) in their SFHs. This is also supported by our finding that the galaxies with higher gas masses have formed more mass in the last Gyr of their SFH. Gas-rich minor mergers therefore may contribute to the diversity of molecular gas observed in high-redshift quiescent galaxies.
While the current sample is small, increasing the sample size of galaxies with both detailed SFH measurements and cold ISM constraints will increase our understanding of the quenching process. Studying the SFHs of galaxies in upcoming deep photometric and spectroscopic surveys with the James Webb Space Telescope will provide an excellent sample to make follow-up molecular gas observations with ALMA. Finally, we note that although CO observations even with ALMA are relatively expensive, future progress to understand the dust to molecular gas ratios in high-redshift quiescent galaxies (e.g. Whitaker et al. 2021b) would more robustly enable efficient constraints on the cold molecular gas in quiescent galaxies based only on dust emission.
We thank the referee for a thorough report that greatly strengthened this work. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 194-A.2005 (The LEGA-C Public Spectroscopy Survey). CW is supported by the National Science Foundation through the Graduate Research Fellowship Program funded by Grant Award No. DGE-1746060. C.C.W. and M.R. were supported by the National Aeronautics and Space Administration (NASA) Contract NAS50210 to the University of Arizona.
This material is based upon High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the UArizona Research Technologies department.
We respectfully acknowledge that the University of Arizona is on the land and territories of Indigenous peoples. Today, Arizona is home to 22 federally recognized tribes, with Tucson being home to the O'odham and the Yaqui.