Evolution of gas flows along the starburst to post-starburst to quiescent galaxy sequence Open Access

ABSTRACT

We measure velocity offsets in the |$\rm{Na {}\rm {\small I}}$| λλ5890, 5896 (⁠|$\rm{Na {}\rm {\small D}}$|⁠) interstellar medium absorption lines to track how neutral galactic winds change as their host galaxies evolve. Our sample of ∼80 000 SDSS spectra at 0.010 < z < 0.325 includes starburst, post-starburst, and quiescent galaxies, forming an evolutionary sequence of declining star formation rate (SFR). We detect bulk flows across this sequence, mostly at higher host stellar masses (log(M_*/M_⊙)) > 10). Along this sequence, the fraction of outflows decreases (76 ± 2 per cent to 65 ± 4 per cent to a 3σ upper limit of 34 per cent), and the mean velocity offset changes from outflowing to inflowing (−84.6 ± 5.9 to −71.6 ± 11.4 to |$76.6\pm 2.3\, \rm km s^{-1}$|⁠). Even within the post-starburst sample, wind speed decreases with time elapsed since the starburst ended. These results reveal that outflows diminish as galaxies age. For post-starbursts, there is evidence for an AGN contribution, especially to the speediest outflows: (1) SFR declines faster in time than outflow velocity, a decoupling arguing against massive stellar feedback; (2) of the few outflows strong enough to escape the interstellar medium (9/105), three of the four hosts with measured emission lines are Seyfert galaxies. For discy starburst galaxies, however, the trends suggest flows out of the stellar disc plane (with outflow 1/2-opening angle > 45°) instead of from the nucleus: the wind velocity decreases as the disc becomes more edge-on, and the outflow fraction, constant at |$\sim 90~{{\ \rm per\ cent}}$| for disc inclinations i < 45°, steadily decreases from |$\sim 90~{{\ \rm per\ cent}}$| to 20 per cent for i > 45°.

1 INTRODUCTION

In the local universe, observational surveys suggest there are two main types of galaxies: star-forming and quiescent. The current understanding of this bimodality is that star-forming galaxies will eventually stop forming stars and evolve to the quiescent phase. Transitioning galaxies evolving between star-forming and quiescent are found in the ‘green valley’ of the galaxy colour–magnitude diagram (Baldry et al. 2004). Among transitioning galaxies are post-starburst galaxies, which experienced a burst of star formation that ended in the past ∼1 Gyr (Dressler & Gunn 1983; Couch & Sharples 1987). Evidence for this star formation history (SFH) comes from their optical spectra, which show strong Balmer absorption lines that indicate a substantial population of A-type stars and thus a recent starburst but lack the nebular emission line signatures of current star formation. Even though star formation must decline along the evolutionary sequence from starburst to post-starburst to quiescent, the reasons why are still unknown.

One likely mechanism for this quenching is feedback. For current cosmological simulations of galaxy formation and evolution, feedback is required to reproduce observed galaxy baryonic mass functions and star formation efficiency (Kereš et al. 2009; Hopkins et al. 2014; Somerville & Davé 2015; Naab & Ostriker 2017, and references therein). Massive stars, supernovae, and active galactic nuclei (AGN) may be sources of negative feedback through thermal or kinetic energy transferred to the interstellar medium (ISM) in galaxies. If negative feedback processes remove the cold gas in winds (see Veilleux, Cecil & Bland-Hawthorn 2005; King & Pounds 2015; Veilleux et al. 2020 for a comprehensive galactic wind review), star formation could rapidly end due to the removal of the fuel for star formation.

Previous work has found observational evidence of outflows in low-redshift galaxies, generally focusing on AGN-host galaxies (e.g. Crenshaw et al. 1999; Nesvadba et al. 2008; Rubin et al. 2011; Cicone et al. 2014; Harrison et al. 2014; Kang & Woo 2018; Baron & Netzer 2019) or on pure star-forming or starburst galaxies (e.g. Heckman, Armus & Miley 1990; Heckman et al. 2000; Rupke, Veilleux & Sanders 2005a, b; Chen et al. 2010; Rubin et al. 2014; Zschaechner et al. 2018; Perrotta et al. 2021). For post-starburst galaxies, outflows have been detected in hosts with (Baron et al. 2021) and without AGN (Coil et al. 2011). The evolution of such outflows along the starburst to post-starburst to quiescent sequence, however, remains unprobed.

Outflows are detected in different gas phases, for example molecular gas (Cicone et al. 2014; Fluetsch et al. 2018), neutral gas (Chen et al. 2010; Bae & Woo 2018; Concas et al. 2019), or ionized gas (Harrison et al. 2014; Concas et al. 2017; Bae & Woo 2018). The complex connections between multiphase outflows and star formation or AGN activity have been well-studied over the past decades, but a unified picture is elusive. For instance, Concas et al. (2017) analysed a complete spectroscopic sample of galaxies in the SDSS DR7 survey and found that the flux percentage of the outflow component of the ionized gas outflows tracer [|$\rm{O {}{\small III}}$|] λ5007 are higher in low-ionization nuclear emission-line regions (LINERs) and Seyfert galaxies than in pure star-forming galaxies, suggesting that such winds are powered by the nucleus. Alternatively, Chen et al. (2010) systematically analysed the stacked |$\rm{Na {}\rm {\small D}}$| lines of star-forming galaxies in the SDSS DR7 survey, finding that significant ISM |$\rm{Na {}\rm {\small D}}$| absorption lines (equivalent widths |$\text{EW}\gt 0.8~\rm{\mathring{\rm A}}$|⁠) are only prevalent in disc galaxies with high star formation rate (SFR), SFR per unit area (Σ_SFR), or stellar mass (M_⋆).

In their multiphase outflow study of post-starburst galaxies, Baron et al. (2021) used both model-fitting and Machine-learning methods to select post-starburst galaxies with strong recent starbursts (⁠|$H\delta ~\text{EW}\gt 5~\rm{\mathring{\rm A}}$|⁠), AGN features (AGN-like narrow-line flux ratios), and ionized outflows (multiple broad emission lines). They found 144 eligible post-starburst galaxy candidates, 40 per cent of which were detected with |$\rm{Na {}\rm {\small D}}$|-traced neutral gas outflows. While they found a connection between the detection of |$\rm{Na {}\rm {\small D}}$| and star-formation (SF) luminosity, the correlations of the neutral outflow velocity, mass outflow rate, and kinetic power with SF luminosity or AGN luminosity were not significant. The ionized outflow velocity is correlated significantly with both AGN luminosity and SF luminosity but has a stronger correlation with AGN luminosity.

The detection of neutral gas outflows is correlated with the disc inclination in starbursting galaxies. Heckman et al. (2000) found ∼70 per cent of such outflows are detected in discs with inclinations (i) of less than 60°. Later works similarly found that most of the neutral gas outflows detected in absorption are in galaxies with low inclination (Chen et al. 2010; Rubin et al. 2014; Roberts-Borsani & Saintonge 2018). Bae & Woo (2018) confirmed the inclination dependence of neutral outflow velocity and velocity dispersion in both pure star-forming and AGN-host galaxies, interpreting this result as support for a perpendicular outflow geometry that could only arise from star formation-driven winds. Concas et al. (2019) observed that the velocities of neutral gas outflows decrease with inclination for their starburst galaxy sample.

Even though many studies have identified connections between SFR and wind properties, none have tested how winds evolve along the galaxy evolutionary sequence where star formation is declining. Based on previous findings that higher outflow detection rates and larger |$\rm{Na {}\rm {\small D}}$| EWs prefer higher SFR and M_⋆ host galaxies, we expect that there should be a decline in winds over this sequence. Therefore, in this paper, we build a sequence of declining star formation, from starburst to post-starburst to quiescent galaxies in the local universe (0.010 < z < 0.325). We then explore the evolution of neutral gas winds along that sequence. Within our post-starburst sample, we know even more: the time elapsed since the starburst ended (French et al. 2018) (hereafter F18). These ‘post-burst ages’ allow us to trace the wind evolution directly over time, from roughly 0 to 2 Gyr after the end of the burst, with an age resolution of about ± 50 Myr.

We use the spectroscopic SDSS DR12 survey (Alam et al. 2015) and analyse the |$\rm{Na {}{\small I}}$| λλ5890, 5896 line to trace the bulk motions of the neutral gas in the starburst, post-starburst, and quiescent samples. We test how the detection of winds and the wind velocities change along the whole sequence and with post-burst age. We further investigate whether wind properties are related to host galaxy properties such as stellar mass and disc inclination. Finally, we consider the potential sources of the winds.

This paper is organized as follows: we describe the selection of our three galaxy samples in Section 2. Our method of subtracting the stellar component from each SDSS spectrum, fitting the profile of the residual (ISM) |$\rm{Na {}\rm {\small D}}$| absorption line, and deriving the direction (inflow versus outflow) and line-of-sight velocity of the flow is described in Section 3. We present the main results in Section 4 and discuss them further in Section 5. Finally, we summarize our findings in Section 6.

Throughout this paper, we assume a standard ΛCDM universe whose cosmological parameters are H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_m = 0.3.

2 SAMPLE SELECTION AND DATA

This study is based on three different galaxy samples: a post-starburst galaxy sample and two comparison samples of starburst and quiescent galaxies. The starburst sample is chosen to be consistent with the evolutionary precursors of the post-starburst sample. The quiescent sample is selected to be consistent with evolving from galaxies like those in the post-starburst sample.

2.1 Post-starburst galaxy sample

F18 drew 532 post-starburst galaxies from the spectroscopic galaxy catalogue of the Sloan Digital Sky Survey Data Release 8 (DR8, Aihara et al. 2011), using galaxy properties from the MPA-JHU Catalogues (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). From the plane of the Lick Hδ absorption line index (Hδ_A, LICKHDA) versus the H α emission line equivalent width (HALPHAEQW) (the left panel of Fig. 1), F18 identified the post-starburst sample from the lower right corner. They selected for low H|$\alpha ~\text{EW}\lt 3~\rm{\mathring{\rm A}}$| in emission, representing a lack of current star formation, and strong Balmer absorption line index |${\rm H}\delta _{\rm A}-\sigma ({\rm H}\delta _{\rm A}) \gt 4~\rm{\mathring{\rm A}}$|⁠, where σ(Hδ_A) is the measurement error, representing a large population of A-type stars and thus a substantial burst of star formation in the past ∼1-1.5 Gyr.

Through a detailed modelling of the SFHs of the post-starburst galaxies, F18 parameterized each SFH with the post-burst age, the fraction of stellar mass produced in the burst (m_burst), and the burst duration. In their age-dating technique, they used Lick indices measured from the SDSS spectra. They excluded galaxies with signal-to-noise (S/N) less than 10 pixel⁻¹ in the continuum to avoid unreliable measurements. The best constraints on the age-dating parameters require ultraviolet (UV) photometry. By matching with Galaxy Evolution Explorer (GALEX, Martin et al. 2005) far-UV (FUV) and near-UV (NUV) data, they finalized a sample of 532 post-starburst galaxies with >3σ detections in both the NUV and FUV bands.

After F18 selected their post-starburst sample, the MPA-JHU Catalogues were updated, and some line measurements were revised. As a result, nine of the original F18 galaxies no longer satisfy the post-starburst Hδ_A-H α selection criteria, and another galaxy has an unreliable global stellar mass measurement (LGMTOTP50). We also find the H α flux measurements of six galaxies are bad (HALPHACHISQ=0). Throughout this paper, after removing those 16 galaxies, we use the remaining 516 post-starburst galaxies from F18. Their redshift range is 0.010 < z < 0.325.

Measuring current SFRs for post-starburst galaxies is complicated by the uncertainties due to dust geometry, the long time-scale of star formation traced by some indicators, and contamination from possible AGN (see recent review in French 2021). H α is a relatively good SFR indicator in the post-starburst phase, because it traces a short star formation time-scale (i.e. that associated with short-lived massive stars) and dust attenuation and AGN effects can be modelled and removed.

Here we convert the H α fluxes from the MPA-JHU catalogues to SFR following French et al. (2015) and Li et al. (2019).¹ We use the L_H α–SFR relation from Kennicutt, Tamblyn & Congdon (1994). We estimate the amount of internal dust extinction from the observed Balmer decrement, H α/H β. Assuming that the hydrogen nebular emission follows Case B recombination, the intrinsic Balmer flux ratio (H α/H β)₀ = 2.86 for T_e = 10⁴ K. We adopt the reddening curve of O’Donnell (1994). When the H β line flux is negative in the MPA-JHU catalogues, we assume the mean value of E(B − V) from the other post-starburst galaxies. The mean attenuation is A_V = 0.92 mag.

Following the methodology from Wild, Heckman & Charlot (2010), we correct for a potential underlying AGN contribution to the H α flux and thus to the derived SFR. We calculate the emission-line ratios [|$\rm{O {}{\small III}}$|]/H β and [|$\rm{N {}{\small II}}$|]/H α of our post-starburst galaxies, place them on the BPT diagram, and determine the likely AGN contributions to their H α luminosity. For 230 post-starbursts with [|$\rm{O {}{\small III}}$|], [|$\rm{N {}{\small II}}$|], and H α S/N < 3, we do not have a reliable AGN correction, so we do not measure their SFRs. We do not set a S/N constraint on the H β line, as it is not usually well detected in our post-starburst galaxies. For the case of negative H β line flux, we use its 1σ flux upper limit to determine the corresponding 1σ upper limit for the AGN correction to the SFR. When a negative correction factor is derived, if the corresponding 1σ upper limit is positive, we use this upper limit to determine the SFR 1σ upper limit; otherwise, we do not measure the SFR.

We then correct for the SDSS fiber aperture using the ratio of the global SFR to the SFR within the SDSS fiber (SFRTOTP50/SFRFIBP50) from the MPA-JHU catalogues. In the end, we obtain SFR measurements for 171 of the 516 post-starburst galaxies, including 79 post-starburst galaxies for which only SFR upper limit measurements are available. The bottom panel of Fig. 2 shows that the SFRs of our post-starburst galaxies decline over time. We discuss the time-scale of this decline in subsection 5.1.

2.2 Comparison samples

We build a broader evolutionary sequence of declining star formation by constructing samples of likely post-starburst progenitors (starburst galaxies) and descendants (quiescent galaxies). We start with a parent sample of 593 807 galaxies drawn from the SDSS DR12 survey (Alam et al. 2015). These galaxies lie within the same redshift range as our post-starburst galaxies and have properties and spectral-line flux measurements from the MPA-JHU catalogues. All galaxies are selected to have continuum S/N ≥10.

2.2.1 Starburst galaxy sample

As the progenitors of post-starburst galaxies, starburst galaxies should have strong ongoing star formation by definition. Because the post-starburst sample from F18 was selected against strong AGN, we also require their progenitors to be without strong AGN features. While there are SDSS starburst galaxies with powerful AGN that could also evolve into galaxies like our post-starburst galaxies, we exclude them here to avoid having to decouple the time-scale of AGN decline from that of the SFR.

We first select an appropriate parent sample of 157 018 star-forming galaxies without strong AGN, using Kewley et al. (2006)’s classification scheme

with S/N ≥3 in the H α, H β, [|$\rm{O {}{\small III}}$|], and [|$\rm{N {}{\small II}}$|] emission line fluxes. Any remaining galaxies that are classified as quasars by the SDSS pipeline (SpecClass = QSO) or have invalid flux values around the |$\rm{Na {}\rm {\small D}}$| region are excluded as well.

To select which of these star-forming galaxies have SFRs characteristic of starburst galaxies, we consider the relation between H α-based specific SFRs (sSFRs) and M_* from the MPA-JHU catalogues (Fig. 1, right panel). The solid line indicates the star-forming main sequence from Elbaz et al. (2007).² The 2–4×SFR levels are often used to distinguish between the main sequence and starburst galaxies (Elbaz et al. 2011; Rodighiero et al. 2011; Sargent et al. 2014). There are 3723 galaxies located above the 5 × SFR level (dashed line); we define these as starburst galaxies and further justify that choice below.

To ensure that the starburst sample is the progenitor of the post-starburst sample, we consider the growth of the stellar mass based on the mass fraction of new stars created in the starburst m_burst. For each post-starburst galaxy, we determine the mass range of its potential progenitors as [M_*/M_⊙(1 − m_burst), M_*/M_⊙] and choose all starburst galaxies within it. F18 found some post-starburst galaxies have poorly constrained burst mass fractions (m_burst > 0.5 and with large uncertainties). Therefore, we set these to 0.5, the 1σ upper limit on the average burst mass fraction for the better-determined values. In the end, we select 3306 galaxies as possible starburst progenitors of the post-starburst sample. The right panel of Fig. 1 confirms that these starburst galaxies have sSFRs comparable to the maximum sSFRs of our post-starburst galaxies during their starburst phase (F18).

2.2.2 Quiescent galaxy sample

Quiescent galaxies are selected to have |${\rm H}\delta _{\rm A}+\delta ({\rm H}\delta _{\rm A})\lt 1~\rm{\mathring{\rm A}}$| and H α EW|$\lt 3~\rm{\mathring{\rm A}}$|⁠, as shown in the left panel of Fig. 1, with the intention of selecting galaxies with no current and recent star formation. These quiescent galaxies follow the example evolutionary tracks of post-starburst galaxies with different burst durations and burst mass fractions provided by F18 (adapted as blue lines in the left panel of Fig. 1). These selection criteria result in 260 526 galaxies from the parent sample. Again, we remove galaxies with significant emission lines to exclude ongoing SF or AGN activity. In this step, we remove 176 125 galaxies with S/N≥3 in at least one of the emission lines H β, [|$\rm{O {}{\small III}}$|] λ5007, H α, [|$\rm{N {}{\small II}}$|] λ6584, and [|$\rm{S {}{\small II}}$|] λλ6717, 6731. Also, 1208 quiescent galaxies with incomplete spectra around the |$\rm{Na {}\rm {\small D}}$| region are excluded. Our quiescent sample is finalized at 72 251 galaxies.

Under the simple assumption that post-starburst galaxies evolve with their current residual SFRs (typically ∼0.1 M_⊙ yr⁻¹) throughout the post-starburst phase, 10⁸ M_⊙ mass growth is expected over 1 Gyr. More generously, a maximum mass growth of 10⁹ M_⊙ can be expected over a Hubble time. Considering these mass growth bounds, we exclude 10 942 quiescent galaxies with stellar masses less than the minimum stellar mass (M_* < 10^8.29 M_⊙) of our post-starburst sample or greater than the maximum allowed stellar mass. Our quiescent sample is finalized at 72 251 galaxies.

For consistency with the starburst and post-starburst samples, we use the H α luminosity to estimate SFRs for the quiescent galaxies, rather than the D4000-based SFRs provided by the MPA-JHU catalogues. Our selection for low levels of star formation in the quiescents requires low H α values. Furthermore, we only measure SFR for those quiescents where H α has S/N > 1. We additionally require that all BPT emission lines for the quiescents have S/N < 3, which prevents us from correcting for any AGN contribution. Due to these restrictions, the SFRs for the quiescents are upper limits. In total, we obtain SFR measurements for 36 718 of 72 251 quiescent galaxies. The mean SFR of the quiescent sample (top panel of Fig. 2) is thus a rough upper limit.

To summarize, by following the earlier steps, we select starburst, post-starburst, and quiescent galaxy samples to build an evolutionary sequence of declining star formation. The top panel of Fig. 2 shows the mean global SFR of starbursts (log₁₀(SFR) = 0.82 ± 0.01 M_⊙ yr⁻¹), young (−1.05 ± 0.08 M_⊙ yr⁻¹), and old (−1.49 ± 0.10 M_⊙ yr⁻¹) post-starbursts (young and old post-starbursts will be defined in subsection 4.4.1), and quiescents (<−1.49 yr⁻¹). It confirms that the global SFR is declining along the sequence that we have constructed. Similarly, the bottom panel of Fig. 2 shows that the global SFR declines with post-burst age for the post-starburst sample.

3 ANALYSIS

In this section, we describe how we fit the |$\rm{Na {}\rm {\small D}}$| doublet absorption line profiles to derive flow properties. Recently, several works studied |$\rm{Na {}\rm {\small D}}$| absorption properties using the ‘stacked’ SDSS spectra of large galaxy samples (Chen et al. 2010; Cicone, Maiolino & Marconi 2016; Concas et al. 2017; Roberts-Borsani & Saintonge 2018); these analyses produce high S/N spectra, but preclude exploration of the |$\rm{Na {}\rm {\small D}}$| properties of individual galaxies. Therefore, in this paper, we focus on the |$\rm{Na {}\rm {\small D}}$| properties of each galaxy in the three galaxy samples from Section 2 to learn how those properties evolve in time. Our analysis consists of three steps (1) modelling, fitting, and subtracting the stellar continuum, (2) measuring the EW of the residual |$\rm{Na {}\rm {\small D}}$| absorption, that is, due to the neutral ISM, as well as the EW of the total |$\rm{Na {}\rm {\small D}}$| absorption before continuum-subtraction, that is, the sum of the stellar and neutral ISM contributions, and (3) determining the offset of the residual |$\rm{Na {}\rm {\small D}}$| from the redshift of the galaxy.

3.1 Stellar continuum modelling

We model the stellar continuum of each galaxy spectrum over the wavelength range of 3800–7000|$~\rm{\mathring{\rm A}}$| in the rest frame, using the penalized pixel fitting (pPXF, Cappellari & Emsellem 2004; Cappellari 2017) Python code. We mask some significant emission and absorption lines, for example Balmer lines, |$\rm{N {}{\small III}}$| λ4640, |$\rm{He {}{\small II}}$| λ4686, [|$\rm{O {}{\small III}}$|] λλ4959, 5007, |$\rm{He {}{\small I}}$| λ5875, [|$\rm{O {}{\small I}}$|] λ6300, [|$\rm{N {}{\small II}}$|] λλ6548, 6584, and [|$\rm{S {}{\small II}}$|] λλ6717, 6731, as they may affect the best fits of the stellar continuum models. The |$\rm{Na {}\rm {\small D}}$| lines are also masked because they include both stellar absorption and ISM contributions. A polynomial with degree=4 is added to avoid mismatches between galaxy spectra and stellar templates and ensure the accuracy of fitted kinematic parameters.

We use MILES (Vazdekis et al. 2010) simple stellar population (SSP) models as stellar templates, adopting the SSP model spectra with four different metallicities: [M/H] = −0.71, −0.40 (sub-solar), 0.00 (solar), and +0.22 (super-solar). For each galaxy, we convert its stellar mass into a range of metallicities using the mass–metallicity relation from Gallazzi et al. (2005) (see Fig. 8 in this paper):

• log(M_*/M_⊙) < 10.0: [M/H] = –0.71, –0.40, 0.00;

• 10.0 ≤ log(M_*/M_⊙) < 10.5: [M/H] = –0.71, –0.40, 0.00, +0.22;

• 10.5 ≤ log(M_*/M_⊙) < 11.0: [M/H] = –0.40, 0.00, +0.22;

• log(M_*/M_⊙) ≥ 11.0: [M/H] = 0.00, +0.22.

To determine the impact of the metallicity choice on our results, we also test the assumption of constant solar metallicity for all galaxies, finding that the following results and conclusions do not change significantly. Therefore, we use the results based on the four different metallicities throughout this paper, because this choice is more physically motivated.

pPXF provides the best fit redshift and its error from the stellar velocity measurement. We use redshifts from the MPA-JHU catalogues as an initial guess and run iterations to obtain the best fit results from pPXF, although the two redshifts are similar. If the difference between the redshifts of two adjacent iterations is less than |$3~\rm km s^{-1}$|⁠, we consider the results as converged and stop iterating; otherwise, the iteration continues until the total number of iterations reaches 10. Finally, we shift the original and the residual (after subtracting the best fit stellar continuum model) spectra to the rest-frame wavelength given by the pPXF redshift.

The left columns of Figs 3–5 show examples of galaxy spectra and the best fit stellar continuum models (red lines) for post-starburst, starburst, and quiescent galaxies, respectively. The region immediately around |$\rm{Na {}\rm {\small D}}$| and the best fit model are shown in the middle columns. The masked wavelength regions are marked as grey areas. Given that pPXF does not provide flux errors for the best fit stellar models, we assume that the residual flux error is the same as that of the original flux. Systematic errors arising from the assumed stellar population models may thus affect our estimations of the ISM properties.

3.2 EW of Excess Na D due to Neutral ISM

The total |$\rm{Na {}\rm {\small D}}$| absorption in the original galaxy spectrum arises from a combination of stellar and neutral ISM components. After stellar continuum subtraction, the residual |$\rm{Na {}\rm {\small D}}$| line indicates the amount of absorption due to the neutral ISM. To quantify the amount of absorption, we measure the total and residual EWs, EW_{NaD, tot} and EW_{NaD, exc}, respectively. In this work, a positive EW indicates absorption. To avoid contamination from the |$\rm{He {}{\small I}}$| λ5875 emission line frequently seen in star-forming galaxy spectra, we mask 5|$~\rm{\mathring{\rm A}}$| around 5875|$~\rm{\mathring{\rm A}}$|⁠. The EW measurement is performed over the wavelength range 5880–5910|$~\rm{\mathring{\rm A}}$|⁠, and the continuum level is defined as the average of the fluxes at 5850–5870|$~\rm{\mathring{\rm A}}$| and 5920–5940|$~\rm{\mathring{\rm A}}$| in the best fit stellar continuum model. We perform a Monte Carlo simulation, generating 1000 mock spectra based on the observed flux errors, to derive the measurement errors of EW_{NaD, tot} and EW_{NaD, exc}.

3.3 Velocity offset of Na D Excess

The next step is to fit the profile of the |$\rm{Na {}\rm {\small D}}$| excess and determine its velocity offset (ΔV) from the galaxy redshift determined by pPXF. Examples of the |$\rm{Na {}\rm {\small D}}$| excess fitting for three galaxy sample groups are shown in the right panels of Fig. 3–5. We apply a double Gaussian profile to fit the |$\rm{Na {}\rm {\small D}}$| residual spectra. Several previous studies (e.g. Rupke et al. 2005a; Chen et al. 2010; Lehnert et al. 2011; Sarzi et al. 2016) used a profile derived by a physical model to fit more physical parameters (the optical depth and the covering factor) and multiple components to distinguish different kinematic parts. Our fitting method is more straightforward to obtain the mean velocities of bulk flows, where a positive ΔV is an inflow and negative ΔV is an outflow.

The double Gaussian profiles that we use for fitting the |$\rm{Na {}\rm {\small D}}$| excess are

Here we assume that the widths of the two |$\rm{Na {}\rm {\small D}}$| lines are the same (σ) and that their height ratio (p₂/p₁) ranges from 1 to 2, corresponding to optically thick and thin scenarios, respectively. We also restrict the line widths to be no smaller than the instrumental resolution of the SDSS spectra (70 |$\rm km s^{-1}$|⁠). In some cases, the |$\rm{Na {}\rm {\small D}}$| lines are resolved (e.g. the first row of Fig. 3), while some are blended (e.g. the fifth row of Fig. 3). The wavelength shift of the |$\rm{Na {}\rm {\small D}}$| excess μ is varied over the range [−20, 20]|$+ 5890~\rm{\mathring{\rm A}}$|⁠, which corresponds to ΔV ∼ ± 1000 |$\rm km s^{-1}$|⁠. We fit the double Gaussian profiles to the |$\rm{Na {}\rm {\small D}}$| excess at 5850–5940|$~\rm{\mathring{\rm A}}$| using the PYTHON PYCMPFIT package.³

Depending on the geometry, it is possible to see |$\rm{Na {}\rm {\small D}}$| emission when absorbed photons are re-emitted along the line of sight. This scattering process can produce P-Cygni profiles of |$\rm{Na {}\rm {\small D}}$| (see the profiles of |$\rm{Mg {}{\small II}}$| λλ2796, 2803 simulated by the cold gas wind models in Prochaska, Kasen & Rubin 2011). We find that some galaxies show a clear P-Cygni profile in their residual spectra (the two bottom rows of Figs 3–5). Given that blueshifted absorption is accompanied by redshifted emission, the overall profile should be fit by two components.

To identify whether the |$\rm{Na {}\rm {\small D}}$| has a pure absorption (or emission) profile or a P-Cygni profile, we use two individual double Gaussian profiles to fit the absorption and emission components separately. Then, we classify the profiles as follows:

• Pure Absorption: the peak of only the absorption component has S/N_ab ≥3;

• Pure Emission: the peak of only the emission component has S/N_em ≥3;

• P-Cygni: both S/N_ab and S/N_em ≥3;

• Nothing: both S/N_ab and S/N_em <3.

The numbers of galaxies identified as ‘Pure Absorption,’ ‘Pure Emission,’ ‘P-Cygni,’, or ‘Nothing’ of the three galaxy samples are shown in Table 1. From here on, we exclude galaxies labelled as ‘Nothing’ or ‘Pure Emission’ to focus on deriving the kinematics properties of the neutral gas from the |$\rm{Na {}\rm {\small D}}$| absorption feature.

To identify the subset of galaxies with reliable fits to the excess |$\rm{Na {}\rm {\small D}}$| profile and the bulk velocity, we first require that, for ‘Pure Absorption’ cases, the |$\rm{Na {}\rm {\small D}}$| excess be significant (EW_{NaD, exc} − 3σ(EW_{NaD, exc}) > 0). We also include all galaxies with ‘P-Cygni’ profiles. We then employ Monte Carlo simulations again, sampling the error distribution of the residual |$\rm{Na {}\rm {\small D}}$| profile for each galaxy, and fitting each mock residual spectrum. For ‘P-Cygni’ cases, we use two double Gaussian profiles to fit the emission and absorption together. If more than half are well fit, that is, the error for the ΔV is determined and the parameter values are not at the boundaries, we consider the flow detectable in that galaxy. As mentioned in subsection 3.1, we assign errors to the residual spectra from the original flux uncertainties and ignore any arising from the stellar population model. This choice does not affect our conclusions; even if we relax the restriction from 3σ to 1σ, we obtain the same main results. In Table 1, we summarize the number of these ‘flow-detectable galaxies’ in the post-starburst, starburst, and quiescent samples.

Table 1 shows that the fraction of flow-detectable galaxies is higher for the post-starburst sample (31 per cent) than for the starburst (11 per cent) and quiescent (6 per cent) samples. The lower flow-detectable rate is driven in part by higher fractions of low S/N spectra in the starburst and quiescent samples, making a higher fraction of the residual |$\rm{Na {}\rm {\small D}}$| line profiles difficult to fit.⁴ The flow-detectable fraction also depends on M_*, which we discuss further in subsection 4.2.

For the flow-detectable galaxies, we measure ΔV and its measurement uncertainty based on the Monte Carlo simulation results of the |$\rm{Na {}\rm {\small D}}$| excess profile fitting. For the ‘P-Cygni’ cases, ΔV is measured from the wavelength shift of the absorption component relative to the galaxy’s systemic velocity. The final uncertainty in ΔV is propagated from the dispersion of the line shift measurements from the Monte Carlo simulations and the error in the pPXF redshift.

4 RESULTS

4.1 ISM Contribution to Na D: EW_{NaD, exc}/EW_{NaD, tot}

The three different galaxy samples in our evolutionary sequence of declining star formation have different EW_{NaD, tot} distributions (panel (a) of Fig. 6). The median values of EW_{NaD, tot} from starburst to post-starburst to quiescent galaxies increase from 0.68|$~\rm{\mathring{\rm A}}$| to 1.99|$~\rm{\mathring{\rm A}}$| to 3.27|$~\rm{\mathring{\rm A}}$| as the dominant stellar population changes from early-type stars to late-type stars and results in stronger |$\rm{Na {}\rm {\small D}}$| stellar absorption. The median fractional error in EW_{NaD, tot} is 23 per cent for the starburst sample, 6 per cent for the post-starburst sample, and 4 per cent for the quiescent sample (panel (b) of Fig. 6).

Panel (c) of Fig. 6 shows the distributions of the EW_{NaD, exc} for the three galaxy samples. All three histograms have peaks at around 0|$~\rm{\mathring{\rm A}}$|⁠. The median EW_{NaD, exc} is -0.2|$~\rm{\mathring{\rm A}}$|⁠, 0.01|$~\rm{\mathring{\rm A}}$|⁠, and 0.07|$~\rm{\mathring{\rm A}}$| for the starburst, post-starburst, and quiescent samples, respectively. The histograms for the starburst and the post-starburst samples have tails towards positive EW_{NaD, exc}, but not for the quiescent sample. The median fractional error on EW_{NaD, exc} is 47 per cent, 25 per cent, and 56 per cent for the starburst, post-starburst, and quiescent samples, respectively, as shown in panel (d) of Fig. 6.

Panel (e) of Fig. 6 shows the distributions of the ratio EW_{NaD, exc}/EW_{NaD, tot} for the three galaxy samples. The EW_{NaD, exc}/EW_{NaD, tot} ratio traces the ISM contribution to the total |$\rm{Na {}\rm {\small D}}$| absorption line. In this panel, we plot only the galaxies with positive EW_{NaD, tot} and EW_{NaD, exc}: 1139 starbursts, 262 post-starbursts, and 41 422 quiescent galaxies. Most (97 per cent) of the quiescent galaxies have EW_{NaD, exc}/EW_{NaD, tot} smaller than 20 per cent, compared 39 per cent for starburst galaxies and 41 per cent for post-starburst galaxies. The median EW_{NaD, exc}/EW_{NaD, tot} is 8 per cent for the quiescent galaxies, 27 per cent for the starbursts, and 25 per cent for the post-starbursts. Given the fact that the median EW_{NaD, tot} increases from starburst to quiescent (see panel (a)), but the peaks of the EW_{NaD, exc} distributions for all three galaxy samples are similar (panel (c)), the low EW_{NaD, exc}/EW_{NaD, tot} ratios typical of the quiescents arise from a low fraction of neutral gas to stars.

4.2 Neutral gas detection versus stellar mass

We test the dependencies of the flow-detectable sample in Table 1 on M_*. The right column of Fig. 7 shows that, for all three types of galaxies, the rate of flow-detectables increases with stellar mass, with most galaxies having log(M_*/M_⊙) > 10. Even controlling for the effects of S/N (see subsection 4.4), the rates are higher in higher mass galaxies.

This result is consistent with the findings of Roberts-Borsani & Saintonge (2018), who found neither inflows nor outflows in their low mass galaxies. This trend arises because |$\rm{Na {}\rm {\small D}}$| is very fragile – its low-ionization potential (5.14 eV) requires it to be shielded from hydrogen-ionizing photons by a large amount of dust. In more massive galaxies, we would expect to detect the |$\rm{Na {}\rm {\small D}}$| ISM feature, because the metal-richness and dust lead to better shielding of |$\rm{Na {}\rm {\small D}}$|⁠. Roberts-Borsani & Saintonge (2018) also found that the |$\rm{Na {}\rm {\small D}}$| residual profiles of their low-mass non-AGN host galaxies (log(M_*/M_⊙) ≤ 10) are typically in emission, but we exclude such ‘Pure Emission’ galaxies from our flow-detectable sample (see subsection 3.3).

The left column of Fig. 7 shows the distributions of redshift for the three galaxy samples; the flow-detectable rate is generally independent of redshift.

4.3 Significant inflows versus outflows

Hereafter we focus only on the 163 post-starburst galaxies, 370 starburst galaxies, and 4598 quiescent galaxies that are flow-detectable. Assuming a possible evolutionary sequence from starburst to post-starburst to quiescent, Fig. 8 shows the histograms of ΔV for the three samples defining this sequence (grey histograms). The histograms for the starburst and post-starburst samples are significantly skewed with a tail at high negative velocity offsets. The mean velocity offsets are |$-84.6\pm 5.9~\rm km s^{-1}$| and |$-71.4\pm 11.5~\rm km s^{-1}$| for the starburst and post-starburst samples, respectively. Alternatively, the velocity offset distribution of the quiescent sample is skewed to the positive side with the mean velocity offset of |$76.6\pm 2.3~\rm km s^{-1}$|⁠. The uncertainty of the mean velocity offset is calculated from the overall dispersion of velocity offset in each galaxy sample.

To test whether these differences are statistically significant, we compare the distribution of observed ΔV with the distribution expected for galaxies with zero offsets, broadened by the 1σ measurement errors (yellow histograms in Fig. 8). For the starburst and post-starburst samples, the ΔV distributions show negative tails that exceed the distributions expected in the absence of outflows. This difference is confirmed by a Kolmogorov–Smirnov (K-S) test. Alternatively, the distribution of ΔV for the quiescent sample shows a significant excess of positive velocity offsets, indicating inflows. To test the effect of underestimating the ΔV errors, we also broaden the zero flow distribution assuming the 1.5σ and the 2σ errors; the excess tails remain and are still significant.

Assuming that the |$\rm{Na {}\rm {\small D}}$| excess represents the average properties from the foreground ISM in the sightline, negative and positive velocities represent ISM bulk flows in opposite directions: negative velocities indicate outflows, while positive velocities indicate inflows⁵ Therefore, the excess negative velocity tails of the starburst and post-starburst samples indicate significant outflow populations, while the positive tail of the quiescent sample suggests a significant inflow population. The numbers of outflows and inflows are shown in Table 1.

4.4 Trends with evolutionary sequence

4.4.1 Velocity offset and outflow fraction from starburst to quiescent

In this section, we test how gas flow properties (i.e. ΔV and the fraction of galaxies with outflows) evolve along the assumed evolutionary sequence of starburst to post-starburst to quiescent. For the post-starburst sample, we further divide it into two groups: young post-starburst galaxies (post-burst age |$\lt 400\, \mathrm{Myr}$|⁠) and old post-starburst galaxies (post-burst age |$\ge 400\, \mathrm{Myr}$|⁠).

As mentioned in subsection 4.3 and shown in Fig. 9, the mean velocity offset changes from negative to positive along the evolutionary sequence, suggesting a transition from outflows to inflows with diminishing SFR and age. The same trend is seen with post-burst age within the post-starburst sample. The value of the old post-starburst sample (⁠|$-27.6\pm 13.3~\rm km s^{-1}$|⁠) is lower than that of the young sample. While we average positive and negative velocities in computing the means, the observed trend is driven by outflows declining along the evolutionary sequence, as shown by a comparison of the panels in Fig. 8. (For the quiescent sample, there is also a contribution from real inflows.) Hereafter we use the terms ‘mean velocity offset’ and ‘mean outflow velocity’ to denote overall flow velocity (in any direction or static) and outflow velocity (negative only), respectively.

To test the effects of the stellar mass and spectral S/N on the mean velocity offset trend in Fig. 9, we further divide our samples into four bins (1) low M_* and low S/N; (2) low M_* and high S/N; (3) high M_* and low S/N; (4) high M_* and high S/N. The stellar mass subgroups are separated using median values, which are log(M_*/M_⊙) = 10.7 for the starburst, young post-starburst, and old post-starburst samples, and log(M_*/M_⊙) = 11.0 for the quiescent sample. The low and high S/N subgroups are separated using the median value S/N = 20.

The trends for the four groups are shown in the left panel of Fig. 10. In all of them, the mean velocity offset changes from negative to positive along the sequence. For the starburst and post-starburst samples, the behaviour of the low and high S/N groups is consistent. For the quiescent sample, the low S/N groups lie closer to zero mean velocity offset; this might be due to the general lack of an ISM in quiescent galaxies, which makes it difficult to perform |$\rm{Na {}\rm {\small D}}$| excess fitting and obtain accurate velocity measurements, especially at low S/N. Or, the less positive mean ΔV could arise physically, from diminished inflows or more outflows.

The groups divided by stellar mass show similar trends to the overall trend in Fig. 10, except that the high M_* groups shift upward at fixed S/N. This shift is larger for the high S/N group, suggesting that it is physical. We recover a similar trend for individual post-starburst galaxies in subsection 4.4.2.

The outflow fraction is defined as the number of galaxies with outflows as a percentage of the number of flow-detectable galaxies. The outflow fractions for the starburst, post-starburst, and quiescent samples are |$76\pm 2~{{\ \rm per\ cent}}$|⁠, |$64\pm 6~{{\ \rm per\ cent}}$|⁠, and 33 per cent (3σ upper limit), respectively. Considering that outflows in the quiescent galaxies are hard to distinguish from static flows (panel c of Fig. 8), the outflow fraction for the quiescent sample is an upper limit.

Fig. 9 displays the overall trend of outflow fraction along the sequence, with the post-starburst sample divided into the young and old subsamples as before. The outflow fraction of the young group is |$70\pm 5~{{\ \rm per\ cent}}$| and that of the old one is |$58\pm 3~{{\ \rm per\ cent}}$|⁠. In a manner similar to the change in the mean velocity offset, the outflow fraction decreases along the sequence. As before, we check the stellar mass and S/N dependencies, which are shown in the right panel of Fig. 10. For all the stellar mass and S/N groups, the outflow fraction shows a declining trend similar to that in Fig. 9. The lower S/N groups have higher (upper limit) outflow fractions for the quiescent sample. As discussed earlier, this effect may arise because more fake ‘outflows’ are detected in quiescents with low S/N.

4.4.2 Velocity offset versus age for post-starbursts

In this section, we focus on the post-starburst phase to consider the correlation between the velocity offset and the post-burst age for individual galaxies. We use the PYTHON Pingouin package.⁶ to calculate the Spearman partial-correlation coefficient that describes the correlation with the effect of controlling variables removed. The controlling variables are carefully considered to ensure that the correlation between velocity offset and post-burst age is not driven by other factors.

Fig. 11 shows the ΔV of the post-starburst sample as a function of stellar mass, post-burst age, redshift, and burst mass fraction. The corresponding Spearman correlation coefficients⁷ are listed in the caption. The figure indicates that the mean ΔV becomes less negative with increasing post-burst age (ρ = 0.22,  2.80σ), suggesting that outflows slow in time. The mean ΔV becomes more negative with increasing redshift (ρ = −0.33,  4.25σ) and marginally more negative with increasing stellar mass (ρ = −0.15,  1.87σ). There is no observed dependence of mean ΔV with burst mass fraction (ρ = 0.03,  0.34σ).

For the partial-correlation coefficient, only factors correlated with both post-burst age and mean ΔV need to be considered as controlling variables. Therefore, we do not consider the burst mass fraction because it is ∼constant with ΔV. We then test the correlations among stellar mass, redshift, and post-burst age, which are shown in Fig. 12. The post-burst age has significant anticorrelations with stellar mass (ρ = −0.35,  4.56σ) and redshift (ρ = −0.29,  3.78σ). There is also a correlation between stellar mass and redshift (ρ = 0.64,  9.27σ), which is caused by the magnitude limit of the SDSS sample, where galaxies at higher redshift will have higher average stellar masses. The strong correlation between stellar mass and redshift suggests that the correlation of post-burst age with stellar mass, and the correlation with redshift, are highly overlapping. Therefore, we only consider stellar mass, which is more likely to be physically correlated with galaxy internal properties, as a controlling variable.

The Spearman partial-correlation coefficient is applied to characterize the correlation between ΔV and post-burst age after controlling for stellar mass. The final partial-correlation coefficient is 0.18, with 2.29σ significance. Therefore, there is a significant correlation between ΔV and post-burst age, even after controlling for stellar mass. This residual correlation implies that outflows diminish during the evolution of the post-starburst phase.

4.5 Trends with disc inclination for starbursting galaxies

In this section, we test the inclination dependency of the mean velocity offset and the outflow fraction to explore the possible mechanism for outflows in starburst galaxies, the only one of our three galaxy samples with a large number of measured disc inclinations. Galactic outflows usually show a biconical structure, that is, the alignment of outflows is perpendicular to the stellar or AGN accretion disc, and the outflow can only be detected within the opening angle. Therefore, the measured wind velocity along the line-of-sight depends on the viewing angle. The fastest velocity would be measured when the viewing angle is zero, that is, aligned with the outflow. The velocity decreases with the viewing angle due to the projection effect, until the viewing angle exceeds the outflow opening angle and the outflow can no longer be detected.

How the observed velocity offset and outflow fraction vary with galactic disc inclination depends on the origin of the winds (e.g. Bae & Woo 2018; Concas et al. 2019). These papers find that the velocities of ionized outflows from AGN do not depend on host galaxy inclination, whereas star formation-driven winds are observed to lie principally perpendicular to the disc. Star formation-driven outflows, which are associated with supernovae and stellar winds, are initially isotropic and then preferentially escape along the galaxy’s minor axis, the path of least resistance. In contrast, outflows driven by an AGN jet are perpendicular to the black hole’s accretion disc, which can be randomly oriented relative to the host galaxy’s disc (Veilleux et al. 2005). Assuming that an AGN-driven outflow is not dramatically affected by the disc and retains much of its initial orientation and velocity (Tanner & Weaver 2022; but also see Hartwig, Volonteri & Dashyan 2018 and Nelson et al. 2019), there will be no dependence of outflow velocity on disc inclination for a statistical sample.

The SDSS photometric pipeline (Lupton et al. 2001) fits each galaxy image with a 2D model of a de Vaucouleurs and an exponential profile. The axial ratios (b/a) from the two models are derived by this fitting procedure, but we only use the exponential profile axial ratio below. The pipeline also calculates the best linear combination of the exponential and de Vaucouleurs models as a parameter called fracDeV. Because inclination is only meaningful for disc galaxies, we require fracDeV < 0.8, which corresponds to the disc component being more than 20 per cent of the total light. In this choice, we follow Padilla & Strauss (2008) and Chen et al. (2010). The b/a is converted to the inclination using

where q = 0.13 is the intrinsic axis ratio (Giovanelli et al. 1994). The typical b/a error is <0.1, bringing the error in inclination to about a few degrees. Due to our requirements here that a galaxy has a measured b/a and fracDeV, and then a confirmed disc, we are only able to measure inclinations for 255 starburst and 18 post-starburst galaxies from our samples.

The correlation between ΔV and inclination for the starburst sample is shown in the left panel of Fig. 13. The mean velocity offset of the starburst sample decreases with inclination and becomes zero when nearly edge-on (i ∼ 80°). We also do linear regression on this trend using PYTHON LINMIX package⁸, which confirms that the declining slope is significantly(>6σ) different than zero.

The declining velocity trend with inclination implies that the alignment of outflows detected in the starburst galaxies is more out of the disc, favouring the star formation-driven wind scenario. It is not surprising that star formation drives winds in our starburst sample – we select for star-forming galaxies in the BPT diagram and against even weak AGN features. Concas et al. (2019) also explore the inclination dependence of outflow velocity, but for galaxies with higher SFRs than ours: a lower-limit of 12.5 versus an average of 6.6 M_⊙ yr⁻¹ (see Fig. 2), respectively. (Some of their galaxies include AGN.) Another difference is that they measure their wind velocities from stacked spectra, instead of for individual galaxies as we do. The slope of their best wind velocity-inclination relation is nevertheless consistent with ours.

For our post-starburst galaxy sample, there are relatively fewer discs, and those with measured inclinations are concentrated at high values. As a result, we cannot determine their ΔV-inclination relation and test for the dominant wind-driving mechanism via wind geometry. We note, however, that the discy post-starbursts in our sample do not obviously deviate from the relation defined by our starburst sample.

The outflow fraction trend with inclination for the starburst sample is displayed in the right panel of Fig. 13. Considering the inclination error, we use an inclination bin size =10°. At i ≤ 45°, the outflow fraction is constant (∼90 per cent), and outflows are ubiquitous in starburst galaxies. For i > 45°, the outflow fraction steadily decreases. Because randomly distributed outflow alignments would allow us to detect outflows at any inclination with equal probability, the alignment of neutral outflows seen here is more likely to be along the minor axis.

The observed outflow fraction trend with inclination also provides some constraints on the distribution of outflow half-opening angles. First, we assume simplistically that the outflow alignments are exactly perpendicular to the disc. We then assume that the distribution of the half-opening angles is uniform and not related to the inclination. Because we observe that outflows are ubiquitous for i ≤ 45°, there can be no winds with half-opening angles less than 45°, else the outflow fraction would decrease before 45°. Thus, the differences in outflow fraction between adjacent inclination intervals when i > 45° represent the proportions of outflows with inclinations in each interval. The predicted distribution of half-opening angles is shown in Fig. 14. As expected, there are no galactic neutral gas outflows with half-opening angles less than 45°. Propagating the outflow fraction uncertainty to the half-opening angle distribution, the typical uncertainty in each bin is ∼0.1(10 per cent). Therefore, for the distribution of broader half-opening angles (>45°), we do not observe any significant patterns.

5 DISCUSSION

5.1 Outflow versus SFR evolution in post-starbursts

The results in subsection 4.4.2 show that winds are declining over the star formation declining sequence and with post-burst age. To further understand the connection between the evolution of outflows and star formation activity, that is, how they affect each other, we fit the outflow and SFR trends for the post-starburst galaxy sample. Here we focus on the outflow-detected post-starburst sample, given that only outflows contribute to the expulsion of gas from the galaxy and could thereby reduce star formation. We fit the post-starburst age–ln (outflow velocity) relation using PYTHON LINMIX package⁹ in Fig. 15. The best fit exponential time-scale for the wind decline is 550 Myr, with an uncertainty of about 200 Myr.

The best fit time-scale for the decline in SFR over the same period is 200 Myr, with an uncertainty of ∼50 Myr. To test the uncertainty introduced by the AGN contribution correction to the SFR measurements, we refit for the SFR time-scale without the AGN correction and obtain consistent results. French et al. (2023) found that the decline in SFR for their post-starburst galaxies occurs over two different time-scales: early (65 Myr) and late (480 Myr), with a turning point at 77 Myr. However, their sample included many galaxies younger than ours here – our youngest post-starburst galaxy with a SFR measurement is 83 Myr, older than French et al. (2023)’s turning point age. Compared with the SFR decline time-scale, the outflow velocity time-scale is longer, suggesting that the outflows are decoupled from the current (low levels of) star formation. This decoupling suggests that outflows in post-starburst galaxies are not driven primarily by stellar winds or supernovae.

5.2 Fate of post-starburst outflows

In this section, we test whether the |$\rm{Na {}\rm {\small D}}$| outflows in the post-starburst galaxy sample have the capability of removing gas from the ISM of the galaxy, which is roughly the scale probed by the SDSS fiber over this redshift range. We compare wind velocity to ISM galaxy escape velocity in Fig. 16. The ISM escape velocity is estimated as V_esc ≈ (2GM_fib/r_fib)^1/2, where M_fib is the stellar mass within the 3” SDSS fiber from the MPA-JHU catalogue, and R_fib is the physical scale of the fiber at the galaxy’s redshift (0.3–7 kpc). We separately estimate the escape velocity necessary to leave the galaxy’s halo (i.e. the circumgalactic medium (CGM) escape velocity) as V_{h, esc} ≈ 3V_{h, vir} (e.g. Weiner et al. 2009), where V_{h, vir} is the halo virial velocity. We convert stellar mass into the halo mass using the relation from Behroozi et al. (2019) and then calculate the halo virial velocity using the relation from Maller & Bullock (2004).

Only 9/105 (9 per cent) of post-starburst galaxies have fast enough outflows to escape the ISM, that is, have wind velocities 3σ higher than the ISM escape velocity. Only 3/9 ISM-escaping outflows can also escape from the halo. To consider the sightline projection effects, which may cause us to underestimate the actual wind velocity, we randomly select projected angles θ from 0° to 90° from a uniform distribution over 1000 trials. Then we perform a projection correction (ΔV_real = ΔV/cos θ). We recalculate the percentage of ISM-escaping winds, finding |$\sim 30~{{\ \rm per\ cent}}$|⁠, with an upper limit of 50 per cent. This result confirms statistically that most of the |$\rm{Na {}\rm {\small D}}$| outflows detected in our post-starburst sample cannot expel gas from the galaxy. This gas would fall back into the disc, perhaps in ‘galactic fountains’.

Even if many winds cannot efficiently remove gas, they might still affect any gas reservoirs within the galaxy. Smercina et al. (2022) found high internal turbulent pressure, revealed by compact CO gas reservoirs and high-velocity dispersions, in six gas- and dust-rich post-starburst galaxies, and French et al. (2023) detected low molecular gas excitation, suggesting that lower gas densities may suppress star formation during the post-starburst phase. Modest outflows within the ISM could both trigger high turbulence and prevent the cold gas from collapsing to form new stars. While these two phenomena were observed for molecular gas, our neutral outflows could also be associated with them, as multiphase outflows are very common in post-starburst galaxies (Baron et al. 2021; Luo et al. 2022).

5.3 Origins of post-starburst outflows

5.3.1 Associated with star formation

As we mentioned in subsection 5.1, the outflow velocities in post-starburst galaxies diminish over a time-scale that is longer than that over which the SFR declines. To further test if this decoupling of winds and SFR disfavours a stellar origin for the winds, we compare the correlation between outflow velocity and SFR to the predictions of the simple analytic model for supernova explosion momentum-driven outflows introduced by Xu et al. (2022). This model describes the behaviours of clouds under a combination of outward momentum from the starburst and the inward force of gravity. It also assumes that the gas we see is swept up, ambient gas at the surface of the expanding wind bubble driven by the wind momentum flux (see subsection 6.2 of Xu et al. 2022 for details). The predicted scaling relation between outflow velocity and SFR is then log (outflow velocity) ∝ 0.25log (SFR). Xu et al. (2022) detected ionized gas outflows in a sample of 45 low-redshift (z < 0.2) star-forming galaxies; their best fit slope to the ionized outflow velocity versus SFR relation (0.224 in the logarithmic scale) is consistent with the model.

If the detected winds in our post-starburst galaxies are driven by supernovae, and cold neutral clouds are coupled to hot ionized winds, we would expect a similar outflow velocity-SFR relation. However, Fig. 17 shows that the model slope is somewhat steeper than for our starburst (slope =0.05 ± 0.04) and post-starburst (0.10 ± 0.09) samples, which have slopes consistent with zero. The apparent discrepancy with the model might arise if the neutral phase of star formation-driven outflows has a weaker SFR dependence than the ionized phase. Another possibility is that the outflows arise mostly from a mechanism other than supernovae. At least for the starburst galaxies, the inclination dependence discussed in subsection 4.5 still implies that star formation-driven winds play a role.

There is another interpretation to consider. Even if the outflows detected in post-starbursts are unrelated to current star formation, the cool neutral gas in the flow may be long lived, representing a relic outflow (e.g. Lochhaas et al. 2018) launched near the peak starburst time that has lingered in the ISM. In this scenario, the covering factor of the outflowing cold clouds as seen by the galaxy decreases in time. However, we find that EW_{NaD, exc}, which is proportional to the covering factor of the neutral gas, is not correlated with post-burst age for our sample, arguing against this interpretation.

5.3.2 Associated with nuclear activity

AGN are likely to contribute to the outflows seen in post-starburst galaxies. Baron et al. (2021) detected neutral outflows in 40 per cent of their 144 post-starburst galaxies with AGN and ionized outflows. The median velocity of neutral outflows in their post-starburst galaxies is |$-633~\rm km s^{-1}$|⁠, much higher than the mean ΔV of our post-starburst sample (⁠|$-71.4~\rm km s^{-1}$|⁠). This difference is likely due to our selection against significant AGN spectral features and to our averaging over all bulk motions, including inflows. If the difference arises from sample selection, it suggests that strong AGN (included in the Baron et al. (2021) sample) produce faster outflows than those driven by weaker AGN or by stellar winds alone (our sample).

Even though our post-starburst sample is biased against strong AGN, nuclear activity may have been stronger at earlier times. Wild et al. (2010) explored the growth of black holes in 400 local galactic bulges that have had a starburst in the past 600 Myr and found that the AGN activity peaks 250 Myr after the starburst¹⁰ On the simulation side, Byrne-Mamahit et al. (2023) used IllustrisTNG to test the supermassive black hole (SMBH) accretion rate and SFR enhancement in post-merger galaxies – post-starburst galaxies often result from mergers – and found that accretion rate enhancements persist for up to 2 Gyr after coalescence. Keel et al. (2012) found that the time-scale for an AGN to vary from quasar-like to LINER-like is much faster than that of the post-starburst phase. Therefore, while the relative time-scales among the starburst, nuclear activity, and outflows are not well constrained, outflows in post-starburst galaxies could be the relics of AGN in the starburst or early post-starburst phase.

To further test the hypothesis that AGN is responsible for the post-starburst winds, we place our post-starburst galaxies on the |$[\rm{S {}{\small II}}]$|-BPT diagram (Veilleux & Osterbrock 1987) in Fig. 18.¹¹The BPT emission lines come from the MPA-JHU catalogues (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). Of the 516 post-starburst galaxies, 336 have positive BPT emission line measurements from the MPA-JHU catalogues, including 61 outflow-detected galaxies.

Although we select against galaxies with strong emission lines, many post-starburst galaxies in our sample still have weak AGN or LINER features after stellar spectral continuum subtraction. Given the high Seyfert/LINER fraction of outflow-detected galaxies (53/61, 87 per cent), it is natural to assume that the AGN and outflows are connected. However, to fully test this point requires that we check whether the chance of a Seyfert/LINER post-starbursts having wind is significantly higher than expected from the fraction of all post-starbursts with winds (whether we can place them on the BPT diagram or not). For our sample, this assumption yields that the outflow detection rate in Seyfert/LINER post-starbursts is not statistically higher than the rate in all post-starbursts; both rates are roughly 20 per cent. So, the high AGN and high outflow fractions in our sample do not in themselves lend support to the picture of AGN-driven winds.

If we instead focus on the four with ISM-escaping outflows (see subsection 5.2), we find that three lie in the Seyfert region, including both galaxies with CGM-escaping outflows. (The remaining galaxy with an ISM-escaping outflow is a LINER.) Compared to the Seyfert-only fractions of the outflow-detected post-starburst sample (24/61, 39 per cent) and of the entire post-starburst sample (124/336, 37 per cent), the Seyfert fraction of post-starbursts with ISM-escaping outflows (75 per cent) is higher. Together with the observed differences in outflow velocity between the Baron et al. (2021) sample and ours, this result, while not statistically significant, suggests that AGN could drive the most effective outflows in post-starburst galaxies.

Further indirect evidence for AGN-driven outflows in post-starbursts is from the correlation between ΔV and M_* (subsection 4.4.2). This correlation may be a secondary correlation caused by the real ΔV-post-burst age relation (panel (c) of Fig. 11) and M_*-post-burst age relation (panel (a) of Fig. 12). It is also possible that the ΔV–M_* relation is driven by trends with black hole mass, that is, more massive galaxies have more massive black holes (BHs) that drive faster outflows. To test this second possibility, we check whether ΔV correlates with bulge mass (M_bulge) as it does with stellar mass. We obtain M_bulge values from Mendel et al. (2014). Fig. 19 shows that the ΔV–M_bulge correlation has similar strength and significance (ρ = −0.14,  1.78σ) to the ΔV–M_* relation (panel (a) in Fig. 11). Given that bulge mass is strongly correlated with black hole mass (M_BH), the observed ΔV-M_* relation could be caused by the ΔV–M_BH relation, a finding consistent with the AGN-driven outflow scenario.

In addition to the possibility of AGN outflows, the high tidal disruption event (TDE) rate in post-starburst galaxies (French, Arcavi & Zabludoff 2016; French et al. 2020) and the likely resulting energy input (Mockler & Ramirez-Ruiz 2021) could generate low-level nuclear winds that contribute to the outflows we see here. Smercina et al. (2022) makes a similar argument for TDEs as the source of the turbulent molecular gas reservoirs they observed in post-starburst galaxies.

6 CONCLUSIONS

In this paper, we use the |$\rm{Na {}\rm {\small D}}$| absorption doublet to analyse bulk neutral gas motions in the ISM of ∼80 000 0.010 < z < 0.325 galaxies from the SDSS DR12 survey. We first construct a galaxy evolutionary sequence in declining SFR, from starburst to post-starburst to quiescence. To construct this sequence, we find likely progenitors (starburst galaxies) and descendants (quiescent galaxies) of the post-starburst galaxies in F18 by constraining their stellar masses and SFRs. To isolate any gas flow contribution, we subtract the best fit stellar continuum from the spectrum to obtain a residual |$\rm{Na {}\rm {\small D}}$| line that arises from the ISM. We then measure the positive (inflow) or negative (outflow) velocity shift ΔV of this line relative to the systemic velocity of the stars.

Our main findings are

• bulk flows – in or out of the galaxy – are detected along the entire evolutionary sequence, especially at higher host stellar masses, that is, log(M_*/M_⊙) > 10. Both the starburst and post-starburst samples have ΔV distributions with statistically significant negative tails, implying that we have detected many real outflows. The quiescent sample has a positive ΔV tail that reveals significant inflows. These inflows are associated with smaller EW residual |$\rm{Na {}\rm {\small D}}$| lines than the starburst and post-starburst outflows, indicating lower mass loading; this result may explain why, despite inflows being common in our ‘red, dead’ galaxies, they do not trigger significant star formation.

• We see further evidence for outflows in the |$\rm{Na {}\rm {\small D}}$| line profiles. Among those hosts with bulk flows, 9/370 (2.4 per cent) of starbursts, 11/163 (6.7 per cent) of post-starbursts, and 13/4598 (0.3 per cent) of quiescents have P-Cygni profiles that indicate outflows. As a consistency check on our fitting methodology, we confirm that all of these systems are classified as having outflows, that is, negative velocity offsets.

• For all three galaxy samples, the mean outflow velocities at high M_* (log(M_*/M_⊙) ≥ 11) are faster than at low M_* (log(M_*/M_⊙) < 11). For the post-starbursts, outflow velocity increases significantly with M_*. It is possible that this trend is tied instead to bulge mass and thus black hole mass, as we find a similarly strong ΔV-M_bulge relation for post-starbursts and independent bulge mass estimates.

• Outflows diminish as galaxies age and star formation declines. The fraction of hosts with outflows decreases along the galaxy evolution sequence from starbursts (⁠|$76\pm 2~{{\ \rm per\ cent}}$|⁠) to post-starbursts (⁠|$64\pm 6~{{\ \rm per\ cent}}$|⁠) to quiescents (33 per cent, 3σ upper limit). Similarly, the mean velocity offset of the gas changes from −84.6 ± 5.9 to −71.4 ± 11.5 to 76.6 ± 2.3, that is, from outflowing to inflowing. Even within the post-starburst sample, and after controlling for host stellar mass, the mean flow speed decreases with time elapsed since the starburst ended.

• For the post-starburst sample, the SFR declines more quickly with post-burst age than the outflow velocity does, suggesting that stellar winds are not the principal driver of the outflows.

• For the post-starbursts, only a small fraction of measured outflows significantly exceed the escape velocity of the ISM defined by the SDSS fiber aperture (9/105) or the circumgalactic medium defined by virial radius of the dark matter halo (3/105). Although we measure only the line-of-sight component of the wind, most (⁠|$\sim 70~{{\ \rm per\ cent}}$|⁠) winds are still statistically likely to remain within the host galaxy.

• For the post-starbursts, 3/4 (75 per cent) with ISM-escaping outflows, including both of the galaxies with CGM-escaping outflows, that can be placed on the BPT diagram lie in the Seyfert region (despite our post-starburst sample selection against H α emission). In comparison, the Seyfert fraction of the entire post-starburst sample is 124/336 (37 per cent) and that of the subsample with outflows of any kind is 24/61 (39 per cent).

• For the starburst galaxies with discs, the outflows are dependent on disc inclination, which favours a wind model in which gas is ejected perpendicularly from the disc plane. The wind velocity decreases as the disc becomes more edge-on, and the outflow fraction, constant at |$\sim 90~{{\ \rm per\ cent}}$| for disc inclinations i < 45°, steadily decreases from |$\sim 90~{{\ \rm per\ cent}}$| to 20 per cent for i > 45°. Our analysis limits the outflow 1/2-opening angle to roughly >45°.

There are conflicting clues as to the nature of the winds. For the post-starburst galaxies, the high Seyfert fraction associated with the strongest outflows suggests an AGN connection to gas removal from these galaxies. Another point favouring a nuclear over stellar origin for post-starburst outflows is the different time-scales over which the SFR and wind speed decline. The post-starburst ΔV-M_bulge relation observed here is also consistent with the AGN-driven picture if M_bulge is a proxy for black hole mass. Furthermore, the ΔV-SFR relations for starbursts and post-starbursts are marginally inconsistent with the prediction of a simple supernovae momentum-driven wind model. Alternatively, for starburst galaxies with discs, the trends of inclination versus ΔV and outflow fraction support a perpendicular wind model, a geometry most likely to arise from star formation throughout the disc.

Critically, the existence of outflows does not demonstrate ‘feedback’, as it is not yet clear how (or even if) these winds regulate star formation in their hosts. The intriguing drop observed here in outflow speed and outflow fraction with SFR along the starburst/ post-starburst/quiescent sequence does not imply causation. Given their slow speeds, many of these outflows are likely to remain bound to their host, unable to remove its gas. Yet, outflows may reduce the SFR by generating turbulence or re-distributing the ISM.

To better understand the whole picture of outflows – how they affect and are affected by SFH – a similar analysis is required within individual galaxies. Spatially resolved spectroscopy, for example from the SAMI IFS data (Croom et al. 2021) and the MaNGA IFU data (Bundy et al. 2015), across the galaxy evolution sequence defined here could reveal the source of the outflows by localizing them to the disc or to the galactic center. Resolving the spatial extent of the outflows will also constrain the mass-/energy-loss (Roberts-Borsani et al. 2020; Avery et al. 2021, 2022), quantifying the winds’ impact on the ISM and CGM, a true test of feedback processes.

ACKNOWLEDGEMENTS

We thank the anonymous referee for a thoughtful report that improved this paper. We also thank George Rieke and Dennis Zaritsky for helpful discussions. YS and AIZ acknowledge support from NASA ADAP grant number 80NSSC21K0988. AIZ also thanks the hospitality of the Columbia Astrophysics Laboratory at Columbia University, where some of this work was completed. YY is supported by the Basic Science Research Programme through the National Research Foundation of Korea funded by the Ministry of Science, ICT & Future Planning (2019R1A2C4069803).

DATA AVAILABILITY

The data underlying this article were accessed from the publicly available Data Release 12 (DR12; Alam et al. 2015) of the Sloan Digital Sky Survey (SDSS-III; Eisenstein et al. 2011). Specifically, we make use of the SDSS optical spectra. See the SDSS web page (https://dr12.sdss.org/) for more details of the DR12 spectra. All spectra can be accessed using the traditional SDSS Science Archive Server (SAS) at https://data.sdss.org/sas/dr12/sdss/spectro/.

We have tabulated the measured properties of our starburst, post-starburst, and quiescent galaxy samples, including the |$\rm{Na {}\rm {\small D}}$| EW and velocity offsets derived from the optical spectra. We also include the inclinations of starburst galaxies and H α-based SFR measurements for the post-starburst and quiescent galaxies. We also tabulate the redshifts, stellar masses, and SFRs of starburst galaxies from the MPA-JHU catalogues, and the post-burst ages and burst mass fractions of the post-starbursts from French et al. (2018). These data will be made publicly available through the online supplementary material upon the acceptance of the manuscript.

Footnotes

References