Signatures of Inflowing Gas in Red Geyser Galaxies Hosting Radio Active Galactic Nuclei

We use data primarily from the recently completed SDSS-IV MaNGA survey (Blanton et al. 2017; Bundy et al. 2015; Drory et al. 2015; Law et al. 2015; Yan et al. 2016; Albareti et al. 2017). MaNGA is an integral field spectroscopic survey that provides spatially resolved spectroscopy for nearby galaxies (z ∼ 0.03) with an effective spatial resolution of 2.5'' (FWHM). The MaNGA survey uses the SDSS 2.5 m telescope in spectroscopic mode (Gunn et al. 2006) and the two dual-channel BOSS spectrographs (Smee et al. 2013), which provide continuous wavelength coverage from the near-UV to the near-IR: 3600–10000 Å. The spectral resolution varies from R ∼ 1400 at 4000 Å to R ∼ 2600 at 9000 Å. An r-band signal-to-noise ratio (S/N) of 4–8 Å⁻¹ is achieved in the outskirts (i.e., 1–2 R_e) of target galaxies with an integration time of approximately 3 hr. MaNGA has observed more than 10,000 galaxies with $\mathrm{log}\,({M}_{* }/{M}_{\odot })\gtrsim 9$ across ∼2700 deg² over its 6 yr duration. In order to balance radial coverage versus spatial resolution, MaNGA observes two-thirds of its galaxy sample to ∼1.5 R_e and one-third to 2.5 R_e. The MaNGA target selection is described in detail in Wake et al. (2017).

The raw data are processed with the MaNGA Data Reduction Pipeline (Law et al. 2016). An individual row-by-row algorithm is used to extract the fiber flux and derive inverse variance spectra from each exposure, which are then wavelength-calibrated, flat-fielded, and sky-subtracted. We use a MaNGA sample and data products drawn from MaNGA Product Launch 9 (MPL-9) and Data Release 16 (Ahumada et al. 2020). We use spectral measurements and other analyses carried out by the MaNGA Data Analysis Pipeline (DAP), specifically version 2.3.0. The data we use here are based on DAP analysis of each spaxel in the MaNGA data cubes. The DAP first fits the stellar continuum of each spaxel to determine the stellar kinematics using the penalized pixel-fitting algorithm pPXF (Cappellari & Emsellem 2004; Cappellari 2017) and templates based on the MILES stellar library (Falcón-Barroso et al. 2011). The templates are a hierarchically clustered distillation of the full MILES stellar library into 49 templates. This small set of templates provide statistically equivalent fits to those that use the full library of 985 spectra in the MILES stellar library. The emission-line regions are masked during this fit. The DAP then subtracts the result of the stellar continuum modeling to provide a (nearly) continuum-free spectrum that is used to fit the nebular emission lines. This version of the DAP treats each line independently, fitting each for its flux, Doppler shift, and width, assuming a Gaussian profile shape. The final output from the DAP is gas and stellar kinematics, emission-line properties, and stellar absorption indices. All spatially resolved 2D maps shown in the paper are outputs from the DAP with a hybrid binning scheme. An overview of the DAP used for Data Release 15 and its products is presented by Westfall et al. (2019), and assessments of its emission-line fitting approach are described by Belfiore et al. (2019).

We use ancillary data drawn from the NASA–Sloan Atlas ¹⁶ (NSA) catalog, which reanalyzes images and derives morphological parameters for local galaxies observed in SDSS imaging. It compiles spectroscopic redshifts, UV photometry (from the Galaxy Evolution Explorer (GALEX); Martin et al. 2005), stellar masses, and structural parameters. We have specifically used spectroscopic redshifts and stellar masses from the NSA catalog.

Red geysers are early-type galaxies lying in the red sequence (NUV − r > 5) that have little or no ongoing star formation activity. Figure 1 shows the star formation rates (SFRs) of the red geyser sample (red circles) compared to all MaNGA galaxies (gray) from the Salim et al. (2016) catalog. This catalog derived SFRs from UV/optical spectral energy distribution fitting with additional constraints from mid-IR data using GALEX+SDSS+WISE. The typical value of the SFR for the red geysers is ∼0.01 M_⊙ yr⁻¹. The sample of red geysers used here is derived from MPL-9 and consists of 140 galaxies, which account for ≈6%–8% of the local quiescent galaxy population observed by MaNGA. The red geyser sample is visually selected based on the characteristic features described in detail in Cheung et al. (2016), Roy et al. (2018), and Roy et al. (2021). They are briefly outlined below:

• 1.  
  Spheroidal galaxies with no disk component as observed by SDSS and with low SFRs
• 2.  
  Bisymmetric feature in spatially resolved EW map of strong emission lines like Hα, [N ii], and [O iii]
• 3.  
  Rough alignment of the bisymmetric feature with the ionized gas kinematic axis, but misalignment with the stellar kinematic axis
• 4.  
  High spatially resolved gas velocity values (typically reaching a maximum of ±300 km s⁻¹) that are greater than the stellar velocity values by at least a factor of 4–5
• 5.  
  High gas velocity dispersion values, reaching ∼220–250 km s⁻¹ in distinct parts of the galaxy
• 6.  
  Predominantly showing LINER or Seyfert-type line ratios in spatially resolved BPT diagrams

Four example red geysers are shown in Figure 2. Optical images (left panels) from SDSS show the spheroidal morphology typical of these galaxies. The middle panels show the characteristic bisymmetric feature in the H EW map. The right panels show the NaD EW map obtained using the MaNGA DAP (see Section 10 in Westfall et al. 2019 for details). The DAP performs spectral index measurements after first subtracting the best-fit emission-line model from the observed spectrum. It does not perform any stellar continuum subtraction on the NaD feature; hence the observed EW shows the presence of both the stellar and the ISM component. However the NaD maps show a highly anisotropic and clumpy spatial morphology, indicating that a major portion of the observed NaD absorption is likely contributed by ISM gas. We also find that the regions with large NaD EWs are spatially offset from the bisymmetric feature observed in the H map suggesting the presence of multiphase gas components tracing different activities in these galaxies.

Zoom In Zoom Out Reset image size

Roy et al. (2018) cross-matched the red geyser sample with the VLA-FIRST catalog to identify the radio-detected red geysers. For the latest sample from MPL-9, we have an updated list of radio-detected red geysers consisting of 30 galaxies.

We construct a control sample of spheroidal quiescent galaxies that are matched in global properties, namely stellar mass, color, redshift, and axis ratio, but do not show resolved geyser-like features as described previously. For each red geyser, we match up to four quiescent galaxies (having NUV − r > 5) with the following criteria:

• 1.  
  log M_{⋆,red geyser}/M_⋆,control < 0.2 dex
• 2.  
  z_{red geyser} – z_control < 0.01
• 3.  
  b/a_{red geyser} – b/a_control < 0.1

where M_⋆ is the stellar mass, z is the spectroscopic redshift, and b/a is the axis ratio from the NSA catalog. This technique results in 458 unique control galaxies with 65 of them detected in FIRST. These constitute the radio-detected control sample.

The sodium doublet feature can arise from photospheric transitions in evolved cool stars, as well as from cool neutral gas in the ISM. Since our samples consist of red and quiescent galaxies with evolved stellar populations, the K–M-type stars throughout the galaxy are likely to make up a major fraction of the spectral NaD profile (Jacoby et al. 1984; Tyson & Rich 1991) with the residual absorption attributed to the ISM component. On the other hand, since the primary contribution in the magnesium triplet Mg ii λλλ5167, 5173, 5184 (Mg b) absorption is from the stars, studies of star-forming galaxies have found a linear relation between the EWs of Mg b and NaD from the stellar component alone (Alatalo et al. 2016). In spite of the observed scatter around the empirical relation in Alatalo et al. (2016), any enhancement in NaD as compared to Mg b beyond the relation likely indicates an ISM contribution.

Figure 3 (upper panel) compares the EWs of the NaD and Mg b of red geysers (cyan and magenta squares, representing radio-detected and non-radio-detected galaxies from the FIRST survey, respectively) and the control sample (gray circles), obtained by averaging the EWs measured by the MaNGA DAP (which include NaD absorption from both the stellar and ISM components) over all spaxels within one effective radius. The figure shows that a large fraction of the galaxies in the radio-detected red geyser sample have a considerable excess in NaD lying above the Alatalo et al. (2016) empirical relation (blue dashed line).

Zoom In Zoom Out Reset image size

This is further confirmed by the lower panel in Figure 3, which shows the probability density function of the deviation of the NaD EW distribution with respect to the Alatalo et al. (2016) relation. We find that a large fraction of radio-detected red geysers show an excess amount of NaD, as indicated by the clear positive wing in the distribution, as compared to the other two samples. In fact, 75% of the radio red geyser sample lie above the relation, compared to 51% of the non-radio-detected red geysers and only 32% of the control sample. This indicates an elevated abundance of cool gas within radio-detected red geysers and a possible correlation with the presence of radio-detected AGNs, which we will return to below.

It is to be noted from Figure 3 that there exists a notable amount of NaD deficit at low Mg b particularly for the control galaxies. This is likely caused by the infilling of the resonant NaD emission line, which has been previously detected in integrated, stacked, and spatially resolved NaD spectra (Alatalo et al. 2016; Chen et al. 2010; Rupke et al. 2021). The non-radio-detected red geysers, on the other hand, show a roughly equal fraction of galaxies with NaD deficit (at low Mg b) and excess (at high Mg b) with a rather symmetric distribution overall.

The next step is to accurately model out the stellar continuum contribution to the NaD absorption in order to extract the residual gaseous component. First, we mask spaxels with low or no spectral coverage, low S/N, or foreground stars, unreliable fits, or noisy observations as flagged by MaNGA under the mask-bit names NOCOV, LOWCOV, FORESTAR, UNRELIABLE, and DONOTUSE. Then, following K. Rubin et al. (2021 in preparation), we bin and stack the spectra in 2'' × 2'' spatial bins (or spaxels) to acquire a good continuum S/N (>10 per dispersion element) per bin throughout the galaxy.

We then fit the stellar continuum for each spatial bin with pPXF (Cappellari & Emsellem 2004; Cappellari 2017), which fits a linear combination of simple stellar population (SSP) model templates. Here, we make use of the MIUSCAT SSP (Vazdekis et al. 2012), although a detailed comparison with other SSP models is presented later. Strong emission line regions are masked during the fit. We also mask the region around the NaD transition, as we assume that the doublet is a result of stellar plus ISM contributions and hence should not be fit by the model. We also mask the red half of the He i emission line at 5875.67 Å, which is close enough to the NaD line that it could affect the residual profile. An example spectrum and its continuum fit are shown in Figure 4. The stellar continuum model shows a satisfactory fit, with minimal residuals around Mg b and an excess around NaD, which hint at the presence of ISM gas.

Zoom In Zoom Out Reset image size

To remove the stellar continuum and isolate the ISM component, the spectrum in each bin is divided by the best corresponding continuum fit. We additionally fit a first-order polynomial to the continuum-normalized spectrum in the wavelength range immediately (20 Å) blueward and redward of the NaD profile and divide the residual by the polynomial fit. This is to account for any systematic errors in the continuum fitting that might give rise to artificial residuals.

We now examine the robustness of our stellar continuum modeling to the assumed stellar library. We repeat the above continuum fitting process for 10 galaxies in each of the red geyser and control samples across every spatial bin using two additional libraries, namely MaStar (Yan et al. 2019), an empirical library constructed using the MaNGA instrumentation, and BPASS (Eldridge et al. 2017), a theoretical template library. The BPASS SSP models are drawn from the v2.2.1 release, which adopted the default BPASS initial mass function. ¹⁷ SSPs over a range of log ages (in gigayears) of −1.2 to 1.1, with stellar metallicity mass fractions Z = 0.008, 0.02, and 0.03, are included. The BPASS fitting will be described in full in K. Parker et al. (2021, in preparation).

An example of the variation of the stellar contribution owing to the different stellar population models is shown in the left panel of Figure 5. The different colors indicate the SSPs used to construct the stellar model. The right panel shows the continuum-normalized NaD spectra using the same SSPs. We find that the fractional differences in NaD EW of the continuum-normalized spectra, computed using the MaStar and MIUSCAT models, range between 5% and 15% in the 10 red geyser galaxies studied. However, the BPASS model can change the estimated EW by about 25%–30%, as compared to the other two models. BPASS is the result of combining theoretical model spectra with as few empirical inputs as possible. Hence it differs substantially from the two other stellar models (MIUSCAT and MaStar), which are based on empirical libraries. Preliminary work has revealed that the BPASS continuum model may be underpredicting the amount of NaD coming from stellar atmospheres (K. Parker et al. 2021, in preparation). On the other hand, the stellar models from MIUSCAT and the MaStar models, drawing on empirical libraries, likely contain some interstellar Na I absorption from the Milky Way, overpredicting NaD absorption from stars. Thus, the three models are different and are shown to explore a broad range of parameter space in stellar continuum modeling and their effect in our analyses.

Zoom In Zoom Out Reset image size

Although the EW value of the continuum-normalized NaD spectra changes owing to the different stellar models, the velocity difference is negligible (<5 km s⁻¹), as shown by the example in Figure 5. Among the three libraries we test, we use the MIUSCAT SSPs in what follows because this particular set of models yields the most conservative (lowest) value of inferred NaD ISM contribution.

Once the spectrum is continuum-normalized, the strength of the NaD absorption arising from the ISM is quantified by the EW. The EW is measured as follows:

$$\begin{eqnarray}&&W={h}_{\lambda }\displaystyle \sum _{j=1}^{M}\left[1-\frac{{F}_{j}}{{F}_{{\rm{c}}}}\right]\end{eqnarray} \tag{ 1 }$$

where h_λ denotes the dispersion in wavelength in angstroms per pixel; M is the number of individual pixels that together constitute the desired wavelength interval, which includes the NaD absorption feature and 20 Å immediately blueward and redward of it; and F _j /F_c is the continuum-normalized flux for each of those pixels.

After continuum-normalizing the spectra as discussed in Section 3.1, we want to determine the kinematics of the neutral gas present in our galaxy samples. We perform a "detection" test for every spaxel to determine the significance of the ISM component in the NaD spectra with the following criteria:

• 1.  
  The depth of the absorption trough in the residual spectra must be at least ∼10% of the continuum level.
• 2.  
  The EW measured from the residual ISM component must be >0.3 Å.
• 3.  
  The S/N of each spectrum must be >10 per dispersion element.

The primary reason behind taking such restrictive criteria is to prevent biases from low-S/N spectra with extremely small inferred residuals, which would provide inaccurate kinematic estimates from erroneous model fits. As discussed later, an input spectrum with S/N < 10 gives an error of up to ∼20 km s⁻¹ in measured velocities, which can be marginally close to the observed velocity amplitudes we wish to investigate. Requiring spectra with S/N > 10 ensures that the velocity errors remain within acceptable limits.

The NaD spectral features in the spaxels that satisfy all the requirements are then modeled with an analytical function (Rupke et al. 2005), which takes the following form:

$$\begin{eqnarray}&&I(\lambda )=1-{C}_{{\rm{f}}}+{C}_{{\rm{f}}}\times {e}^{-{\tau }_{{\rm{B}}}(\lambda )-{\tau }_{{\rm{R}}}(\lambda )}\end{eqnarray} \tag{ 2 }$$

where C_f is the velocity-independent covering factor, and τ_B(λ) and τ_R(λ) are the optical depths of the Na i λ5891 and Na i λ5897 (vacuum wavelength) lines, respectively. The optical depth of the line, τ(λ), can be expressed as

$$\begin{eqnarray}&&\tau (\lambda )={\tau }_{0}\times {e}^{-{\left(\lambda -{\lambda }_{0}+{\rm{\Delta }}{\lambda }_{\mathrm{offset}}\right)}^{2}/{\left(({\lambda }_{0}+{\rm{\Delta }}{\lambda }_{\mathrm{offset}}){b}_{{\rm{D}}}/c\right)}^{2}}\end{eqnarray} \tag{ 3 }$$

where τ₀ and λ₀ are the central optical depth and central wavelength of each line component, respectively, b_D is the Doppler line width, and c is the speed of light. The wavelength offset is converted from a velocity offset, given by Δλ_offset = Δv λ₀/c. For NaD, τ_0,B/τ_0,R = 2. The optical depth parameter can be derived from the column density of sodium described as

$$\begin{eqnarray}&&N(\mathrm{Na}\,{\rm\small{i}})=\frac{{\tau }_{0}b}{1.497\times {10}^{-15}{\lambda }_{0}f}\end{eqnarray} \tag{ 4 }$$

where λ₀ and f are the rest-frame wavelength (vacuum) and oscillator strength, respectively. Throughout this study we assume λ₀ = 5897.55 Å and f = 0.318 (Morton 1991).

The absorption feature is fitted here using a single kinematic component. Fitting the neutral gas absorption with one component is certainly an oversimplification in a scenario where contributions from more than one gas cloud along the line of sight are embedded within a single line profile. However, this approach allows us to characterize the global kinematics obtained from the data. The absorption model is convolved with MaNGA's instrumental resolution (∼65 km s⁻¹) before the fit is performed.

In order to explore possible degeneracies, obtain unbiased fits, and estimate the errors, we wrap our fitting procedure in a Bayesian inference approach using the dynamical nested sampling algorithm (Higson et al. 2019) as implemented in the Python package dynesty (Speagle 2020). Nested sampling (Skilling 2006, 2004) estimates both the Bayesian evidence and posterior distributions in an iterative fashion until the convergence criteria are met. The method has the flexibility to sample from complex, multimodal distributions using adaptive sampling to maximize the accuracy and efficiency of the fitting process. The greater flexibility and the ability to compute Bayesian evidence (Z) on top of posteriors are advantages over traditional Markov Chain Monte Carlo (MCMC) methods. The default stopping criterion for dynesty is used in our analyses with no limitation on the maximum number of iterations. We have set reasonable priors in our MCMC fitting routine: the line width is constrained within ±100 km s⁻¹ of the stellar line width obtained from the continuum model and the limits on the allowed velocity are set at ±500 km s⁻¹. The covering fraction takes values between 0 and 1. No constraints are placed on the column density parameter. Figure 6 shows two examples of model fits applied to two different spatial bins of a red geyser (MaNGA ID: 1-217022). The residuals vary from 2% to 5%, indicating an acceptable model fit.

Zoom In Zoom Out Reset image size

For each input spectrum, our procedure delivers estimates of the column density, velocity, line width, and covering fraction. The best-fit velocity with respect to the systemic velocity is assigned to the kinematics of the neutral ISM at that location.

As a key aspect of this paper is a characterization of the global NaD kinematics in our red geyser sample, it is important to test the robustness of our velocity measurements to systematic biases. Such biases might be expected: we are fitting an asymmetric absorption doublet model to the spectral residual obtained after subtracting a broadened and possibly offset version (from the stellar continuum) of a similar doublet model.

We begin by generating synthetic stellar NaD absorption spectra assuming a Gaussian stellar velocity profile of the amplitude and width observed from the stellar continuum model of our example galaxy with MaNGA ID 1-217022. The mean of the Gaussian is set to the systemic velocity of the galaxy. We simulate the ISM component using the doublet absorption model of Rupke et al. (2005) with average values set for each parameter, namely, column density = 2.9 × 10¹³ cm⁻², line width (σ) = 150 km s⁻¹, and covering fraction = 0.2. These parameter values are obtained by fitting the observed NaD ISM spectra for the same galaxy with the same Rupke et al. (2005) model. Mock spectra are created for five assumed input velocity offsets (−100 km s⁻¹, −50 km s⁻¹, 0 km s⁻¹, 50 km s⁻¹, and 100 km s⁻¹) relative to the systemic velocity. We convolve the resultant profile with the instrumental resolution of MaNGA and then add normally distributed noise to the simulated spectra. We perform the above process for different S/N values ranging from 1 to 20 at each input velocity. We fit the resultant spectra and retrieve the model parameters according to the technique described in this subsection and in Section 3.1. Error estimates on the retrieved velocities are obtained by bootstrapping the simulated "data."

The differences between the retrieved and input velocities for different values of S/N are shown in Figure 7 via solid lines. Each colored line represents a different initial velocity offset value. The plot reveals systematic biases that result from our model fitting and velocity retrieval process. We see that a modest positive (redshifted) bias exists, but for S/N > 10 per dispersion element (gray shaded region), similar to the noise level in the NaD spectra used in our analyses, the bias in all cases is less than ∼5 km s⁻¹. As we show here, this systematic error is roughly an order of magnitude smaller than the average velocity offsets we detect in our red geyser samples.

Zoom In Zoom Out Reset image size

Such a bias probably arises from the combination of the inherent shapes of the subtracted stellar component and the interstellar component of the NaD feature. To check that, we produce similar mock NaD spectra, but this time we reverse the shape of the doublet feature, while keeping the same stellar absorption model. We perform a similar velocity retrieval process and the results are shown as dashed lines in Figure 7. We see that the bias has now switched to the negative side for input velocities <0, and find a much smaller positive bias compared to that in the previous exercise for positive velocities. Varying one or a combination of other parameters like column density, covering fraction, or line width has no impact on the amplitude of the bias.

In what follows, we are interested in the average behavior of the NaD ISM component across various subsamples as a check on our model fits of spatially binned spectra in individual galaxies. To compute spatially integrated spectra per galaxy, we shift the spectra from each spaxel to the rest frame of the galaxy and continuum-normalize them before co-addition to remove the stellar contribution by the method discussed in Section 3.1. After generating pure-ISM spectra from each spaxel, we construct a circular radius of 8'' in each galaxy (which is roughly the angular scale corresponding to the average effective radii of our targets) and include only those spaxels within the radius that satisfy the detection criteria. This is done to remove low-S/N spaxels as well as spaxels near the extreme edge of the IFU. We co-add the spectra and obtain weighted averages, where the weights are the S/N of each spectrum used in the stacking. The spatially integrated velocity thus obtained represents the average kinematics from the particular galaxy. Once we obtain an integrated spectrum for each galaxy, we also compute a stacked spectrum across each galaxy sample.

We begin by evaluating the fractional contribution of the ISM component to the total NaD absorption in our sample of 140 red geysers. After fitting the stellar continuum according to the method described in Section 3.1 we calculate the mean EW of the continuum-normalized NaD ISM spectra for each red geyser galaxy by averaging over those spaxels that satisfy the detection test described in Section 3.3. We then compute the ratio of the NaD EW from the ISM component to the total EW derived from the original spectra. The distribution of this ratio for the entire red geyser sample is presented in Figure 8.

Zoom In Zoom Out Reset image size

Since passive early-type galaxies are dominated by an evolved stellar population, which is the major contributor to the NaD absorption feature, the fraction of ISM component in the NaD EW is generally observed to be low. The fraction is <20% in a majority (80%) of the control galaxies. On the other hand, we find that in 50 out of 140 red geysers (i.e., 35% of the sample), the fractional ISM contribution is >40%, i.e., greater than twice the typical fraction observed in the control galaxies. Seventy-seven red geysers (i.e., 55% of the sample) have ISM fractions >30%. If we consider the spatially resolved map for each galaxy, roughly 56% of the red geyser sample show an appreciable presence of ISM gas (i.e., the spaxels satisfy the detection criteria in Section 3.3) over a projected area of ≥5 × 5 kpc². The fraction increases to 78% if we consider only the radio-detected red geyser sample. This contrasts with the mere 25% for the entire control sample and with 41% for the radio-detected control galaxies.

Having shown that detections of substantial ISM components are 2–3 times more common among red geysers, we now investigate whether this implies that the total amount of NaD-associated gas, as seen in absorption, is also greater. After removing the stellar contribution, we calculate the column-density-weighted mean NaD EWs across all galaxies in the FIRST-detected radio-RG and nonradio-RG samples and compare them to those of the radio-CS and nonradio-CS galaxies.

The comparison is shown as cumulative distribution functions in Figure 9. The distribution of the NaD EW from the radio-RG galaxies appears to be different from those of both the radio-detected and non-radio-detected control samples. The nonradio-RG sample, however, show a very similar distribution of NaD EWs to the radio-CS sample. We perform a Kolmogorov–Smirnov (K-S) test to measure the statistical significance of the differences of the radio-RG sample from the radio-CS and nonradio-CS samples and find that the null hypothesis that the said distributions are similar is rejected at a level <1% with p values = 0.008 and 7 × 10⁻⁵, respectively. Additionally, the radio red geysers have the highest mean with $\mathrm{EW}={0.9}_{-0.30}^{+0.31}$ Å, greater by about 0.45 Å than that of the radio-CS sample, which has a mean $\mathrm{EW}={0.46}_{-0.28}^{+0.29}$ Å, and greater than that of the nonradio-CS sample ($\mathrm{EW}={0.27}_{-0.26}^{+0.41}$ Å). The nonradio-RG sample also has a larger NaD EW than both of the control samples with a mean $\mathrm{EW}={0.51}_{-0.33}^{+0.35}$ Å. This indicates a greater level of NaD-absorbing material in the ISM of the red geysers generally, but particularly so for the radio-detected sample.

Zoom In Zoom Out Reset image size

Figure 10 shows the spatially resolved NaD EW maps from the ISM components in four representative red geysers (MaNGA ID: 1-217022, 1-273933, 1-575742, and 1-94168) with ISM contributions $\tfrac{{\mathrm{EW}}_{\mathrm{ISM}}}{{\mathrm{EW}}_{\mathrm{tot}}}\gt 30 \% \mbox{--}40 \%$. While all four of them exhibit an ISM NaD enhancement, the EWs in two of them even exceed 2.5 Å in certain regions of these galaxies.

Zoom In Zoom Out Reset image size

Figure 10 illustrates several common features in the spatial morphology of the ISM NaD distributions observed within the red geysers. These components typically have irregular shapes and are not smooth and symmetric, unlike what is typically expected from the stellar distribution (as in Cazzoli et al. 2014). These asymmetric spatial distributions reinforce the ISM origin of this enhanced NaD absorption.

We have shown that ISM-associated NaD absorption is more common and globally stronger in red geysers, but we can also characterize the typical density distributions of this absorption to distinguish between dense but relatively concentrated components and those that are more extended and diffuse. To quantify the spatial morphology and extent of NaD absorption in each galaxy, we convert the number of spaxels comprising the NaD enhancements in each galaxy to the on-sky projected area using the following conversions:

$$\begin{eqnarray}&&{A}_{\mathrm{flow}}={A}_{\mathrm{NaD}}=N\times {\left({S}_{{\rm{p}}}\right)}^{2}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{A}_{\mathrm{gal}}=\pi \,{R}_{{\rm{e}}}^{2}\,({\text{}}b/{\text{}}a)\end{eqnarray} \tag{ 6 }$$

where N is the number of spaxels; S_p is the pixel scale, which is 05 for MaNGA; R_e indicates the effective radius of the galaxy; and b/a gives the axis ratio. Hence the quantity A_flow/A_gal gives the fractional area covered by the NaD ISM absorption as compared to the on-sky projected area of the entire galaxy. Figure 11 shows the cumulative distributions of this fraction for radio-RG, nonradio-RG, radio-CS, and nonradio-CS galaxies. We see that when ISM NaD absorption above our threshold is detected, it covers a wider fractional area (${0.28}_{-0.07}^{+0.09}$) within the average radio-detected red geyser sample compared to its areal fractions—${0.09}_{-0.09}^{+0.15}$ and ${0.1}_{-0.09}^{+0.1}$—in the radio-detected and non-radio-detected control samples, respectively. This is consistent with Rupke et al. (2021), where the percentage of galaxy disk area showing inflows/outflows via NaD signatures is found to be within 10%–25% for a sample of Seyfert galaxies. A K-S test reveals that the distribution of fractional area from the radio-RG sample is similar to that from the nonradio-RG sample (p value = 0.161) but differs significantly from those from the radio-CS (p = 5 × 10⁻⁴) and nonradio-CS (p = 1 × 10⁻⁵) galaxies. Here, a p value < 0.01 indicates the null hypothesis of the distributions being similar is rejected at a level <1%. If we define a column density threshold of log N(Na i)/cm⁻² > 12.3, which corresponds to the average NaD column density across the entire red geyser and control samples, we find that 63% of the radio-RG sample show NaD absorption that covers at least 10% of the on-sky projected area in the galaxy. This contrasts with 42% in the nonradio-RG sample and with 27% and 21% in the radio-CS and nonradio-CS galaxies, respectively. Thus, the spatial extent of the cool gas present in the FIRST-detected radio-enhanced red geysers far exceeds that in the other samples of galaxies.

Zoom In Zoom Out Reset image size

By a number of measures, the previous sections have illustrated that red geysers, in particular radio-detected red geysers, feature more prominent levels of ISM-associated cool material as traced by NaD. We now turn to the kinematics of this material as a way to probe its physical origin and evolution. We first calculate the spatially resolved kinematics of the cool neutral ISM in the red geyser and control samples by measuring the Doppler shift of the ISM component of the NaD absorption line with respect to the systemic velocity for each spaxel in every galaxy (see Section 4.2).

We calculate the mean Doppler velocity offset within each galaxy after weighting by the column density, similar to the computation of the mean EW. A comparison of the mean velocities in the radio-detected and non-radio-detected red geysers with those in the control samples (split by radio detection) is shown in Figure 12. The radio-RG sample shows an average positive (or redshifted) velocity V = 41 ± 14 km s⁻¹, with respect to the systemic velocity. The nonradio-RG sample also shows a positive velocity, though of lesser amplitude, with a mean value of 13 ± 9 km s⁻¹. The radio-detected and non-radio-detected control galaxies, on the other hand, have average V = 8 ± 12 km s⁻¹ and 4 ± 11 km s⁻¹, respectively. Thus the radio red geyser sample shows a more significantly redshifted velocity on average. Looking at the statistics across the samples, we find that 86% of the radio-RG sample has V ≥ 0, while only 24% and 13% of the radio-CS and nonradio-CS galaxies show a redshifted velocity. The radio-RG sample shows a statistically different distribution from the radio-detected and non-radio-detected control samples from a K-S test with p values = 0.009 and 6 × 10⁻⁴, respectively. The nonradio-RG sample however shows a similar distribution to radio-CS with p > 0.01, but a significant difference from nonradio-CS (p value = 1 × 10⁻³).

Zoom In Zoom Out Reset image size

Figure 13 shows spatially resolved velocity maps of the NaD ISM absorption spectra of the four red geysers in Figure 10 (MaNGA ID: 1-217022, 1-273933, 1-575742, and 1-94168). Except for the galaxy in the lower left panel, all are radio-detected. The radio-detected red geysers show significantly redshifted spaxels, generally clustered on one side of the galaxy.

Zoom In Zoom Out Reset image size

We have found a clear, redshifted signature (positive velocity, V = 41 km s⁻¹) in the average of the spatially resolved NaD offset velocity in the radio-detected red geyser galaxies. However, to estimate these average velocities, we perform Bayesian fitting assuming a specific absorption model to each spaxel (Section 3.3) and compute the mean velocities after weighting by column density. We have shown in Figure 7 that a small but positive bias (≤5 km s⁻¹ at S/N > 10 per dispersion element) can emerge from the fitting algorithm used for the kinematic analyses. Thus to ensure that the velocities we get do not arise from artifacts of model fitting or averaging, we compute the kinematics of stacked NaD profiles from all four samples of galaxies to check if we will see similar trends.

We only use spaxels that meet the criteria described in Section 3.3. Additionally we impose another constraint that the considered spaxels should have a measured EW of the total NaD feature (including both stellar and ISM contributions) of >2 Å. This picks out spaxels with an appreciable NaD signal in each galaxy and provides confidence in the measured EW from the spectra. We then perform a spaxel-wise stacking of the continuum-normalized spectra around the NaD feature and over each sample. The stacked spectra for all four samples are shown in Figure 14. The red geyser samples feature a positive velocity offset. The radio-detected red geysers show a higher redshifted velocity (+53 ± 6 km s⁻¹) compared to the non-radio-detected geysers (31 ± 5 km s⁻¹). The control samples do not show any significant redshift signature relative to the systemic velocity and show values within ±10 km s⁻¹ with large errors (see Table 1, Column 5). The errors quoted in Table 1 indicate the mean 1σ uncertainty over the whole sample. The stacked velocity values agree well with the average resolved kinematic measurements obtained in Section 4.4.

Zoom In Zoom Out Reset image size

Table 1. Summary of the NaD ISM Properties Obtained for the Four Galaxy Samples

Sample	Mean EW (Å)	Fractional Area On-sky	Mean Resolved V (km s−1)	Stacked V (km s−1)
Radio red geysers	${0.90}_{-0.30}^{+0.31}$	${0.28}_{-0.07}^{+0.09}$	40.86 ± 14.07	52.8 ± 6.27
Non-radio red geysers	${0.51}_{-0.33}^{+0.35}$	${0.22}_{-0.13}^{+0.15}$	13.56 ± 9.21	31.0 ± 5.02
Radio control sample	${0.46}_{-0.28}^{+0.29}$	${0.09}_{-0.09}^{+0.15}$	8.13 ± 12.10	10.2 ± 7.78
Non-radio control sample	${0.27}_{-0.26}^{0.41}$	${0.1}_{-0.09}^{+0.1}$	4.00 ± 11.97	−28.1 ± 25.3

Download table as:  ASCIITypeset image

In order to get a rough estimate of the implied mass inflow rate in red geysers, we adopt the following equation from Roberts-Borsani & Saintonge (2019):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}={\rm{\Omega }}\,\mu \,{m}_{{\rm{H}}}\,N({\rm{H}})\,v\,r\end{eqnarray} \tag{ 7 }$$

where Ω is the solid angle subtended by the wind at its origin, m_H is the mean atomic weight (with a μ = 1.4 correction for the relative He abundance), N(H) is the column density of hydrogen along the line of sight, v is the Doppler velocity, and r is the extent of NaD absorption in the galaxy with respect to the center.

Making the same assumptions as Rupke et al. (2005), the column density of hydrogen can be written as

$$\begin{eqnarray}&&N({\rm{H}})=\displaystyle \frac{N(\mathrm{Na}\ {\rm{I}})}{\chi (\mathrm{Na}\ {\rm{I}})\,d(\mathrm{Na}\ {\rm{I}})\,Z(\mathrm{Na}\ {\rm{I}})}\end{eqnarray} \tag{ 8 }$$

where N(Na i) is the sodium column density, χ(Na i) = N(Na i)/N(Na) is the assumed ionization fraction, d(Na i) is the fraction depletion onto dust, and Z(Na i) is the Na abundance. We assume a 90% ionization fraction (χ(Na i) = 0.1), a Galactic value for the depletion onto dust (log d(Na i) = −0.95), and solar metallicity (Z_solar(Na i) = log[N(Na)/N(H)] = −5.69) similar to that assumed in Rupke et al. (2005). Putting the log N(Na i) [cm⁻²] values (which range from 11.6 to 13.7) obtained from modeling the NaD absorption feature in the radio-detected red geyser sample, we obtain an average log N(H)/cm⁻² ∼ 19.50–21.30. This is roughly in agreement with the log N(H)/cm⁻² ∼ 21 estimated according to Bohlin et al. (1978) from the average dust extinction in similar galaxies.

With these assumptions and setting Ωvr in Equation (7) equal to C_f A_NaD/t_acc, where A_NaD is the on-sky projected area of the galaxy comprising NaD absorption from Section 4.3 and t_acc is the accretion timescale $\sim \tfrac{{R}_{{\rm{e}}}}{v}$, we get

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}=\mu \,{m}_{{\rm{H}}}\,{C}_{{\rm{f}}}\,N({\rm{H}})\,{A}_{\mathrm{NaD}}\,\frac{v}{{R}_{{\rm{e}}}}.\end{eqnarray} \tag{ 9 }$$

The total mass of the gas can be estimated by the assumption ${M}_{\mathrm{tot}}={\dot{M}}_{\mathrm{out}}\times {t}_{\mathrm{acc}}$:

$$\begin{eqnarray}&&{M}_{\mathrm{tot}}=\mu \,{m}_{{\rm{H}}}\,{C}_{{\rm{f}}}\,N({\rm{H}})\,{A}_{\mathrm{NaD}}.\end{eqnarray} \tag{ 10 }$$

The cloud covering factor C_f is taken to be 0.4, the average value obtained from our model fits. For the radio-detected red geysers, the NaD velocity (v) varies from 10 km s⁻¹ to 100 km s⁻¹ (Figure 12) with a mean of ∼45 km s⁻¹ (Table 1, Columns 4 and 5). The estimated total mass ranges from M_cool ∼ 0.1 × 10⁷ M_⊙ to 5 × 10⁸ M_⊙ within the radio-detected red geyser sample while the mass inflow rate spans a large range ∼0.02–5 M_⊙ yr⁻¹. The estimated hydrogen column density, the mass flow rates, and the mass of the gas for each galaxy in the radio-RG sample are given in Table 2 along with their 1σ uncertainties. The 1σ errors in the column density are derived from MCMC fitting of the absorption line model (Section 3.3) to the observed NaD spectrum in each spaxel satisfying the detection criteria within each galaxy and adding them in quadrature. The corresponding errors in flow rate and gas mass are obtained from simple error propagation.

A few galaxies show negative velocity (or blueshift) indicating possible outflowing gas and the outflow rate, also estimated from Equation (9), is given with a negative sign in the table. As expected, the majority of radio-detected red geysers show inflowing gas. Those galaxies with no NaD ISM component detected in any spaxels (following Section 3.3) are tabulated as having zero inflow rate and M_tot although there may be cool gas present below our detection level.

Infalling gas in red geysers at the estimated accretion rates can come from a variety of sources:

• 1.  
  Accretion from outside the galaxy and its halo. Galaxy mergers, specifically minor mergers, present an obvious mechanism (Sanders et al. 1988; Kaviraj 2014; Pace & Salim 2014; Weston et al. 2017; Martin et al. 2018).
• 2.  
  A completely internal source, such as gas injected into the ISM during stellar mass loss (Conroy et al. 2015) that undergoes cooling and condensation, perhaps eventually feeding an AGN (Ciotti & Ostriker 1997, 2001, 2007).
• 3.  
  A fountain scenario, where warm ionized outflowing gas, driven out by the central AGN from earlier geyser episodes, rises into the halo, cools and condenses as it expands, and then falls back onto the galaxy (e.g., Shapiro & Field 1976; Bregman 1980; Norman & Ikeuchi 1989; Tremblay et al. 2016, 2018; Voit et al. 2017).

Numerous observations have measured the merger rates in nearby galaxies (see Ellison et al. 2015 and references therein) and several authors have emphasized widespread merger activity when observational evidence sensitive to long timescales and low mass ratios is scrutinized (e.g., van Dokkum 2005; Sancisi et al. 2008). For galaxies with stellar masses between 10^10.5 and 10¹² M_⊙, the average predicted merger rate per galaxy (assuming typical mass ratios of minor mergers ∼1:10) from semiempirical models and hydrodynamical simulations is roughly 0.1 Gyr⁻¹ (Hopkins et al. 2010; O'Leary et al. 2021). Thus galaxy mergers with a 1:10 (or greater) mass ratio can generate an M_* accretion rate of $\tfrac{0.1\times 0.1\times 1.\times {10}^{11}\,{M}_{\odot }}{{10}^{9}\,\mathrm{yr}}=1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ on average. If 10% of this value lies in cold gas, the estimated accretion rate is roughly 0.1 M_⊙ yr⁻¹. This is similar to the accretion rates quoted in Capelo et al. (2015), who performed detailed simulations by exploring a wide merger parameter space including different initial mass ratios, orbital configurations, and gas fractions with high spatial and temporal resolution. The average cool gas inflow rate we estimate from NaD in red geysers is similar to the accretion rate implied by minor mergers. Merging gas-rich dwarfs may therefore be important sources of infalling gas. The apparently higher detection rate of NaD-associated inflows in red geysers compared to that of control sample galaxies may be telling us that red geysers are a signpost of accretion activity. Indeed, some red geysers in our sample with accretion rates ∼ 5 M_⊙ yr⁻¹ may be signaling a peak in the instantaneous accretion rate during a merger event.

In addition to mergers, Ciotti et al. (1991) and Ciotti & Ostriker (2007) predicted that a significant amount of gas, ∼20%–30% of the initial stellar mass of the galaxy, is donated throughout the galaxy over the course of stellar evolution and supernova explosions. This gas can be accreted onto an AGN. The radiative cooling time of this warm gas in typical elliptical galaxies is roughly 10^7.5−8 yr. Ciotti & Ostriker (2007) quoted a mass return rate of ∼5–10 M_⊙ yr⁻¹ from the stellar populations over the cooling timescale, which exceeds the cool gas inflow rates we observe in red geysers. Moreover, some fraction of warm ionized gas (M ∼ 10⁵⁻⁶ M_⊙) ejected by red geyser winds may accumulate within the ISM or the circumgalactic medium of the galaxies over time. After cooling and condensation, this material might also become a source for the observed inflows.

Considering the different sources of infalling gas discussed above, we have found from visual inspection that most of the radio-detected red geysers show a spatially irregular and clumpy NaD distribution (see Figure 10). Additionally, we see from Figure 11 that the average fractional area occupied by NaD absorption is around 30% of the projected area of the entire galaxy for the radio red geyser sample, which implies a concentrated, lumpy distribution of gas confined in a small region. The recycled gas from the fountain scenario might be expected to present a more diffuse and smooth component than is observed here.

Cooling ISM gas originating from stellar winds and supernovae might also be expected to cover a larger area of the whole galaxy. Moreover, cool gas originating from stars might show greater correlations with the stellar kinematics, which are not observed here.

It is also interesting to consider the relatively short timescales of dust-enshrouded cool gas. As discussed in Rowlands et al. (2012), dust in early-type galaxies gets destroyed on relatively short timescales, ranging from 50 to 500 Myr, due to thermal sputtering from the hot X-ray halo surrounding the ISM or from shocks generated from Type Ia supernovae. The signatures of clumpy, dust-associated NaD gas are likely short-lived, suggesting stochastic inflow episodes from merging systems over a more continuous "rain" of recycled material.

Finally, we note that the detected cool gas in red geysers, with masses ranging from 10⁷ to 10⁸ M_⊙ and average accretion rates ∼ 0.1−5 M_⊙ yr⁻¹, might be expected to trigger an SFR of at least 1 M_⊙ yr⁻¹. However, the observed average SFR (∼10⁻² M_⊙ yr⁻¹) is several orders of magnitude lower, as evident from Figure 1. This again is suggestive of active feedback mechanisms, perhaps related to the red geyser winds themselves, that are effective in suppressing star formation in these galaxies.