SQuIGG$\vec{L}$E: Buried star formation cannot explain the rapidly fading CO(2-1) luminosity in massive, $z\sim0.7$ post-starburst galaxies

TL;DR

This study uses deep ALMA observations to investigate CO luminosity in post-starburst galaxies at z~0.7, revealing rapid gas depletion not explained by dust-obscured star formation, challenging previous assumptions about galaxy quenching.

Contribution

It provides new evidence that CO luminosity declines rapidly after starbursts, suggesting a shorter gas-removal timescale and questioning the evolutionary sequence of post-starburst galaxies.

Findings

CO luminosity correlates with starburst age

Buried star formation cannot fully explain CO decline

Gas-rich galaxies may rejuvenate soon

Abstract

Observational and theoretical studies have long held that rapid gas consumption in starbursts is responsible for the formation of quiescent galaxies. However, studies of recently quenched ``post-starburst" galaxies have discovered that a number of them are surprisingly luminous in CO, challenging this assumption. We present deep ALMA CO(2-1) observations of 50 massive ($\log(M_\star/M_\odot)\sim11.2$) post-starburst galaxies from the SQuIGG$\vec{L}$E sample at $z\sim0.7$. We detect a large fraction (27/50) of the galaxies in CO(2-1). Furthermore, we find that the CO luminosity correlates with the age of the recent starburst, suggesting a gas-removal timescale of $\lesssim140$ Myr, an order of magnitude shorter than is implied by their rest optical star formation rates. We perform new spectral energy distribution fits incorporating mid- and far-IR photometry to test whether dust-obscured star formation can explain this trend. We find that while allowing for buried star formation can raise star formation rates by $\sim0.5$ dex, for almost all galaxies it is neither required to fit the observed IR SED, nor is it sufficient to explain the observed depletion trend. Even the combination of significant buried star formation and ULIRG-like $\alpha_{CO}$ is not enough to explain this decay in CO luminosity. Furthermore, there is no strong evidence to support either of those modifications to the depletion time. Therefore, it remains a distinct possibility that the age-CO luminosity trend should not be interpreted as an evolutionary sequence, and that gas-rich SQuIGG$\vec{L}$E galaxies will soon rejuvenate.

Figures

Tables

Table 1. Table 1: Physical properties of SQuIGG L → \vec{L} E  galaxies. Only the three galaxies included in Figure 6 , 7 , and 8 are presented, but the full table is available in a machine readable form online.

ID	$\log \left(M_{\star }/M_{\odot }\right)$	$\mathrm{SFR}$aaMeasured in the 10 Myr before observation	$\log \left(\mathrm{sSFR}\left[{\mathrm{yr}}^{-1}\right]\right)$	${\tau }_{BC}$bbdust1 and dust2 in fsps.	${\tau }_{ISM}$bbdust1 and dust2 in fsps.	$\mathrm{Dust}\mathrm{Index}$	$t_{\mathrm{PSB}}$ccDefined as the time at which the galaxy formed 99% of the total stellar mass in the Gyr before observation.	$L_{IR}$ddDefined as the total luminosity of all dust components in our modeling.
		[$M_{\odot }{\mathrm{yr}}^{-1}$]					[Gyr]	[${10}^{11}{L}_{\odot }$]
No Buried Star Formation:
J1157+0132	${11.41}_{-0.04}^{+0.04}$	${17.57}_{-6.32}^{+9.55}$	$-{10.16}_{-0.19}^{+0.19}$	${0.59}_{-0.20}^{+0.26}$	${0.55}_{-0.10}^{+0.08}$	$-{0.05}_{-0.12}^{+0.11}$	${0.04}_{-0.01}^{+0.03}$	${6.25}_{-0.90}^{+0.96}$
J0910+0218	${11.13}_{-0.08}^{+0.06}$	${14.86}_{-3.46}^{+6.34}$	$-{9.97}_{-0.13}^{+0.22}$	${0.80}_{-0.13}^{+0.18}$	${0.82}_{-0.07}^{+0.09}$	$-{0.55}_{-0.08}^{+0.09}$	${0.02}_{-0.00}^{+0.00}$	${5.86}_{-0.72}^{+1.20}$
J1141-0109	${11.40}_{-0.03}^{+0.02}$	${0.00}_{-0.00}^{+0.08}$	$-{14.89}_{-1.96}^{+2.38}$	${0.23}_{-0.12}^{+0.10}$	${0.19}_{-0.05}^{+0.07}$	$-{0.45}_{-0.22}^{+0.25}$	${0.23}_{-0.03}^{+0.05}$	${0.64}_{-0.16}^{+0.18}$
…	…	…	…	…	…	…	…	…
Buried Star Formation Allowed:
J1157+0132	${11.40}_{-0.05}^{+0.04}$	${75.68}_{-24.17}^{+25.60}$	$-{9.52}_{-0.17}^{+0.13}$	${2.06}_{-0.64}^{+0.41}$	${0.48}_{-0.14}^{+0.09}$	$-{0.07}_{-0.17}^{+0.08}$	${0.02}_{-0.01}^{+0.00}$	${8.67}_{-1.23}^{+1.44}$
J0910+0218	${11.12}_{-0.05}^{+0.06}$	${26.58}_{-12.31}^{+15.50}$	$-{9.71}_{-0.29}^{+0.22}$	${1.26}_{-0.48}^{+0.41}$	${0.81}_{-0.07}^{+0.08}$	$-{0.56}_{-0.07}^{+0.10}$	${0.02}_{-0.00}^{+0.02}$	${6.47}_{-1.05}^{+1.25}$
J1141-0109	${11.44}_{-0.05}^{+0.03}$	${0.00}_{-0.00}^{+0.13}$	$-{14.47}_{-1.66}^{+2.17}$	${1.22}_{-0.74}^{+1.06}$	${0.21}_{-0.05}^{+0.06}$	$-{0.26}_{-0.19}^{+0.19}$	${0.25}_{-0.04}^{+0.06}$	${0.61}_{-0.14}^{+0.14}$
…	…	…	…	…	…	…	…	…

Table 2. Table 2: Properties of the ALMA observations.

ID	RA	Dec	z	Program IDaaSome sources were observed in multiple cycles due to not being flagged as duplicates when they were observed between submission of proposals and the beginning of the new cycle. For these sources, we list both program IDs and combine all data our imaging.	Observing Date	Integration Time	Angular Resolution
	[deg]	[deg]				[s]	[′′]
J0131+0034	22.885654	$+$0.577446	0.618268	2021.1.01535.S	3/26/22	2999.8	2.1
J0153-0207	28.265102	$-$2.126556	0.66585	2021.1.01535.S/2022.1.00604.S	12/22/2022, 4/30/2022	7801.9	1.2
J0221-0646	35.363172	$-$6.778846	0.661296	2021.1.01535.S	3/22/22	2993.8	2.2
J0224+0015	36.178065	$+$0.253695	0.684893	2021.1.01535.S	3/26/22	2963.5	2.1
J0224-0034	36.126895	$-$0.571105	0.743342	2022.1.00604.S	12/24/22	2993.8	1.5
J0224-0630	36.039585	$-$6.511067	0.751483	2021.1.01535.S	1/28/22	2993.8	1.4
J0227-0225	36.99331	$-$2.430015	0.61386	2021.1.00988.S	1/28/22	5927.0	1.4
J0232-0331	38.209972	$-$3.52757	0.738828	2022.1.00604.S	1/26/22	5927.0	1.4
J0237+0123	39.265248	$+$1.390005	0.709438	2021.1.01535.S/2022.1.00604.S	9/11/2022, 12/27/2022	5999.6	1.2
J0240-0148	40.054477	$-$1.81099	0.689513	2021.1.01535.S	1/27/22	2993.8	1.4
J0851+0244	132.76778	$+$2.742179	0.659483	2022.1.00604.S	12/17/22	2999.8	1.4
J0907+0423	136.99116	$+$4.384259	0.663469	2021.1.00988.S	4/21/22	5806.1	1.5
J0909-0108	137.47026	$-$1.134891	0.702088	2021.1.01535.S	1/27/22	2993.8	1.4
J0910+0218	137.61906	$+$2.309311	0.769447	2022.1.00604.S	12/31/22	2993.8	1.6
J0940-0008	145.1991	$-$0.138624	0.618503	2022.1.00604.S	12/28/22	2999.8	1.4
J1017-0003	154.30202	$-$0.061552	0.596369	2022.1.00604.S	12/27/22	2993.8	1.3
J1040+0223	160.04608	$+$2.391646	0.680527	2022.1.00604.S	1/3/23	2993.8	1.4
J1042+0500	160.71635	$+$5.010805	0.626593	2022.1.00604.S	12/31/22	2993.8	1.4
J1046+0123	161.7372	$+$1.399393	0.764902	2022.1.00604.S	1/19/23	2993.8	1.0
J1114+0115	168.60707	$+$1.265318	0.748912	2022.1.00604.S	1/19/23	2993.8	1.0
J1141-0109	175.38035	$-$1.151463	0.657525	2022.1.00604.S	1/2/23	2993.8	1.5
J1142+0006	175.53723	$+$0.105861	0.593473	2021.1.00988.S	5/8/22	5745.6	1.3
J1148+0204	177.1627	$+$2.069523	0.596461	2022.1.00604.S	12/31/22	2993.8	1.4
J1157+0132	179.48879	$+$1.537532	0.755935	2021.1.00988.S	5/9/22	2993.8	1.5
J1211+0240	182.90725	$+$2.673695	0.589034	2022.1.00604.S	10/16/22	2993.8	1.7
J1222+0342	185.65243	$+$3.711715	0.679837	2021.1.01535.S	6/3/22	2993.8	1.0
J1240-0057	190.18837	$-$0.951925	0.795872	2022.1.00604.S	1/19/23	2993.8	1.1
J1244+0248	191.00193	$+$2.807904	0.625788	2021.1.01535.S	8/30/22	2993.8	0.9
J1332+0256	203.11942	$+$2.933643	0.698313	2021.1.01535.S/2022.1.00604.S	8/28/2022, 12/27/2022	5987.5	1.2
J1416+0255	214.00084	$+$2.924476	0.760227	2022.1.00604.S	12/27/22	2993.8	1.6
J1436+0447	219.16637	$+$4.795968	0.633855	2021.1.00761.S	8/28/22	5987.5	0.9
J1437+0311	219.43569	$+$3.19169	0.667036	2021.1.01535.S/2022.1.00604.S	8/28/22	5987.5	1.2
J1444-0006	221.03227	$-$0.106326	0.700658	2021.1.01535.S	4/6/22	2999.8	2.3
J1449+0206	222.33699	$+$2.101912	0.743967	2022.1.00604.S	12/26/22	2999.8	1.7
J1455-0048	223.93174	$-$0.8092	0.625287	2021.1.01535.S	4/6/22	2999.8	2.1
J2213-0050	333.45469	$-$0.833433	0.66203	2021.1.01535.S/2022.1.00604.S	4/06/2022, 12/20/2022	4814.2	1.5
J2232+0007	338.16968	$+$0.131622	0.722158	2021.1.01535.S/2022.1.00604.S	9/10/2022, 12/29/2022	5987.5	1.2
J2310-0047	347.59742	$-$0.798654	0.737807	2021.1.00988.S	1/25/22	5999.6	1.4

Table 3. Table 3: CO(2-1) and continuum measurements in a 2 ′′ aperture. All non-detections are listed as 3 σ \sigma upper limits. For sources in the Herschel/SPIRE 250 250 μ \mu m source catalog, we also list their flux density.

ID	SdvCO(2–1)	LCO(2-1)	$M_{H_2}$aaAssuming $r_{21}=1.0$ and ${\alpha }_{CO}=4.0$	$f_{2mm}$	$f_{250\mu m}$
	[J km $s^{-1}$]	[${10}^9$ J km $s^{-1}{\mathrm{pc}}^2$]	[${10}^{10}{M}_{\odot }$]	[$\mu$Jy]	[mJy]
J0131+0034	$0.05\pm 0.01$	$0.26\pm 0.07$	$0.11\pm 0.03$		–
J0153-0207	$0.38\pm 0.04$	$2.37\pm 0.23$	$0.95\pm 0.09$		–
J0221-0646					–
J0224+0015					–
J0224-0034	$0.51\pm 0.06$	$3.99\pm 0.45$	$1.59\pm 0.18$	$146.14\pm 15.17$	–
J0224-0630	$0.63\pm 0.06$	$5.0\pm 0.51$	$2.0\pm 0.2$		–
J0227-0225	$1.09\pm 0.04$	$5.79\pm 0.22$	$2.32\pm 0.09$		–
J0232-0331					–
J0237+0123	$0.62\pm 0.06$	$4.38\pm 0.44$	$1.75\pm 0.17$		–
J0240-0148					–
J0851+0244	$0.61\pm 0.06$	$3.77\pm 0.39$	$1.51\pm 0.16$	$65.49\pm 18.65$	–
J0907+0423	$2.32\pm 0.03$	$14.41\pm 0.17$	$5.76\pm 0.07$	$93.69\pm 10.64$	–
J0909-0108	$1.0\pm 0.14$	$6.98\pm 1.01$	$2.79\pm 0.4$		–
J0910+0218	$1.82\pm 0.04$	$15.28\pm 0.37$	$6.11\pm 0.15$	$42.91\pm 11.76$	$50.3\pm 13.8$
J0940-0008					–
J1017-0003					–
J1040+0223					–
J1042+0500	$0.2\pm 0.05$	$1.12\pm 0.29$	$0.45\pm 0.11$		–
J1046+0123					–
J1114+0115	$0.5\pm 0.14$	$3.99\pm 1.09$	$1.6\pm 0.44$		–
J1141-0109					–
J1142+0006	$2.27\pm 0.05$	$11.22\pm 0.24$	$4.49\pm 0.1$	$77.11\pm 12.65$	$43.2\pm 13.1$
J1148+0204					–
J1157+0132	$3.49\pm 0.05$	$28.26\pm 0.4$	$11.3\pm 0.16$	$102.37\pm 13.71$	$84.3\pm 19.6$
J1211+0240	$0.24\pm 0.06$	$1.18\pm 0.27$	$0.47\pm 0.11$		–
J1222+0342	$0.3\pm 0.08$	$1.96\pm 0.54$	$0.78\pm 0.21$		–
J1240-0057					–
J1244+0248	$0.86\pm 0.18$	$4.73\pm 0.97$	$1.89\pm 0.39$		–
J1332+0256	$0.23\pm 0.06$	$1.61\pm 0.4$	$0.64\pm 0.16$	$47.93\pm 12.06$	–
J1416+0255					–
J1436+0447	$1.23\pm 0.07$	$6.94\pm 0.42$	$2.78\pm 0.17$		–
J1437+0311					–
J1444-0006	$0.3\pm 0.04$	$2.06\pm 0.27$	$0.82\pm 0.11$		–
J1449+0206	$0.24\pm 0.06$	$1.87\pm 0.51$	$0.75\pm 0.2$		–
J1455-0048					–
J2213-0050					–
J2232+0007					–
J2310-0047	$1.13\pm 0.04$	$8.7\pm 0.34$	$3.48\pm 0.14$		–

Table 4. Table 4: Alternate measures of the age of the starburst for these galaxies. The full table is available in a machine readable form.

ID	$t_{\mathrm{PSB},90}$	$t_{q,{10}^{-10}{\mathrm{yr}}^{-1}}$
	[Gyr]	[Gyr]
No Buried Star Formation:
J1157+0132	0.117${\pm }_{0.018}^{0.03}$	0.01${\pm }_{0.009}^{0.041}$
J0910+0218	0.084${\pm }_{0.017}^{0.017}$	0.001${\pm }_{0.0}^{0.009}$
J1141-0109	0.272${\pm }_{0.023}^{0.028}$	0.225${\pm }_{0.034}^{0.052}$
…	…
Buried Star Formation Allowed:
J1157+0132	0.168${\pm }_{0.049}^{0.053}$	0.001${\pm }_{0.0}^{0.0}$
J0910+0218	0.082${\pm }_{0.014}^{0.019}$	0.001${\pm }_{0.0}^{0.0}$
J1141-0109	0.284${\pm }_{0.024}^{0.048}$	0.251${\pm }_{0.038}^{0.06}$
…	…

Full text

††journal: ApJ††facilities: ALMA, SDSS, Subaru, WISE, Herschel††software: Astropy (Astropy Collaboration et al., 2013, 2018, 2022), Matplotlib (Hunter, 2007), Flexible Stellar Population Synthesis (Conroy et al., 2009; Conroy & Gunn, 2010), SEDPy (Johnson, 2019), Prospector (Johnson & Leja, 2017; Leja et al., 2017; Johnson et al., 2021)††thanks: Email: davidsetton@princeton.edu††thanks: Brinson Prize Fellow

1 Introduction

One of the largest unsolved questions in galaxy evolution is the physical cause of the cessation of star formation–or “quenching”–in the most massive galaxies. The observed bi-modality in galaxy color (e.g,. Shen et al., 2003) and morphology (e.g., van der Wel et al., 2014) suggests that blue, star-forming disks and red, quiescent elliptical galaxies constitute different populations. This implies that the transformation of morphology must have accompanied the cessation of star formation. Additionally, this implies that once galaxies are “red and dead,” they must continue to suppress star formation to remain so, likely due to the influence of active galactic nucleus (AGN) feedback (Hopkins et al., 2006; Croton et al., 2006).

A clear feature that has emerged in the study of large galaxy populations is the presence of a so-called “green valley,” a gap between star-forming galaxies and the quiescent populations that emerges in color-magnitude diagrams and the star forming main sequence (SFMS). The sparse population of this region, especially at high masses, indicates that galaxies spend very little time at intermediate age and star formation rates, and instead quench rapidly (e.g., Schawinski et al., 2014). This has led to paradigm of two quenching pathways, with a slow mode dominating at low-redshift and a rapid mode dominating in the early universe (e.g., Belli et al., 2019). Far-IR evidence of a hidden population of dust-obscured star-forming galaxies in the valley can allow for smoother evolution that diminishes the importance of the rapid quenching track (e.g., Eales et al., 2018). However, the implied short formation timescale of the most massive quiescent systems from abundance and star formation history measurements, both in the oldest galaxies in the local universe (e.g., Thomas et al., 2005; Greene et al., 2013; McDermid et al., 2015; Ferré-Mateu et al., 2017) and in evolved systems at and above cosmic noon (e.g., Kriek et al., 2016; Glazebrook et al., 2017; Carnall et al., 2019; Man et al., 2021; Antwi-Danso et al., 2025; Carnall et al., 2023a; Glazebrook et al., 2024; Beverage et al., 2024; de Graaff et al., 2025; Carnall et al., 2024; Beverage et al., 2025; Zhang et al., 2025; McConachie et al., 2025), does suggest that a rapid quenching pathway is important in the formation of massive galaxies at high-redshift.

Observationally, one can study the rapid quenching process via its immediate descendants, post-starburst galaxies (Dressler & Gunn, 1983; Zabludoff et al., 1996). These galaxies exhibit deep Balmer absorption features, indicating that that their stellar populations are dominated by $\sim$a few hundred Myr old stellar populations as the result of a rapid decline in their star formation rate. Post-starburst galaxies are exceedingly rare in the local universe (e.g., Wild et al., 2009; Pattarakijwanich et al., 2016). However, their rising number density with redshift in conjunction with the falling number density of older quiescent galaxies (e.g., Whitaker et al., 2012a; Wild et al., 2016; Rowlands et al., 2018; Belli et al., 2019; Park et al., 2023; Setton et al., 2023; Park et al., 2024; Adscheid et al., 2025) indeed suggests that they represent the evolutionary link between rapid quenching and massive elliptical galaxies. As the direct products of rapid quenching, studying post-starburst galaxies can illuminate the physical causes of the rapid quenching process.

Perhaps the most important physical pieces of the puzzle of rapid quenching is the state of the molecular gas that fuels star formation. Massive, evolved quiescent galaxies are predominantly gas poor, and have very low molecular gas fractions ($M_{\mathrm{gas}}/M_{\star}$) relative to co-eval star forming systems, whether in CO and dust continuum (the best extragalactic tracers of the molecular gas mass) surveys and stacking experiments in the local volume (e.g., Saintonge et al., 2011; Young et al., 2011; Davis et al., 2016; Michałowski et al., 2019, 2024) or at high-z (Sargent et al., 2015; Spilker et al., 2018; Williams et al., 2021; Whitaker et al., 2021a; Adscheid et al., 2025). While there is indirect evidence from stacking analysis of dust continuum and $A_{V}$ measurements that typical quiescent systems at high-z may be more ISM-rich than their low-z counterparts (e.g., Gobat et al., 2018; Martis et al., 2019; Akhshik et al., 2023; Setton et al., 2024; Lee et al., 2024; Siegel et al., 2025), the bimodality in gas content across cosmic time implies that depletion of molecular gas reservoirs is a prerequisite to quenching. As such, theoretical models have largely sought to explain the mechanisms by which galaxies halt additional accretion from the cosmic web (e.g., Kereš et al., 2005; Dekel & Birnboim, 2006; Feldmann & Mayer, 2015). In these models, the majority of the molecular gas that fueled the burst of star formation prior to quenching would be consumed in the burst and driven out by the combination of quasar- and star-formation-driven outflows, and by the time galaxies are observed as post-starburst, they should already be gas poor (e.g., Hopkins et al., 2008).

It has therefore been surprising that follow-up observations of young quiescent galaxies, in contrast with more evolved quiescent systems, have found that most are CO luminous, implying molecular gas fractions as high as 20% in systems that appear to lie below the star forming main sequence (French et al., 2015; Rowlands et al., 2015; Alatalo et al., 2016a; Suess et al., 2017; Yesuf et al., 2017; Smercina et al., 2018; French et al., 2018; Belli et al., 2021; Bezanson et al., 2022; Woodrum et al., 2022; Otter et al., 2022; Wu et al., 2023; Baron et al., 2023; Michałowski et al., 2024; Umehata et al., 2025; Suess et al., 2025). Perhaps even more surprising, the molecular gas fraction is strongly correlated with the post-burst age, such that the youngest post-starburst galaxies are the most CO-luminous. This implies that these gas reservoirs are being removed on short ($\sim$hundreds of Myr) timescales in a way that is incommensurate with the star formation rate, at least as measured by rest-optical tracers (French et al., 2015; Bezanson et al., 2022).

There have been numerous proposed explanations for the state of this molecular gas, and for the ultimate evolutionary pathway, of these post-starburst systems. Ionized (Tremonti et al., 2007; Sell et al., 2014; Baron et al., 2018, 2020; Diamond-Stanic et al., 2021; Fodor et al., 2025), neutral (Alatalo et al., 2016b; Baron et al., 2020, 2022; Luo et al., 2022; Belli et al., 2024; Park et al., 2024; Sun et al., 2025; Wu, 2025; Valentino et al., 2025), and molecular (Geach et al., 2014, 2018; Spilker et al., 2020a, b) outflows have been identified in recent starburst and post-starburst systems, but it is unclear whether these starburst or AGN driven winds can remove the bulk of the molecular gas seen in these systems. The CO reservoirs of post-starburst galaxies have been found to be morphologically disturbed (Smercina et al., 2022; Sun & Egami, 2022), and in some systems as much at 50% of the CO luminosity has been found on extended scales of tens of kiloparsecs along tidal features (Spilker et al., 2022; D’Onofrio et al., 2025), suggesting that merger driven disruption may play a role in removing gas and suppressing star formation. Together, the injection of energy by the supernovae, AGN, and mergers that drive this removal and the dynamical effects of the compact post-starburst cores (Almaini et al., 2017; Maltby et al., 2018; Wu et al., 2020; Diamond-Stanic et al., 2021; Setton et al., 2022; Zhang et al., 2024) may play a role in morphologically quenching the remaining gas and stopping it from forming new stars (e.g., Martig et al., 2009).

However, the finding that many gas-rich post-starburst galaxies are mid-IR luminous (e.g., Alatalo et al., 2017; Yesuf et al., 2017; Smercina et al., 2018; Baron et al., 2022; Suess et al., 2022a) has also promoted the potential interpretation that these CO-luminous galaxies that present in the rest-optical as post-starburst are actually still in the midst of their bursts, hiding deeply obscured starbursts behind optically thick dust (e.g., Smail et al., 1999; Poggianti & Wu, 2000; Baron et al., 2023). In this physical picture, the CO luminosity can be explained as the fuel for these active starbursts, and the depletion time of the molecular gas can be commensurate with the observed fading reservoirs on timescales of a few hundred Myr (French et al., 2015; Bezanson et al., 2022).

In this work, we set out to test the aforementioned claims, using the SQuIGG$\vec{L}$E  (Studying Quenching in Intermediate-z Galaxies: Gas, angu$\vec{L}$ar momentum, and Evolution, see Suess et al., 2022a) sample, the largest spectroscopic Sample of massive post-starburst galaxies at $z=0.5-0.9$. By combining new and ancillary CO(2-1) and dust continuum observations, quadrupling the total sample size, with rest-optical spectroscopy and rest-optical-to-far-infrared photometry from the Sloan Digital Sky Survey DR14 (SDSS, Abolfathi et al., 2018), the Wide Infrared Survey Explorer (WISE, Wright et al., 2010), and the Herschel Space Observatory (Pilbratt et al., 2010; Poglitsch et al., 2010; Griffin et al., 2010), we perform new spectrophotometric fitting to quantify the level of allowable buried star formation in these systems and to test whether it can explain their over-luminous CO reservoirs. This work is laid out as follows. In Section 2, we present the SQuIGG$\vec{L}$E  sample and describe our new Atacama Large Millimeter Array (ALMA) CO(2-1) and dust continuum observations to assemble the largest sample of molecular gas measurements of post-starburst galaxies to date. In Section 3, we contextualize our CO(2-1) measurements and present new spectrophotometric fits, building on the detailed Prospector (Johnson & Leja, 2017; Leja et al., 2017; Johnson et al., 2021) modeling presented Suess et al. (2022a) via the inclusion of mid- and far-IR constraints. In Section 4, we analyze how these new star formation histories impact our conclusions about the star formation rate, the IR luminosity, and the depletion time of these systems. In Section 5, we discuss the implications of these star formation histories on the depletion time in detail. Finally, in Section 6, we summarize our results. In an accompanying research note (Kumar et al. submitted), we also characterize the serendipitously CO-detected “buddy” galaxies that are likely satellites of the SQuIGG$\vec{L}$E  host galaxies to constrain the environments of these post-starburst galaxies.

Throughout this work, we adopt the best-fit cosmological parameters from the WMAP 9 year results (Hinshaw et al., 2013): $H_{0}=69.32\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$, $\Omega_{m}=0.2865$, and $\Omega_{\Lambda}=0.7135$, utilize a Chabrier initial mass function (Chabrier, 2003), and quote AB magnitudes. For convenience, throughout the work we quote $\alpha_{CO}$ without its units of $M_{\odot}\ \mathrm{(K\ km\ s^{-1}\ pc^{-2})^{-1}}$.

2 Data

2.1 The SQuIGG$\vec{L}$E  Sample

The SQuIGG$\vec{L}$E  Sample was selected from the Sloan Digital Sky Survey DR14 (SDSS, Abolfathi et al., 2018) to constitute a pure sample of post-starburst galaxies at $z=0.5-0.9$ based on their spectroscopic shapes. In particular, they were selected using $U_{m}$, $B_{m}$, and $V_{m}$ rest-frame colors defined in Kriek et al. (2010) to have strong Balmer breaks, but also to have blue slopes redward of the break to select against old or dusty populations. A full description of the sample and detailed rest-optical spectrophotometric fitting can be found in Suess et al. (2022a). In total, the sample consists of 1318 galaxies, with a median redshift of 0.7 and SDSS $i_{AB}=19.5$, corresponding to $\log(M_{\star}/M_{\odot})\sim 11.2$.

Initially, twelve of the most luminous and massive SQuIGG$\vec{L}$E  galaxies (see Figure 1) were targeted for follow-up. Spatially resolved rest-optical spectroscopy of these systems with Gemini/GMOS revealed a range of ordered motion and flat age gradients as traced by $H\delta_{A}$ equivalent width (Hunt et al., 2018; Setton et al., 2020). Additionally, deep ($\sim 2$ hour/galaxy) ALMA Band 4 observations of these systems revealed that many of them are surprisingly luminous in CO(2-1), suggesting significant residual molecular gas despite their apparent quiescence as determined by rest-optical spectrophotometric fitting (Suess et al., 2017; Bezanson et al., 2022). While much of this CO(2-1) appears bound to the galaxies, a significant amount (as high as $\sim 50\%$) of the luminosity is found to emit at distances of 10s of kiloparsecs along tidal tails that are revealed in HST imaging and high resolution ALMA observations (Spilker et al., 2022; D’Onofrio et al., 2025). All 12 of the CO(2-1)-observed galaxies were targeted for rest-optical spectroscopy of $H\alpha$, with weak emission confirming their low star formation rates, albeit with systematic uncertainty in the dust correction due to the inability to correct for Balmer decrement (Zhu et al., 2025).

There is an additional subset of SQuIGG$\vec{L}$E  that is distinguished from the rest of the sample by its overlap with the Hyper-Suprime Cam imaging survey (HSC Aihara et al., 2018), enabling the study of their structures. Setton et al. (2022) performed Sérsic  fitting of these 145 galaxies and found that they lie systematically below the mass-size relation, even relative to other quiescent galaxies, suggesting that their proto-elliptical structure is in place at the time of quenching, and Verrico et al. (2023) showed that they host a very high incidence of tidal features, especially in the youngest systems, suggesting a link between their burst and quenching and major mergers.

2.2 ALMA observations

In this work, we compile a range of Band 4 observations of SQuIGG$\vec{L}$E  galaxies to expand on previous analysis of CO(2-1) and dust continuum measurements, quadrupling our sample relative to Bezanson et al. (2022). Furthermore, this new is far more representative sample of SQuIGG$\vec{L}$E  , as galaxies were selected to span a range of key parameters (provided they overlapped with HST imaging, see Setton et al., 2022), in contrast with the Bezanson et al. (2022) sample that drew from the brightest sources. Here, we describe those programs, all of which were observed in compact configurations as detection experiments with beam sizes ranging between 0.7-2.2*′′*.

The majority of the new data presented in this work comes from Program 2021.1.01535.S and its re-submission, 2022.1.00604.S (PI: D. Setton). These observations, which were taken at various times throughout 2022 and 2023, were part of a survey program that targeted a sample of 31 galaxies spanning a range in stellar mass, star formation rate, and age (as determined by the fits in Suess et al., 2022a), in addition to redshift and rest-optical size (as measured in Setton et al., 2022; Verrico et al., 2023). Observations in this program were targeted for 1 hour, although a small number of targets that were observed between the Cycle 9 submission deadlines and the end of Cycle 8 were targeted twice, resulting in a total integration time closer to two hours.

An additional six galaxies were targeted as a part of Cycle 8 Program 2021.1.00988.S (PI: D. Setton). These galaxies were specifically targeted as young mergers (via their tidal features, see Verrico et al., 2023), similar to those presented in (Spilker et al., 2022; D’Onofrio et al., 2025), with the goal of identifying tidally stripped gas. As such, galaxies in this program were observed for $\sim 2$ hours/galaxy to mimic the depths of the galaxies presented in Suess et al. (2017) and Bezanson et al. (2022) where the spatially-extended CO was first identified. Because these galaxies were specifically targeted to select for CO luminous sources, they were detected at a high rate (5/6) and are among the most CO luminous sources in our sample. Finally, one galaxy was observed as a part of the partially completed Cycle 8 Program 2021.1.00761.S (PI: K. Suess), which targeted AGN candidates identified in Greene et al. (2020) that overlapped with the HSC footprint. Details of the targets and observing strategy are presented in Table 2.

In Figure 1, we show the distribution of these new observations in $\log(M_{\star}/M_{\odot})$ versus redshift as black points, as compared to the full SQuIGG$\vec{L}$E  sample (grey contours) and the (Bezanson et al., 2022) sample (blue diamonds). In contrast with the previous observations, which targeted a subset that was heavily biased toward lower-redshift massive sources, this new sample does a much better job at spanning the full range of SQuIGG$\vec{L}$E  galaxies. While we cannot purport that this sample is “complete” (as the SQuIGG$\vec{L}$E  sample itself is a heavily biased subset of SDSS galaxies, largely drawn from the BOSS CMASS sample, but also spanning a wide range of other selection criteria), it will at least span a wide range in the properties of galaxies that meet the rest-frame selection criteria that Suess et al. (2022a) demonstrated result in a pure sample of galaxies that just underwent a rapid starburst, followed by a precipitous drop in star formation rate. In the next two sections, we describe our reduction of the ALMA data, the synthesis of CO(2-1) and continuum cubes, and our extraction of line and continuum fluxes. All imaging was performed with CASA version 6.6.3-22 (CASA et al., 2022).

2.2.1 CO maps and flux measurements

We image all galaxies in the channel centered on CO(2-1) with 0.08*′′* cells, using natural weighting and combining all available UV data. We utilize the tclean algorithm with the auto-multithresh algorithm, setting pblimit=-0.1 and nsigma=3. For each galaxy, we generate maps at 50, 100, and 200 km/s, centered at the systemic velocity of our galaxy. We then extract 1D spectra using a 2*′′* aperture at the location of our galaxy. We measure the total CO flux by integrating at $\pm$500 km/s, estimating uncertainties using the rms of the spectrum outside this region. All measurements are first made on the 50 km/s resolution cube; if the galaxy is not detected at the 3$\sigma$ level, we then iteratively use the 100 km/s and 200 km/s cubes, marking a galaxy as a non-detection if it is not significantly detected in any of the 3 cubes. The measured CO flux densities and luminosities are presented in Table 3. In Figure 2 and 3, we show the Moment 0 maps in addition to the collapsed CO(2-1) spectra, with the velocity channel used for integration indicated on the spectra.

For illustrative purposes, we also generate Moment 1 maps for our four most CO-luminous sources, which we obtain by performing a velocity-weighted average over the CO luminosity within $\pm 500$ km/s. All four of these galaxies show some degree of velocity gradients, with disrupted kinematics that often mirror the tidal features seen in the rest-optical HSC imaging. However, as our resolution is low and the sources are only marginally resolved, we defer any detailed kinematic and structural analysis to future work. In Figure 3, we show the rest of the new galaxies observed in this program, restricting to only the HSC image, CO(2-1) moment 0 map, and extracted CO spectrum.

2.2.2 2 mm continuum measurements

We image the continuum with identical tclean parameters to the previous section, using the three spectral windows that do not contain CO(2-1) emission and specmode=‘mfs’. We also measure fluxes within a $2^{\prime\prime}$ aperture, estimating uncertainties with the RMS of the image under the assumption that our sources are unresolved. In Figure 2, we show the 2mm continuum images for our four most luminous sources from the new sample, all of which are detected. However, only 11/50 galaxies in the full SQuIGG$\vec{L}$E  galaxies are continuum detected at the 3$\sigma$ level. The measured continuum flux densities are presented in Table 3.

2.2.3 Herschel photometry

In order to further supplement our FIR constraints near the peak of the dust SED, we search for archival Herschel Space Observatory (Pilbratt et al., 2010; Poglitsch et al., 2010; Griffin et al., 2010) photometry of the SQuIGG$\vec{L}$E  targets in our survey. Specifically, we query both the PACS and SPIRE point source catalogs (Marton et al., 2017; Schulz et al., 2017) at the locations of our sources, searching all bands for matches where the centroid in the point source catalog is within 3*′′* of our source. There are no PACS 70$\mu$m, 100$\mu$m, or 160$\mu$m matches for any of our sources that do not also appear in the Rejected Source List. Similarly, no sources appear in the SPIRE 350 $\mu$m or 500 $\mu$m catalogs. However, we find three sources (J1157+0132, J0910+0218, and J1142+0006) with matches in the 250 $\mu$m catalog, all only barely detected at 3-4$\sigma$ significance. We include this photometry in our analysis.

While we do not include non-detections in our analysis (as non-detections at the typical depth of Herschel are not particularly constraining for sources in our redshift range), we do note that while we only found that 3 galaxies were securely detected, 22/50 sources are either detected or have an entry in the point source catalog within 60*′′* (suggesting that they are in the archival imaging footprint of Herschel). As such, we conclude that while a few exceptional sources were detected, the lack of matches is not primarily due to the lack of coverage, and instead it seems that the typical SQuIGG$\vec{L}$E  galaxy is not IR-luminous enough to be detected in Herschel imaging.

3 Analysis

In this section, we seek to contextualize the observed CO properties of the SQuIGG$\vec{L}$E  post-starburst sample within the evolutionary trajectory that can transform star forming galaxies into red-and-dead quiescent galaxies. Doing so relies heavily on our understanding of the stellar populations of these galaxies. Suess et al. (2022a) presented detailed spectrophotometric fits to the rest-optical SED of the full SQuIGG$\vec{L}$E  sample, finding that these sources are characteristically massive ($\log(M_{\star}/M_{\odot})\sim 11.2$) and quiescent (log(sSFR $[\mathrm{yr^{-1}}])\lesssim-11$). In the first sections, we utilize these fits, in conjunction with our new and previously published (Suess et al., 2017; Bezanson et al., 2022) CO(2-1) measurements, to place these sources into context under the standard set of assumptions in that previous work. Then, we shift our focus to the stellar populations of SQuIGG$\vec{L}$E  themselves, addressing how changing those assumptions and incorporating data from the far-IR SED might change the inferences about the quiescence of these systems.

3.1 The molecular content of post-starburst galaxies

First, we evaluate how our the molecular gas mass fraction of the SQuIGG$\vec{L}$E  sample compares to the expectation for a population that has finished its primary epoch of star formation and is entering quiescence. To do so, in Figure 4 we compare our sample to scaling relations for the time evolution of the molecular gas fraction for star forming massive ($\log(M_{\star}/M_{\odot})=11$) galaxies from Tacconi et al. (2018). The molecular gas fraction evolves with cosmic time, as the entire population of galaxies has become less gas-rich and less star forming as the universe has evolved in the wake of cosmic noon. However, across cosmic time, massive quiescent galaxies (shown in red) lie well below this relation, indicating that a large part of why they are no longer star forming is because they no longer host sufficient fuel for substantial new star formation.

Post-starburst galaxies (shown in green), however, which potentially represent the evolutionary link between these populations, show much more mixed properties. While many post-starburst galaxies, especially at high-z, are not detected in CO or dust continuum, finding limits that are consistent with the gas-poor quiescent population, many sources are detected with gas fractions that place them on (or at low-z, above, predominantly in the SPOG sample, see Alatalo et al. 2016a) the star forming main sequence at their redshift. Our compiled SQuIGG$\vec{L}$E  galaxies reflect this trend near cosmic noon. Under the standard assumption of $\alpha_{CO}=4.0$ (Bolatto et al., 2013) and $r_{21}=1.0$, we find that our most strongly detected galaxies have molecular gas fractions as high as 30%, while non-detections are constrained to have fractions as low as $0.3\%$. This finding reflects a broad consensus in the literature that many optically selected post-starburst galaxies at and below cosmic noon can host molecular gas fractions that are typical, if slightly low, relative to co-eval star forming galaxies (Suess et al., 2017; Spilker et al., 2018; Bezanson et al., 2022; Zanella et al., 2023; Wu et al., 2023). However, there also exists a significant population of post-starburst systems, largely probed by SQuIGG$\vec{L}$E, that fall significantly below the main sequence, with molecular gas fractions that are commensurate with the limits from older high-z quiescent systems (e.g., Sargent et al., 2015; Magdis et al., 2021; Whitaker et al., 2021a; Williams et al., 2021).

In the context of SQuIGG$\vec{L}$E, the high molecular gas fractions found in $\sim$half our sources is surprising given that these galaxies are, by all rest-optical tracers, quiescent (Suess et al., 2017, 2022a; Bezanson et al., 2022; Zhu et al., 2025). To illustrate this, in Figure 5, we plot the star formation rate versus the molecular gas mass, with literature scaling relations illustrated as contours (Saintonge et al., 2011; Tacconi et al., 2018; Freundlich et al., 2019). The star formation rates of SQuIGG$\vec{L}$E  galaxies measured in Suess et al. (2022a) place the full sample below this relation (black points, see also Bezanson et al., 2022). However, in contrast, literature post-starburst samples at similar redshift (Belli et al., 2021; Woodrum et al., 2022; Wu et al., 2023, shown as green diamonds) all appear to have star formation rates that are largely within the scatter on star forming systems with similar $M_{H_{2}}$.

A key difference may lie in the way the star formation rates were inferred in these samples. All the aforementioned works measured their star formation rates using SED fits that incorporated rest-MIR SED information via Spitzer/MIPS 24 $\mu$m photometry, enforcing that the IR luminosity from attenuated stellar populations could produce that observation. In contrast, the Suess et al. (2022a) fits only utilized data from the rest-optical SED out to rest-frame $\sim 2$ $\mu$m, and the fits to the SQuIGG$\vec{L}$E  sample presented in that work were shown to under-predict WISE W3 and W4 photometry in sources that were detected. As such, it is possible that much of the tension between the inferred rest-optical star formation rates and the high molecular gas mass can be resolved by accounting for buried star formation that does not contribute in the rest-optical but produces significant IR luminosity. In the following section, we update the SED fits for the full SQuIGG$\vec{L}$E-CO(2-1) sample, addressing the IR SED in addition to a number of other modifications to test whether the previous fits were under-estimating the star formation rate of SQuIGG$\vec{L}$E  post-starburst galaxies.

3.2 Updated panchromatic spectrophotometric fitting

The SED fits presented in Suess et al. (2022a), using Prospector (Johnson & Leja, 2017; Leja et al., 2017; Johnson et al., 2021), did not attempt to fit redward of W2, and therefore remained agnostic to the information contained at mid- and far-IR wavelengths. Prospector handles the dust SED by assuming energy balance, where the luminosity that is attenuated by dust is re-radiated in the mid- and far-IR with an SED that is described by the Draine et al. (2007) templates. The Suess et al. (2022a) fits acknowledged that the mid-IR SED of WISE-detected SQuIGG$\vec{L}$E  galaxies could not be produced by the a set of fixed Draine et al. (2007) parameters they chose. However, the choice to fix these parameters, and a number of other choices related to the treatment of dust in the fitting, make the Suess et al. (2022a) fits ill-suited for answering questions about the presence of buried star formation, which has been suggested as a possible explanation for the gas-rich nature of many post-starburst galaxies (e.g., Baron et al., 2023).

For example, the Suess et al. (2022a) fits fixed the PAH-mass-fraction, $q_{PAH}$, to 2%, limiting the ability of over-luminous PAHs that have been seen in post-starburst systems (see Smercina et al., 2018) to boost the mid-IR by shifting energy from cold dust to PAHs. The Suess et al. (2022a) fits similarly fixed the $U_{min}$ and $\gamma$ parameters that govern the shape of the dust SED to values typical of very cold dust, leaving no freedom for deviation from this distribution toward hotter dust. Additionally, the Suess et al. (2022a) fits did not include a mid-IR AGN torus, which can be included in Prospector using the Nenkova et al. (2008) torus templates (Leja et al., 2018), with normalization that is independent of the energy balance of the SEDs, providing a potential boost to the mid-IR flux without requiring any additional source of attenuated luminosity.

Finally, the Suess et al. (2022a) fits adopt a two-screen Charlot & Fall (2000) dust model, where all stars see dust with optical depth $\tau_{ISM}$, and young ($t<10^{7}$ yr) stars see an additional “birth cloud” dust with optical depth $\tau_{BC}$. However, motivated by Wild et al. (2016) and by observations that show a relation between stellar- and Balmer decrement- $A_{V}$ (Price et al., 2014), tracing these two quantities, the Suess et al. (2022a) fits fixed $\tau_{BC}=\tau_{ISM}$. This assumption is quite conservative relative to other SED codes where the two parameters are left independent of one another (for example, MAGPHYS allows the terms to vary independently and places a fairly weak prior on the total attenuation, see da Cunha et al., 2008). Because the SQuIGG$\vec{L}$E  post-starburst galaxies are by definition selected to look like low-dust A-star populations that lack emission, this prior effectively did not allow for any star formation to be buried behind optically thick dust.

Given the previous SED fitting, it is not actually clear on a per-galaxy basis: 1) if any additional star formation is necessary to explain the mid-IR luminosity, given that young stellar populations subjected to modest dust attenuation can still produce significant IR luminosity (e.g., Wild et al., 2025), and, 2) how much star formation the mid-IR detections (or upper limits) can tolerate given that there must be significant IR contribution from dust heated by attenuated A stars. Here, we address those questions head on by fully folding our full set of constraints on the IR spectral energy distribution into our analysis and re-fitting the SQuIGG$\vec{L}$E-CO sample with IR dust templates under the assumption of energy balance. We make the following key modifications from the Prospector fitting procedure outlined in Suess et al. (2022a):

We include WISE W3 and W4 photometry presented (but excluded from the fits) in Suess et al. (2022a). We exclude SDSS u-band imaging, as most of our sources are at the edge of the SDSS detection limit ($u_{AB}\gtrsim 22$), and we have found that the uncertainties on supposed detections appear to be underestimated when inspecting these images. The exclusion of this photometry does not affect our fits. In total, our fits include the SDSS $griz$ photometry and WISE W1/W2/W3/W4 photometry, spanning $\lambda_{\mathrm{rest}}=0.3-13$ $\mu$m. Additionally, for the 3 galaxies with photometry in the Herschel/SPIRE point source catalog at 250 $\mu$m ($\lambda_{rest}\sim 150$ $\mu$m), we include those measurements in our fits. 2. 2.

The Draine et al. (2007) parameters that govern the reprocessed attenuated luminosity are left free, with $q_{PAH}$ (the PAH mass fraction) $\in[0.1,10]$, and $\gamma_{e}$ and $U_{min}$ (which, broadly speaking, control the temperature of the dust SED) $\in[0.0001,1]$ and $[0.1,25]$ respectively. We sample the PAH mass fraction linearly (favoring the high PAH luminosity contribution seen in local post-starburst galaxies, see Smercina et al., 2018), and sample in log space for $U_{min}$ and $\gamma_{e}$ to allow for maximally flexible dust SEDs (as the dust SED changes logarithmically with these parameters). This model is scaled to the total attenuated luminosity for a given star formation history/dust law draw during the fitting. 3. 3.

The Nenkova et al. (2008) AGN torus templates are also included in the fit, with a free AGN luminosity (parameterized as a fraction of the total bolometric luminosity of the galaxy) $\log(f_{\mathrm{AGN}})\in[-3,0.48]$ and $\log(\tau)\in[0.7,2.2]$ following Leja et al. (2018). The luminosity of this AGN torus is not tied to any observed properties of the galaxy (e.g., rest-optical AGN indicators), although we note that as many as 20% of the youngest SQuIGG$\vec{L}$E  galaxies may host an AGN, as traced by the [OIII]/H$\beta$ line ratio (Greene et al., 2020). 4. 4.

We modify the Suess et al. (2022b) post-starburst star formation history, adding one final bin that covers exactly the most recent 10 Myr of star formation. We do so because the two-screen Charlot & Fall (2000) dust model assumed in fsps by default allows young ($t<10$ Myr) stars to see an additional “birth cloud” dust screen that is added to the “ISM” dust screen (Conroy et al., 2009). Forcing this bin of star formation to be independent from the rest of the star formation history allows dust-obscured star formation to be independent from the rest of the star formation in the final bin. The total length of time of the final two bins of star formation, $t_{last}$, is fit as a free parameter drawn from a flat (in log-space) distribution with a minimum time of 20 Myr and a maximum time of 15% of the age of the universe at the redshift of observation. The final two bins have lengths of 10 Myr and $t_{last}-10$ Myr, respectively. 5. 5.

Finally, a number of priors from the Suess et al. (2022a) fits are relaxed. First, we do not include a prior on the star formation history based on the average star formation history of a mass- and redshift-matched galaxy in the Universe Machine (Behroozi et al., 2019); instead, we place a Student’s t prior on the logarithmic ratio of star formation rate between neighboring bins (with $\sigma=0.3$ and $\nu=1$ for all but the final two bins, which have $\sigma=0.5$ to allow for a more rapid rise or fall). We also relax the mass-metallicity relation prior (based on the Gallazzi et al., 2005) relation that was used in Suess et al. (2022a), allowing the metallicity to vary freely with $\log(Z/Z_{\odot})\in[-0.4,0.19]$.

We perform two sets of fits, meant to test two assumptions about the relationship between the ISM dust and the birth cloud dust. In the first model, we mimic the Wild et al. (2020) (motivated by Price et al., 2014) prior choice that $\tau_{BC}\sim\tau_{ISM}$. However, in contrast with Suess et al. (2022a), rather than simply fixing the parameters to have the same value, we instead place a Gaussian prior on the fraction between these two parameters, with $\mu=1$ and $\sigma=0.25$, leaving $\tau_{ISM}$ free $\in[0.0,2.5]$ with a flat prior. In a second set of fits, we allow for significantly more freedom in $\tau_{BC}$ to allow for buried star formation. To do so, we adopt a clipped Gaussian prior on $\tau_{BC}$, with $\mu=0.5$ and $\sigma=2.5$ (though it is not the same functional form, this prior is very similar to the one adopted in MAGPHYS, see da Cunha et al., 2008). The shape of the dust law for the ISM dust is parameterized with a dust index as in Kriek & Conroy (2013), while the birth cloud dust law shape is fixed as a power law with a slope of -1. Throughout this work, we refer to these fits as “No Buried SF” and “Buried SF Allowed”, respectively.

Other than the aforementioned changes, our remaining modeling choices are largely similar to the Suess et al. (2022a) fits. We enforce a signal-to-noise ceiling of 20 on all photometry We use the dynesty dynamic nested sampling package (Speagle, 2020), the Flexible Stellar Population Synthesis (FSPS) stellar population synthesis models (Conroy et al., 2009; Conroy & Gunn, 2010), the MILES spectral library (Sánchez-Blázquez et al., 2006; Falcón-Barroso et al., 2011), and the MIST isochrones (Choi et al., 2016; Dotter, 2016). We assume a Chabrier (2003) Initial Mass Function and fix the model redshift to the spectroscopic redshift as measured by SDSS. Before fitting, we convolve all models to the SDSS spectroscopic resolution, and we fit with an additional velocity dispersion term $\in[0,400]$ km/s. Nebular continuum and emission is turned on, and the nebular ionization parameter $\log(U)$ is left free $\in[-4,-1]$. However, we mask the [OII] and [OIII] lines as they are highly susceptible to ionization by shocks and AGN, both of which are not accounted for in our modeling framework. Finally, we utilize the PolySpecModel procedure with a 2nd order polynomial, and we include both a spectroscopic jitter term and the Prospector pixel outlier model. In Table 1, we show some of the derived parameters for two galaxies under both sets of fits; full tables of the entire SQuIGG$\vec{L}$E  sample are available online.

4 Results

4.1 Star formation histories

In the previous section, we performed two sets of fits to the panchromatic SED and spectra of SQuIGG$\vec{L}$E  post-starburst galaxies, one which used a commonly adopted prior where the ISM dust and the birth cloud dust are tied together (No Buried SF) and one which relaxed that prior to allow significantly more freedom in the birth cloud dust (Buried SF Allowed). Here, we explore the different conclusions about the star formation histories of post-starburst galaxies from the application of those priors.

In Figures 6, 7, and 8, we present example SED fits to three galaxies. In each figure, we first show the observed SDSS spectrum, with the best fitting model spectra for the No Buried SF (pink) and Buried SF Allowed (grey) models shown as dashed lines. In the bottom two rows, we show the median fitting (with 68% confidence interval) star formation history (middle, showing only the 25% of the galaxy’s star formation history probed by the flexible bins and omitting earlier star formation for clarity) and SED (bottom, including photometric residuals normalized by the uncertainty) for each of the models. In the SED fit, we break down the spectrum into three components, active star formation (red) which shows the stellar and dust contributions from stellar populations with $t<10$ Myr, the contribution from older stars that formed $t>10$ Myr ago (green), and the mid-IR Nenkova et al. (2008) CLUMPY torus (blue, though in all of these sources, an AGN torus is not needed to fit the mid-IR).

In Figure 6, we highlight J1157+0132, our most CO-luminous and Herschel-luminous source. For this source, the two sets of priors yield remarkably different star formation histories. The fit which does not favor buried star formation essentially reproduces the class of models presented in Suess et al. (2022a), where the star formation histories of SQuIGG$\vec{L}$E  galaxies indicate a large burst with SFR $\sim 200\ M_{\odot}\ \mathrm{yr^{-1}}$ and a sharp decline, in this case to a modest star formation rate (in the final 10 Myr bin) of $18\pm^{10}_{6}$ $M_{\odot}$/yr ($\log(\mathrm{sSFR\ [yr^{-1}]})=-10.16\pm^{0.19}_{0.19}$). However, in contrast with the Suess et al. (2022a) fits that did not allow for any freedom in the shape of the FIR SED, our more flexible dust modeling has no problems producing a mid-IR SED that is in agreement with the data without the need for any additional buried star formation. The star formation rate constrained in this fit is quite similar to the original Suess et al. (2022a) fits, indicating that the problem producing the mid-IR photometry was not due to a lack of sufficient IR luminosity. Instead, it was simply the rigidity of the dust SED assumptions (for example, the Suess et al. (2022a) fits fixed $q_{PAH}=2$, whereas the fit we present here requires $q_{PAH}=6.4\pm^{2.0}_{2.0}$) that drove the inability to produce the red IR continuum. However, while the predicted IR luminosity from the attenuation of the stars that formed in the $\sim 100$ Myr old burst is nearly enough to account for all the mid-IR luminosity, this model undershoots the Herschel/PACs detection by a factor of 2.7, indicating that additional IR luminosity is needed to fully explain the FIR.

In our Buried SF Allowed modeling framework, this additional IR luminosity is provided by additional star formation behind optically thick dust, which limits its influence in the rest optical. The Buried SF Allowed fits can tolerate significantly higher star formation rates, with a highly uncertain best fitting SFR of $76\pm^{26}_{24}$ $M_{\odot}$/yr. In this class of solutions, the active star formation contributes comparably to the IR luminosity relative to the older stars, with the optically thick dust mirroring a sharp drop in star formation that can produce the deep Balmer absorption seen in the spectrum. Interestingly, despite a 0.4 dex discrepancy between the model and observed 250 $\mu$m flux in the No Buried SF model, the Buried SF Allowed model is able to reproduce the full SED with only a 0.14 dex shift in the total IR luminosity. This IR luminosity boost can be achieved even with star formation rates that are contained in the posterior of the No Buried SF framework, as evidenced by the overlapping posteriors. As such, while we conclude that in our modeling framework, buried star formation is required to produce this shift, the degree of buried star formation required is currently highly uncertain. J1157+0132 could indeed still be in the midst of its starburst, with optically thick dust attenuating its emission, but it could also be in a period of rapid decline from a recent period of much higher star formation rate during a $\sim 100$ Myr old burst.

In Figure 7, we highlight J0910+0218, our second most CO-luminous and Herschel-luminous source. In contrast with J1157+0132, here both sets of fits do an equally good job at describing the mid-IR SED, as evidenced by the practically identical constraints on for the IR luminosity between the two models despite the allowance for a factor of $\sim 2\times$ more star formation in the Buried SF Allowed fit. As such, in this source, which which is far more representative of our CO-detected sample than J1157+0132, we definitively cannot state that the data itself exhibits any strong preference for any level of buried star formation; the luminosity from the attenuated $\sim$100 Myr stellar populations that dominate the rest-optical continuum is more than enough to produce the observed FIR SED.

Finally, in Figure 8, we highlight J1141+0109, a galaxy which is not detected in any IR band or in CO(2-1). In contrast with J1157+0132 and J0910+0218, the two sets of fits to J1141-0109 return essentially identical star formation histories and physical properties. Because the galaxy is not detected in the infrared, the combination of the W3/W4 non-detections place a hard cap on the maximum IR luminosity that is tolerable in the fits. As such, there is very little room to add an appreciable contribution from buried star formation even when the priors are more permissive to hiding emission in the rest-optical. This is reflected in the overlapping posteriors for the star formation rate and the IR luminosity, even when the best fitting $\tau_{BC}$ is significantly greater than $\tau_{ISM}$.

In Figure 9, we compare the full set of galaxies in the inferred total IR luminosity (defined as the sum of the luminosity of all dust components included in our fits, including an AGN torus). The vast majority of the galaxies in SQuIGG$\vec{L}$E  show essentially no deviation in their total IR luminosity between the two fits, indicating not only that producing the observed IR luminosity is possible with fits that do not require buried star formation, but that the attenuated luminosity of the moderately dusty A-type populations required to fit these post-starburst spectra is already so high that there is little room in the dust SED to hide a contribution from a starburst that is invisible in the rest-optical. As such, the median shift in $\log(L_{IR})$ from the No Buried SF fit to the Buried SF Allowed fit is only 0.02 dex, with a scatter that encompasses 0, and even in the most extreme cases, J1157+0132 (the most IR luminous galaxy in the sample), there is not room for more than a factor of $\sim 2$ more IR luminosity than is implied in the fits without buried star formation.

This manifests in star formation rates that are boosted only moderately in the Buried SF Allowed fits, with a median shift of 0.5 dex in sources with a star formation rate $>1M_{\odot}\ \mathrm{yr^{-1}}$ where constraints from this kind of fitting are meaningful (see Suess et al., 2022b). This factor of $\sim 3$ shift, even under the assumption that there is a maximal amount of star formation behind optically thick dust, is quite small given the freedom of the fits, and in all cases, the posterior for the star formation rate measured in the Buried SF Allowed fit overlaps significantly with posterior for the the star formation rate measured under the No Buried SF priors.

We caution that the tight constraints on the IR luminosity in these fits, and subsequently, the star formation rate, in many ways rely on the validity of the assumption that the Draine et al. (2007) dust templates span the full range of possible dust temperature and grain size distributions such that our extrapolations from our rest-frame 10 $\mu$m (and, occasionally, a single point at $\sim$150 $\mu$m) observations are constraining. Dust that emits with a distribution that maximizes output at warm temperatures where we lack constraints could allow for more IR luminosity, and therefore more star formation. Future work, ideally with better FIR data from an observing facility like PRIMA (Moullet et al., 2023) or with more robust dust-obscured star formation rate tracers, could better allow for these degeneracies to be broken. For now, we proceed forward under the assumption that our two sets of measured star formation histories span a reasonable range of the possible present day star formation activity in SQuIGG$\vec{L}$E  post-starburst galaxies. In the next section, we investigate how the changes in the inferred star formation rate would affect the interpretation of where these galaxies lie relative to the star forming main sequence.

4.2 The star forming main sequence

Many works have found that many post-starburst galaxies lie systematically below the star forming main sequence and the SFR-$M_{H_{2}}$ relation for normal star forming galaxies for a variety of rest-optical SFR tracers (e.g., French et al., 2015; Suess et al., 2017; Belli et al., 2021; Bezanson et al., 2022), though with considerable disagreement between different star formation rate tracers. Baron et al. (2023) argues that the high IR luminosities of post-starburst galaxies imply star formation rates that are well above the main sequence, resolving all tension related to the depletion time but requiring that significant star formation is occuring behind optically thick dust. However, it has also been argued that direct conversions from the IR luminosity to star formation rate in these systems can significantly overestimate the true instantaneous star formation rate, as heating from more evolved stars (especially young A type stars formed on timescales of $\sim 100$ Myr before observation) can contribute significantly to the IR luminosity (Utomo et al., 2014; Hayward et al., 2014; Leja et al., 2019; Wild et al., 2025).

In the previous section, we showed that we can measure instantaneous (10 Myr) star formation rates under two sets of assumptions about the level of allowed buried star formation that simultaneously produce the UV and the IR for our galaxies. In Figure 10, we compare the star formation rates we infer under these two sets of assumptions against the star forming main sequence at $z=0.7$ from Leja et al. (2022), with the No Buried SF fit shown in pink on the left and the Buried SF Allowed fit in grey on the right. For clarity, galaxies with a median SFR $<0.1$ $M_{\odot}$/yr are represented by their $3\sigma$ limits.

As stated previously, the median shift in the measured SFR under a much more permissive prior for the dust around birth clouds is only $0.5$ dex. While this is enough to boost some galaxies onto the star forming main sequence (SFMS), especially when upper limits are considered, the vast majority of galaxies are still consistent with lying below this relation, and only one galaxy (J1157+0132, see Figure 6) is able to hide enough IR luminosity that it can be consistent with being in the midst of an active starburst. As such, we conclude that the SQuIGG$\vec{L}$E  post-starburst galaxies can at most be hosting star formation rates that are consistent with lying on the Leja et al. (2022) main sequence, which is notably lower than other published sequences such as the Whitaker et al. (2012b) relation (also shown, in purple) at high masses due to the accounting for the heating from old stars in the conversion from IR luminosity to the star formation rate.

Even with buried star formation allowed, the SQuIGG$\vec{L}$E  galaxies still lie below the Whitaker et al. (2012b) relation, which assumes that the entirety of the IR luminosity is produced by active star formation. This suggests that while SQuIGG$\vec{L}$E  galaxies may still be forming star stars at moderate levels, they are not still in the midst of active starbursts, and an appreciable part of their IR luminosity must come from evolved stellar populations. And, once again, the data provides no strong evidence for the higher star formation rates measured with more permissive priors; moderate star formation rates simply cannot be ruled out with existing data.

4.3 The molecular gas content as a function of age

Previous analysis of SQuIGG$\vec{L}$E  CO(2-1) observations in Bezanson et al. (2022) found a clear bi-modality in the observed properties of galaxies, with the youngest (time since quenching $\lesssim$200 Myr, see Suess et al., 2022a) exhibiting highly luminous CO and the older galaxies remaining undetected. Similar age trends are seen in French et al. (2018) as a function of the post-burst age.

Here, we further investigate that trend with our full sample. Rather than using the time since quenching, which requires a sharp drop in measured star formation rate (which does not occur in all acceptable models under the Buried SF Allowed priors, see Figure 6), we instead define $t_{\mathrm{PSB}}$ as the lookback time at which a galaxy formed 99% of its stars in the Gyr before observation. This quantity has the advantage of being largely insensitive to the choice of prior; the median scatter in $t_{\mathrm{PSB}}$ between the two sets of priors is 0, with a scatter of 0.1 dex. In Appendix B, we show the comparison between this quantity and the “time since quenching” measured in Suess et al. (2022a), finding very good agreement.

In Figure 11, we show the trend between the molecular gas mass and molecular gas fraction (assuming $r_{21}=1.0$ and $\alpha_{CO}=4.0$) and $t_{\mathrm{PSB}}$. We find that the bimodality presented in Bezanson et al. (2022) holds for our larger and more representative SQuIGG$\vec{L}$E  sample; essentially all CO(2-1) luminous galaxies have $t_{\mathrm{PSB}}<200$ Myr, and a stack of all non-detections is only detected at the 2$\sigma$ level with 1-2 dex lower CO(2-1) luminosity than our detected galaxies. We fit an exponential decay to the gas fraction as a function of age, estimating uncertainty by running 1000 fits to draws from the stellar mass and age posteriors for the ensemble of galaxies and incorporating the uncertainty in the gas mass in the gas fraction at each stage. We find a characteristic timescale of decay of $70\pm 20$ Myr. In Appendix B, we show that this qualitative trend is robust to different definitions of age, and the decay timescale is robust within a factor of 2. If the SQuIGG$\vec{L}$E  sample is indeed tracing a single population drawn from a distribution of post-burst age, this implies galaxies must fully deplete their gas content within $\sim 250$ Myr of quenching, and our detected galaxies must deplete their remaining fuel within $\sim 140$ Myr (two e-folding times). In the following section, we explore the possible physical scenarios that might explain this trend.

5 Why does the CO(2-1) luminosity rapidly fade in the wake of a starburst?

In the previous section, we demonstrated that the previously identified trend between CO luminosity and age wherein essentially all CO-luminous post-starburst galaxies have ages $\lesssim 200$ Myr (French et al., 2018; Bezanson et al., 2022) holds for the full SQuIGG$\vec{L}$E  sample. This, along with the low luminosity of stacked non-detections, indicates that if the trend in Figure 11 is taken at face value and these galaxies represent an evolutionary sequence, that the CO luminosity must drop by an order of magnitude on a timescale of $\lesssim 140$ Myr.

Here, we confront this trend by asking whether the depletion times of SQuIGG$\vec{L}$E  galaxies (the time it would take to remove their gas reservoirs if they maintain their current star formation rate) can be commensurate with this implied $\sim 140$ Myr timescale. Essentially all studies of the molecular gas content in post-starburst galaxies bake two key assumptions into their analysis (e.g., French et al., 2015; Suess et al., 2017; Smercina et al., 2018; French et al., 2018; Belli et al., 2021; Bezanson et al., 2022; Woodrum et al., 2022; Zanella et al., 2023):

There is no significant buried star formation, and the low attenuation measured from stellar emission can be extrapolated using a relation similar to Price et al. (2014) to infer the attenuation around birth clouds. 2. 2.

The Milky Way CO-to-$H_{2}$ conversion holds in these systems.

Using the resolved Kennicutt-Schmidt relation (SFR vs $M_{H_{2}}$, Schmidt, 1959; Kennicutt, 1998), in this section we systematically test how varying these assumptions effects the measured depletion time ($M_{H_{2}}$/SFR) of our systems, which should be similar to 100 Myr if SQuIGG$\vec{L}$E  galaxies can passively evolve along the $L_{CO(2-1)}$ versus age relation. In Section 5.1, we first present our fiducial model under the standard assumptions of no buried star formation and a Milky Way $\alpha_{CO}$. In Section 5.2, we explore the effect of buried star formation by using the inferred star formation rates from the Buried SF Allowed fits. In Section 5.3, we explore the effect of relaxing the assumption of a Milky Way $\alpha_{CO}$, and also discuss the potential effects of merger-driven tidal stripping and/or outflows. In Section 5.4, we discuss whether there is evidence in our data or in the literature to support either of these modifications to the standard assumptions about the state of the post-starburst ISM. Finally, in Section 5.5, we discuss the possibility that the trend in $L_{CO(2-1)}$ versus $M_{H_{2}}$ is not a trend at all and that SQuIGG$\vec{L}$E  post-starburst galaxies may rejuvenate before re-entering the quiescent sequence.

5.1 Can SQuIGG$\vec{L}$E  galaxies deplete their gas under standard assumptions?

We begin by exploring the standard assumptions of a Milky Way-like $\alpha_{CO}$ and star formation rates that are inferred from the rest optical under the assumption that there is no additional component behind optically thick dust. In the top left panel of Figure 12, we show the SQuIGG$\vec{L}$E  sample under the typical No Buried Star Formation and $\alpha_{CO}=4$ assumptions, replicating the finding from Bezanson et al. (2022) that these galaxies lie below the relation with depletion times that are much larger than the characteristic CO-fading timescale and that post-starburst galaxies are offset from the population of typical star forming galaxies (Saintonge et al., 2011; Tacconi et al., 2018; Freundlich et al., 2019). It is not trivial that the entire sample still lies below the resolved Kennicutt Schmidt relation in the current fits, since unlike Suess et al. (2022a), we explicitly include the mid- and far- infrared fluxes that are luminous in some of sources. In other words, there is sufficient heating from A-type stars to satisfy the long wavelength constraints in our No Buried Star Formation model without adding significant star formation under the assumption that $\tau_{ISM}\sim\tau_{BC}$, though we cannot rule out buried star formation from this fact alone.

Given that we recover the same location in the Keniccut Schmidt relation for SQuIGG$\vec{L}$E  galaxies as Bezanson et al. (2022), we also reach the same conclusion with regard to their depletion times. Under standard assumptions, we measure depletion times $>1$ Gyr, with the bulk of the population lying closer to $\sim 5$ Gyr, far longer than the $\sim 140$ Myr timescale necessary to explain the falling CO luminosity. Therefore, if SQuIGG$\vec{L}$E  galaxies do indeed evolve from CO luminous to CO faint on the observed timescale, there must be some modification to either their inferred star formation rate or the total gas mass consumption that is necessary in that time period. In the following sections, we explore both those possibilities.

5.2 Buried star formation?

We now confront the possibility of deeply buried star formation in our systems. As we demonstrated in Section 4, there is some latitude to raise the star formation rates of SQuIGG$\vec{L}$E  post-starburst galaxies by relaxing assumptions about the relationship between the ambient ISM dust and the birth cloud dust, with the potential to raise a typical galaxy’s star formation rate by $\sim 0.5$ dex. We stress again that this preference is not strongly favored by the data, as the vast majority of our galaxies are not detected in the mid-IR, and even those that are detected in W3/W4 can be well-modeled under the No Buried SF priors. Additionally, studies of low-z post-starburst galaxies have found good agreement between observed $H\alpha$ star formation rates and IR tracers like (Smercina et al., 2018; Luo et al., 2022; French et al., 2023), providing evidence against dust obscured star formation beyond what is expected from standard prescriptions for the dust geometry. However, we cannot rule out the possibility that this buried star formation truly is there, resulting in a significant reduction in the depletion time.

In Figure 12b, we again show the star formation rate versus the molecular gas mass, this time using the star formation rate measured from the Buried SF Allowed fits (and indicating the typical shift from the previous panel with a bold grey vector). Perhaps predictably, a factor of $\sim 3$ boost in the typical star formation rate is not enough to resolve the order-of-magnitude tension between the observed and predicted CO depletion times, and the majority of the sources are consistent with depletion times closer to $\sim 1$ Gyr than $\sim 140$ Myr, especially in the most CO-luminous sources.

Therefore, we conclude that even if there is significant star formation occurring behind optically thick dust in our galaxies, the implied star formation still could not deplete the gas quickly enough. Depletion of $\sim$a few $\times 10^{10}$ $M_{\odot}$ of molecular gas in $\sim 100$ Myr requires star formation rates $\gtrsim 100$ $M_{\odot}$/yr. Such high star formation rates cannot coexist with the strong Balmer absorption, the lack of strong emission lines, and the mid- and far-IR limits on the total IR luminosity within our modeling framework.

5.3 Systematically lower $H_{2}$ masses?

After testing the potential impact of buried star formation, we now turn to the other lever arm in our depletion time calculation–the total $H_{2}$ mass that must form into stars to explain the rapidly dropping CO luminosity. Many works have shown a preference for mergers among post-starburst samples (e.g., Alatalo et al., 2016a; Sazonova et al., 2021; Wilkinson et al., 2022; Ellison et al., 2024), and that same preference manifests in a very high merger fraction within SQuIGG$\vec{L}$E  ($\sim 70\%$ in our youngest galaxies, see Verrico et al., 2023). Gas rich mergers have the potential to systematically lower the CO-to-$H_{2}$ conversion by allowing more CO luminosity per-unit-$H_{2}$ to escape a turbulent, more diffuse post-merger ISM (e.g., Narayanan et al., 2011; Bournaud et al., 2015; Renaud et al., 2019), and we investigate the impact of the assumed $\alpha_{CO}$ on our conclusions. In order to do so, we adopt a ULIRG-like $\alpha_{CO}=0.8$ (Bolatto et al., 2013). In Figure 12c and 12d, we illustrate the effect of this shift CO-detected SQuIGG$\vec{L}$E  galaxies, with No Buried SF models on the left and Buried SF Allowed models on the right.

A 0.7 dex shift in the inferred $H_{2}$ mass from an altered $\alpha_{CO}$ is not enough to place the majority of the galaxies at depletion times where the $L_{CO(2-1)}$ versus age trend can be fully explained by the low-lying star formation in the No Buried SF models. Even combining this effect with buried star formation cannot resolve the tension; while many sources can reach depletion times of $\sim 250$ Myr, the observed $\sim 140$ Myr target remains elusive as none of the SQuIGG$\vec{L}$E  galaxies are consistent with still being in the midst of active starbursts.

We also address one other possibility related to lowering the $H_{2}$ mass in these post-starburst galaxies: the mechanical removal of star forming material via outflows and mergers. Previous works leveraging higher resolution CO(2-1) observations have shown that as much as half the CO luminosity in luminous SQuIGG$\vec{L}$E  sources lies in tidal features on scales of tens of kiloparsecs (Spilker et al., 2022; D’Onofrio et al., 2025), as well as in neighboring galaxies that may have recently interacted with SQuIGG$\vec{L}$E  systems (Kumar et al. 2025 submitted). Even in the low-resolution data we present in this work, similar extended CO can be seen most clearly in J1157+0132 (see Figure 2), our most CO-luminous system, where extended CO along the tidal features is kinematically distinct from the core of the galaxy. We also note that such extended CO luminosity is seen in co-eval massive starbursts and may be partially driven by outflows (Geach et al., 2018); similar evidence for outflows playing a large role in ISM removal has been found in high-redshift quiescent systems (Belli et al., 2024; Sun et al., 2025), which could potentially drive the extended CO seen in the aforementioned works. However, in the absence of higher resolution CO maps, we can only speculate, and at least at present, the CO centroid and velocity are well-aligned with the rest optical locations of our galaxies, suggesting that the majority of the molecular gas in SQuIGG$\vec{L}$E  post-starbursts will not be fully removed by any present outflow.

5.4 Do the data support modifications to the star formation rate or molecular gas content?

In the previous sections, we demonstrated that while even the combination of buried star formation and a and a reduction in the CO-to-$H_{2}$ conversion cannot totally explain the trend of post-starburst age versus gas fraction, the depletion times can be lowered by an order of magnitude if both these effects are at play in our systems. Here, we speculate briefly on the observational consequences of those modifications, with an eye toward clear tests that could validate or dispute their importance in explaining our observed trends.

Perhaps the most straightforward of these possibilities to test is the presence of buried star formation. From the perspective of our SED fitting, only J1157+0132 is better fit in the IR by the set of priors that favor buried star formation (and even then, only a $\sim 0.14$ dex boost to the IR luminosity is needed to fit the 250 $\mu$m data). In all remaining galaxies, there is very little appreciable improvement in the fits, and only the choice of prior impacts the star formation rate. Outside of SED fitting, direct measurements of the star formation rate in SQuIGG$\vec{L}$E  galaxies have come from rest-optical tracers, with the most rigorous constraints coming from weak $H\alpha$ emission (Zhu et al., 2025). However, $H\alpha$ is highly susceptible to dust attenuation. If the effective $A_{V}$ in $H_{\mathrm{II}}$ regions is $\gtrsim 3$, as allowed in our Buried SF Allowed fits, the true $H\alpha$ luminosity, and therefore star formation rates, could be higher than is inferred in that work by a factor of 10. This degeneracy can be broken with longer wavelength observations of instantaneous star formation rate tracers like Paschen and Brackett series emission, as the rest-NIR is far less susceptible to the effects of dust attenuation than the rest optical. Such observations of low-z post-starbursts found that there is little evidence from instantaneous SFR tracers of highly buried star formation (Luo et al., 2022; French et al., 2023). However, no such observations have been presented for massive, CO-luminous post-starburst systems like those in SQuIGG$\vec{L}$E.

Addressing the conversion from CO luminosity to $H_{2}$ mass remains important given that varying assumptions can lead to a factor of 5 difference in the inferred star forming material (Bolatto et al., 2013). One potential way forward would be instead estimating the molecular gas mass from the dust continuum, but different choices in dust emissivity (e.g., Dunne et al., 2003; Draine et al., 2007) can easily induce factor-of-2 changes in the dust mass, as can shifts of $\sim 20$ K in dust temperature. Furthermore, the gas-to-dust ratio, which would underpin this conversion, is itself found to have order-of-magnitude level uncertainty in quiescent systems in both simulations (Whitaker et al., 2021b) and observations (Spilker et al., 2025), making any comparison from dust to gas just as systematically uncertain. That said, the application of some very standard choices of $T_{dust}=25$ K, Dunne et al. (2003) emissivity, and the $\delta_{GDR}=100$ (gas-to-dust ratio) to our four most luminous CO systems (all of which are detected in the continuum at 2mm) suggests that $\alpha_{CO}=2-3.5$, nowhere near the factor-of-5 reduction that would be obtained by invoking a ULIRG-like $\alpha_{CO}=0.8$. This same preference for Milky Way-like $\alpha_{CO}$ was found for the two SQuIGG$\vec{L}$E  systems studied in D’Onofrio et al. (2025) under these same assumptions. However, all comparisons under standard $\delta_{GDR}$ assumptions are currently highly uncertain given that quiescent galaxies clearly depart from the typical relations between these quantities (Spilker et al., 2025).

Addressing the importance of mechanical gas removal is also fairly straightforward. Spilker et al. (2022) and D’Onofrio et al. (2025) have already shown that in many CO-luminous galaxies, as much as 50% of the CO luminosity lies at distances $r>10$ kpc, and that gas, in the harsh environment of the circumgalactic medium, could potentially be disrupted without the need to form additional stars. Future high-resolution studies could better understand the bulk kinematics and spatial distribution of potential star forming gas. Absorption line studies could also search for evidence of outflowing material, which could have been launched by supernovae or AGN activity during the starburst and may play a role in the removal of the last of the ISM (e.g., Geach et al., 2018; Belli et al., 2024; Sun et al., 2025).

At present, there is little data-driven pressure toward any modification to the star formation rate (from buried star formation) or the $H_{2}$ mass (from a systematic lowering of $\alpha_{CO}$ or mechanical removal of star forming material) in the depletion time calculations for SQuIGG$\vec{L}$E  post-starburst galaxies. However, with current data, it is also not possible to rule out either of these modifications. Given that without them, there is no way to interpret the trend in Figure 11 as an evolutionary trend, future observations that can test these assumptions are crucial to understanding the future evolution of these perplexing CO-luminous systems.

5.5 A final possibility: rejuvenation removes the need for rapid depletion

Finally, we acknowledge that the assumption that underpins the aforementioned discussion–that the trend in Figure 11 should be interpreted as an evolutionary sequence of CO luminosity with time since quenching–may be wrong. It remains a distinct possibility, regardless of exactly what the present day star formation rate of SQuIGG$\vec{L}$E  galaxies is, that the young, gas rich galaxies in our sample are simply in a temporary lull where the combination of dust and a bursty star formation history post-merger make them temporarily appear quiescent (an effect that is seen in simulations, see Cenci et al. in preparation). A small number ($\sim 25\%$) of the young ($t_{\mathrm{PSB}}<200$ Myr) SQuIGG$\vec{L}$E  galaxies are not detected in CO(2-1), and even among CO-detected galaxies, there is an order of magnitude spread in implied molecular gas mass. The CO-poor (or undetected) SQuIGG$\vec{L}$E  galaxies may represent “true” post-starburst galaxies that are entering into quiescence, while the gas rich galaxies have not yet finished their life as star forming galaxies and will experience one or more additional periods of star formation before quenching in actuality.

There is some qualitative evidence that this might be the case if we take the ansatz that all massive SQuIGG$\vec{L}$E  post-starburst galaxies have their starburst and quenching induced by a gas rich major merger, and that the absence of clear tidal features signals that the merger occurred further in the past. Of the most CO-luminous galaxies in our sample, nine present very bright tidal features that were flagged in Verrico et al. (2023), and the remaining one (J0909+0108) was flagged as a close pair–though it too exhibits a clear tidal feature to the east. These galaxies could be those where a first burst has occurred, imparting enough feedback (stellar and AGN) into the remaining material to lower the star formation efficiency temporarily. These galaxies will eventually begin bursting again, and only after this will they enter the gas-poor post-starburst population, with considerably fainter tidal features owing to the additional time that has passed since the initial merger. In low-z systems, studies have found that AGN activity occurs almost immediately after coalescence, while an elevated post-starburst fraction occurs 160-480 Myr after the merger (Ellison et al., 2024, 2025). A similar sequence could occurring in the SQuIGG$\vec{L}$E  sample, and our CO luminous galaxies may be fated to become starbursts again on very short ($\sim$100 Myr) timescales.

Quantitatively addressing the impact of rejuvenation within the SQuIGG$\vec{L}$E  sample is more difficult, because the completeness of our selection and the SDSS targeting as a function of age, burst mass, etc. is not well quantified. Self-consistently characterizing the timescales of quiescence would require a full census of star-forming and quiescent populations above a mass limit. With such a sample, one could imagine a detailed accounting of the number density of quiescent galaxies as a function of age to test whether galaxies indeed experience multiple periods where they appear quiescent before re-joining the star forming population. Such a signal would manifest in a pile up of sources at young ($t<200$ Myr) age, where sources like the gas rich galaxies in SQuIGG$\vec{L}$E  could live once after an initial burst and again after they exhaust their final gas supply. Future surveys like the Prime Focus Spectrograph Galaxy Evolution Survey (Greene et al., 2022) may be able to make significant strides here, with mass complete spectroscopic samples of galaxies at cosmic noon where rapid quenching must be happening in significant numbers. However, future detailed studies of the ISM in SQuIGG$\vec{L}$E  galaxies may also offer clues about whether rejuvenation is likely in individual galaxies.

6 Conclusions

In this work, we perform a systematic study of CO(2-1) in 50 massive ($\log(M_{\star}/M_{\odot})\sim 11.2$) post-starburst galaxies at $z\sim 0.7$, finding that a large fraction (27/50) are detected at the 3$\sigma$ level. We then set out to understand why these galaxies still host large molecular gas reservoirs despite appearing to be quiescent in the rest-optical (Suess et al., 2022a; Zhu et al., 2025).

First, we address the question of whether these galaxies are still forming stars at an elevated rate behind optically thick dust. To do so, we perform new spectrophotometric fitting incorporating WISE W3/W4 photometry that was not included in previous fits, in addition to Herschel/SPIRE 250 $\mu$m photometry for a small number of sources that are detected. While we find that the full panchromatic SEDs of these galaxies can be well replicated without the need for buried star formation (see Figures 6, 8, and 9), allowing for optically thick birth cloud dust can increase inferred star formation rates by $\sim 0.5$ dex. This shift can place some of these galaxies onto the star forming main sequence as measured in (Leja et al., 2022, see Figure 10), but there is very little room for any of them to be in the midst of active starbursts in our modeling framework.

We then address the time evolution of the CO luminosity. Previous findings indicate the CO luminosity drops by an order of magnitude in the few hundred Myr after quenching (French et al., 2015; Bezanson et al., 2022, see Figure 11). We find that this implies depletion times of $\lesssim 100$ Myr, and we set out to understand the conditions that could lower their CO luminosities by an order of magnitude on that timescale. We find that under the standard assumptions made in the literature about post-starburst galaxies (that there is little to no star formation behind optically thick dust and that Milky Way assumptions about the CO thermalization and CO-to-$H_{2}$ conversions apply, see French et al., 2015; Suess et al., 2017; Smercina et al., 2018; French et al., 2018; Belli et al., 2021; Bezanson et al., 2022; Zanella et al., 2023), the depletion time is an order of magnitude longer than the necessary $\sim 140$ Myr (see Figure 12a).

We also find that the $\sim 0.5$ dex boost from buried star formation is also not sufficient to match the observed $L_{CO(2-1)}$ versus age trend (see Figure 12b), and the assumption of a lower starburst $\alpha_{CO}$ alone similarly can achieve the necessary depletion times (see Figure 12c). Together, some amount of contribution from either of these mechanisms can lower the depletion times by more than an order of magnitude, but even still, the target depletion time of $\sim 100$ Myr remains out of reach and cleanly explaining the age versus gas fraction trend as an evolution from gas-rich to gas-poor remains difficult.

At present, the evolutionary state of CO-luminous SQuIGG$\vec{L}$E  galaxies remains unclear. Is there moderate star formation occurring behind optically thick dust, or are even the CO-luminous galaxies consistent with lying well below the star forming main sequence? Are they host to $>10^{10}$ $M_{\odot}$ molecular gas reservoirs, or is their molecular gas mass considerably overestimated due to the application of a Milky Way $\alpha_{CO}$ value? If the former is the case, what is stabilizing these large molecular gas reservoirs against collapse? Is significant gas being removed in tidal driven stripping or outflows? And in any of the aforementioned scenarios, will these galaxies directly evolve directly into a quiescent population, or will they rejoin the star forming population before quenching permanently?

Given the significant implications of all these possibilities to the evolutionary pathway that is thought to produce red-and-dead quiescent galaxies, further study of the ISM in these systems is needed. The formation of massive quiescent galaxies via rapid quenching process is empirically necessary, especially at high-z where it is essentially impossible to model observed spectra without extreme star formation that shuts off on short timescales (e.g., Carnall et al., 2023b; Glazebrook et al., 2024; de Graaff et al., 2025; Beverage et al., 2025). Furthermore, the presence of sub-millimeter galaxies with strong Balmer breaks has demonstrated that the physical process that is acting in SQuIGG$\vec{L}$E  galaxies is also likely at work at cosmic noon (e.g., Cooper et al., 2025). It is crucial that we use these lower-z analogues where detailed studies of gas and stars across the panchromatic SED are possible to understand the physics that drive the rapid cessation of star formation after starbursts.

DJS acknowledges Joel Leja for many helpful conversations regarding the implementation of dust emission in Prospector, Jared Siegel for helping to hack out modifications to star formation histories, and Dalya Baron for helpful discussions about FIR SED parameterization and Prospector prior choices. Support for this work was provided by The Brinson Foundation through a Brinson Prize Fellowship grant. DJS gratefully acknowldges support from NRAO Student Observing Support from NRAO-SOSPA8-008. KAS, JSS, and VRD gratefully acknowledge support from NSF-AAG#2407954 and 2407955, and NRAO-SOSPA11-006. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2016.1.01126.S, ADS/JAO.ALMA #2017.1.01109.S, ADS/JAO.ALMA #2021.1.01535.S, ADS/JAO.ALMA #2021.1.00988.S,ADS/JAO.ALMA #2021.1.00761.S, and ADS/JAO.ALMA #2022.1.00604.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss4.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics — Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, 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. The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory. Based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.

Appendix A Full tables of observing parameters

Here, we present Table 2 and Table 3, which contain the properties of our ALMA observations and the measured CO(2-1) and continuum properties of our sources.

Appendix B Comparison of “time since quenching” and “post-starburst time”, and alternate definitions

Here, we compare our measured post-starburst times, which once again is defined as the lookback time where the galaxy formed 99% of the total stellar mass it formed in the Gyr before observation, to the time since quenching measured in Suess et al. (2022a). The Suess et al. (2022a) time since quenching relied on taking a derivative of the star formation history, and is therefore sensitive to the parameterization and also whether the galaxy dropped enough to meet the specific definition. This makes the measurement impossible for galaxies that do not fully quench, which becomes increasingly more likely the star formation histories inferred under our Buried SF Allowed priors.

In Figure 13, we compare these metrics for the star formation histories derived in our No Buried SF and Buried SF Allowed priors, finding very good agreement with the Suess et al. (2022a) time since quenching in both cases. The largest outlier, which is $\sim 150$ Myr old in the Suess et al. (2022a) fits but is closer to 15 Myr old by our new metric, is J0910+0218 (see Figure 7). This difference in age reflects an actual change in the star formation history of the source relative to the Suess et al. (2022a) fits; the inclusion of mid- and far-IR photometry necessitated a younger age to properly capture the full SED in this source. Thus, we propose that this $t_{\mathrm{PSB}}$ metric, which has the added advantage of being very easy to compare to simulations, can serve as a good proxy for the time since quenching.

That said, our definition of $t_{PSB}$ is not without problems. As burst gets older, the measure becomes increasingly susceptible to being pulled to younger ages by low lying levels of star formation. As such, it may also be desirable to compute other moments of the star formation history that are more robust to this effect by tracing star formation on longer timescales. Thus, we compute a measure of the post-starburst time, $t_{\mathrm{PSB},90}$, that is identical in definition to our fiducial $t_{\mathrm{PSB}}$, but instead measuring the lookback time where 90% of the stars formed in the last Gyr. This metric tends to measure a value that is closer to the peak of the burst, especially in younger galaxies.

Finally, we seek to measure a quantity which is as close as possible to the quenching time defined in Suess et al. (2022a) but which is still possible to calculate for galaxies which do not drop significantly enough to be classified by a derivative. For this, we adopt a metric $t_{q,10^{-10}\ \mathrm{yr}^{-1}}$, which we define as the lookback time where the star formation rate first exceeds $M_{\star}\times 10^{-10}\ \mathrm{yr}^{-1}$. If the galaxy is already measured to be instantaneously forming stars above this rate, the metric is measured as 1 Myr. This metric has the disadvantage of being scale variant, insofar as we assume an arbitrary cutoff in the specific star formation rate that should in all likelihood evolve with redshift. However, is still serves as a useful measure that can trace the last time the galaxy was forming stars at anywhere near main sequence star formation rates.

In Figure 14, we present the same gas fraction vs age plot as in Figure 11, but instead using these different age metrics measured from our Buried SF Allowed fits. We use a similar fitting procedure to measure the decay timescale, and find we measure a timescale for both these parameters that are 1$\sigma$ consistent with the measurements using our fiducial $t_{\mathrm{PSB}}$. As such, the finding of decay that is too rapid to be produced by even the combination of buried star formation and a significantly lower $H_{2}$ mass is robust to the definition of age.