The EDGE-CALIFA survey: exploring the role of the molecular gas on the galaxy star formation quenching

Understanding how galaxies cease to form stars represents an outstanding challenge for galaxy evolution theories. This process of"star formation quenching"has been related to various causes, including Active Galactic Nuclei (AGN) activity, the influence of large-scale dynamics, and the environment in which galaxies live. In this paper, we present the first results from a follow-up of CALIFA survey galaxies with observations of molecular gas obtained with the APEX telescope. Together with EDGE survey CARMA observations, we collect $^{12}$CO observations that cover approximately one effective radius in 472 CALIFA galaxies. We observe that the deficit of galaxy star formation with respect to the star formation main sequence (SFMS) increases with the absence of molecular gas and with a reduced efficiency of conversion of molecular gas into stars, in line with results of other integrated studies. However, by dividing the sample into galaxies dominated by star formation and galaxies quenched in their centres (as indicated by the average value of the H$\alpha$ equivalent width), we find that this deficit increases sharply once a certain level of gas consumption is reached, indicating that different mechanisms drive separation from the SFMS in star-forming and quenched galaxies. Our results indicate that differences in the amount of molecular gas at a fixed stellar mass are the primary driver for the dispersion in the SFMS, and the most likely explanation for the start of star-formation quenching. However, once a galaxy is quenched, changes in star formation efficiency drive how much a retired galaxy separates in star formation rate from star-forming ones of similar masses. In other words, once a paucity of molecular gas has significantly reduced star formation, changes in the star formation efficiency are what drives a galaxy deeper into the red cloud, retiring it.


Introduction
The appearance of galaxies in the nearby Universe is largely shaped by their star formation activity. The cessation of star formation that accompanies the transformation of a blue, spiral galaxy into a "red-and-dead" elliptical is usually called "star formation quenching" (e.g., Faber et al. 2007). Quenching is generally associated with the shortage of the raw fuel that feeds the star formation: the cold gas, and in particular its molecular phase. In low-mass galaxies, the gas, which is weakly bound due to their shallow potential wells, can be promptly removed by stellar feedback (e.g., Dekel & Silk 1986). High-mass galaxies, instead, might require a more powerful way to disperse the gas, such as Active Galactic Nuclei (AGN) outflows (e.g., Lacerda et al. 2020). The AGN activity, by heating up the gas, can also block the accretion from the intergalactic medium causing "quenching by starvation" (e.g., Cicone et al. 2014). Likewise, dcolombo@mpifr-bonn.mpg.de the suppression of cold gas accretion can result from shockheating in dark matter halos with mass > 10 12 M ("halo quenching"; e.g., Dekel & Birnboim 2006). Small galaxies falling towards a galaxy cluster can have their gas removed by tidal stripping (e.g., Abadi et al. 1999), or through the interaction with hot intra-cluster medium (e.g., Moore et al. 1996) that prevents further accretion from the intergalactic medium ("environmental quenching"). Alternatively, galaxies can stop forming stars efficiently even if a substantial amount of gas is present. In this "morphological quenching" scenario the development of a bulge or a spheroid, together with dispersive forces such as shear, stabilises the galactic gaseous disk against collapse, preventing star formation (Martig et al. 2009). To fully disentangle the dominant quenching mechanisms, their time-scales, and their parameter dependencies requires the analysis of cold gas conditions for a statistically significant sample of galaxies.
Major nearby galaxy cold gas mapping surveys (Regan et al. 2001, Wilson et al. 2009, Rahman et al. 2011, Leroy et al. 2009, Donovan Meyer et al. 2013 Article number, page 1 of 12 arXiv:2009.08383v1 [astro-ph.GA] 17 Sep 2020 A&A proofs: manuscript no. Colombo_APEX-EDGE_SFM et al. 2019, Sun et al. 2018 have focused on observations of the molecular gas (through CO lines). Despite a few notable exceptions (e.g., Alatalo et al. 2013;Saintonge et al. 2017), these surveys observed mainly spiral or infrared-bright galaxies (i.e. galaxies with significant star formation) and have emphasized understanding of how star formation happens, rather than how it stops. This boils down to quantifying the relation between molecular gas and star formation rate (SFR), which appears nearly linear in nearby disks (Kennicutt 1998;Bigiel et al. 2008;Leroy et al. 2013;Lin et al. 2019). This relationship is often parametrized via the ratio between the SFR and the molecular gas mass, called "molecular star formation efficiency" (SFE=SFR/M mol = 1/τ dep ), where the inverse of the SFE is the "depletion time," τ dep . The depletion time indicates how much time is necessary to convert all the available molecular gas into stars at the current star formation rate. On kpc scales and in the disks of nearby, star-forming galaxies, τ dep is approximately constant around 1-2 Gyr (Bigiel et al. 2011;Rahman et al. 2012;Leroy et al. 2013;Utomo et al. 2017), and it appears to weakly correlate with many galactic properties such as stellar mass surface density or environmental hydrostatic pressure Rahman et al. 2012). Nevertheless, small but important deviations for a constant SFE have been noticed, which can be the first hints of star formation quenching. In some galaxies, the depletion time in the centres appear shorter Utomo et al. 2017) or longer ) with respect to their disks. These differences may correlate with the presence of a bar or with galaxy mergers ; see also Muraoka et al. 2019) and do not seem to be related to unaccounted variation in the CO-to-H 2 conversion factor Utomo et al. 2017). Spiral arm streaming motions have also been observed to lengthen depletion times (Meidt et al. 2013;Leroy et al. 2015).
Besides variation of the SFE within galaxies, differences in global SFE between galaxies have been explored more widely in the nearby Universe, thanks especially to molecular gas galaxyintegrated studies (see Saintonge et al. 2011, 2017, andreferences therein). Those studies are less expensive in terms of exposure time compared to resolved mapping studies and provide the opportunity to collect data for larger galactic samples. In particular, integrated samples provide access to the molecular gas content of galaxies below the "star formation main sequence" (SFMS or MS, i.e the locus of the star-forming galaxies in the SFR-stellar mass diagram e.g., Brinchmann et al. 2004;Whitaker et al. 2012;Renzini & Peng 2015;Cano-Díaz et al. 2016), that is, galaxies that are slowly shutting down their star formation (located in the so-called "green valley"; Salim et al. 2007), down to passive galaxies (or "red sequence" galaxies, as defined on a colour-magnitude diagram). In general, galaxies on the main sequence have longer depletion times than similar stellar mass galaxies located below the main sequence (Saintonge et al. 2016(Saintonge et al. , 2017. The reasons for this are unclear and might be due to a combination of effects. Barred and interacting galaxies show generally shorter global depletion times compared to other systems (Saintonge et al. 2012). Early-type galaxies show lower SFE compared to late-type objects (Davis et al. 2014) and some of the lowest values of SFE are observed in bulge-dominated galaxies (Saintonge et al. 2012). As star formation quenching is more often seen to happen inside-out (e.g., González Delgado et al. 2016), galaxy morphology and structural properties seem to play a role in modifying the SFE. On average, τ dep,mol appears to decrease moving from early-to latetype systems (Colombo et al. 2018), following the decrease in shear (Davis et al. 2014;Colombo et al. 2018) as described by the "morphological quenching" scenario (but see Koyama et al. 2019). The molecular depletion time is also seen to decrease with stellar surface density and increase with molecular gas velocity dispersion (Dey et al. 2019). This might indicate that τ dep,mol is longer for bulged systems, and where gas is less gravitationally bound (as the gas boundness is proportional to Σ mol /σ CO , i.e. the ratio between the molecular gas mass surface density and the CO velocity dispersion; see Leroy et al. 2015). The environment in which a galaxy lives might also be important: galaxies in clusters appear to have longer depletion time than group galaxies (e.g., Mok et al. 2016), possibly due to the turbulent pressure and additional heating induced by the cluster itself.
Nevertheless, most theories attribute star formation quenching to the absence of molecular gas, rather than to a less efficient conversion from gas to stars. This has been explored observationally mostly by integrated surveys (Genzel et al. 2015;Saintonge et al. 2016;Tacconi et al. 2018;Lin et al. 2017), which usually parameterised the shortage of gas through the molecular gas fraction, f mol = M mol /M * . This quantity seems to drop drastically for galaxies with M * > 10 10.5−11 M (Saintonge et al. 2017;Bolatto et al. 2017) and for redder objects (Saintonge et al. 2011). It appears also reduced in barred galaxies . This absence appears tentatively connected to the presence of an AGN in a galaxy (Saintonge et al. 2017); negative feedback due to the AGN may provide an efficient gas-removal mechanism, which is suggested by the fact that AGN-hosting galaxies are almost exclusively observed in the green valley (see Lacerda et al. 2020). Nevertheless, this is still matter of debate (e.g., Kirkpatrick et al. 2014;Rosario et al. 2018).
Despite several studies exploring the methods for quenching, there is not a clear conclusion as to whether the quenching is driven by the reduction in molecular gas content, a change in the star formation efficiency of the molecular gas or both effects. Furthermore, these studies have not assessed how these effects may change throughout the quenching process. To address these shortcomings, in this work, rather than examining causes for SFE and f mol variations, we use the 12 CO(1-0) maps from EDGE  in combination with new 12 CO(2-1) single-dish measurements to investigate whether SFE or f mol changes are the main cause of star formation quenching in the centre of more than 470 CALIFA (Sánchez et al. 2016a) galaxies at different quenching stages. The paper is structured as follows. Section 2 exposes the data used in this paper, while Section 3 describes the quantities derived from these data such as star formation rates (SFRs) and molecular gas masses (M mol ). The results of the analyses are shown in Section 4. Those results are discussed and summarised in Section 5. For the derived quantities we assume here a cosmology H 0 = 71 km s −1 Mpc −1 , Ω m =0.27, Ω Λ =0.73.

Sample and data
We collect a homogenised compilation of 472 galaxies with molecular gas measured in an aperture of diameter 26.3" (corresponding to the APEX beam at 230 GHz). This compilation comprises our new observations using APEX together with a reanalysis of CARMA observations acquired by the EDGE collaboration ). All galaxies were already covered by spatially resolved IFS observations by the CALIFA survey (Sánchez et al. 2012). The sample of 472 galaxies covers a wide range of galaxy parameters in terms of morphology (from E to Sm, including a few irregular galaxies), stellar masses (10 7.2 − 10 11.9 M ), and star-formation rates (10 −4.0 − 10 1.3 M yr −1 ). Thus, it is one of the first explorations of a large sample not systematically biased by selection. The CARMA sample is generally made up of galaxies concentrated along the SFMS (as it has been assembled considering 22µm bright WISE galaxies), while the APEX sample targets cover more uniformly the so-called "green valley" and "red sequence", as we can see in Fig. 1. Together, the CARMA and APEX samples provide good coverage of the full extended CALIFA sample. The APEX and CARMA samples do not overlap, meaning that objects observed by CARMA in 12 CO(1-0) have not been re-observed with APEX in 12 CO(2-1). Fig. 2 gives an example of the quality of the APEX data and the rich variety of the targets in our dataset. Spectra of APEX 12 CO(2-1) observations are illustrated in the first row. The centres of all the galaxies in this example are well detected in CO, but the continuum and Hα equivalent width (W Hα ) maps show that the objects are in three different phases of their evolution. On the left middle panel, the continuum map indicates that the galaxy (NGC0873) is entirely blue, therefore dominated by star formation. Along with colour, W Hα is also a faithful proxy for the star formation properties of the galaxies (e.g., Sánchez 2020). In particular, W Hα > 6Å is found in galaxy areas dominated by HII regions, while where W Hα < 6Å the galaxy is quenched, or dominated by other effects than recent star formation. Indeed, the bottom left panel of Fig 2 indicates that W Hα > 6Å almost everywhere and closely follows the continuum map information. In NGC0170 (middle column), instead, the star formation in the centre is fully quenched, as indicated by the median value of W Hα within the APEX beam aperture, W Hα,b < 6Å (see Section 3 for further details), and by the yellow colour of the continuum map in the central region. Nevertheless, on the outskirts, stars are forming, and this galaxy appears globally dominated by star formation (as suggested by the median W Hα across the full galaxy W Hα,g > 6Å). The last galaxy displayed in the figure, NGC7550, is fully retired, as suggested by W Hα < 6Å basically everywhere and by the yellow colour of the whole continuum map. The galaxy, however, still possesses a measurable amount of molecular gas. In the following we assume as "centrally starforming" galaxies the objects that show a median W Hα > 6Å within the APEX beam, and "centrally retired, quenched or quiescent" the targets where median W Hα < 6Å within the APEX beam.

CALIFA data
CALIFA is an integral field spectroscopy (IFS) optical survey that imaged more than 1000 galaxies (667 included in the data release 3, and 416 in the extended sample) using the PMAS/PPak integral field unit instrument mounted on the 3.5m telescope of the Calar Alto Observatory (Sánchez et al. 2012(Sánchez et al. , 2016aLacerda et al. 2020). The CALIFA sample is drawn from the Sloan Digital Sky Survey (SDSS, York et al. 2000) to reflect the present-day galaxy population (0.005< z <0.03) in a statistically meaningful manner (log(M * /[M ])=9.4-11.4; E to Sd morphologies, including irregulars, interacting, and mergers; Walcher et al. 2014;Barrera-Ballesteros et al. 2015). Here we consider the galaxies observed with the low-resolution (V500) setup which covers between 3745-7500 Å with a spectral resolution FWHM=6 Å. CALIFA datacubes possess a spatial resolution of FWHM∼ 2.5 arcsec (García-Benito et al. 2015). Given the limits on the redshift, CALIFA allows the study of galaxies on kpc-scale. Additionally, the maps extend beyond 2.5 R eff , covering most of the optical disks. Ionized gas and stellar continuum map properties have been obtained through the Pipe3D pipeline (Sánchez et al. 2016c,b). Pipe3D analyses the stellar population applying the GSD156 simple stellar populations (SSP) library (Cid Fernandes et al. 2013). A stellar population fit is performed to the spatially rebinned V-band datacubes in order to estimate a spaxelwise stellar population model. This model is used to calculate the stellar mass density value within each spaxel. The ionised gas datacube is then generated by subtracting the stellar population model from the original cube. Each of the 52 sets of emission line maps is performed, calculating flux intensity, centroid velocity, velocity dispersion, and equivalent width for every single spectrum.

APEX observations and survey goal
We observed the 12 CO(2-1) emission (rest frequency, ν12 CO(2˘1) = 230.538 GHz) from 296 galaxy centres and 39 off-centre positions. In this paper, we present only the centre observations. The project was carried out with the APEX 12 m sub-millimetre telescope (Güsten et al. 2006) in ON-OFF mode using the wobbler (which ensures stable baselines), and the PI230 receiver which operates in the 1.3 mm atmospheric window. Galaxies have been observed across two projects M9518A_130 and M9504A_104 (PI: D. Colombo) which allocated 180 and 205 hours in the summer and winter semesters 2019, respectively, for a total of approximately 385 hours which include calibrations, additional overheads, and further test observations. All galaxies have been drawn from the CALIFA extended sample, with the only requirement to be accessible by APEX, i.e. all galaxies in the sample have declination ≤ 30 • .
The APEX resolution at 230 GHz is 26.3 arcsec. The median ratio of the beam radius to the effective radius of the full sample of galaxies is 1.12, with an inter-quartile range of 0.60. These do not change much if we consider only the face-on targets (with an inclination less than 65 • ). This means that, on average, the APEX beam covers roughly half of the radial extent of the CAL-IFA maps (see Section 2.1 and Fig. 2).
Article number, page 3 of 12 APEX 12 CO(2-1) spectra of three observed galaxies in different quenching phases. Top row panels show the spectra for each galaxy in green, where the dotted line represents the observation σ RMS . In the top row panel titles, the name of the galaxy (as in the CALIFA database) and its morphology, CO signal-to-noise ratio (SNR), and the logarithmic ratio of the galaxy global star formation rate to the global star formation main sequence (∆SFMS g ) are shown. Middle row panels show continuum RGB images extracted from the CALIFA datacubes using u− (blue), g− (green) and r− (red) bands. In the bottom row, the W Hα maps are displayed in diverging red colours where W Hα < 6Å, while the diverging blues illustrate the part of the map where W Hα > 6Å. The colour maps are centred at W Hα = 6Å in logarithmic units, and this value is indicated as a black vertical line in the colour-bars. Black contours mark the 25, 50, and 75 percentiles of the log(F Hα ) distribution, previously masked at 3σ RMS . W Hα map is also masked below 3σ RMS of Hα flux map. In panel legends, the median Hα equivalent width across the whole map ( W Hα,g ) and within the beam aperture ( W Hα,b ) are presented. In the panels of the two bottom rows, the green circle shows the APEX beam (FWHM=26.2 arcsec at 230 GHz), while white ellipsoids indicate 1 and 2 R eff .
The survey is designed to reach a uniform rms of 2 mK (70 mJy) per δv = 30 km s −1 wide-channels. For several targets for which we have detections but low signal-to-noise ratio (SNR<3), we integrate longer to achieve a rms of 1 mK. These requirements allowed us to attain pointed observations 10× more sensitive than CARMA (in term of achievable minimum molecular gas mass surface density; see also Section 2.3) and thereby detect the CO line in even the most gas-poor galaxies, which constituted the main goal of our APEX observations. In particu-Article number, page 4 of 12 D. Colombo, S. Sanchez, A. Bolatto et al.: Molecular gas and star formation quenching in CALIFA galaxies lar, we obtain 207 CO detections (SNR≥ 3) for a detection rate of 70%, with 50% of the observed galaxies detected at 5σ RMS .
Data calibration and reduction of the APEX data have been performed using the "Grenoble Image and Line Data Analysis Software" (GILDAS 1 ) and "Continuum and Line Analysis Single-dish Software" (CLASS) package with which we fit and remove a linear baseline to each spectrum outside a window of 600 km s −1 centred on the the galaxy V LSR . Afterwards, we smooth the data to a common spectral resolution δv = 23 km s −1 . The final median rms from the full sample at 23 km s −1 is 2.2 mK, which corresponds to a median, 3σ RMS , M mol = 2.8 × 10 8 M (Σ mol ∼ 2.8 M pc −2 ) at the median distance of the sample of ∼ 67 Mpc and using a constant α CO(2−1) = 6.23 M (K km s −1 pc 2 ) −1 . Full details about survey specifics and data reduction will be presented in an upcoming survey paper.

CARMA data
In addition to the APEX data, in this paper, we also make use of the CARMA database which constitute the first 12 CO(1-0) (and 13 CO(1-0)) follow-up of CALIFA galaxies undertaken by the EDGE collaboration. A full description of the CARMA CO data is given in Bolatto et al. (2017); here we provide a brief summary. The EDGE-CALIFA collaboration originally mapped 177 infrared-bright CALIFA galaxies with CARMA Econfiguration. A sub-sample of 126 higher signal-to-noise galaxies were observed also in D-configuration; subsequently, for these galaxies, D+E combined cubes were produced. In this paper, we use the 12 CO(1-0) data of the 126 D+E galaxies and the remaining 51 E-configuration galaxies that were not followed-up with the D-array of CARMA. The data cubes were smoothed to a spectral resolution of 20 km s −1 . The final D+E galaxies have 4.5 arcsec resolution, while the E-configuration only galaxies show a 9 arcsec resolution. The average rms noise in the D+E galaxies is 38 mK per 20 km s −1 channel width. This dataset was extensively exploited in Utomo et al. (2017) (2020), exploring different aspects of the interconnection between the ionised and cold molecular gas phases in galaxies.

Derived quantities
In this paper, we use spatially unresolved (i.e. single beam) data from APEX as well as resolved data from CARMA and CALIFA. In order to make these data comparable and simulate the effect the APEX beam would have on CARMA and CALIFA maps, we introduce a "tapering" function W T , i.e. a bidimensional Gaussian, centred on the centre of the galaxy, with unitary amplitude and FWHM θ= θ 2 APEX − θ 2 CARMA,CALIFA , where the APEX beam FWHM θ APEX = 26.3 arcsec, while the CARMA beam FWHM θ CARMA = 4.5 arcsec or θ CARMA = 9 arcsec, for the D+E or E-only configuration, respectively; and θ CALIFA = 2.5 arcsec. Integrated quantities within the APEX beam aperture will be indicated with the sub-script "b" and are calculated by co-adding the pixels within the CARMA or CAL-IFA maps, previously multiplied by the Gaussian filter, W T . Average quantities within the APEX beam aperture are obtained using the weighted median of pixels values in the resolved maps 1 http://www.iram.fr/IRAMFR/GILDAS where the weights for each pixel are given by the "tapering" function, W T . This operation is equivalent to the convolution of the CALIFA property maps to the APEX beam size and sampling the result at the pointing centre of the APEX beam. Globally integrated quantities, denoted with the subscript "g", are measured by summing up all the pixels in a given map, without applying the Gaussian taper. Where no subscript is indicated, the quantity is computed spaxel or pixel-wise for CALIFA and CARMA maps, respectively.

CO luminosity from ON-OFF APEX observations
We derive the 12 CO(2-1) flux within a spectral window of 400 km s −1 centred on the systemic velocity of the galaxy which is derived from the stellar redshift. The CO line velocityintegrated flux is expressed by the following equation: where δv = 23 km s −1 is the data final channel width, using T mb = T * A /η mb , with η mb = 0.78 for the APEX beam efficiency at 230 GHz. S CO,b is converted to Jy km s −1 through a conversion factor between Kelvin and Jansky 2 of Jy/K=37.
The statistical error for the flux is given by: where σ RMS is the standard deviation of the flux variations measured in the first and last 20 line-free channels of each spectrum: i.e., measured in two spectral windows of 20 channels before and after the velocity range used to measured emission line intensities. Finally, W 50 is the full width at half maximum derived as W 50 = 8 log(2)σ v , where σ v is the second moment calculated in the spectral window selected to measure the emission line (i.e,. the 400 km s −1 range centred on the systemic velocity of a galaxy). The values of W 50 obtained in this way coincides well with the ones derived using the "two slopes method" presented by Springob et al. (2005) and used elsewhere (Saintonge et al. 2011(Saintonge et al. , 2017Cicone et al. 2017). For non-detected galaxies, or galaxies showing SNR< 3 we make use of CO,b to provide an upper limit for flux given by 3 CO,b , where we assume a constant W 50 = 200 km s −1 . This value was derived using the procedure described before for W 50 resulting from stacking the central APEX spectra of all the APEX galaxies. This CO flux upper limit serves as a basis for calculating all other CO properties of non-detected galaxies presented here.
From the CO flux, we derived the 12 CO(2-1) luminosity using equation 3 of Solomon et al. (1997): where D L is the luminosity distance in Mpc (derived from the stellar redshift, z), ν obs is the observed frequency of the emission line in the rest frame in GHz, and S CO,b is the CO velocityintegrated flux derived using equation 1 (but in Jy km s −1 ).

CO luminosity from CARMA data cubes
The 12 CO(1-0) CARMA observations in the original EDGE database comes as position-position-velocity data cubes. Therefore, we use the tapering function, W T , here. For the CARMA data, the CO velocity-integrated intensity within the APEX beam aperture is given by: where the summation runs on the bi-dimensional pixels of the integrated intensity map I CO of the whole galaxy. After this, the CO luminosity within the APEX beam aperture is calculated by: where δx, δy are the pixel sizes in pc, and z is the galaxy redshift. As in Section 3.1, we provide an upper limit for the CO flux-related quantities of the non-detected galaxies (SNR< 3) as 3 CO given by equation 2, where we assume a constant W 50 =200 km s −1 . For detected targets, the full width at half maximum, W 50 , of the line that enters in the calculation of the flux statistical error, CO , is obtained for the CARMA galaxies from the CO spectrum built using by the integrated flux in each channel map.

Molecular gas mass and α CO
The molecular gas mass follows from the CO luminosity by assuming a CO-to-H 2 conversion factor, α CO,b : For the CARMA CO(1 − 0) data we use α CO,b ≡ α CO(1−0),b ). However, for the APEX observations, an additional correction factor is required, with α CO,b = α CO(1−0),b /R 21 , where R 21 is the CO(2-1)/CO(1-0) ratio. This value is determined to be ∼ 0.7 in nearby galaxies Saintonge et al. 2017).
Given that our sample consists of a variety of galaxies often quite far from the star formation main sequence, here we assume a variable α CO(1−0) based on Bolatto et al. (2013), equation 31: where Z is the gas-phase metallicity relative to solar metallicity and Σ * is the stellar mass surface density measured at each pixel in the CALIFA data and γ = 0.5 where Σ * > 100 M pc −2 or γ = 0 otherwise. Unlike Bolatto et al. (2013), to avoid iterative solving here we simply assume that Σ total ≡ Σ * , since for our sample galaxies the gas mass surface density is generally one order of magnitude lower than the stellar mass surface density. Also, Σ 100 GMC , the Giant Molecular Cloud (GMC) molecular gas mass surface density in units of 100 M pc −2 , does not appear in equation 7 as we assume Σ 100 GMC = 1 here, considering that GMC molecular gas mass surface density inner regions of nearby galaxies and Milky Way is largely consistent with 100 M pc −2 (see Sun et al. 2018;Colombo et al. 2019). For galaxies where optical emission lines remain undetected we assume, we assume Z = 1. Details on the Z and Σ * calculation are given in Section 3.4. In equation 7, α CO is calculated across the whole CALIFA map; within the APEX beam aperture, we used the weighted median from the CO-to-H 2 conversion factor map, α CO,b , where the weights for each pixel are given by the "tapering" function, W T .
Using this method we derive a median µ α CO(2−1),b = 3.93 M (K km s −1 pc 2 ) −1 with an inter-quartile range σ α CO(2−1),b = 2.04 M (K km s −1 pc 2 ) −1 for the new dataset observed with APEX. In contrast, for the CARMA dataset we obtain µ α CO(1−0),b = 2.76 M (K km s −1 pc 2 ) −1 and σ α CO(1−0),b = 1.23 M (K km s −1 pc 2 ) −1 . Those values are few times lower than the canonical α CO(1−0) = 4.35 M (K km s −1 pc 2 ) −1 and α CO(2−1) = 6.21 M (K km s −1 pc 2 ) −1 of the Milky Way, possibly due to the fact that the stellar mass surface density in the centre of our sample galaxies (which extend to massive red sequence galaxies) is typically higher than in the Milky Way or generally star-forming galaxies, while the gas-phase metallicity of most of our galaxies is close to solar.

IFS-derived parameters
For the purpose of this work, we make use of both nebular lines as well as stellar continuum derived maps provided by CALIFA data. In particular we use Hα, Hβ, [OIII] λ5007, [NII] λ6583 flux maps, F Hα , F Hβ , F [OIII] , F [NII] , respectively; the Hα equivalent width, W Hα , and the stellar mass surface density maps.
We calculate the extinction-corrected star formation rate (SFR) spaxel-by-spaxel using the nebular extinction based on the Balmer decrement: where the coefficients K Hα = 2.53 and K Hβ = 3.61 follow the Cardelli et al. (1989) extinction curve (see also Catalán-Torrecilla et al. 2015).
The SFR is then computed as: as indicated by Kennicutt (1998) which assumes the Salpeter initial mass function (Salpeter 1955). To obtain the integrated SFR within the APEX beam aperture we coadded all spaxels where W Hα > 6Å. In regions where W Hα < 6Å the Hα flux is not due to recent star formation but is dominated by the old-stellar population or other effects Espinosa-Ponce et al. 2020). Additionally, we remove from the summation all spaxels where the ionisation is due to AGN, i.e. all the spaxels that fall above the BPT diagram (Baldwin et al. 1981) demarcation line given by Kewley et al. (2001) in their equation 5. We distinguish between the beam SFR, SFR b = i SFR i × W T,i , as the integrated SFR within the APEX beam aperture, and the global SFR, SFR g = i SFR i , as the co-addition of pixels over the whole map.
The beam stellar mass, M * ,b = i M * ,i × W T,i , is obtained in the same way by the summation of the stellar masses from the spaxels within the APEX beam aperture, while the global stellar mass, M * ,glob = i M * ,i , is given by the summation over the whole map. In this formulae, the index i runs over the spaxels of the whole maps, and the quantities without sub-script indicate the respective CALIFA data-derived maps. D. Colombo, S. Sanchez, A. Bolatto et al.: Molecular gas and star formation quenching in CALIFA galaxies To calculate α CO within the APEX beam aperture we measure the gas-phase metallicity over the CALIFA maps using the O3N2 method (Marino et al. 2013).
As before the median gas-phase metallicity within the APEX beam aperture is calculated assuming the weights given by the Gaussian filter, W T . Using this method we obtain a median of 8.43 dex from the full galaxy sample with an inter-quartile range of 0.07 dex. The median metallicity with respect to the Solar metallicity (8.69; Allende Prieto et al. 2001 Lastly, we calculate the star formation efficiency (SFE) and the molecular gas mass fraction (with respect to the stellar mass) within the APEX beam aperture, respectively, using:

The SFR -M * diagram
Fig. 3 displays the distribution of galaxies across the SFR-M * diagram using SFR and M * measurements integrated over the entire CALIFA field of view (FoV), SFR g and M * ,g , respectively.
Each panel illustrates a different quantity calculated over the APEX beam aperture. Quantities shown include: a) equivalent width of Hα ( W Hα,b ), b) molecular gas mass (M mol,b ), c) star formation efficiency (SFE b ), and d) the ratio of molecular gas mass to stellar mass (molecular gas mass fraction f mol,b ), respectively. On average, the CALIFA maps extend to approximately 2 R eff , while the APEX beam covers the inner 1 R eff of the galaxies (see Fig. 2). The measurements within the APEX beam aperture can be considered, therefore, as inner galaxy measurements.
In addition, we include the SFMS fit derived by Cano-Díaz et al.
(2016) (log(S FR) = (0.81 ± 0.02) log(M * ) − (8.34 ± 0.19)), to illustrate the location of the star-forming galaxies across this diagram. In the following, we will provide the Pearson and Spearman correlation coefficients r as r p and r s , respectively. For all correlation discussed in the paragraphs we obtain extremely low p−value<< 10 −16 .
Global star-forming galaxies can be largely separated from galaxies on the way to quenching using a threshold given by W Hα,b = 6 Å, i.e. by considering the quenching stage of their centre. Panel a of Fig. 3, where the data points are colourencoded by W Hα,b , shows that most galaxies where this value is above 6 Å are tightly distributed along the SFMS locus defined by the Cano-Díaz et al. (2016) fit. We consider the lower boundary of the SFMS to be 3σ (∼ 0.6 dex) below the fit, which encompasses ∼ 99% of the galaxies dominated by star-formation in their centre (which show similar medians W Hα,b and W Hα,g around 13Å). However, we still observe ∼ 25% of the centrally retired galaxy sub-sample (having W Hα,b < 6 Å) above this line. Those galaxies have median W Hα,b ∼ 4.3 Å and median W Hα,g ∼ 5.3 Å, quite close to our threshold of 6 Å and are mostly star-forming in their outskirts but not in their centres.
It is interesting to note that quenched objects still possess a significant amount of molecular gas in their centres (Fig. 3,  panel b). In particular, we see objects with high or average values of this quantity well within the "green valley" and the "red sequence"; for example 22 galaxies with M mol,b values above the 75 th percentile of its distribution are below 3σ from the SFMS fit. At the opposite extreme, 46 objects with M mol,b values below the 25 th percentile of the M mol,b distribution are above the dashed green line defining galaxies within 3σ of the SFMS fit. Fig. 3  (panel b) shows that M mol,b is directly related to the global SFR and M * , but correlates more tightly with SFR g (given r p = 0.65 and r s = 0.73) than with M * ,g (r p = 0.46 and r s = 0.33).
The SFE in the centre of the galaxies varies widely across the SFR-M * diagram. However, it appears roughly constant along the SFMS, but it drops sharply right below it (panel c). Quantitatively, µ SFE b = 8.85 × 10 −10 yr −1 and σ SFE b = 0.10 × 10 −10 yr −1 above the "green valley" boundary, while below it µ SFE b = 1.80 × 10 −10 yr −1 and σ SFE b = 0.62 × 10 −10 yr −1 . In other words, galaxies below the main sequence show, on average, depletion times a factor 5 longer than galaxies across the SFMS in their central regions. Additionally, we have a few galaxies that do not appear to form stars in their centre (R < R eff ), resulting in SFE∼ 0 in the aperture defined by the APEX beam. Those are galaxies with Hα flux below the detection limit, or their Hα map spaxels are masked since W Hα < 6 Å. Nevertheless, most of these galaxies do have a few star-forming regions in the outskirts, typically outside the circle defined by their R eff (such as NGC0171 in Fig. 2). Some of these galaxies are non-detected in CO (9 targets); therefore the absence of recent star formation could be directly attributed to the absence of molecular gas (in the matched aperture). However, some of the SFE∼ 0 galaxies are detected in CO (5 targets), meaning that the star formation quenching is the centre is due to causes other than a simple shortage of raw fuel. As for M mol,b , we also observe a few galaxies with SFE b much lower than the sample average very close to the main sequence (5 targets with SFE b below the 25 th percentile of the SFE distribution are above the 3σ line). SFE b decreases with stellar mass (as observed also in Colombo et al. 2018 with a limited sample of EDGE galaxies; see their Fig. 3, panel d). The SFE in the galaxy centres appear quite strongly correlated with SFR g (given r p = 0.73 and r s = 0.46), but only moderately with M * ,g (r p = −0.38 and r s = −0.53). The SFE b reaches very low values, however, those values are generally SFE b lower limits, being mostly driven by CO non-detected galaxies, for which we use an M mol,b upper limit to calculate SFE b .
The ratio of the molecular gas mass to stellar mass f mol,b shows a general behaviour across the SFR-M * diagram quite similar to SFE b (panel d). As for SFE b , f mol,b is largely constant along the SFMS and sharply decreases below from the SFMS. Quantitatively, µ f mol,b = 2 × 10 −2 and σ f mol,b = 2 × 10 −2 above the "green valley" boundary, and µ f mol,b = 5 × 10 −3 and σ f mol,b = 9 × 10 −3 below it. We observe galaxies along and below the main sequence with f mol,b values much lower and much higher than the sample average. In particular we have 14 targets with f mol,b below the 25 th percentile of the distribution above the 3σ line (dashed green line) from the SFMS fit, while 5 objects with f mol,b values above the 75 th are located below this line. Nevertheless, the ratio of molecular gas mass to stellar mass appears only moderately correlated with SFR g (r p = 0.50 and r s = 0.50) and M * ,g (r p = −0.44 and r s = −0.50).

What quenches galaxies: variable SFE or shortage of molecular gas?
Star formation quenching can be parameterised using the logarithmic difference between the observed SFR and the SFR expected from the best fit to the SFMS, ∆SFMS (e.g., Genzel et al. 2015;Tacconi et al. 2018;Ellison et al. 2018;Thorp et al. 2019). In panels a and b of Fig. 4 we plot ∆SFMS g with respect to SFE b and f mol,b , respectively, in order to understand whether the star formation quenching is more tightly connected to variations in SFE or to the absence of molecular gas in galaxy centres. This is basically a reorganisation of the star formation-mass diagram presented in Fig. 3, removing the dependence of SFR on M * (i.e., the SFMS trend). As before, SFE b and f mol,b are measured within the APEX beam aperture, while ∆SFMS g uses the SFR and M * measured over the entire CALIFA map. Following the arguments discussed for Fig. 3 (panel a), we divide the sample into two sub-samples based on the average W Hα within the APEX beam aperture ( W Hα,b ): galaxies largely quenched in the centre ( W Hα,b < 6Å) and galaxies dominated by star formation in their centres( W Hα,b > 6Å). The two sub-samples are well balanced in terms of target size. Centrally star-forming galaxies number 256, i.e. ∼ 54% of the full sample, while the centrally retired galaxies constitute ∼ 46% of the sample, i.e. 216 objects. The behaviour of ∆SFMS g versus SFE b is somehow similar for the two sub-samples. The ∆SFMS g -SFE b relationship measured using the principal component analysis (PCA, see Colombo et al. 2018) shows that the slope from the confidence ellipsoids between the two sub-samples is on the same order ∼ 0.3 for centrally star-forming galaxies and ∼ 0.6 for centrally quenched galaxies (see Table 1, where also Spearman correlation coefficient r s for the two sub-samples are reported).
Nevertheless, data points for star-forming galaxies are tightly concentrated close to the SFMS and have SFE b values between 10 −10 −10 −8 yr −1 (i.e. τ dep = 0.1−10 Gyr). By contrast, quenched galaxies cover a much larger parameter space in both ∆SFMS g and SFE b , in particular, they span 6 orders of magnitudes in SFE. Additionally, ∆SFMS g and SFE b appear strongly correlated, showing a Pearson r p = 0.9. However, this tight correlation is mostly driven by the centrally quenched galaxies, for which r p = 0.8, while for the star-forming targets r p = 0.2, which indicates that ∆SFMS g and SFE b are basically uncorrelated for this kind of object and the calculated slope of the relationship is meaningless. SFE in our centrally star-forming galaxies is quite constant, in line with results from several other resolved and unresolved studies of nearby, star-forming galaxies. Note also that the values of the correlation coefficients do not change significantly for the ∆SFMS g -SFE b and ∆SFMS g f mol,b relationships if only the detected targets are considered (Table 1).
On the other hand, the slopes of relationship between ∆SFMS g and f mol,b are starkly different if we consider the starforming and the quenched targets separately. Galaxies largely quenched in the centre span a few orders of magnitude in ∆SFMS g as well as in f mol,b . Indeed, the PCA shows that in quenched galaxies there is a steep relationship between ∆SFMS g and f mol,b (with a slope of ∼ 3.15). However, this correlation is shallower for galaxies with central star-formation activity (slope ∼ 0.64). They have f mol,b approximately one order of magnitude larger than for centrally quenched galaxies ( log( f mol,b ) = −1.77 for centrally star-forming and log( f mol,b ) = −2.51 for centrally retired objects; see Table 1).
Nevertheless, the Pearson correlation coefficients are lower with respect to the SFE case: for the full sample, we observe r p = 0.7, while for the two sub-samples separately we observe a similar r p ∼ 0.5. This indicates that, in contrast to the SFE case, for both centrally star-forming and quenched targets ∆SFMS g and f mol,b are moderately correlated. Those conclusions do not change significantly if only the CO detected galaxies are considered (see Fig. 4 and Table 1).
It is worth noting that the SFE b exhibits a bimodal distribution similar to the one found for ∆SFMS g , i.e. values of SFE b below 10 −10 yr −1 are almost exclusively associated with quenched targets. By contrast, a bimodal distribution is not evident for f mol,b , for which the difference in f mol,b between centrally star-forming and retired objects is not as sharp. Additionally, the bimodality in SFE b and ∆SFMS g is driven by the same group of galaxies. In other words, centrally star-forming and retired galaxy sub-groups are equally well separated in ∆SFMS g and in SFE b . Thus, SFE in the galaxy centres (in particular retired centres) is a better predictor of the separation between the two groups than the respective f mol,b .
Panels c, d, and e show the histograms of ∆SFMS g , SFE b , and f mol,b colour-coded by W Hα,b in a given bin. Generally,  Table 1. Summary of the PCA fit of the ∆SFMS g -log(SFE b ) and ∆SFMS g -log( f mol,b ) relationships for centrally star-forming ( W Hα,b > 6 Å) and centrally retired ( W Hα,b < 6 Å) (or globally star-forming, W Hα,g > 6 Å, and globally retired, W Hα,g < 6 Å) for different galaxy sub-samples or quantity calculations. "Best" indicates the whole galaxy datasets with molecular gas masses calculated using a variable α CO , which results are shown in Fig. 4; for "SNR>3" only the detections are considered; "Const. α CO , represents the whole galaxy sample, where a constant α CO is used to convert CO luminosities into molecular gas masses; "Only APEX" marks the fit results when only the APEX data are used. In the columns, m and q indicate slope and intercept of the relations, respectively, inferred from PCA; ∆SFMS g , log(SFE b ) , and log( f mol,b ) , show the medians of global ∆SFMS , beam SFE, and beam f mol of the distributions, respectively; r p and r s are the Pearson and Spearman correlation coefficients, respectively. For most of the correlation realizations we obtain p−values largely below 10 −5 from both correlation tests, except for ∆SFMS g -SFE b relations for the centrally star-forming galaxies, for which we measure Pearson p−values of the order of 10 −2 and Spearman p−values of the order of 10 −1 . The uncertainties are obtained by 1000 bootstrap iterations of the PCA fit, and are provided as 75 th -50 th percentiles and 50 th -25 th percentiles of the m and q distributions. the median of global ∆SFMS g (−0.3) and beam SFE b (∼ 4.4 × 10 −9 yr) are close to the values of these parameters that separate star-forming and quenched galaxies (i.e., ∆SFMS(W Hα = 6Å) and SFE(W Hα = 6Å)). In particular, the median SFE corresponds to a τ dep = 2.3 Gyr, which is equivalent to the value measured from kpc-resolved EDGE objects (see Utomo et al. 2017;Colombo et al. 2018) and other nearby spiral galaxies ). However, the median of the beam f mol,b (∼ 10 −2 ) distribution is shifted towards the retired sub-sample as this value is slightly below the demarcation f mol,b that separates centrally star-forming and quiescent galaxies ( f mol,b (W Hα = 6Å)=10 −1.95 ).

Discussion and conclusions
In this paper, we use 472 galaxies to test whether the star formation quenching of CALIFA galaxies is mostly due to changes in the SFE or to the absence of molecular gas (as described by the ratio between the molecular and stellar gas masses, f mol ) in their centres. We observe that for galaxies dominated by star formation activity in their centre, distance from the main sequence correlates better with the molecular to stellar mass ratio. For centrally quiescent galaxies, instead, distance from the main sequence cor-relates better with SFE. This suggests a scenario where the progressive loss of the cold gas reservoir is what causes galaxies to move out of the main sequence. Once this happens, the star formation efficiency in the remaining cold gas reservoir is what modulates their retirement, with lower efficiencies corresponding to more quiescent galaxies. In this scenario both amount of (molecular) gas and SFE matter, but they have different roles. In particular, the stabilisation of the molecular gas reservoir plays a role once the galaxy enters the green valley, but it is less important than the size of the reservoir to move the galaxy out of the main sequence. Furthermore, this quenching happens from the inside-out, with centrally quenched galaxies leading the path towards totally quenched ones.
Those results do not change significantly if we consider a constant α CO instead of our preferred α CO from Eq. 7 in converting CO luminosity to molecular gas mass, or if we divide the full sample using the value of the Hα equivalent width obtained over the full maps, or if only the APEX sub-sample of galaxies is used (see Table 1).
The importance of the absence of molecular gas for understanding why some galaxies are located far from the SFMS has been acknowledged in the past from other (integrated, but aperture limited) studies that used a direct molecular gas tracer as notice that changes in the total gas fraction (calculated including the contribution of the atomic gas) are more significant than the total SFE in explaining distance from the resolved SFMS for star-forming galaxies, as we observe here using integrated measurements.
Integrated CO surveys cannot reach the level of detail regarding the molecular gas organisation achieved by kpc-resolved studies. But they do provide the ability to obtain samples that are several times larger and to detect galaxies with much less molecular gas using much less observing time. In this paper, we give the first presentation of a new integrated CO survey using APEX that follows up CALIFA targets. Once completed, this survey will give 12 CO(2-1) (and possibly also 13 CO(2-1) and C 18 O(2-1) for the brightest targets) observations of 450 CALIFA galaxy centres and a few off-centre detections. Thus, the size of this survey is similar to that of the most recent explorations at redshift ∼ 0, like xCOLDGASS (532 galaxies), but with aperturematched optical spectroscopic data (not restricted to the central 3", which could cause several issues in the classification of the ionisation stages, Sánchez 2020), and for a much narrower range of cosmological distances (i.e., with less bias introduced by possible cosmological evolution). The survey is unbiased by construction, having as its only requirement that the targets are observable by APEX (δ < 30 • ). Together with CARMA data, we will collect CO data for ∼ 630 galaxies fully covered by highresolution IFS information, providing the largest CO database of any major IFS survey to date.
Nonetheless, to fully exploit the IFS information and take the next step in understanding the mechanisms that drive these galaxy changes requires high-resolution interferometric gas imaging of the sample, something that needs to be strongly supported by proper time allocation.
ALMA would be particularly appropriate for this scope. Within our centrally quenched sample, we measure a median molecular gas mass upper limit M mol ∼ 10 8 M , within a 26.3 arcsec APEX beam. A short (∼ 20 min) integration on CO(2-1) emission over the full disk of a CALIFA galaxy with ALMA 12m array would provide a 1σ sensitivity of ∼ 8.7 mJy in 1 km s −1 channel. This is equivalent to Σ mol ∼ 2 M pc −2 for 30 km s −1 -wide lines, which would correspond to a M mol 10 6.2 M for the median distance to CALIFA galaxies, in a beam that matches the CALIFA resolution (∼ 3 arcsec). Therefore a short integration with ALMA has a 5σ limit of M mol ∼ 10 7 M (depending on the precise distance and line-width of the target), and would be able to improve significantly on our limits and potentially resolve the faint CO emission of a quenched galaxy. At these integration times, large surveys are possible, so about 300 CALIFA targets could be done in approximately 100 hours. This would provide a representative, invaluable high-resolution sample of galaxies to study the star formation quenching process in the local Universe.