The Fate of the Interstellar Medium in Early-type Galaxies. III. The Mechanism of Interstellar Medium Removal and the Quenching of Star Formation

Understanding how galaxies quench their star formation is crucial for studies of galaxy evolution. Quenching is related to a decrease of cold gas. In the first paper we showed that the dust removal timescale in early-type galaxies (ETGs) is about 2.5 Gyr. Here we present carbon monoxide and 21 cm hydrogen line observations of these galaxies and measure the timescale of removal of the cold interstellar medium (ISM). We find that all the cold ISM components (dust and molecular and atomic gas) decline at similar rates. This allows us to rule out a wide range of potential ISM-removal mechanisms (including starburst-driven outflows, astration, or a decline in the number of asymptotic giant branch stars), and artificial effects like the stellar mass–age correlation, environmental influence, mergers, and selection bias, leaving ionization by evolved low-mass stars and ionization/outflows by Type Ia supernovae or active galactic nuclei as viable mechanisms. We also provide evidence for an internal origin of the detected ISMs. Moreover, we find that the quenching of star formation in these galaxies cannot be explained by a reduction in the gas amount alone, because the star formation rates (SFRs) decrease faster (on a timescale of about 1.8 Gyr) than the amount of cold gas. Furthermore, the star formation efficiency (SFE) of the ETGs ( SFE≡SFR/MH2 ) is lower than that of star-forming galaxies, whereas their gas mass fractions ( fH2≡MH2/M* ) are normal. This may be explained by the stabilization of gas against fragmentation, for example due to morphological quenching, turbulence, or magnetic fields.


introduction
In order to have a full picture of galaxy evolution, we need to understand how galaxies become passive, i.e., how they stop forming stars, the process called quenching.Star formation ceases either when gas is removed from a galaxy or is made unable to form stars.A galaxy can run out of cold gas when it is used for star formation (astration; Schawinski et al. 2014;Peng et al. 2015).On the other hand, gas can be expelled or ionized by either supernovae (SNe; Dekel & Silk 1986;Ceverino & Klypin 2009;Muratov et al. 2015;Hopkins et al. 2018a;Li et al. 2020) or evolved low-mass stars (Binette et al. 1994;Conroy et al. 2015;Herpich et al. 2018;Hopkins et al. 2018b).
The existence of gas in many passive galaxies may contradict the interstellar medium (ISM) removal as the only mechanism of quenching.
The origin of the ISM in ETGs remains unclear.One possibility is that it is of internal origin, e.g., the left-overs from past star formation or released by low-mass stars (Knapp et al. 1992;Rowlands et al. 2012;Micha lowski et al. 2019).If gas is brought in from the outside, then the orientation of its rotation is expected to be random with respect to the kinematic axes of a galaxy.However, analysing the kinematic misalignment of stellar and gas components in ETGs, Davis & Bureau (2016) found that the paucity of counter-rotating gas disks implies very short gas depletion rates and unrealistically high merger rates (in order to match the gas detection rate which would otherwise be low for short depletion times).Alternatively, a very long gas relaxation timescale must be invoked, which is consistent with cosmological simulations (van de Voort et al. 2015).The alignment of the gas/dust disk and stellar component has also been used to argue against an external origin of the ISM in ETGs (Bassett et al. 2017;Sansom et al. 2019;Richtler et al. 2020).Moreover, Griffith et al. (2019) found that the stellar and gas metallicities of ETGs are similar, suggesting an internal origin of gas.Finally, Babyk et al. (2019) found that the molecular gas mass in ETGs is correlated with their hot gas mass, also suggesting an internal origin.In a similar vein, cold gas in simulated ETGs comes from cooling from the hot halo (Lagos et al. 2014).
The other possibility is an external origin of ISM.In this scenario the ISM is acquired by ETGs by mergers with gasrich dwarf galaxies or gas inflows.Davis & Young (2019) found that only 7% of ETGs have a lower gas metallicity than stellar metallicity (a clear signature of an external origin of gas), but given very short enrichment timescale (and hence short visibility of low-metallicity features), they estimated that for at least a third of ETGs the gas is of external origin.Moreover, ETGs with dust lanes contain cold ISM, which has also been shown to be brought in by minor mergers (Davis et al. 2015).An external source of the ISM has been claimed for around half of ETGs, which have misaligned stellar and gas disks (Davis et al. 2011;Barrera-Ballesteros et al. 2014, 2015;Jin et al. 2016;Bryant et al. 2019).This has also been supported by other works (Young et al. 2014;Woodrum et al. 2022;Cao et al. 2022;Lee et al. 2023).This topic has been investigated in simulations by Lagos et al. (2015) who found that a misalignment of the rotation axis of the gas and stellar components is mostly the consequence of cold gas flows.This is the third paper in a series in which we analyze ISM removal from ETGs.We use the term ETG for galaxies which are morphologically classified as ellipticals, lenticulars (S0), or early-type spirals (Sa and SBa).In Micha lowski et al. (2019, Paper I, M19 hereafter) we presented the decline of dust mass as a function of stellar age, measurement of the dust removal timescale and the origin of dust in these galaxies.In Leśniewska et al. (2023) and Ryzhov et al. (in prep.)we present an expanded analysis of 2 000 of these galaxies, allowing us to analyse which galaxy properties influence the dust decline.In Nadolny et al. (in prep) we found similar galaxies in simulations, providing a physical insight into the mechanism of this process.The objectives of the present paper are: i) to determine the mechanism of the ISM decline in ETGs, and ii) to constrain the mechanism of quenching.
The paper is structured as follows.In Section 2 we present the ETG sample and in Section 3 we describe our new CO and H I data.Section 4 describes the numerical galaxy evolution model used to interpret the data.We present the results in Section 5. We discuss the implication of the gas removal in ETGs on quenching of star formation in Section 6.In Section 7 we discuss possible mechanisms for this gas removal and in Section 8 we discuss the source of energy needed for this process.We close with a summary of our results in Section 9. We use a cosmological model with H 0 = 70 km s −1 Mpc −1 , Ω Λ = 0.7, and Ω m = 0.3.We also assume a Chabrier (2003) initial mass function (IMF), to which all star formation rates (SFRs) and stellar masses were converted (by dividing by 1.6) if given originally assuming the Salpeter (1955) IMF.Errors are given as 1σ.

sample
As in M19, we use the sample of dusty ETGs from Rowlands et al. (2012).This sample includes all galaxies with elliptical/lenticular morphology and red spirals from the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) ∼ 14 deg 2 Science Demonstration Field (Ibar et al. 2010;Pascale et al. 2011;Rigby et al. 2011;Smith et al. 2011) that are detected at 250 µm.The specific selection criteria used by Rowlands et al. (2012) are: 2. Matched to an optical Sloan Digital Sky Survey (SDSS) source with a spectroscopic redshift in a range 0.01 < z spec < 0.32 within a 10 ′′ radius and with a match reliability greater than 0.8 .
4. Visually classified by Rowlands et al. (2012) as early-type (elliptical or S0), or red spiral with nearultraviolet (NUV) to r-band color of NUV−r > 4.5.The color selection has not been applied to ellipticals or lenticulars.
The sample consists of 61 galaxies, including 42 ellipticals or lenticulars, and 19 red spirals (mostly Sa or SBa).These galaxies will here be referred to as ETGs.Recently, Zhou et al. (2021) found that such red spirals have similar star formation histories as ellipticals, so they are treated collectively.We used the galaxy properties derived by Rowlands et al. (2012) based on the spectral energy distribution modelling using the data from the Galaxy And Mass Assembly (GAMA) survey (Driver et al. 2011(Driver et al. , 2016;;Hill et al. 2011;Robotham et al. 2010;Baldry et al. 2010).In order to assess the AGN power, we also used the [O III] data from single fiber spectroscopy in the GAMA survey (Gordon et al. 2017) 16 .Twelve of our galaxies are located within their respective star formation main sequence (Micha lowski et al. 2019).

data
We observed the CO and H I lines of 13 galaxies from the sample (nine ellipticals and four red spirals).We randomly selected them to have even coverage of the entire stellar age range.They also cover the entire stellar mass range from 10 9.7 to 10 11 M ⊙ .

IRAM30m/EMIR: CO lines
We performed observations with the IRAM 30-m telescope17 using the Eight MIxer Receiver18 (EMIR; Carter et al. 2012).We implemented the wobbler switching mode (with the offset to the reference positions of 60 ′′ ), which provides stable and flat baselines and optimizes the total observing time.We centered one intermediate frequency (IF) at the frequency of the CO(1-0) line and the other at the frequency of the CO(2-1) line (the latter was not possible for all sources, see Table 1).We used the Fourier Transform Spectrometers 200 (FTS-200) backend providing 195 kHz spectral resolution and 16 GHz bandwidth in each linear polarization.The observations were divided into 6-min scans, each consisting of 12 scans 30 s long.The pointing was verified every 1-2 hr on a nearby quasar 0823+033.The observing log is presented in Table 1.We reduced the data using the Continuum and Line Analysis Single Dish Software (Class) package within the Grenoble Image and Line Data Analysis Software19 (Gildas; Pety 2005).Each spectrum for a given galaxy was calibrated, and corrected for baseline shape.Then all spectra were averaged.
The IRAM30m/EMIR spectra were binned to 30 km s −1 channels.A Gaussian was fitted to the binned data, and in the case of detections the 2σ width of the Gaussian was adopted to integrate the line flux.In cases of non-detections, a [−200, 200] km s −1 width was adopted.The error per spectral channel was calculated using the ranges [−900, −400] km s −1 and [400, 900] km s −1 .Then the uncertainty of the flux estimation was calculated by a Monte Carlo simulation.We calculated the line luminosities based on eq. 3 in Solomon et al. (1997).The molecular masses were calculated using α CO = 5 M ⊙ (K km s −1 pc 2 ) −1 (the conversion includes helium).The choice of the Galactic CO-to-H 2 conversion factor is justified by the fact that most ETGs have solar metallicity (Conroy et al. 2014;Davis & Young 2019).The widths of the lines, integrated fluxes, luminosities and the resulting molecular gas masses are presented in Table 2.In addition to the results from the CO(1-0) line, the table shows those from the CO(2-1) line if such tuning was possible.The CO spectra are shown in the left and middle columns of Fig. A10 in the Appendix.Out of 13 targets for CO(1-0) we detected nine, and out of nine targets for CO(2-1) we detected five.

GBT: H I line
We performed observations with the Green Bank Telescope (GBT)20 using the Versatile GBT Astronomical Spectrometer (VEGAS).We used the mode with 61 kHz spectral resolution (corresponding to 14 km s −1 at the H I frequency).The observations were divided into scans lasting 3 or 5 min.The flux calibration was done using observations of 3C286, whereas pointing and focus was verified using observations of 0744-0629 (radio source 4C -06.18) every three hours.The observing log is presented in Table 3.We used the GBTIDL package21 to reduce the data.We calibrated each spectrum individually and then averaged them.
We processed the H I GBT/VEGAS spectra in a similar way as for CO.The data were not usable for J085934.1+003629 and J091205.8+002656due to strong radio-frequency interference (RFI).We calculated the atomic gas masses based on eq. 2 in Devereux & Young (1990).The widths of the lines, integrated fluxes, luminosities and the resulting atomic gas masses are presented in Table 4.The H I spectra are shown in the right column in Fig. A10 in the Appendix.Out of eight targets with usable data we detected seven.The features at 500 km s −1 for J085828.5+003814 and J090551.5+010752 are likely due to RFI, because they are present only in a fraction of the data.

numerical model
To aid the interpretation of the observed data we use the chemical dust evolution model from Gall et al. (2011a, Table 1 Observing log for the IRAM30m/EMIR observations with integration times and 225 GHz atmospheric opacity.

Galaxy
Obs. date t int τ 225 GHz (hr) J085828.5+0038142015 Jul 28 302 ± 10 1.83 ± 0.52 11.59 ± 3.32 58.0 ± 16.6 254 ± 16 3.13 ± 0.33 19.97 ± 2.09 99.9 ± 10.5 Note. -The columns show the Gaussian FWHM of the line, the integrated line flux within the dotted lines on Fig. A10, the line luminosity and the molecular gas mass assuming α CO = 5 M ⊙ (K km s −1 pc 2 ) −1 .The left set is for CO(1-0) and the right set is CO(2-1) for galaxies for which simultaneous setting for this line was possible.Note.-The columns show the Gaussian FWHM of the line, the integrated line flux within the dotted lines on Fig. A10, the line luminosity and the atomic gas mass.
hereafter GAH11), allowing us to follow the dust and gas removal with stellar age.The model considers SN and asymptotic giant branch (AGB) star dust production at different efficiencies (Gall et al. 2011b,c).It assumes that all material is recycled, and hence available for star formation, instantaneously after the death of the stars.Furthermore, ISM dust destruction through SN shocks is either turned off or set at a moderate rate.We calculate a suite of models for which we either switch on or off SN and AGB star dust and gas recycling.Additionally we investigate gas and dust removal through outflows.
The evolution of dust and gas in a galaxy is: (1) (2) Here, E d,SN (t) and E d,AGB (t) are the SN and AGB dust production (injection) rates as defined in GAH11 while E g,SN (t) and E g,AGB (t) are the rates for the recycled gas phase elements (Eq.6 in GAH11).The variable η d (t) = M d (t)/M ISM (t) can be understood as the 'dustto-gas mass ratio' or the fraction of dust present in the ISM.The evolution of the SFR, ψ(t), here is chosen to represent the measured SFRs (Fig. 1 in M19 and eq. 8 below) as where ψ ini is the initial SFR, T gal (t) is the age of the galaxy and τ SFR = 1.8 Gyr is the star formation decline timescale (duration of the star formation episode).We note that this SFR evolution is not coupled to the total amount of gas in the galaxy as it would be by a Schmidt-Kennicutt law (Kennicutt 1998) and as defined in GAH11.Instead, the formulation here is based on observations (Fig. 1 in M19) of the galaxies in question.An additional cold gas removal is considered, either in a form of gas heating or physical ISM removal.It is modeled as where M ISM,out is the amount of ISM material to be removed over a timescale τ ISM .The variable ξ SN (t) = M cl R SN (t) defines how much ISM mass, M cl , is completely cleared of dust by core-collapse SNe (CCSNe; Dwek et al. 2007), with R SN (t) being the CCSN rate.We have been testing zero to moderate dust destruction with M cl of 0, 50, 100 or 500 M ⊙ .Both the CCSNe rate and AGB rate are calculated as which assumes that only single stars form.R AGB,SN (t) is regulated by the IMF, φ(m), and ψ(t − τ ) is the lifetime of a star with a given zero-age main sequence mass m.The model parameters are summarized in Table A6.

results
Our aim is to investigate the time evolution of various galaxy properties, using the stellar age as a time proxy, so first we verify that indeed the ages of these galaxies increase linearly with time.In principle this is not the case for all galaxies, especially for highly star-forming ones for which the mean stellar age may even decrease with time.Strictly speaking the stellar age increases linearly with time when star formation rate is zero.Our galaxies have very low levels of star formation and large stellar masses, so they should be close to this approximation, because the mean stellar age is not affected significantly by the presence of new populations of stars and therefore should grow linearly with time.
We checked the age evolution with time using 2 000 simulated analogs of our galaxies selected in a very similar way: similar redshift distribution and stellar masses, dust masses, so that the 250 µm flux would be detectable as for the real sample, early-type morphological classification, and low star-formation activity placing them below the main-sequence (Nadolny et al., in prep.).These are Millennium Simulations with LGalaxies semi-analytical models with ∼ 100 3 Mpc 3 box (Lemson & Virgo Consortium 2006;Springel et al. 2005b;Henriques et al. 2020).Unlike for observed galaxies, for each of the simulated galaxies we have a complete knowledge of the light-weighted stellar age as a function of the Universe age (cosmic time).
In order to analyse the time evolution of ages for the entire sample, we need to choose a common time zeropoint.Galaxies observed at similar ages of the Universe have a range of stellar ages, so in order to assess how the stellar age evolve with time for the entire sample, we shifted horizontally (in time) the age-time tracks, so that at the time of observation the cosmic time is equal to the measured stellar age.With such a choice of the time zeropoint, galaxies with higher levels of current star formation activity than in the past should be located below the age = time line (ages lower than the cosmic time value) due to a significant number of young stars.On the other hand, galaxies with long periods of constant star formation rate should be located above this line (ages higher than the cosmic time value), because their mean stellar age does not change in time.
Fig. 1 shows the age-time evolution for simulated dusty ETGs.It is clear that neither of the non-linear scenarios described above applies to them.Their mean stellar ages increase linearly with time at least during the last 8 Gyr, which is the range of interest of our study.Very low scatter around the age = time line indicates that very few galaxies exhibit levels of star formation strong enough to break the linearity of stellar ages with time.
Having demonstrated that the stellar ages of dusty ETGs increase linearly with time, at least for the dusty ETGs selected in the way we do here, we come back to the observed sample.In order to compare galaxies with different masses, we normalize the gas and dust masses by their stellar masses.We show the molecular and atomic gas-to-stellar mass ratios (M H2 /M * , M HI /M * ) as a function of stellar age in Fig. 2. We detect a decline of the gas-to-stellar mass ratios with age, similar to the evolution of the dust-to-stellar ratio (M dust /M * ) from M19 and Leśniewska et al. (2023).For the M H2 /M * -age and  et al., in prep.).The black solid line is the running average for given time, whereas the grey shaded region represent the 1σ range.The dashed line indicates the linear evolution when the stellar age is equal to the value of cosmic time.The time reaches negative values because of how the moment when time is equal zero was set, in order to have the time equal to the stellar age at the moment of observation.This figure shows that stellar ages increase linearly with time for dusty ETGs, so can be used as a time proxy.
The decline of the gas amount in ETGs with time was seen in models (Calura et al. 2017), but has not been measured directly before.On the other hand, Smercina et al. (2018) detected a dust decline with age for post-starburst galaxies with ages up to 1 Gyr, and a very weak gas decline.They interpret this as the effect of sputtering of dust grains in hot gas.With a larger sample of post-starburst galaxies, French et al. (2018a) detected a gas decline on the timescale of 100-200 Myr.Similarly, molecular gas was found in post-starburst galaxies at z ∼ 0.6 only for those which were quenched less than 150 Myr ago, indicating a rapid gas removal (Bezanson et al. 2022).The timescale we measure is comparable to that inferred by Gobat et al. (2020), based on the low gas fraction of low-z ETGs, and higher than the quenching timescale due to environmental influence (Gobat et al. 2015).
The gas-to-dust ratios are either independent of or weakly declining with age (Fig. 3).There is a hint of a decline of the M H2 /M dust ratio with age at a ∼ 2.9σ level (probability of the null hypothesis of no correlation of 0.0037 and the Spearman rank correlation coefficient of −0.85).There is no indication of the decline of the M HI /M dust ratio (the Spearman rank correlation coefficient of −0.09, the probability of the null hypothesis of 0.87).We note that the gas-to-dust ratios we found are similar to those of star-forming galaxies.Hence, these ETGs do not belong to the class of passive galaxies at 0 < z < 3 with extremely high ratios, recently identified in simulations (Whitaker et al. 2021;Donevski et al. 2023).
The ETGs in our sample follow the evolutionary sequence found for other galaxies based on ISM proper-ties.In Fig. 4 we show the dust-to-baryon mass ratio [M dust /(M * + M HI )] as a function of atomic gas fraction [M HI /(M * + M HI )] together with dust-selected galaxies (Clark et al. 2015).Modelling shows that galaxies evolve from left (high gas fractions) to right (low gas fractions), increasing their dust-to-baryon ratios initially and then following a decline (Clark et al. 2015;De Vis et al. 2017a,b, 2019;Donevski et al. 2020;Nanni et al. 2020).The ETGs in our sample have low gas fractions, consistently with being among the oldest galaxies in this plot.
The only galaxy with high gas fraction is J085915.7+002329,having the lowest stellar mass in our sample [log(M * /M ⊙ ) = 8.62].It is indeed one of the youngest as well with an age of 10 9.3 yr.
We note that the exponential decline parametrization for SFR in the above equation and in the model (eq.3) is justified for these old galaxies.First, this function fits the data (Fig. 1 of M19).Second, the SED fits from Rowlands et al. (2012) resulted in no recent starbursts in this sample with the time since the last starburst > 500 Myr for all galaxies, > 1 Gyr for 80% of the sample, and > 2 Gyr for 56% of the sample.The existence of these bursts does not affect our analysis because in eq. 8 we quantify the evolution of the average SFR after these episodes.Similarly, the assumption of the exponential decline of the SFR in the model has almost no impact on the resulting properties, because the SFRs are low, so the astration is weak, independent of what parametrization is adopted, as long as it is consistent with the low measured values.
Galaxies in our sample have low SFRs for their molecular gas masses, as demonstrated in Fig. 5.They are below the Schmidt-Kennicutt law by a factor of several.
These efficiencies are similarly low as for the sample of ETGs with dust lanes (Davis et al. 2015) and some poststarburst galaxies (Alatalo et al. 2015b;French et al. 2015French et al. , 2023;;Rowlands et al. 2015;Smercina et al. 2018Smercina et al. , 2022;;Bezanson et al. 2022;Luo et al. 2022); at the lower limit for high-z ETGs (Magdis et al. 2021); and lower than the full ETG ATLAS 3D sample (Davis et al. 2014;Kokusho et al. 2017).In Fig. 6 we compare the SFEs of our sample with those of the general star formation population.
For molecular and atomic gas we used equation 1 of Micha lowski et al. ( 2018) and equation 1 of Micha lowski et al. (2015), respectively (see similar estimates in Sargent et al. 2014;Saintonge et al. 2017;Liu et al. 2019;Magdis et al. 2021).Again, the molecular SFEs of galaxies in our sample are lower than those of star-forming galaxies.The difference is stronger for older  -Ratios of molecular (large red squares and arrows for upper limits), atomic (large blue stars) gas, and dust (black circles; M19) masses to stellar masses as a function of luminosity-weighted stellar age of ETGs galaxies detected by Herschel from the sample of Rowlands et al. (2012).The exponential fits to the gas-to-stellar and dust-to-stellar mass ratios (eq.6 and 7) are shown as solid lines, colored as the datapoints.The red and blue dotted lines denote the curve for the dust-to-stellar ratio shifted upwards by the measured average gas-to-dust ratios of log(M H 2 /M dust ) = 2.144 ± 0.018 and log(M HI /M dust ) = 2.610 ± 0.011.Small red squares and blue stars mark the galaxies which we observed at CO and H I, respectively.Open circles denote galaxies which are within the main sequence (see Fig. 1 in M19).
All ISM components decline at a similar rate.ETGs are below the relation for star-forming galaxies, indicating that their gas reservoir is ceasing to be able to form stars. ages, so these galaxies keep decreasing their efficiency with time.However, the atomic SFEs are similar to those of star-forming galaxies.This could be because the atomic gas in ETGs is only indirectly connected with star formation.

ism removal and quenching
We have provided the a measurement of the timescale of ISM removal in ETGs of ∼ 2.3 Gyr (eq.6).This slow ISM removal suggests that either quenching in these galaxies is a slow process or that the main reason for quenching is not exhaustion of the gas supply, but rather gas stabilization that prevents further star formation.Indeed, we have measured a slow SFR decline (eq.8).Similarly, the timescale of quenching was measured to be of the order of several Gyr in observed (Peng et al. 2015;Trussler et al. 2020;Kipper et al. 2021;Tacchella et al. 2022;Noirot et al. 2022;Bravo et al. 2023;Donevski et al. 2023) and simulated (Trayford et al. 2016;Nelson et al. 2018;Wright et al. 2019;Park et al. 2022;Walters et al. 2022, Nadolny et al., in prep.)galaxies, with the exception of cluster members in which it is shorter (Muzzin et al. 2014;Socolovsky et al. 2018;Davis et al. 2019b;Zavala et al. 2019).
As detailed below, the galaxies in our sample are shutting down their star formation not only because they are running out of gas, but because their SFRs decline faster than the gas amount and they have low SFEs and normal gas fractions.
Even in samples of post-merger and post-starburst galaxies the gas supply was actually higher than in other galaxies with similar masses (Rowlands et al. 2015;French et al. 2015;Alatalo et al. 2016;Suess et al. 2017;Ellison et al. 2018), so quenching in these galaxies is likely related to turbulence, not exhaustion or expulsion of gas.Morphological quenching, the influence of the bulge making the gas resilient against fragmentation (Martig et al. 2009(Martig et al. , 2013;;Bluck et al. 2014Bluck et al. , 2020a;;Bitsakis et al. 2019;Lin et al. 2019;Gensior et al. 2020) or the influence of turbulence and magnetic fields (Padoan & Nordlund 2002;Federrath & Klessen 2012) can also be responsible for decreasing SFRs.Indeed galaxies in our sample evolve upwards from the SFR-M dust relation (having higher M dust than what their SFRs would imply), as predicted for quenching that is not caused by the removal of gas (Hjorth et al. 2014;Leśniewska et al. 2023, M19).However, in simulations morphological quenching is found to be effective at gas fractions below a few % (Martig et al. 2013;Gensior et al. 2020), much lower than for our sample (Fig. 7).
The very low SFE of our ETGs (Fig. 6) indicate that star formation is suppressed even in comparison with other ETGs.This means that the process shutting down the SFRs in these galaxies is not due to physical gas removal, but to its inability to form stars.This again supports the internal quenching scenario, either morphological or connected with turbulence or magnetic fields.
For the ETGs in the ATLAS 3D sample, Kokusho et al. (2017)  black circles denote molecular gas, atomic gas, and total gas, respectively (the latter only for galaxies with both CO and H I measurements).The solid lines represent the SFE of star forming galaxies with SFR = 0.01, 0.1, and 1M ⊙ yr −1 (from bottom to top) using eq. 1 in Micha lowski et al. (2015, atomic gas) and eq. 1 in Micha lowski et al. (2018, molecular gas) colored in the same way as datapoints.The molecular SFEs decline with age and are lower than those of star-forming galaxies, whereas the atomic SFEs stay consistent with those of star-forming galaxies.).The red and blue symbols correspond to averages of mass-selected galaxies and only those on the main sequence, respectively (Saintonge et al. 2017).Gas fractions of ETGs are comparable or higher than those of star-forming galaxies, indicating that the ETGs do not stop forming stars due to lack of gas.
forming galaxies.Hence, together with the decline in the gas fraction with age, this was interpreted as quenching driven by a decrease of the gas reservoir.In contrast, we do see a decline in the SFE in our sample (Fig. 6) and they are below the SFR-M H2 relation (low SFRs for their molecular gas masses; Fig. 5), again pointing at quenching being connected with the decreased ability of gas to form stars, not with a lack of gas.
Moreover, the high molecular gas fractions also point to declining gas amount not being the main reason for declining SFRs.Fig. 7 shows the molecular gas fraction (f H2 ≡ M H2 /M * ) as a function of stellar mass compared with mass-selected low-redshift galaxies (Saintonge et al. 2017).The ETGs in our sample on the main sequence (lower ages) exhibit larger gas mass fractions than these star-forming galaxies, whereas the ETGs below the main sequence have comparable gas fractions to star-forming galaxies.Hence the lack of gas is not the main reason for the ETGs to become passive.This is different from greenvalley galaxies, for which both SFE and f H2 were found to be suppressed (Brownson et al. 2020;Lin et al. 2022).
We note that our galaxies may not be representative of the entire ETG population, because they were selected based on Herschel 250 µm detections.This corresponds to the dust mass detection threshold of 10 5.2 M ⊙ at z = 0.05 and 10 6.7 M ⊙ at z = 0.3 (M19).These are not very high limits, especially given the high stellar masses of our galaxies, but most ETGs were shown to have less dust.In particular, in the ETG parent sample of Rowlands et al. (2012), only 5.5% were detected by Herschel.The re-maining dust-poor ETGs may follow different evolutionary tracks than those presented here.

mechanism of ism removal
Here we discuss the possible physical mechanisms explaining the trend depicted in Fig. 2. For each mechanism we list the predictions which can then be compared with existing and future datasets.The predictions are summarized in Table 5.Only two mechanisms are fully consistent with all the current data: removal of the entire cold ISM (Section 7.1) and outflows (Section 7.2).They are not mutually exclusive, so it could be that they operate together.
In Section 8 we also discuss the energy source required for these mechanisms.

Removal of the entire cold ISM
The destruction of all components of the cold ISM involves destroying dust, molecular and atomic gas.We will consider mechanisms having a stronger effect on dust in Section 7.3.There may be several physical mechanisms responsible for the removal of the entire ISM (planetary nebulae, SNe Type Ia, AGN, hot gas, cosmic rays).We will return to discussing these energy sources in Section 8.
The gas which is ceasing to be in the cold (molecular or atomic) phase is transformed into ionized hot gas.In our sample this requires heating a few times 10 9 M ⊙ of gas.The hot gas in elliptical galaxies has a mass corresponding to 1% of stellar mass (or 10-20% for the most massive ones; Mathews & Brighenti 2003; Sparke & Gallagher bias against old ISM-rich galaxies X bias against young ISM-poor galaxies X Note. -: the prediction is consistent with the data, X: the prediction is inconsistent with the data, ?: the data needed to test this prediction is not available yet.R12: Rowlands et al. (2012) 2006).For our galaxies with stellar masses of M * = 10 10−11 M ⊙ , this corresponds to 10 8−9 M ⊙ of hot gas (or more for the most massive galaxies in our sample).Hence, the amount of gas needed to be ionized by this mechanism is similar to the typical hot gas reservoirs in such galaxies, taking into account that some of the gas will be expelled or used for star formation.
This scenario predicts the following.
1. M gas /M dust should be constant, as both dust and gas are destroyed, or decreasing with age if gas particles are destroyed faster due to their more diffuse distribution than dust, which has a more clumpy distribution.
2. M * should be constant or slightly increasing with age, given the low SFRs of galaxies in our sample.
3. Gas and dust should be relatively uniformly distributed with possible either central concentration, reflecting the initial distribution before the quenching, if the process removing the cold ISM is not violent and operates throughout galaxies, or central deficit, if the source of the energy is in the galaxy center (AGN).
The M gas /M dust ratio indeed declines slightly with age (Fig. 3) and M * is slightly increasing (Fig. 1 in M19).
Our numerical model (Section 4, Fig. A12, Table A6) shows that in order to explain the data, around 3-5 × 10 10 M ⊙ of gas (including both the atomic and molecular phase) needs to be removed on the timescale of 10 Gyr (the average rate of 3-5 M ⊙ yr −1 , which includes the phase of being a normal star-forming galaxy).We note that this is much higher than the measured SFRs for higher ages, at which the gas decline is the strongest, so astration is unlikely to dominate the ISM mechanism (see Section 7.4).

Outflows
Gas and dust can be expelled from galaxies either by AGN-or stellar-induced winds.
This scenario predicts the following.
1. M gas /M dust should be constant with age, as both components are being removed, or decreasing with age if dust distribution is more clumpy (as in Section 7.1).
2. M * should be constant or slightly increasing with age, given low SFRs of galaxies in our sample.
3. Some molecular gas should be located away from optical centers of the galaxies, possibly in filamentary or plume-like structures with scales of around 1-2 kpc (Walter et al. 2002;Feruglio et al. 2010Feruglio et al. , 2015;;Alatalo et al. 2011;Morganti et al. 2015).
4. Multiple velocity peaks with locations incompatible with a rotating disk, or broad line wings should be present (e.g.Cicone et al. 2014).
Predictions 1 and 2 are consistent with the data.We do not have spatial information on the distribution of gas to test prediction 3.In principle we can look for high-velocity wings in the spectra (Fig. A10, prediction 4), but the expected signal is at the level of 10% of the main CO line peak (e.g.Cicone et al. 2014).Our data are not of enough signal-to-noise ratio to test this.
We also assessed how the required outflow rate compares with an empirical calibration.In our sample the molecular gas mass goes down from ∼ 10 10 at an age of 1 Gyr to a few times 10 9 M ⊙ at 9 Gyr.Hence, the average outflow rate during this time must be around 1M ⊙ yr −1 .
We estimated the actual outflow rate from the empirical calibration of Fluetsch et al. (2019, their eq. 5) based on SFR, M * and AGN luminosity (L AGN ).We estimated the AGN luminosity from the [O III] luminosity (L [OIII] ), as done in Fluetsch et al. (2019) for cases with no X-ray data: (Heckman et al. 2004).Only five galaxies in our sample have > 3σ detection of the [O III] line, so we first calculated the molecular outflow rates setting the AGN luminosity to zero.This is shown as circles in Fig. 8. Then we calculated the maximum outflow rate by setting the AGN luminosity to the 2σ upper limit allowed by the [O III] data (arrows on this figure).It is clear that the measured outflow rate are at the level of the required average rate only at the beginning of the evolution.Hence, the outflows calculated in this way are not powerful enough to explain the gas mass decline throughout the entire period.
To demonstrate this we fitted an exponential function to the outflow rate as a function of age and obtained the functional form of the outflow rate in M ⊙ yr −1 of (3.9±1.0)exp{−age/[(1.36±0.24)Gyr]} (the half-life time is 0.41 ± 0.07 Gyr, a dashed line on Fig. 8).Then starting at the molecular mass of 10 10 M ⊙ at an age of 1 Gyr we show in Fig. A11 in the Appendix how much gas and dust is removed with the outflow of such strength.It is clear that the evolution is nearly flat, as the outflow rate is too low (especially at older ages) to explain the gas decline.

Dust grain destruction
Dust grains can be destroyed e.g. by SN shocks, AGN or by sputtering in the hot gas.This scenario predicts the following.
1. M gas /M dust should increase with age as gas is not destroyed, just dust.
2. M * should be constant or slightly increasing with age, given low SFRs of galaxies in our sample.
This mechanism is ruled out by prediction 1.The gasto-dust ratio is not increasing with age (Fig. 3), as directly predicted by this mechanism.The gas masses are also decreasing with time (Fig. 2).Arrows denote the molecular outflow rates with the maximum (2σ) contribution of AGN allowed by the [O III] data.The dashed line is an exponential fit to the data with the form (outflow rate / M ⊙ yr −1 ) = (3.9± 1.0) exp{−age/[(1.36± 0.24) Gyr]}.These outflow rates are too low to explain the ISM decline (see Fig. A11).

Astration/strangulation
Gas and dust are incorporated into newly formed stars (astration), so if the SFR is higher than the gas inflow rate, or if the gas inflow is stopped (a process called strangulation; Peng et al. 2015;Trussler et al. 2020), then the gas reservoir will be depleted This scenario predicts the following.
1. M gas /M dust should be constant with age, as both components are being removed, or slightly increasing if dust distribution is more clumpy and this is where star formation is occurring.
2. M * should be constant or slightly increasing with age, given low SFRs of galaxies in our sample.
3. Given the decreasing SFR with age (Fig. 1 in M19), the dust and gas removal should be weaker at older ages, so the drop of M dust /M * and M gas /M * should flatten at older ages (see the M dust panel in Fig. 1 in M19).
This mechanism is ruled out by prediction 3 and, with a lower statistical significance, by prediction 1.First, Fig. 2 clearly shows steepening of the decline of the M dust /M * ratio, so the dust removal is stronger at high ages.Hence this cannot be connected with the SFR, which is the lowest at high ages.Indeed, M19 showed that the effect of astration with measured SFRs on the dust amount is minor, not able to explain the dust decline.
In order to quantify what effect the astration might have on gas evolution we made calculations similar to those presented in M19 for dust evolution.Fig. 1 there presents an exponential fit to the SFR evolution (see eq, 8 above).We used this, starting with 10 10 M ⊙ of gas to see the effect of astration.This is shown in Fig. A11 in the Appendix.The evolution is nearly flat, showing that the SFRs are not high enough, so astration cannot be responsible for most of the gas and dust decline.Hence, we conclude that, while astration is happening, the effect is too weak to fully explain the data.This is consistent with the conclusion drawn in M19.
Similarly, Fig. A12 in the Appendix shows that models with only astration are inconsistent with the data because the gas removal is the weakest for high ages.Only incorporating additional gas removal makes the models agree with the data.
Finally, the gas-to-dust ratio is not constant (or increasing) with age (Fig. 3).The Spearman correlation coefficient for that plot is −0.85 with a probability of a null hypothesis (no decline) of ∼ 0.0037 (∼ 3σ).

A decline in the number of dust-producing AGB stars
It could be that the decline of M dust /M * and M gas /M * is due to a decline in the number of dust-producing AGB stars at older ages, which also release gas to the ISM.This is only possible if the dust and gas produced by these AGB stars dominate the current dust and gas budget in these galaxies.Otherwise, if their contribution is minor, then the decline in their numbers cannot explain the decline in the dust or gas content.Rowlands et al. (2012) and M19 showed that indeed AGB stars do not contribute significantly to dust production in these galaxies.Given their stellar masses, the number of AGB stars is too low and each of these stars would need to have produced a too high amount of dust.
We now examine whether the gas masses we detected in the ETGs could have had a large contribution from the ejecta of AGB stars.We compare the gas masses with the expected numbers of AGB stars calculated from the stellar masses, using a similar argument as presented in Micha lowski et al. (2010a,b), Micha lowski (2015), and Leśniewska & Micha lowski (2019) (see also Morgan & Edmunds 2003;Dwek et al. 2007;Dwek & Cherchneff 2011;Rowlands et al. 2014).A stellar mass of 10 11 M ⊙ implies around 10 10 AGB stars (1.5 − 8 M ⊙ ), assuming the Chabrier ( 2003) IMF.A gas mass of ∼ 10 10 M ⊙ for the galaxies in our sample at the lowest ages implies that each AGB star would need to release approximately 1 M ⊙ of gas.This is the order of magof the AGB total ejecta mass (Morgan & Edmunds 2003;Ferrarotti & Gail 2006).This calculation neglects past gas consumption due to star formation and outflows (which would increase the required total ejecta mass per star) and inflows (which would decrease the required ejecta mass).It seems that AGB stars are numerous enough to explain the observed gas masses.However, as mentioned above, they do not contribute significantly to dust production, so it is likely that the gas decline is connected with the same mechanism as the dust decline, and not a decreasing number of AGB stars.On the other hand, this calculation supports the internal origin of the ISM in these ETGs, as claimed in M19.
This leaves the question on why the majority of ETGs selected in a similar way to our sample do not have detectable ISM content, given that they have similar stellar masses, so a similar number of AGB stars.We speculate that this may be due to the non-detected galaxies being older, or exhibiting stronger ISM destruction mechanism (e.g. more powerful AGN), or having a higher amount of very hot gas, leading to immediate ionization of all the gas released by AGB stars.

Dust cooling
It could be that with time dust is not destroyed, but cools down, so becomes invisible for Herschel.This scenario predicts the following.

As a result of this extra cold dust component,
longer-wavelength observations should result in high flux levels above the extrapolation from the Herschel wavelengths.

Dust temperatures measured with the
Herschel data should decrease with age.
3. Gas-to-dust ratios should increase, because the gas mass measurement is unaffected.
The data are inconsistent with all these predictions.Our JCMT/SCUBA-2 data rule out this mechanism, because we do not see the flux excess at 850 µm (Fig. 5 in M19).Moreover, we do not see any trend of the dust temperature with age and the range of the temperatures is relatively narrow (Fig. 1 in M19).Finally, the gas content also decreases (Fig. 2) and hence, the gas-to-dust ratios are constant or slightly decreasing (Fig. 3).

Dust heating
If dust is progressively heated to high temperatures (but not destroyed), then it could become invisible for Herschel/SPIRE.This scenario predicts the following.
1.The SEDs of these galaxies should have the peak shifted towards shorter wavelengths.
2. Dust temperature measured with Herschel data should increase with age.
3. Gas-to-dust ratios should increase, because the gas mass measurement is unaffected.
The data are inconsistent with all these predictions.The SEDs of the galaxies in our sample do not show any evidence of a shift of the peak towards shorter wavelengths (Rowlands et al. 2012).Moreover, we do not see any trend of the dust temperature with age (Fig. 1 in M19).Finally, the gas-to-dust ratios do not increase with age (Fig. 3).
An additional warm component with a temperature of 30 K could only contribute a few percent to the cold dust mass in order not to overshoot the 100 µm detections or limits.The contribution of even hotter dust could be even lower because even a small amount of hot dust is too bright.Hence, the possible contribution of hot dust is too small to explain the decrease of the M dust /M * ratio.

Environmental influence
In principle it is possible that the observed trend is the reflection of environmental influence.Galaxies in richer environments are quenched quicker (so the stellar ages are higher) and the ISM is being removed quicker due to environmental effects like interactions or ram-pressure stripping.
This scenario predicts the following.
1.There should be a trend of the environmental density with age (and therefore with the M dust /M * ratio).
2. These galaxies should live in rich environments in which the influence on their ISM is significant.
None of these predictions are consistent with the data.As shown in Fig. 1 in M19, the galaxy density does not depend on the age.Moreover, these galaxies do not reside in rich environments.The projected galaxy densities are below 10 Mpc −2 , below those of galaxy groups.We also do not see any dependence of the dust decline on environment in the extended analysis of 2 000 dusty ETGs (Leśniewska et al. 2023).7.9.Mergers with gas-rich galaxies It could be that the trend cannot be interpreted as a time evolution, but results from mergers of passive gaspoor galaxies with gas-rich star forming galaxies of different masses, or that each passive galaxy experienced different number of mergers.If the merging star-forming galaxy was relatively large (or a passive galaxy had experienced more mergers), then the resulting mean age would be low (because of many young stars brought in during mergers) and the resulting M dust /M * would be high (because more dust is brought).A merger with a small starforming galaxy (or fewer number of mergers) would result in a much higher derived age and a much lower M dust /M * ratio.This scenario predicts the following.
1. M gas /M dust should be increasing with the derived age, because passive galaxies with high derived ages would need to have merged with smaller galaxies, which have high M gas /M dust ratios (Grossi et al. 2010(Grossi et al. , 2015;;Galametz et al. 2011;Leroy et al. 2011;Cortese et al. 2012b;Hunt et al. 2014;Rémy-Ruyer et al. 2014).
2. M * should be constant or slightly decrease with the derived age because merging dwarf galaxies would not contribute much to M * , and the galaxies in our sample with lower inferred ages would have on average experienced more merging and so may be slightly more massive.
3. M gas /M * should be decreasing with the derived age, for the same reason why M dust /M * is decreasing (higher-mass galaxies bring more gas).
4. M dust should only be weakly correlated with M * , because the gas-rich companions would not bring significant amount of stars, so the final M * should only weakly depend on the number of merger events.
5. Sizes should be increasing with decreasing age, because in this scenario a low derived age means more merging and galaxies grow in size as a result of mergers (Naab et al. 2009;Hopkins et al. 2010;Trujillo et al. 2011;Furlong et al. 2017).
This mechanism is ruled out by predictions 1, 2, 4, and 5.It has also been ruled out by Rowlands et al. (2012) and M19 in the context of the analysis of the source of dust in these galaxies (see also Donevski et al. 2023).
The gas-to-dust ratio is not increasing with age (Fig. 3), inconsistent with prediction 1. Fig. 1 in M19 shows that the stellar mass is increasing slightly with age, which should not be the case if the apparent dust decline with age was a result of merging with smaller (or fewer) galaxies at high ages and larger (or more numerous) galaxies at low ages (prediction 2).Moreover, Fig 3 in M19 shows that M dust is correlated with M * , inconsistent with prediction 4. The Spearman rank correlation coefficient is 0.5 with a very small probability (∼ 3 × 10 −5 , ∼ 4σ) of the null hypothesis (no correlation) being acceptable.Fig. 2 in M19 also shows that there is no correlation of the size of these galaxies with age, as would be expected if mergers are the mechanism responsible for the trends we observe (prediction 5).
Conversely, Davis et al. (2015) concluded that the gas in most of ETGs with dust lanes is of external origin, because of the large range of gas-to-dust ratio extending to high values (∼ 800), typical for dwarf galaxies, which have apparently merged with gas-poor ETGs (see also -Gas-to-dust ratios as a function of stellar mass.Large symbols denotes galaxies in our ETG sample, and small symbols denote those from HRS. Red squares, blue stars, and black circles denote molecular gas, atomic gas, and total gas, respectively (the latter only for galaxies with both CO and H I measurements).ETGs in our sample have high gas-to-dust ratios, but consistent with the HRS population arguing against external source of the ISM by merging with dwarf galaxies, which have high gas-to-dust ratios.Lianou et al. 2016).For our sample the gas-to-dust ratios are smaller, suggesting high-metallicity ISM, and hence internal origin.We show the gas-to-dust ratio as a function of stellar mass in Fig. 9.We compare our galaxies with spirals from the Herschel Reference Survey (Boselli et al. 2010(Boselli et al. , 2014;;Bendo et al. 2012;Cortese et al. 2012aCortese et al. , 2014;;Ciesla et al. 2012Ciesla et al. , 2014)).Their dust is likely of internal origin (because they are currently star-forming) and the ETGs in our sample fully overlap with these spirals, which supports the internal origin of dust in these ETGs.On the other hand, dwarf galaxies exhibit much higher atomic gas-to-dust ratios of 300-10 000 (Grossi et al. 2010(Grossi et al. , 2016;;Cormier et al. 2014;Hunt et al. 2014Hunt et al. , 2015Hunt et al. , 2017)).

Stellar mass-age correlation
More massive galaxies are on average older (e.g.McDermid et al. 2015), so in principle the M dust /M * -age anti-correlation could be driven by the M * -age correlation.This scenario predicts that this anti-correlation should disappear if only a narrow range of stellar masses is analysed.
The sample analysed in this paper is too small to subdivide it in narrow bins of stellar mass, but this can be done for galaxies selected similarly by Leśniewska et al. (2023) in a much bigger field.Fig. A13 in the appendix shows that the M dust /M * -age anti-correlation persists even for very narrow ranges of stellar masses, which is inconsistent with it being driven by the M * -age correlation.The scatter in the stellar mass bins log(M * /M ⊙ ) = 10.5-10.6 and 10.6-10.7 is 0.39 and 0.36 dex, respectively, so it does not increase compared to the full sample with a scatter of 0.43 dex.

Selection bias
Our sample has been selected based on redshift, elliptical morphology and dust detection and is a 250 µm fluxlimited sample.If the M dust /M * -age decline was the a result of a selection bias, then this would need to imply that 1.The selection criteria should remove old very dustrich galaxies.
2. The selection criteria should remove young galaxies with dust content similar to that detected in older galaxies.
None of these biases are introduced by our selection criteria.If old galaxies with high dust content or young galaxies with dust content equally low as we see for older galaxies existed, we would have detected them.
8. energy source While observations favor a scenario involving the removal of the entire cold ISM or outflows, the question remains what is the source of energy responsible for the ISM removal.This cannot be connected with the current star formation, because the SFRs are low in the galaxies in our sample and the SFR level gets progressively lower, which would result in a flattening evolution of the ISM content, inconsistent with the observations.We consider planetary nebulae (PNe), SNe Type Ia, cosmic rays, hot halo gas, and AGNs.PNe or SNe Type Ia (either directly or through cosmic rays), or AGNs (if their duty cycle is low) are the most likely explanations.They are short lived sources of energy, so consecutive generations of stars need to constantly go through these phases to make these scenarios possible.
For all potential energy sources, we will discuss how the state of the cold gas changes.It can either be removed from the galaxy entirely, transformed into the warm ionized gas (temperature of 10 4−5 K) or hot ionized gas (temperature > 10 6 K).Given the lack of X-ray data, we do not have a constrain on the latter, but we calculated the amount of warm ionized gas from the Hα luminosity using eq. 1 of Pagotto et al. (2021), following their assumption of electron density of n e = 100 cm −3 .The resulting masses of warm ionised gas of 10 4−5 M ⊙ are much lower than the difference between the cold gas masses of the youngest and oldest galaxies analysed here of 10 9 M ⊙ .Hence, cold gas cannot be transformed into the warm ionised phase for these galaxies and can only be heated to much higher temperatures or removed by outflows.

Planetary nebulae
Low-mass stars during the main-sequence or AGB phases are not energetic enough to ionize the gas around them (which is required to explain the decline of both molecular and atomic gas).However, the post-AGB/PN phase is a possibility, usually referred to as a hot low-mass evolved star (HOLMES; Binette et al. 1994;Flores-Fajardo et al. 2011;Herpich et al. 2018).During this phase a star ejects its envelope with an expansion velocity of around 30 km s −1 , and this gas has the temperature of ∼ 10 4−5 K (e.g.Cuisinier et al. 1996;Milingo et al. 2002;Sharpee et al. 2007;Sahai & Chronopoulos 2010;Bohigas et al. 2015;Ali et al. 2015;Ali & Dopita 2019).This is enough to move gas to the warm ionized phase, but is not enough to expel gas from a galaxy or to heat it to very high temperatures of > 10 6 K.However, in environments with high stellar velocity dispersion, the interaction of PNe with the ambient gas can lead to heating to such high temperatures (Conroy et al. 2015).This scenario was invoked by Conroy et al. (2015) as a mechanism responsible for preventing star formation in quiescent galaxies.There are indeed indications that post-AGB stars are responsible for photoionization of gas in ETGs (Binette et al. 1994;Stasińska et al. 2008;Cid Fernandes et al. 2010, 2011;Sarzi et al. 2010;Kehrig et al. 2012;Yan & Blanton 2012;Papaderos et al. 2013;Singh et al. 2013;Gomes et al. 2016;Herpich et al. 2018).
We made simple calculations to estimate whether PNe are numerous enough in the galaxies in our sample to explain the gas decline.A galaxy with a total stellar mass of ∼ 10 11 M ⊙ has ∼ 1.9 × 10 10 stars with masses between 1-8 M ⊙ , which ended their lives as PNe.In our sample the total gas mass goes down from ∼ 10 10 at the age of 1 Gyr to ∼ 2 × 10 9 M ⊙ at 9 Gyr.Hence, one PN would need to remove ∼ 8 × 10 9 M ⊙ /1.9 × 10 10 = 0.4M ⊙ of gas.
This calculations does not take into account that some stars go through the PN phase before the period considered here for gas removal of 8 Gyr.At z = 0.1-0.3there are additional 2-4 Gyr after the Big Bang.If we make the most conservative assumption that all stars were created at the time of the Big Bang, then during this extra time stars with masses of more than 1.9-1.4M ⊙ have already went through the PN phase before our considered ISMremoval phase.This decreases the number of available PNe to 1.1-0.8× 10 10 and therefore increases the required gas removal to 0.7-1 M ⊙ per PN.
The largest (i.e.oldest) PNe or circumstellar envelopes of evolved stars have radii of the order of 1 pc (e.g.O'Dell et al. 2004;Sahai & Chronopoulos 2010;Matthews et al. 2015).For an ISM density of 20-50 cm −3 for the cold atomic medium (Ferrière 2001), this corresponds to ∼ 2-5 M ⊙ of swept-up gas.This is similar to our estimate of the required gas removal per PN, so they are numerous enough to be responsible for the detected gas decline.

Supernovae Type Ia
SNe Ia can ionize gas and destroy dust in the swept up part of the ISM (Li et al. 2020).The swept-up gas can reach temperatures of > 10 6 K (Ceverino & Klypin 2009;Hopkins et al. 2018a;Li et al. 2020), required to move it to the very hot gas phase.
Similarly as for PNe we estimate the required effectiveness of SN Ia.Using eq.19 and Table 1 of Andersen & Hjorth (2018) we estimate the rate of SN Ia for a galaxy with M * = 10 11 M ⊙ and SFR = 0-10 M ⊙ yr −1 to be 0.008-0.014yr −1 .For a period of 8 Gyr this corresponds to ∼ 1.1-6.4× 10 7 SNe Ia.Hence, one SN Ia would need to remove ∼ 730-130 M ⊙ of gas.
SNe Ia release around 10 51 erg of energy (Khokhlov et al. 1993).For a SN with such energy in the simulations of Yepes et al. (1997) ∼ 20-5000 M ⊙ of gas is ionized.The required strength of SN feedback for our sample is consistent with this range.We therefore find this mechanism feasible.
We do not consider core-collapse SNe here, because the lifetimes of their progenitors are very short, so their numbers are closely connected with recent star formation.Therefore their numbers decrease for galaxies with high ages, which makes them unable to explain the detected gas decline, similarly to astration.

Cosmic rays
The interaction with cosmic rays can lead to the destruction of dust particles by thermal evaporation and sputtering (Dwek & Arendt 1992), and also to the heating and ionization of gas (Hayakawa et al. 1961;Spitzer & Tomasko 1968;Ferrière 2001;Indriolo et al. 2009;Padovani et al. 2020;Gabici 2022).
The lifetime of cosmic rays is only 10 6−7 yr (Jokipii & Parker 1969;Garcia-Munoz et al. 1975, 1977;Jokipii 1976;Yanasak et al. 2001;Bisbas et al. 2015Bisbas et al. , 2017;;Gabici 2022).Hence, cosmic rays produced by core-collapse SNe are unlikely to affect the ISM of the galaxies in our sample, due to their low SFRs, resulting in a low number of currently exploding core-collapse SNe unable to ionize the amount of gas we detect.However, cosmic rays can also be produced by SNe Type Ia (Chan et al. 2019).
Simulations show that cosmic rays are capable of dispersing or launching gas clouds on the timescale of the order of 10 Myr (Brüggen & Scannapieco 2020).This is much shorter than the timescale over which the gas con-tent declines in our sample.Therefore, in order for this mechanism to be viable, large amounts of cold dust would need to be shielded in dense clouds against the cosmic ray influence, which would make the timescale much longer and consistent with our measurements.
It has been shown with simulations that cosmic rays can reduce SFRs of galaxies by a factor of five and the density of gas in galaxy centers by an order of magnitude (Hopkins et al. 2021b;Byrne et al. 2023).However, they do not heat gas to temperatures beyond 10 6 K, required here.Hence, the only viable mechanism to explain the decline in cold gas reservoir detected in our sample is the CR-driven outflow.
Summarizing, cosmic rays produced by SN Type Ia can also be responsible for the ISM decline we detect if they drive outflows.

Active galactic nuclei
AGN feedback has been invoked as a mechanism of quenching for massive galaxies and they can heat gas to very high temperatures (see the review of Fabian 2012).
As stated in Section 7.2, only five galaxies in our sample have > 3σ detection of the [O III] line, which is used as an AGN indicator.The median 2σ upper limit on the AGN luminosity, calculated as L AGN = 3500L [OIII] (Heckman et al. 2004) is 2 × 10 42 erg s −1 .This indicates a low level of current AGN activity.This is consistent with our analysis of emission lines of 2 000 dusty ETGs, where we found that only up to 15% exhibit line ratios typical for AGNs (Ryzhov et al. in prep.).
However, the timescale during which a typical AGN is active is much shorter than the Gyr timescale considered here (Novak et al. 2011;Hickox et al. 2014;Padovani et al. 2017).Therefore, the visibility of AGN in a given population depends mainly on the duty cycle (fraction of time during which the AGN is active).For a conservatively low Eddington ratio of 0.01, simulations predict a duty cycle of ∼ 1-10% (Novak et al. 2011, their Fig. 8), so in a sample fully dominated by AGNs, only this fraction is expected to show current AGN activity.For a higher Eddington ratios (stronger AGN activity) the expected detected fraction is even lower.These fractions are similar to what we obtain, so the AGN feedback can be a possible mechanism of removing the ISM in these galaxies.

Hot halo gas
Hot, X-ray-emitting gas is frequently found in elliptical galaxies with stellar masses similar to those in our sample (e.g.Sarzi et al. 2013;Su et al. 2015;Kim et al. 2019;Kokusho et al. 2019).The dust destruction timescales in hot (10 000 K or more) medium are estimated to be 10 4 −10 8 years (Draine & Salpeter 1979;Jones et al. 1994;Jones 2004;Micelotta et al. 2010;Bocchio et al. 2012;Hirashita et al. 2015;Hirashita & Nozawa 2017).This is much shorter than the dust and gas removal timescale measured here, in M19, and in Leśniewska et al. (2023).However, if the cold ISM is partially shielded from the influence of hot gas, then the destruction process would be slower and could be responsible for the decline we detect.Smercina et al. (2018) interpreted the dust decline for post-starburst galaxies as the effect of grain sputtering in a hot medium.However, they did not detect a gas decline, in contrast to the galaxies in our sample.Galliano et al. (2021) interpreted the low dust-to-gas ratios of ETGs as the result of dust grain sputtering in hot gas, but our sample exhibits much higher dust-to-gas ratios (compare Fig. 3 and their Fig. 8).
Finally, dust destruction in hot halo gas applies only to low-density medium (Bocchio et al. 2012).Hence, it is unlikely to be the main mechanism in the galaxies in our sample, given their substantial ISM masses.

conclusions
We present carbon monoxide (CO) and 21 cm hydrogen (H I) line observations of dusty ETGs and measure the removal of the cold interstellar medium (ISM).We find that all the cold ISM components (dust, molecular and atomic gas) decline at similar rates.This allows us to rule out a wide range of potential ISM removal mechanisms (including starburst-driven outflows, astration, a decline in the number of asymptotic giant branch stars), and artificial effects like stellar mass-age correlation, environmental influence, mergers, and selection bias, leaving ionization by evolved low-mass stars or ionization/outflows by supernovae Type Ia or active galactic nuclei as viable mechanisms.We also provide the support of internal origin of the detected ISM.Moreover, we find that the quenching of star formation in these galaxies cannot be explained by a reduction in gas amount alone, because the star formation rates decrease faster (at a time scale of about 1.8 Gyr) than the amount of cold gas (a timescale of 2.3 Gyr).Furthermore, the star formation efficiencyies of the ETGs (SFE ≡ SFR/M H2 ) are lower than those of star-forming galaxies, whereas their gas mass fractions (f H2 ≡ M H2 /M * ) are normal.This may be explained by the stabilization of gas against fragmentation, for example due to morphological quenching, turbulence, or magnetic fields.
M. Micha lowski) as well as support through a VILLUM FONDEN Young Investigator Grant (project number 25501)    -The same as Fig. 2, but showing the effect of outflows and astration.The dotted and dashed lines denote the molecular gas and dust mass evolution assuming that outflows (with a rate as measured in Fig. 8) and astration (with the SFR as measured in Fig. 1 of M19) are the only cause of the gas mass change, respectively.This shows that these two effects are too weak to explain the data.A6.This shows that some additional gas removal mechanism is necessary to explain the data.

Fig. 1 .
Fig.1.-Light-weightedstellar age as a function of cosmic time for simulated dusty ETGs(Nadolny et al., in prep.).The black solid line is the running average for given time, whereas the grey shaded region represent the 1σ range.The dashed line indicates the linear evolution when the stellar age is equal to the value of cosmic time.The time reaches negative values because of how the moment when time is equal zero was set, in order to have the time equal to the stellar age at the moment of observation.This figure shows that stellar ages increase linearly with time for dusty ETGs, so can be used as a time proxy.
dust / M star ) log(M H2 / M star ) log(M HI / M star ) Fit Scaled with dust/gas ratio

Fig
Fig.2.-Ratios of molecular (large red squares and arrows for upper limits), atomic (large blue stars) gas, and dust (black circles; M19) masses to stellar masses as a function of luminosity-weighted stellar age of ETGs galaxies detected by Herschel from the sample ofRowlands et al. (2012).The exponential fits to the gas-to-stellar and dust-to-stellar mass ratios (eq.6 and 7) are shown as solid lines, colored as the datapoints.The red and blue dotted lines denote the curve for the dust-to-stellar ratio shifted upwards by the measured average gas-to-dust ratios of log(M H 2 /M dust ) = 2.144 ± 0.018 and log(M HI /M dust ) = 2.610 ± 0.011.Small red squares and blue stars mark the

Fig
Fig.3.-Gas-to-dustratios as a function of luminosity-weighted stellar age.Red squares, blue stars, and black circles denote molecular gas, atomic gas, and total gas, respectively (the latter only for galaxies with both CO and H I measurements).Only very week trends are present indicating that the ISM components are affected by the same mechanism.

Fig
Fig.5.-Star formation rates as a function of molecular gas masses.Arrows denote upper limits.Open circles denote galaxies which are within the main sequence (see Fig.1 in M19).The solid line denotes the relation for star-forming galaxies (eq. 1 in Micha lowski et al. 2018).ETGs are below the relation for star-forming galaxies, indicating that their gas reservoir is ceasing to be able to form stars.
Fig. 6.-Star formation efficiency (SFE ≡ SFR/Mgas) as a function of luminosity-weighted stellar age.Red squared, blue stars, and Fig.7.-Molecular gas fraction (f H 2 ≡ M H 2 /M * ) as a function of stellar mass for the sample of ETGs in our sample (black circles and arrows).Open circles denote galaxies which are within the main sequence (see Fig.1 in M19).The red and blue symbols correspond to averages of mass-selected galaxies and only those on the main sequence, respectively(Saintonge et al. 2017).Gas fractions of ETGs are comparable or higher than those of star-forming galaxies, indicating that the ETGs do not stop forming stars due to lack of gas.

Fig. 8 .
Fig. 8.-Molecular outflow rate (eq. 5 of Fluetsch et al. 2019) as a function of luminosity-weighted stellar age.Circles denotes the measurement without taking into account the AGN contribution (because most of the galaxies in our sample have not been detected at [O III]).Arrows denote the molecular outflow rates with the maximum (2σ) contribution of AGN allowed by the [O III] data.The dashed line is an Fig.9.-Gas-to-dust ratios as a function of stellar mass.Large symbols denotes galaxies in our ETG sample, and small symbols denote those from HRS. Red squares, blue stars, and black circles denote molecular gas, atomic gas, and total gas, respectively (the latter only for Fig.A10.-CO(1-0)(left), CO(2-1) (middle), and H I (right) spectra from the IRAM30m and GBT observations.The vertical dotted lines show the velocity range over which the spectra were integrated in order to obtain line fluxes.
Fig. A12.-The same as Fig. 2, but showing the results of the models of Gall et al. (2011a) presented in Section 4. The left panel show the models with astration (with the SFR as measured in Fig. 1 of M19) as the only mechanism of ISM removal.In the right panel models with additional cold gas removal (gas heating or outflows) are shown.The model parameters are shown in TableA6.This shows that some additional gas removal mechanism is necessary to explain the data.

Fig
Fig. A13.-The same as Fig. 2, but for the larger sample of Leśniewska et al. (2023) only for a narrow range of stellar masses: log(M * /M ⊙ ) = 10.5-10.6 (left) and 10.6-10.7 (right).The solid black line shows the fit to the full sample of Leśniewska et al. (2023) and the dashed green line shows the fit to the sample in Micha lowski et al. (2019).This shows that the M dust /M * -age anti-correlation is not driven by the M * -age correlation, because it would disappear for narrow ranges of M * .

Table 2
CO fluxes and luminosities from the IRAM30m/EMIR observations.

Table 3
Observing log for the GBT/VEGAS observations.Note.aThe data not usable due to strong radio-frequency interference.

Table 4 H
I fluxes and luminosities from the GBT/VEGAS observations.

Table 5
Predictions of the mechanisms that can explain dust decline in Fig.2and their consistency with the currently available data.
. J.H. was supported by a VILLUM FONDEN Investigator grant (project number 16599).A.-L.T. acknowledges the support of the National Science Centre, Poland through the POLONEZ grant 2015/19/P/ST9/04010.T.T.T. has been supported by the Grant-in-Aid for Scientific Research (No. 17H01110 and 21H01128).P.B. acknowledges funding through the Spanish Government retraining plan 'María Zambrano 2021-2023' at the University of Alicante (ZAMBRANO22-04) This research was funded in whole or in part by National Science Centre, Poland (grant numbers: 2021/41/N/ST9/02662 and 2020/39/D/ST9/03078).For the purpose of Open Access, the author has applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.This article has been supported by the Polish National Agency for Academic Exchange under Grant No. PPI/APM/2018/1/00036/U/001.O.R. acknowledges the support of the National Science Centre, Poland through the grant 2022/01/4/ST9/00037.This work is based on observations carried out under project number 198-14, 62-15, and 174-15 with the IRAM 30m telescope.IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).The research leading to these results has received funding from the European Commission Seventh Framework Programme (FP/2007-2013) under grant agreement No 283393 (Ra-dioNet3).The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.This research has made use of the Tool for OPerations on Catalogues And Tables (TOPCAT; Taylor 2005): www.starlink.ac.uk/topcat/;SAOImage DS9, developed by Smithsonian Astrophysical Observatory (Joye & Mandel 2003); the NASA's Astrophysics Data System Bibliographic Services and the WebPlotDigitizer 22 of Ankit Rohatgi.