Molecular gas during the post-starburst phase: low gas fractions in green-valley Seyfert post-starburst galaxies Free

Abstract

Post-starbursts (PSBs) are candidate for rapidly transitioning from starbursting to quiescent galaxies. We study the molecular gas evolution of PSBs at z ∼ 0.03–0.2. We undertook new CO (2–1) observations of 22 Seyfert PSB candidates using the Arizona Radio Observatory Submillimeter Telescope. This sample complements previous samples of PSBs by including green-valley PSBs with Seyfert-like emission, allowing us to analyse for the first time the molecular gas properties of 116 PSBs with a variety of AGN properties. The distribution of molecular gas to stellar mass fractions in PSBs is significantly different from normal star-forming galaxies in the CO Legacy Database (COLD) GASS survey. The combined samples of PSBs with Seyfert-like emission line ratios have a gas fraction distribution that is even more significantly different and is broader (∼0.03–0.3). Most of them have lower gas fractions than normal star-forming galaxies. We find a highly significant correlation between the WISE 12 and 4.6 μm flux ratios and molecular gas fractions in both PSBs and normal galaxies. We detect molecular gas in 27 per cent of our Seyfert PSBs. Taking into account the upper limits, the mean and the dispersion of the distribution of the gas fraction in our Seyfert PSB sample are much smaller (μ = 0.025, σ = 0.018) than previous samples of Seyfert PSBs or PSBs in general (μ ∼ 0.1–0.2, σ ∼ 0.1–0.2).

1 INTRODUCTION

Despite extensive observational and theoretical studies, the physical mechanisms that regulate the star formation rates (SFRs) of galaxies are still poorly understood. Star formation quenching, by yet unknown mechanisms, causes star-forming galaxies to migrate to the ‘red sequence’ (Gladders et al. 1998; Faber et al. 2007). One likely such formation of mechanism of the red sequence is the transformation of star-forming, disc-dominated, gas-rich galaxies into early types via mergers and their associated feedback (Toomre & Toomre 1972; Hopkins et al. 2006).

Post-starburst (PSB) galaxies are candidate post-merger objects, rapidly transitioning from the blue cloud to the red sequence (e.g. Dressler & Gunn 1983; Couch & Sharples 1987; Zabludoff et al. 1996; Wild et al. 2009; Snyder et al. 2011; Yesuf et al. 2014; Pawlik et al. 2016). In their quiescent phase, their spectra reveal little-to-no current star formation, but a substantial burst of star formation before an abrupt cessation ∼1 Gyr ago, long enough for the ionizing O&B stars to evolve away, but recent enough for A stars to dominate the stellar light. Recent observational efforts have enlarged the traditional definition of PSB to include a more complete and less biased sample of galaxies with ongoing star formation or/and AGN activity (Wild, Heckman & Charlot 2010; Alatalo et al. 2014; Yesuf et al. 2014). With the identification of this sample, which spans the entire starburst-to-quenched PSB evolutionary sequence, better constraints on theoretical models of galaxy evolution have started to emerge. One of the firm constraints is the observed time delay between the starburst phase and the AGN activity by about 200 Myr (Wild et al. 2010; Yesuf et al. 2014). This time delay suggests that AGNs are not primarily responsible for the original quenching of starbursts but may be responsible for keeping star formation at a low level by removing gas and dust during the PSB phase. The other emerging but contested phenomenon is molecular gas evolution after the starburst has ended (Schawinski et al. 2009; Rowlands et al. 2015). This work aims to further examine the molecular gas evolution along the starburst sequence using both new and existing data, by adding molecular gas observations from green-valley Seyfert PSBs, which were excluded from other PSB samples.

In simulations of gas-rich mergers, gas is funnelled to galaxy centres, powering intense nuclear starbursts and obscured nuclear AGN activity. At the end of the starburst, the leftover gas and dust are cleared out due to feedback from the AGN (e.g. Sanders et al. 1988; Barnes & Hernquist 1991; Silk & Rees 1998; Di Matteo, Springel & Hernquist 2005; Springel, Di Matteo & Hernquist 2005; Hopkins et al. 2006, 2008; Kaviraj et al. 2007; Wild et al. 2009; Snyder et al. 2011; Cen 2012). For instance, Narayanan et al. (2008) found that galactic winds are a natural consequence of merger-induced star formation and black hole growth. In their simulated galaxies, the galactic winds can entrain molecular gas of ∼10⁸–10⁹ M_⊙, which, the authors showed, should be observable in CO emission. The molecular gas entrained in the winds driven by AGN is predicted to be longer lived than the gas entrained solely in starburst-driven winds. The wind velocities in the simulated galaxies with AGN feedback can reach close to 2.5 times the circular velocity. Thus, making AGN feedback a viable mechanism to get rid of a residual gas and dust at the end of a starburst.

Despite its theoretical appeal, the evidence that connects AGN activity with the end of star formation in galaxies has been elusive, with evidence both for (Schawinski et al. 2009; Alatalo et al. 2011; Cicone et al. 2014; García-Burillo et al. 2014) and against (Fabello et al. 2011; French et al. 2015; Geréb et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016).

Now we review previous works on molecular gas contents of PSBs and AGN. Using Institut de Radioastronomie Millimétrique (IRAM) CO observations, Rowlands et al. (2015) investigated the evolution of molecular gas and dust properties in 11 PSBs on the starburst-to-quenched PSB sequence at z ∼ 0.03. Two of these PSBs are Seyfert galaxies while the rest are either star-forming or composite galaxies of star formation and AGN activity. 10/11 of the PSBs were detected in the CO (1–0) transition and 9/11 of the PSBs were detected in CO (2–1) transition. The gas and dust contents, the star formation efficiency and the gas depletion time of majority of these PSBs are similar to those of local star-forming spiral galaxies (Saintonge et al. 2011; Boselli et al. 2014) and gas-rich elliptical galaxies (Young et al. 2011; Davis et al. 2014). In addition, the authors found a decrease in dust temperature with the starburst age but they did not find evidence for dust heating by AGN at late times.

Similarly, Alatalo et al. (2016) studied 52 PSBs with shock signatures at z = 0.02–0.2 using IRAM and Combined Array for Research in Millimeter-Wave Astronomy (CARMA). About half of these PSBs are at z > 0.1, 14/52 are Seyferts and 90 per cent have CO (1–0) detections. The molecular gas properties of these PSBs are also similar to those of normal star-forming galaxies. More than 80 per cent of the PSBs in Alatalo et al. (2016) and Rowlands et al. (2015) samples are located in the blue cloud.

French et al. (2015) studied 32 PSBs in the green valley at z = 0.01–0.12 using the IRAM 30 m and the Submillimeter Telescope (SMT). Almost all of these PSBs have signatures mimicking low-ionization nuclear emission line regions (LINERs), and 53 per cent of them have CO detections. Those detected in CO have gas masses and gas-to-stellar mass fractions comparable to those of star-forming galaxies while the non-detected PSBs have gas fractions more consistent with those observed in quiescent galaxies.

The three aforementioned studies on the molecular gas contents of PSBs suggested that the end of starburst in these galaxies cannot be ascribed to a complete exhaustion or removal or destruction of molecular gas. The studies also suggested that multiple episodes of starburst or/and AGN activities may be needed for the eventual migration of these galaxies to the red sequence and that a transition time longer than 1 Gyr may be needed for this migration to take place.

Saintonge et al. (2012) found that, among the gas-rich, disc-dominated galaxy population, those that are ongoing mergers or are morphologically disturbed have the shortest molecular gas depletion times. They found no link between the presence of AGN and the long depletion times observed in bulge-dominated galaxies. Even though their AGN sample has lower molecular gas fractions than the control sample matched in NUV − r colour and stellar mass surface density, the depletion times of the two populations are similar. More than 90 per cent of the AGN studied by Saintonge et al. (2012) are not Seyferts, and instead are LINERs.

In contrast, Schawinski et al. (2009) presented evidence that AGN are responsible for the destruction of molecular gas in morphologically early-type galaxies at z ∼ 0.05. Their sample included 10 star-forming galaxies, 10 star formation and AGN composite galaxies and 4 Seyfert galaxies. The galaxies studied by Schawinski et al. (2009) are not PSBs but might have experienced mild recent star formation (see Schawinski et al. 2007). The authors found that the molecular gas mass drops significantly 200 Myr after a recent star formation in the composite galaxies, and none of their Seyferts have CO detections. The authors interpreted their observations as evidence for a destruction of molecular gas and for a suppression of residual star formation by low-luminosity AGN. Likewise, many studies have reported molecular outflows with high mass-outflow rates as evidence for AGN feedback in non-PSB AGN host galaxies (e.g. Feruglio et al. 2010; Fischer et al. 2010; Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Cicone et al. 2014; García-Burillo et al. 2014; Sun et al. 2014).

In this work, we study molecular gas in 22 green-valley, Seyfert PSBs candidates using SMT CO (2–1) observations, in combination with an existing sample of 94 PSBs from the literature. Our observations were motivated by French et al. (2015) and Rowlands et al. (2015) and were designed to be complementary to the samples in these two works. Our sample is also complementary to the recently published sample of shocked PSB galaxies (Alatalo et al. 2016), which contains 14 galaxies with Seyfert-like emission line ratios, which are mainly located in the blue cloud. Our sample is comparable in number to the existing sample of Seyfert PSBs but represents those in the UV–optical green valley. When combined with other samples, our sample is indispensable in sampling PSBs with a variety of AGN properties.

1.1 Simple energetic argument for no feedback or delayed AGN feedback after the molecular gas fraction is ≲4 per cent

In this section, we present the condition required for momentum-driven AGN wind to clear gas in a galactic disc or very close to the galaxy, following Silk & Nusser (2010). These authors argued, from the condition they derived, that AGN cannot supply enough momentum in radiation to unbind gas out of haloes of galaxies. They theoretically estimated the total gas mass within a dark matter halo assuming a gas-to-total mass fraction of 10 per cent. In this work, we measure the molecular gas mass-to-the stellar mass fraction, |$f_{\rm H_2} = \frac{M_{\rm H_2}}{M_{\star }}$|⁠, in the discs of our galaxies (or very close to the vicinities of the galaxies). Our estimated condition is local, while Silk & Nusser (2010) were interested if the AGN-driven winds eventually escape from the haloes of galaxies. We recast their derivation in terms of molecular gas fraction to show that AGN may only impart sufficient momentum to clear galactic discs once the molecular gas fraction is below about ≲4 per cent, perhaps explaining the delayed AGN molecular gas destruction hinted in this work, and the minor role of AGN in quenching starbursts (e.g. Wild et al. 2010; Yesuf et al. 2014). Note that starbursts have molecular gas fractions of 20–30 per cent while normal star-forming galaxies have gas fractions of about 10 per cent.

Let us assume an isothermal sphere galaxy. Within a halo radius r, its total mass is M_gal = 2σ²r/G, where G is Newton’s gravitational constant and σ is the velocity dispersion.

The above assumed values are reasonable for PSB galaxies. Similar assumption of R_max/R = 10 is also made by Rupke, Veilleux & Sanders (2002) in studying galactic winds in local ultraluminous infrared galaxies. The gas measurements are made within R ∼ 10–20 kpc of our galaxies while their dark matter haloes may extend to few 100 kpc. R_max/R may range between 10 and 100. Adopting R_max/R = 100 instead will only change our |$f_{\rm H_2}$| limit by a factor of 1.3. Only few per cent of the AGN luminosity is expected to couple to the gas. Therefore, the above lower limit is very liberal estimate since it assumes 100 per cent feedback efficiency. 5 per cent feedback efficiency is often used in AGN feedback models (e.g. Scannapieco & Oh 2004; Di Matteo et al. 2005; Zubovas & King 2012). Therefore, our estimate implies that that AGN feedback may only be effective in late-stage PSBs after the molecular gas fraction is below ≲4 per cent.

The rest of the paper is organized as follows: Section 2 presents the sample selection. Section 3 presents overview of statistical methods used in the paper. Section 4 presents the main results. Section 5 discusses our sample in comparison with other samples of PSBs. Section 6 summarizes the main findings of this work. Appendix A provides ancillary information on how our sample relates to existing samples of PSBs with molecular gas measurements. We assume |$(\Omega _{\rm m},\Omega _\Lambda ,h) = (0.3,0.7,0.7)$| cosmology.

2 SAMPLE SELECTION AND OBSERVATIONS

2.1 Sample selection

Using the cross-matched catalogue of Sloan Digital Sky Survey (SDSS), GALEX and WISE (Martin et al. 2005; Wright et al. 2010; Aihara et al. 2011; Yesuf et al. 2014), we select a sample of 22 transition PSB Seyfert galaxy candidates based on the evolutionary path that starbursts and PSBs follow in the dust-corrected NUV − g colour, Hδ absorption equivalent width and the 4000 Å break (Yesuf et al. 2014), and based on the BPT line ratio AGN diagnostic (Baldwin, Phillips & Terlevich 1981; Kewley et al. 2001; Kauffmann et al. 2003), as shown in Fig. 1. D_n(4000) probes the average temperature of the stars responsible for the continuum emission and is a good indicator of the mean stellar age (Bruzual 1983; Balogh et al. 1999). The combination NUV − g colour with Hδ and D_n(4000) is useful in identifying late-stage PSB galaxies, which are outliers at intermediate age from loci of galaxies with regular (continuous) star formation histories. As discussed in Yesuf et al. (2014), we acquired the measurements for physical parameters such as stellar masses and spectral indices from the publicly available catalogues on SDSS website¹ (Aihara et al. 2011). Our sample is restricted to galaxies with redshift, z < 0.06, stellar mass, M >10¹⁰ M_⊙, (NUV − g)_dc > 2.2 and D_n(4000)_dc < 1.55 and Hδ > 1 Å. The subscript ‘dc’ denotes dust correction. We use the Hα/Hβ flux ratio with two-component dust attenuation model of Charlot & Fall (2000) to correct for dust attenuation of the nebular emission lines, and the empirical relationship between the emission line and continuum optical depths found in Wild et al. (2011) to correct the continuum fluxes (for more details, see Yesuf et al. 2014). Two candidate PSB galaxies with CO observations are removed from the main analyses because they do not satisfy the above cuts. TPSB12 is removed because it has Hδ < 1 Å while TPSB13 is removed because it has (NUV − g)_dc < 2 and it is also an ongoing merger. Including these two PSBs instead does not change the main conclusions. Furthermore, note that in Fig. 1 we only show galaxies with NUV detections (Yesuf et al. 2014) in previous samples of PSBs (French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016). In later analyses, this restriction is not required.

The redshift cut, z < 0.06, was imposed due to the sensitivity of the SMT to achieve the desired signal-to-noise ratio in 6–8 h. This severely limited the number of Seyferts available for the observation. In addition, in the second year of the observation, the sample was restricted to be above a declination of 35° (away from the Sun avoidance zone for the SMT and the bulge of the Galaxy), further limiting the observable sample. Therefore, we did not impose the Hδ > 4 Å absorption cut that was done in the original sample selection of Yesuf et al. (2014). Most (19/22) of the galaxies in our sample have, within the measurement errors, Hδ absorptions ≳3 Å, which is expected in PSBs with weak or strong bursts or rapidly truncated star-forming galaxies (e.g. Poggianti et al. 1999). Our sample on average has lower Hδ than do the aforementioned previous samples. Our PSBs are consistent with being later stage PSBs or PSBs with weaker starbursts compared to most of the PSBs in the other samples considered here. Note that the Hδ equivalent width is ∼2 Å at 1 Gyr after a starburst with star formation time-scale τ = 0.1 Gyr and burst mass fraction b_f = 20 per cent. For starburst with τ = 0.1 Gyr and b_f = 3 per cent, Hδ ∼ 0 at 1 Gyr (see Fig. A1 in Appendix A). Thus, not imposing Hδ > 4 Å cut helps select old PSBs. However, it remains uncertain whether these galaxies underwent weak or strong starbursts or rapidly truncated star formation.

It should be noted that at z < 0.06, the SDSS fibre covers only ∼3.5 kpc of the central region of a galaxy, and the spectroscopic measurements may not reflect the galaxy-wide values. On the other hand, the NUV − g colour is an integrated galaxy-wide quantity. By combining the spectroscopic measurements with NUV − g colour, we select only (PSB) galaxies that are fading galaxy wide.

2.2 SMT CO observations

The observations were carried out using the Arizona Radio Observatory SMT on Mount Graham, Arizona. The observing runs were performed in 2015 February 25–March 10 and in 2016 March 4–25. We follow the same instrument set-up and observing strategy as French et al. (2015). Namely, we used the 1 mm Atacama Large Millimeter Array (ALMA) Band 6 dual polarization sideband separating SIS (superconductor–insulator–superconductor) receiver and 1 mHz filterbank to measure the CO (2–1) 230.5 GHz emission line. The beam size of the SMT for this line is about 33 arcsec. Beam switching was done with the secondary at 2.5 Hz switching rate and a throw of 120 arcsec, in the BSP (beam switching plus position switching) mode. Calibration using a hot load and the standard chopper wheel method was performed every 6 min. Calibration using a cold load was performed at every tuning. The observing times range between 4 and 9 h. French et al. (2015) observed a subset of 13 PSBs using both SMT and IRAM, which we use to guide comparisons between SMT- and IRAM-observed samples.

The data reduction is done using class, a program within the gildas software package.² The main beam efficiency is calculated using Jupiter in each polarization. A first-order polynomial baseline is subtracted from the spectra using data between [−600, 600] km s⁻¹, excluding the central regions of [−300, 300] km s⁻¹. The spectra are scaled using the main beam efficiency, and are co-added by weighting with the root-mean-square (rms) noise. The spectra are rebinned to 14 km s⁻¹ velocity bins. The typical rms error per bin is 1–2 mK. Thus, we achieve similar sensitivity as previous works (Saintonge et al. 2011; French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016).

To calculate the integrated CO line intensity, I_CO, we fit a Gaussian profile to each line, allowing the peak velocity to differ from the optical systemic velocity by up to 200 km s⁻¹. The statistical uncertainty in the line intensity is calculated following Young et al. (2011) as |$\sigma _{\rm I}= (\Delta v)^2\sigma ^2N_{\rm I} (1+\frac{N_{\rm I}}{N_{\rm b}})$|⁠, where Δv is the channel velocity width, σ is the channel rms noise, N_I is the number of channels used to integrate over the line and N_b is the number of channels used to fit the baseline. When the line is not detected, the upper limits of the line intensity are calculated as three times the statistical uncertainty. Following Solomon et al. (1997), the CO line luminosity in K km s⁻¹ pc² is |$L^\prime _{\rm CO} = 23.5\, \Omega _{s*b} D_{\rm L}^2 I_{\rm CO} (1+z)^{-3}$|⁠, where Ω_s★b is the solid angle of the source convolved with the beam, z is the redshift from the SDSS optical spectrum and D_L is the luminosity distance in Mpc. If the source is much smaller than the beam, then Ω_s★b ≈ Ω_b. We do not know the CO emitting sizes of our galaxies. We adopt this approximation as a simplifying assumption. We do not expect this to affect our conclusion. Previous works also adopt the same approximation. The fact that SMT has larger beam size than IRAM, and CO (2–1) is known to be more centrally concentrated than CO (1–0), makes the effect of this approximation less significant in this work compared to previous works. French et al. (2015) estimated for their sample that the |$L^\prime _{\rm CO}$| may be underestimated by about 1.5 × because of this effect.

The molecular gas mass can be calculated from |$L^\prime _{\rm CO}$| by assuming a CO conversion factor (mass-to-light ratio) α_CO, M(H|$_2)= \alpha _{ \rm CO} L^\prime _{\rm CO}$|⁠. We assume an α_CO = 4.3 M_⊙(K km s⁻¹pc²)⁻¹  in Milky Way disc (Bolatto, Wolfire & Leroy 2013) and |$R_{21}=L^\prime _{\rm CO(2{\rm -}1)}/L^\prime _{\rm CO(1{\rm -}0)}=1$|⁠. As discussed in Section 5.3, choosing a lower value of α_CO strengthens our main conclusion. The ratio R₂₁ is uncertain (Leroy et al. 2013; Sandstrom et al. 2013). Leroy et al. (2013) found R₂₁ ∼ 0.7 ± 0.3 in nearby disc galaxies. We adopt the Galactic α_CO value throughout the paper unless explicitly stated, and all previous measurements are adjusted to our assumed value of α_CO. The observations are summarized in Table 1, the co-added spectra are given in Figs 2 and 3 and the SDSS cut-out images of the Seyfert PSBs are shown in Fig. 4.

3 OVERVIEW OF STATISTICAL METHODS

In this section, we briefly review the statistical methods and tests used in the next section. Readers interested in more details of the methods should refer to the references provided.

In this section, we are concerned with how to analyse data that include both CO detections and CO upper limits, and make inference that is consistent with the incomplete information at hand. In statistics, a data set is called censored when the value of a measurement is only partially known to be above or below a threshold. Survival analysis is an area of statistics that mostly deals with modelling censored data such as time to an event (e.g. Feigelson & Nelson 1985; Klein & Moeschberger 2005; Halsel 2012). We use existing survival analysis methods to properly extract information from our censored CO observations.

Let T denote a positive random variable representing time to an event of interest, say death. The survival function is the probability that an individual survives beyond time t, S(t) = Pr(T ≥ t) = 1 − F(t), where F(t) is the cumulative distribution function. The Kaplan–Meier estimator (Kaplan & Meier 1958) is a non-parametric maximum likelihood estimator of the survival function, even in the presence of censoring. A non-parametric method does not require an assumption that the data follow a specific probability distribution. Instead, the data are ranked from smallest to largest, providing information on the relative positions of each observation. In the case of ties, Kaplan–Meier method assigns the smallest rank to each observation. The censored observations are also used in calculating the ranks. Let there be a set of |$\lbrace t_k\rbrace ^n_{k=1}$| data points and of these let |$t^\prime _1 < t^\prime _2 < t^\prime _3, \ldots , t^\prime _r$| be ranked, distinct uncensored values. At each time point |$t^\prime _j$|⁠, we observe d_j, the number of deaths, c_j, the number of censored observations between the time |$t^\prime _j$| and |$t^\prime _{j-1}$| and, n_j, the number of individuals at risk just prior to the time |$t^\prime _j$|⁠. In other words, n_j is the total sample size minus those who are censored or have died before |$t^\prime _j$|⁠, n_j = n − c_j − d_j = ∑_{l ≥ j}(c_l + d_l). The Kaplan–Meier estimator has a form: |$\hat{S}(t)= \prod _{j:t^\prime _j <t} (1-d_j/n_j)$|⁠. It is a step function with jumps at times |$t^\prime _j$|⁠.

To compute the empirical cumulative distribution function (ECDF), using the Kaplan–Meier method, we use the nada package in r programming language (Halsel 2012; Lee & Lee 2015). We use the function cenfit in the nada package to compute the Kaplan–Meier estimator for the CO data (Halsel 2012; Lee & Lee 2015).

The survival curves or ECDF of two groups can be compared using the log-rank test, which tests the null hypothesis that two groups have the same distribution against the alternative hypothesis that two groups have different distributions (Mantel 1966; Cox 1972; Harrington & Fleming 1982; Martinez 2007). If the data are uncensored, the log-rank test gives similar result to the Mann–Whitney test. We use the routine cendiff in the nada package to do the log-rank test.

The Kendall’s τ is a non-parametric correlation coefficient that can be used for testing trends in both censored and uncensored data (e.g. Halsel 2012). τ ranges between −1 and 1. τ = 0 is no correlation, τ = 1 is a perfect correlation and τ = −1 is a perfect anti-correlation.

In the uncensored case, the slope of a linear regression can be estimated using the Theil–Sen slope (Sen 1968), which is the median of all pairwise slopes between two data points. This slope results in the Kendall’s τ of 0 for the correlation between the residuals and the covariate, X. For censored data, the slope can be estimated using the extended version called Akritas–Theil–Sen estimator (Akritas, Murphy & LaValley 1995). The intercept of the linear regression is the median of the residual. We use the cenken routine in nada package to compute the Kendall’s τ correlation coefficient and fit the ATS regression line to our censored CO data (Halsel 2012; Lee & Lee 2015).

4 RESULTS

Fig. 5 compares the ECDFs of log f_gas  of normal star-forming galaxies (D_n(4000) < 1.6) with that of PSBs while Fig. 6 compares the ECDF of the star-forming galaxies with Seyfert PSBs. The vertical probability axis in both plots tracks the estimated Kaplan–Meier percentiles of galaxies with gas fraction less than a given threshold (Kaplan & Meier 1958; Halsel 2012). The log f_gas  curves for the PSBs and non-PSBs are similar at high gas fractions but are different at lower gas fractions. The log-rank test indicates that ECDFs of all PSBs or Seyfert-only PSBs are significantly different from that of the young non-PSBs. When the ECDF for Seyfert PSBs is compared to that of the young non-PSBs, the test gives a χ² of 24 on 1 degree of freedom corresponding to a p-value of 10⁻⁶. Therefore, the null hypothesis that the EDCFs for these two populations are the same can be rejected at a 4.8σ significance level. Similarly, the test gives a χ² of 11.4 on 1 degree of freedom, corresponding to a p-value of 7 × 10⁻⁴, when EDCFs of all PSBs and non-PSBs with D_n(4000) < 1.6 are compared.

Having established that there is a statistically significant difference between the gas fraction of young star-forming non-PSB galaxies and of PSBs in general, and of Seyfert PSBs in particular, in the rest of this section we explore in more detail where these trends come from. We first demonstrate that the WISE flux ratio between 12 and 4.6 μm is an excellent proxy for log f_gas  using the CO Legacy Database (COLD) GASS sample of non-PSBs (Saintonge et al. 2011) and then use it as a tool to constrain the gas fraction evolution in PSBs. Stellar populations younger than 0.6 Gyr dominate the 12 μm emission, and this ratio is known to correlate well with the specific SFR (sSFR; Donoso et al. 2012). The AGN hot dust emission is mostly between ∼ 3  and  5 μm; some have suggested that this ratio is not appreciably affected by AGN emission and that it can trace the sSFR in AGN (Donley et al. 2012; Donoso et al. 2012). Yesuf et al. (2014) have shown that starbursts and PSBs form a time sequence in this ratio, with starbursts having the highest ratios and quiescent PSBs having ratios similar to those of green-valley and red-sequence galaxies.

Fig. 7 shows the WISE flux ratio log f₁₂/f_4.6  versus log f_gas  for all normal galaxies in the COLD GASS survey. As shown in Fig. A3 of Appendix A, the general population of SDSS galaxies show bimodal log f₁₂/f_4.6  distribution that correlates with D_n(4000). Star-forming galaxies have ratios between log f₁₂/f_4.6  ∼ 0.2 and 1 while quiescent galaxies have ratios between log f₁₂/f_4.6  ∼ −0.8 and 0.0. A tight correlation is observed between log f₁₂/f_4.6  and log f_gas  in the COLD GASS sample, indicating an evolutionary process in normal galaxies in which sSFR declines as gas is used up. The slope of the best-fitting line is β = 1.13 and its intercept is β₀ = −1.8. The Kendall’s τ correlation coefficient for the fit is 0.64. The statistical significance of the correlation is more than 5σ. One can easily show that β and β₀ are related to the logarithm of the molecular gas depletion time at a given sSFR by the relation: |$\log \tau _{\rm dep}= \hat{\beta }\log {\rm sSFR}+ \hat{\beta _0}$|⁠, where |$\hat{\beta }=(\beta /\alpha -1)$|⁠, |$\hat{\beta _0}= -\beta (\alpha _0/\alpha +0.769)+\beta _0$|⁠, where α and α₀ define the relation: log sSFR = α(log f₁₂/f_4.6 +0.769) + α₀. Donoso et al. (2012) computed the relation between sSFR and 4.6 − 12 μm colour in Vega magnitude and 0.769 is subtracted to change the colour to AB magnitude. They found α = 2.5 × 0.775 = 1.938 and α₀ = −12.56 for the SDSS star-forming galaxies. Using these values with our fitted values of β and β₀ gives |$\hat{\beta } = -0.42$| and |$\hat{\beta _0} = 4.7$|⁠. In comparison, Huang & Kauffmann (2014) found |$\hat{\beta }=-0.37 \pm 0.04$| and |$\hat{\beta _0}=5.45 \pm 0.42$| for the COLD GASS sample. The authors used the GALEX FUV and WISE 22 μm to estimate the sSFR. Our inferred depletion times are consistent with Huang & Kauffmann (2014) but there is likely a systematic offset in the two SFR estimates.

Fig. 8 shows the relationship between the WISE flux ratio log f₁₂/f_4.6  and log f_gas  for PSBs. They also show statistically significant (5σ) correlation between these two quantities. The slope of the best-fitting line is β = 1.46 and its intercept is −2.1. Collectively, the broad overall trend suggests same correlation is at work: lower molecular gas fractions mean lower sSFR or dust heating. But the scatter is much larger, and looking in detail we can see that certain subpopulations of PSBs statistically differ. A major insight from this figure is that the sSFR histories of the PSBs and of Seyferts PSBs may not be as simple or well behaved as normal galaxies. Hence, the use of standard age indicators may lead to different results depending on indicator used. This topic will be explored in future papers. The reason for the low log f_gas  of Seyferts in Fig. 6 is consistent with low sSFR upper limits inferred from log f₁₂/f_4.6.

Table 2 summarizes the distribution of molecular gas fraction of our sample and of the previous samples of PSBs. Our sample has a log f_gas  distribution with mean of 0.025 and dispersion 0.018. 95 per cent of our PSBs have gas fractions below 0.07 while 50 per cent have gas fractions below 0.02. The previous samples have means and dispersions of ∼0.1–0.2. It is particularly interesting that the French et al. (2015) sample, which may also contain late-stage PSBs, has log f_gas  distribution with mean of 0.08 and dispersion 0.12. The majority of PSBs from Rowlands et al. (2015) and Alatalo et al. (2016) samples have both high log f₁₂/f_4.6  and high molecular gas fractions in comparison to our sample. The PSBs in these two works may be precursors to the PSBs in our sample or to those in the French et al. (2015) sample.

5 DISCUSSION

In this section, we discuss in detail how gas fractions in our sample compare with those observed in previous samples of PSBs. We will also discuss the general implication of our result and the uncertainty in the hitherto assumed CO conversion factor.

5.1 Comparison of our sample with the previous samples of PSBs

We study the molecular gas evolution of 116 PSB galaxies, of which 22 have new observations undertaken by us while the rest are compiled from the literature (French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016). The galaxies in the new sample are selected using their spectral indices and NUV − g colours, to identify them as green-valley PSBs, and their emission line ratios, to identify them as Seyferts. Our sample represents about two-thirds of the Seyferts PSBs in the combined sample. The Seyfert PSBs from the previous works were serendipitously selected and are mostly in the blue cloud. We show that log f_gas  distribution of the Seyfert PSBs and that of PSBs in general are significantly different from that of normal star-forming galaxies. The molecular gas fraction in the Seyfert PSBs decreases with the WISE ratios log f₁₂/f_4.6 . If there is an evolutionary sequence in Fig. 8, from upper right to lower left, it would be plausible to conclude that our new Seyferts are fairly aged galaxies located nearing the middle of the gas-exhaustion cycle. This conclusion would be consistent with our selection of these galaxies based on integrated UV colour as members of the green valley and well on their way to gas exhaustion. The systematically lower gas content of the new Seyferts relative to the PSBs in Fig. 8 might then be interpreted as the action of late AGNs in them in finally clearing out the gas.

However, we call attention to an important puzzle. In Fig. 8, the French et al. (2015) sample as a whole shows a significant correlation between log f₁₂/f_4.6  and log f_gas; the PSBs with high log f₁₂/f_4.6  have high gas fractions. This is puzzling because the French et al. (2015) galaxies were selected to have very low Hα, which is the gold standard for SFR. Thus, the log f₁₂/f_4.6  and Hα as SFR indicators disagree for the subset of the French et al. (2015) sample with high log f₁₂/f_4.6 . However, A-star heating is significant during the PSB phase, which can artificially boost IR-based SFR indicators (Hayward et al. 2014; Utomo et al. 2014). An analysis of multiple SFR indicators, including extinction-free 1.4 GHz measurements in section 2.5 of French et al. (2015), shows these galaxies are unlikely to have significant dust obscuration. Analyses of the near- and far-IR properties of this sample are the subject of a forthcoming study (Smercina et al., in preparation).

The episodic lifetime of an individual AGN event is currently not well constrained, and it may range between ∼0.3 and 1 times the effective AGN lifetime (Hopkins & Hernquist 2009). The effective AGN lifetime for the Seyfert PSBs is likely to be ∼2 Gyr. This estimate uses the median σ = 110 km s⁻¹ to estimate the median black hole mass, 10^7.4 M_⊙, from the M_BH–σ relation (Kormendy 2013), the median dust-corrected O iii luminosity, 10^7.3 L_⊙, to estimate the Eddington ratio, λ = 0.003 (Lamastra et al. 2009), and fig. 9 of Hopkins & Hernquist (2009) to estimate the effective AGN lifetime from M_BH and λ. The effective AGN lifetime for the Seyfert PSBs is longer than the quenched PSB lifetime of ∼1 Gyr, which suggests that, given the above estimate of the episodic lifetime, the quenched PSBs may experience second episode of AGN activity or no AGN activity at all within their typical lifetime. Perhaps the variation in AGN activity time-scale is linked to a wide range in gas fraction observed in the French et al. (2015) sample, if these PSBs are truly quiescent galaxies, and AGN are linked to the removal of gas in the French et al. (2015) galaxies with low log f₁₂/f_4.6  and in our sample.

We used the log f₁₂/f_4.6  ratio as a crude first-order indicator of star formation activity in both AGN and non-AGN galaxies. Some variation is expected between the ratios of the two classes (Donoso et al. 2012). It is not clear if these variations are due to additional AGN dust heating or preferentially aged stellar populations in AGN host galaxies. Understanding the details of this variation is not important for conclusions of the current work, but future studies of the effect of AGN contamination on log f₁₂/f_4.6  will be useful.

Furthermore, it should be noted that the gas fractions in very nearby early-type galaxies are found to be ≲0.1 per cent (e.g. Young et al. 2011). Many of the PSBs considered here have molecular gas fractions already consistent with detections in these early-type galaxies. A full comparison with early-type galaxies is complicated by the lower redshifts and smaller physical apertures of the Young et al. (2011, ATLAS 3D) sample. Deeper observations of the non-detected PSBs will be necessary to study the full range in their molecular gas properties, and to compare to the deep observations of the ATLAS 3D sample of early-type galaxies.

5.2 Implication for galaxy evolution

It is thought that star formation quenching happens in both slow and fast modes (e.g. Faber et al. 2007; Cheung et al. 2012; Fang et al. 2012; Barro et al. 2013; Dekel & Burkert 2014; Schawinski et al. 2014; Yesuf et al. 2014; Woo et al. 2015; Pandya et al. 2016). The existence of these modes, their relative importance and the physical mechanisms that drive them are still open questions in galaxy evolution. Some works have indicated that PSBs quench rapidly and may exemplify an important mode quenching for the build-up of the red sequence (Kaviraj et al. 2007; Wild et al. 2009, 2016; Wong et al. 2012; Yesuf et al. 2014), while others have questioned their importance (De Lucia et al. 2009; Dressler et al. 2013).

Recent studies found large gas fractions in PSBs complicating the picture of their rapid evolution to the red sequence (French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016). However, these works have been biased against finding PSBs with Seyfert-like emission lines in the UV–optical green valley. In this work, we have observed the molecular gas properties of 22 such galaxies. The combined data presented in this work cannot rule out AGN feedback, and as a whole are consistent with a removal or destruction of molecular gas at later stages by AGN feedback or no AGN feedback at all if different subsamples of PSBs are not evolutionary connected. Better and deeper observations, especially at a later stage of PSB evolution, are needed to firmly and more directly test the merger-induced AGN feedback hypothesis. In the theoretical front, developments in cosmological simulations with AGN feedback are needed to be compared to these future observations.

5.3 The effect of CO conversion factor

The CO luminosity to the molecular hydrogen conversion factor is a widely acknowledged source of uncertainty in studies that use CO as a tracer of molecular hydrogen (Bolatto et al. 2013). Up to this point, we have assumed that all objects including PSBs have a conversion factor similar to the factor observed in the Milky Way disc. There is no study to date that measured the conversion factor in PSBs. As discussed in detail in Bolatto et al. (2013), departure from the Galactic conversion factor is both observed and expected in starburst galaxies. It has been shown that adopting the Galactic conversion factor would cause the inferred molecular gas mass to exceed the dynamical mass for the central region of a starburst galaxy (Solomon et al. 1997; Downes & Solomon 1998). Using galaxy merger hydrodynamic simulations, which incorporate dust and molecular line radiative transfer calculations, Narayanan et al. (2011) have shown that in merger-induced starbursts, the combined effect of increased velocity dispersion and kinetic temperature increases the velocity-integrated CO intensity, and lowers the CO conversion factor from the Galactic value by a factor of ∼2–10. The authors also note that in the PSB phase, it is less trivial to simply relate the conversion factor to the gas velocity dispersion and temperature owing to varying physical conditions in the PSB galaxy. Some of their merger simulations returned to a Galactic conversion factor value quickly after the peak of the starburst, while others remained low. Furthermore, if AGN heat the bulk of molecular gas, the conversion factor could be lower from the Galactic value in PSBs with AGN. Note that a lower CO conversion factor than the assumed Galactic value will lower the gas fractions observed in PSBs compared to normal galaxies.

6 SUMMARY AND CONCLUSIONS

Using the SMT, we undertook new CO (2–1) observations of 22 green-valley Seyfert PSB candidate galaxies at redshift z = 0.02–0.06. The sample was selected using the dust-corrected NUV − g colour, Hδ absorption and 4000 Å break to indicate a PSB signature, and the emission line ratios to indicate Seyfert activity. We analysed our sample with previous samples of 94 PSBs. The combined sample probes a variety of stages in a possible evolutionary sequence, and spans a range of AGN properties. Our main results are as follows.

• We detect molecular gas in only 6 out of 22 Seyfert PSB galaxies we targeted. Using comparable sensitivity limits, this is consistent with the detection in PSB sample by French et al. (2015) using the SMT. However, taking into account our upper limits, the mean and the dispersion of the distribution of the gas fraction in our green-valley Seyfert PSBs (μ = 0.025, σ = 0.018) are much smaller than previous samples of Seyfert PSBs or PSBs in general (μ ∼ 0.1–0.2, σ ∼ 0.1–0.2; French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016). This difference may partly be explained by the evolution of molecular gas fraction with PSB age if our sample probes the late stage of the evolution. However, PSBs in French et al. (2015) also have similar green-valley colours but they show a wider range of gas fraction, despite their low SFR inferred from their Hα emission. The comparison between these samples is complicated by the differing selections on Hδ absorption, as weaker Hδ is consistent with either a weaker starburst or a later PSB age.

• The distribution of gas fraction in PSBs is significantly different from young star-forming galaxies from COLD GASS survey (Saintonge et al. 2011). PSBs are more likely to have lower gas fraction. The Seyfert PSBs have a distribution of gas fraction that is even more significantly different from that of normal star-forming galaxies of similar D_n(4000).

• The WISE flux ratio, log f₁₂/f_4.6 , is an excellent proxy for gas fraction for both PSBs and non-PSBs. We find a statistically significant (5σ) relationship between log f₁₂/f_4.6  and log f_gas  for both PSBs and normal star-forming galaxies.

Acknowledgments

We thank Xavier Prochaska and Ann Zabludoff for their useful comments and suggestions. HY would like to acknowledge support by NSF grants AST-0808133 and AST-1615730. KDF acknowledges support from NSF grant DGE-1143953, PEO, and the ARCS Phoenix Chapter and Burton Family.

The SMT is operated by the Arizona Radio Observatory (ARO), Steward Observatory and University of Arizona. We thank the operators and staff of the ARO.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society and the Higher Education Funding Council for England. The SDSS website is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory and the University of Washington.

REFERENCES

APPENDIX

This appendix presents materials that support the analysis in the main sections of the paper.

A1 Evolution of Hδ equivalent width after a starburst

Fig. A1 shows the evolution of Hδ equivalent width for Bruzual & Charlot (2003) starburst models, with a star formation time-scale, τ = 0.1 Gyr and a burst mass fraction (b_f) 3 or 20 per cent. As detailed in Yesuf et al. (2014), the Hδ is dependent on the magnitude of the burst and there is generally a degeneracy between the burst mass fraction and the burst time-scale. The figure demonstrates, with particular models above, that late-stage PSBs can have Hδ < 3 Å, which is lower than the threshold often used to define PSBs.

A2 A comparison of our PSB sample with previous samples

This section presents additional figures that compare the stellar population of our Seyfert PSBs with previous samples of PSBs (French et al. 2015; Rowlands et al. 2015; Alatalo et al. 2016). Fig. A2 shows the NUV − g versus g − z colours before and after the dust correction. Before the dust correction, star-forming galaxies form a diagonal track that stretches from blue to red colours. The red end of this track is inhabited by dusty galaxies. The quiescent galaxies form a separate concentration above the dusty star-forming galaxies. The locations our Seyfert PSBs in the two-colour diagrams indicate that they are not dusty star-forming galaxies but are transitional galaxies in green valley. The dust-corrected diagram in addition shows that stellar ages of our PSBs and the PSBs in French et al. (2015) are comparable and the PSBs in Rowlands et al. (2015) and Alatalo et al. (2016) are younger than our sample. Fig. A3 depicts log f₁₂/f_4.6  versus log f_4.6/f_3.4  (Wright et al. 2010). Stellar populations younger than 0.6 Gyr dominate the 12 μm emission and log f₁₂/f_4.6  is known to correlate well with sSFR (Donoso et al. 2012) and log f_4.6/f_3.4  is sensitive to hot dust emission from AGN. Normal galaxies form a tight bimodal sequence with some vertical scatter. Star-forming galaxies are located at high log f₁₂/f_4.6 , while quiescent galaxies are located at low log f₁₂/f_4.6 .

Fig. A4 depicts the WISE flux ratio log f_4.6/f_3.4  versus log f_gas  for PSBs. The log f_4.6/f_3.4  >−0.06 criterion can identify hot dust emission from AGN albeit only with ∼50 per cent reliability (Assef et al. 2013). Anti-correlation between log f_4.6/f_3.4  and log f_gas  may naively be expected if AGN feedback destroys the molecular gas by heating it. Such a trend is not seen in the data. There is instead somewhat significant (∼2σ) correlation between log f_4.6/f_3.4  and log f_gas  for the BPT-identified Seyfert PSBs. The current data cannot rule out AGN feedback, and the observed trend is consistent with delayed AGN feedback, which is inefficient at earlier times when the molecular gas fraction is above ∼10 per cent. Future study with large number of PSBs with log f_4.6/f_3.4  >−0.06 may constrain this model better. It should be noted that correlation between log f_4.6/f_3.4  and log f_gas  for all PSBs in the combined sample is highly significant (>5σ) and it is likely tracing stellar dust heating process that declines with the starburst age.

A3 CO luminosities of Seyfert PSBs

Table A1 presents the observed CO luminosities of our new Seyfert PSB sample or 3σ upper limits in the case of no detections. Note that TPSB2 and TPSB23 are the same object, which is observed twice and we have combined the two observations in the analysis presented in the main text of the paper.