On the Existence, Rareness, and Uniqueness of Quenched H i-rich Galaxies in the Local Universe

Using data from ALFALFA, xGASS, H i-MaNGA, and the Sloan Digital Sky Survey (SDSS), we identify a sample of 47 “red but H i-rich” (RR) galaxies with near-UV (NUV) − r > 5 and unusually high H i-to-stellar mass ratios. We compare the optical properties and local environments between the RR galaxies and a control sample of “red and H i-normal” (RN) galaxies that are matched in stellar mass and color. The two samples are similar in the optical properties typical of massive red (quenched) galaxies in the local Universe. The RR sample tends to be associated with slightly lower-density environments and has lower clustering amplitudes and smaller neighbor counts at scales from several hundred kiloparsecs to a few megaparsecs. The results are consistent with the RR galaxies being preferentially located at the center of low-mass halos, with a median halo mass ∼1012 h −1 M ⊙ compared to ∼1012.5 h −1 M ⊙ for the RN sample. This result is confirmed by the SDSS group catalog, which reveals a central fraction of 89% for the RR sample, compared to ∼60% for the RN sample. If assumed to follow the H i size–mass relation of normal galaxies, the RR galaxies have an average H i-to-optical radius ratio of R HI/R 90 ∼ 4, four times the average ratio for the RN sample. We compare our RR sample with similar samples in previous studies, and quantify the population of RR galaxies using the SDSS complete sample. We conclude that the RR galaxies form a unique but rare population, accounting for only a small fraction of the massive quiescent galaxy population. We discuss the formation scenarios of the RR galaxies.


Introduction
In the current theory of structure formation, galaxies form at the center of dark matter halos through gas cooling and condensation (White & Rees 1978;Mo et al. 2010).As the dominant component of cold gas, the atomic hydrogen (H I) reservoir is expected to play a vital role in regulating the rise and fall of star formation in galaxies.In fact, the H I surface mass density in low-redshift galaxies has been found to tightly correlate with the surface density of their star formation rate (SFR), known as the Kennicutt-Schmidt law (Kennicutt 1998).In addition, large surveys of multiband imaging and spectroscopy as well as H I emission at 21 cm have revealed strong correlations between the H I-to-stellar mass ratio (M H I /M * ) and the optical and UV properties of galaxies, with lower H I mass fractions in galaxies with redder colors, lower SFRs, and higher stellar mass densities (e.g., Haynes & Giovanelli 1984;Kannappan 2004;Zhang et al. 2009;Catinella et al. 2010;Huang et al. 2012;Li et al. 2012;Brown et al. 2015;Catinella et al. 2018;Li et al. 2022;Liu et al. 2023).In particular, it has been widely accepted that the cessation of star formation in a galaxy must be associated with the reduction of its H I reservoir.In a recent review, Saintonge & Catinella (2022) concluded that there is no evidence for a significant population of passive galaxies with H I reservoirs comparable to those of star-forming galaxies, based on extensive analyses of the extended GALEX Arecibo Sloan Digital Sky Survey (SDSS) survey (xGASS; Catinella et al. 2010Catinella et al. , 2018)).
On the other hand, given the large scatter of the H I scaling relations found in previous studies, typically 0.2-0.3dex in log 10 (M H I /M * ), one can expect that some galaxies will deviate significantly from the average relations.Studies of galaxies with H I anomalies should in principle provide insights into the physical processes that drive gas-related star formation and quenching, and probably also cosmological models (e.g., Peebles 2022).In this regard, there has been a rich history of studies on the H I content of early-type galaxies (ETGs; e.g., Knapp et al. 1985;Wardle & Knapp 1986;Bregman et al. 1992;Morganti et al. 2006;di Serego Alighieri et al. 2007;Grossi et al. 2009;Serra et al. 2012;Young et al. 2014).The ATLAS 3D H I survey (Serra et al. 2012) obtained resolved H I observations for a volume-limited sample of 166 nearby ETGs selected from ATLAS 3D (Cappellari et al. 2011), detecting significant H I in ∼10% of all of the ETGs inside the Virgo cluster and ∼40% of those outside.In addition, the H I disks of these galaxies are found to have much larger sizes than their optical parts, as well as much lower H I column densities than what was found typically in spiral galaxies (Serra et al. 2012(Serra et al. , 2014)).Although with large H I detection rates, ETGs usually have low H I mass fractions as expected from their relatively low SFRs.Some other attempts have been made to search for H I content in E+A, or post-starburst galaxies (e.g., Chang et al. 2001;Buyle et al. 2006Buyle et al. , 2008;;Zwaan et al. 2013;Pracy et al. 2014), a rare population of galaxies that are believed to have their star formation recently shut down, with little to no Hα Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
emission but strong high-order Balmer absorption indicative of young stars.Zwaan et al. (2013) detected H I emission in a high fraction (6/11) of E+A galaxies selected from the 2dF Galaxy Redshift Survey (Colless et al. 2001) and SDSS (York et al. 2000).A follow-up survey by Pracy et al. (2014) obtained integral field spectroscopy (IFS) for two E+A galaxies with H I emission, finding them to be consistent with a post-starburst population only in the central region, with the outer regions of both galaxies presenting strong optical emission lines.It is thus not surprising that a large fraction of the E+A galaxies in Zwaan et al. (2013) has H I gas extending to large radii.
There have also been many observational studies to search "H I-excess" galaxies (e.g., Wang et al. 2013;Huang et al. 2014;Lee et al. 2014;Lemonias et al. 2014;Wang et al. 2015;Geréb et al. 2016;Lutz et al. 2017;Geréb et al. 2018;Lutz et al. 2018Lutz et al. , 2020;;Randriamampandry et al. 2021;Zasov & Zaitseva 2022), which are usually strongly star-forming galaxies with late-type morphologies, identified as extreme outliers to the scaling relation between H I and stellar content of galaxies.Compared to galaxies with normal H I contents, these galaxies are found to have much larger H I disks and lower H I surface densities, but similar optical sizes and radial distributions of the gas-phase metallicity.Geréb et al. (2018) obtained H I resolved maps and optical spectra for four H I-excess galaxies, and found that the H I disk and stellar bulge are counter-rotating in one of the galaxies, suggesting an external origin of the excessive H I. However, Lutz et al. (2020) found no indication for counter-rotation in any of their H I-excess galaxies.By examining the star formation activity and H I kinematics, Lutz et al. (2018) suggested that these galaxies stabilize their large H I disks by their higher-than-average baryonic specific angular momentum, likely inherited from the high spin of their host dark halos.
Several recent studies have focused on H I-rich galaxies with suppressed SFRs (Lemonias et al. 2014;Parkash et al. 2019;Guo et al. 2020;Wang et al. 2022a;Sharma et al. 2023).Lemonias et al. (2014) selected a parent sample of 258 H I-rich galaxies from the Arecibo Legacy Fast Arecibo L-band Feed Array Survey (ALFALFA; Giovanelli et al. 2005;Haynes et al. 2018), defined as the top 5% fraction in the distribution of the H I mass-to-stellar mass ratio.From this parent sample, 20 massive and low specific star formation rate (sSFR) galaxies with log 10 (M * /M e ) > 10.6, <log sSFR 10.75 10 ( ) and b/a > 0.5 were selected for followup observations with the Very Large Array.Guo et al. (2020) analyzed a sample of 279 massive red spirals with M * >10 10.5 M e and red colors defined by u − r, and found 74/166 in the ALFALFA catalog to have H I detections.Wang et al. (2022a) obtained follow-up H I observations with the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan 2006;Jiang et al. 2020) for the remaining 113 massive red spirals from Guo et al. (2020), detecting H I in 75 of them.Parkash et al. (2019) and Sharma et al. (2023) selected H I-detected galaxies from H I Parkes All-Sky Survey Catalog (Meyer et al. 2004) and the H I-MaNGA program (Stark et al. 2021), respectively, both using a mid-infrared color selection of W2 − W3 < 2.0 to select galaxies with sSFR <10 −10.4 yr −1 .Sharma et al. (2023) further applied a selection of M H I > 10 9.3 M e .
The high detection rates of H I emission in galaxies with suppressed SFRs appear to challenge the well-established correlation of low star formation with low H I gas mass (see the review by Saintonge & Catinella 2022, and references therein).
Based on an analysis of the data from SDSS, xGASS, and ALFALFA, Zhang et al. (2019) claimed that nearly all massive quiescent disk galaxies in the local Universe have a large H I reservoir, similar to star-forming galaxies of similar masses.This claim was attributed by Cortese et al. (2020) to the use of aperture-corrected SFR estimates from the MPA/JHU SDSS catalog (Brinchmann et al. 2004), which do not provide a fair representation of the global SFR of H I-rich galaxies with extended star-forming disks.In a recent study, Li et al. (2023) made careful image analysis and obtained reliable SFR estimates for a sample of nearby galaxies, and they confirmed that the MPA/JHU SDSS catalog indeed significantly overestimates the SFRs for low-SFR galaxies.For massive red spirals, Zhou et al. (2021) and Wang et al. (2022a) pointed out that many of the galaxies selected by optical color u − r in Guo et al. (2020) are actually green or even blue in near-UV (NUV) − r, the NUV-to-optical color index.
In this work, we attempt to identify a sample of "red but H I-rich" (RR) galaxies, in a way that differs from previous studies.First, we make use of three existing H I surveys/ catalogs to maximize our sample size.These include xGASS, ALFALFA, and H I-MaNGA. Second, we use NUV − r as a reliable indicator of the global star formation status of a galaxy, and we use the commonly adopted color criterion of NUV − r > 5 to select truly red, thus fully quenched galaxies.As we will show, this color index works better than both the optical color u − r and the mid-infrared color W2 − W3.Third, we use a complete (volume-limited) sample of galaxies selected from the SDSS as a reference, and we identify our RR galaxies as outliers located outside the 95% contour of the SDSS galaxy distribution in the H I mass fraction versus NUV − r space.For each of the SDSS galaxies, we have estimated an H I mass based on Bayesian inferences of a number of optical properties, a technique developed in our previous work (Li et al. 2022).With the help of the SDSS sample, we are able to quantify the existence and rareness of the RR galaxies as a population.Finally, we use SFRs from both the MPA/JHU SDSS catalog and the GSWLC (Salim et al. 2018).The latter applied SED fitting to photometry from UV to IR, expected to provide more reliable estimates particularly for low-SFR galaxies (e.g., Li et al. 2023).As we will show, the RR sample selected this way is unique compared to all of the previous similar samples.We will compare the optical properties and local environment of the RR galaxies with control samples of H I-normal galaxies of similar mass and color.
The paper is organized as follows.In Section 2 we introduce the data used in this analysis.We present our results in Section 3 and discussion in Section 4. We summarize our main results in Section 5. Throughout this paper, we assume a Λ cold dark matter cosmology with Ω M = 0.3, Ω Λ = 0.7, and h = 0.7.We denote log 10 as log for simplicity.

HI Galaxy Surveys
We use H I observations of nearby galaxies from three surveys: ALFALFA, xGASS, and the H I-MaNGA program.
ALFALFA is a blind survey of H I 21 cm emission in galaxies over 7000 deg 2 of the sky and up to redshift z ∼ 0.06.Two sky areas are covered by ALFALFA: one in the northern Galactic hemisphere and one in the southern Galactic hemisphere.The final ALFALFA catalog, the α.100 catalog, contains ∼31,500 extragalactic H I line sources out to z ∼ 0.06 (Haynes et al. 2018).For each source, a category code is assigned according to data quality, with "code 1" for sources of the highest data quality and "code 2" for sources with lower data quality but a high likelihood to be true.There are 25,434 sources of "code 1" and 6068 sources of "code 2." In this work we use both categories.
The xGASS was also accomplished with the Arecibo telescope for targeted galaxies that are randomly selected from a parent sample of galaxies located in the overlapping region among SDSS data release 6 (Adelman- McCarthy et al. 2008), GALEX Medium Imaging Survey (Martin et al. 2005), and the ALFALFA survey footprint.The whole program was separated into two successive phases resulting in two complementary surveys: GASS and xGASS.The GASS obtained H I observations for a sample of galaxies with redshift 0.025 < z < 0.05 and a flat stellar mass distribution in the range 10 10 M e < M * <10 11.5 M e , down to a uniform limit of ∼1.5% in H I-tostellar mass ratio.The xGASS extends the stellar mass limit of the survey down to M * ∼10 9 M e , by further observing a sample of galaxies with 10 9 M e <M * <10 10.2 M e in the redshift range 0.01 < z < 0.02.In this work we use the xGASS representative sample constructed by Catinella et al. (2018), including 1179 galaxies from GASS and xGASS.
The H I-MaNGA is an H I follow-up program for the SDSS-IV MaNGA survey (Bundy et al. 2014).The MaNGA galaxies with cz < 15,000 km s −1 and located outside the ALFALFA footprint are observed with the Green Bank telescope, down to an rms noise comparable to the ALFALFA survey (around 1.5 mJy at a velocity resolution of 10 km s −1 ).The final H I-MaNGA catalog contains 6358 MaNGA galaxies.

The SDSS Galaxy Sample
To compare our H I-selected samples with the general galaxy population, we have selected a volume-limited sample of galaxies out of the New York University Value-Added Galaxy Catalog (NYU-VAGC), constructed by Blanton et al. (2005) from the SDSS spectroscopic galaxy sample.This sample includes 69,296 galaxies with stellar masses M * >10 9 M e and spectroscopically measured redshifts 0.02 < z < 0.05.

Galaxy Properties
For each galaxy in the H I surveys and the SDSS sample, we have collected or measured the following properties: 1. M * .Stellar mass from the NASA-Sloan Atlas 5 (NSA; Radial gradient in g − r. In addition to the optical properties listed above, we have estimated the H I gas mass for each galaxy using the estimator developed in Li et al. (2022).In short, the estimator uses a linear combination of four photometric parameters to estimate the H I-to- stellar mass ratio (M H I /M * ) of galaxies, with coefficients of the parameters constrained using Bayesian inferences of real H I data from the xGASS.Extensive tests show that the estimator can provide unbiased H I mass for optically selected samples like SDSS.The reader is referred to Li et al. (2022) for more details about the tests and applications of the estimator.The H I mass distribution of the SDSS volume-limited sample provides a predicted reference, to be compared below with samples selected from the H I surveys.

Selection of Red HI-rich galaxies
Figure 1 shows log 10 M H I /M * versus NUV − r, for the three H I surveys and the SDSS volume-limited galaxy sample.The SDSS sample is plotted as gray contours and repeated in every panel.Galaxies in different H I surveys are shown as colored symbols in different panels.In each panel, H I detections are plotted as solid dots, while upper limits are plotted as open circles.For clarity, the H I-MaNGA sample is divided into two panels, one for H I detections and one for upper limits.The anticorrelation between H I mass fraction and color, and the color bimodality can be clearly seen in the SDSS sample.All of the current H I surveys are biased for relatively H I-rich and blue galaxies due to limited survey depths.Although the smallest in sample size, the xGASS is the deepest and thus most complete among all of the surveys considered here.However, even the xGASS lacks all of the red galaxies with M H I /M * 1.5%.Interestingly, the shallower ALFALFA survey has detected H I for a (small) number of lowz galaxies with red colors (NUV − r>5) and low gas fraction because of their relatively high H I masses, covering nearly the same range of M H I /M * as predicted for the red galaxy population in the SDSS.
We select RR galaxies from the three H I surveys by requiring them to have substantially red colors, NUV − r > 5 and , and H I-MaNGA (lower-left for detections and lower-right for upper limits), after removing potentially confused galaxies and galaxies with σ NUV−r > 0.3.Solid and open symbols represent H I detections and upper limits separately.Galaxies included in the final RR and RN sample are shown as red squares and green triangles, respectively.The color is set to NUV − r = 6.5 for a few galaxies that have large uncertainties in NUV (see the text in Section 2.4 for details).The gray contours represent the SDSS volume-limited sample, of which the H I mass of each galaxy is predicted by the estimator developed by Li et al. (2022).The outermost contour includes 95% of the total sample.σ NUV−r < 0.3, as well as outstandingly high H I masses, log 10 (M H I /M * )+0.4 × (NUV − r) − 1.3 > 0. The two criteria are indicated by the vertical and upper-right dashed lines in Figure 1.The latter criterion is manually determined, and it roughly matches the 95th percentile contour of the SDSS galaxy distribution.We consider only H I detections when identifying RR galaxies, to make sure their H I mass fractions meet our requirement.For comparison, we have selected a sample of red and H I-normal (RN) galaxies with the same criteria in color and color error as adopted for RR, but requiring the H I-to-stellar mass ratio to fall below the 75th percentile contour, or specifically log 10 M H I /M * +0.4 × (NUV − r) − 0.9 < 0, as indicated by the lower dashed line plotted to the right of each panel.We consider both H I detections and upper limits when selecting RN galaxies.This is because a galaxy must fall below the lower dashed line if its H I upper limit is already below this line.There are a small number of galaxies with unusually red colors (NUV − r ?6.5) due to their extremely large values of NUV magnitude with large uncertainties.We have visually examined their NUV images and found them to be truly faint in NUV.In order to include them in our analysis, we assign them a fixed color index of NUV − r = 6.5, the reddest color that can be reached by a simple stellar population in current stellar population models (e.g., Bruzual & Charlot 2003).We apply the same criteria as above to select these galaxies into the RR or RN sample.In addition, we require all of the selected RR and RN galaxies to not be confused with neighboring galaxies located in the same beam.We identify galaxies free of confusion by requiring qcrHIconf_flag=0 for those from xGASS, and qcrconf_prob<0.1 for those from H I-MaNGA.For those from ALFALFA, we require them to have no neighboring galaxies within an angular radius of ¢ 3. 5 and a radial velocity separation of 600 km s −1 .Finally, we visually examine the DESI image of all of the selected galaxies, and exclude the edge-on galaxies with obvious dust lanes, which could be intrinsically blue but reddened by dust attenuation.
With these restrictions, we end up with 47 RR galaxies (37 from ALFALFA, four from xGASS, and six from H I-MaNGA), as well as 573 RN galaxies (43 from ALFALFA, 237 from xGASS, and 293 from H I-MaNGA).The basic properties of the RR sample are listed in Table 1.In order to take care of the mass dependence when comparing the RR and RN samples, we further trim the RN sample so that its M * distribution is the same as the RR sample.This yields a final sample of 154 RN galaxies, with 23 from ALFALFA, 80 from xGASS, and 51 from H I-MaNGA.The RR and RN galaxies are highlighted in red and green color in Figure 1, respectively.Note that the RN sample is dominated by upper limits.
The upper-left panel of Figure 2 displays the M H I /M * versus NUV − r relation for the final samples of RR and RN galaxies.In addition to these two samples, we have selected a stellar mass-controlled (SC) sample of 177 galaxies from the xGASS, which has the same M * distribution as the RR and RN samples but with no restrictions on the color and the H I mass fraction.The galaxies in the SC sample are plotted as blue diamonds.The lower-left panel plots the H I mass instead of the H I-to-stellar mass ratio as a function of NUV − r for the same samples.In the middle and right panels, we plot M H I /M * and M H I as functions of SFR, using the SFR estimates from GSWLC and the MPA/JHU catalog, respectively.
The RR and RN galaxies are roughly separated at M H I ∼ 10 9.5 M e , and both samples are dominated by galaxies with low SFRs (log 10 SFR  −1) consistent with their red colors, for SFRs given by both GSWLC and MPA/JHU.When the GSWLC SFRs are used and at the low-SFR end (log 10 SFR  −1), the RN and SC samples appear to follow roughly the overall trend of the SDSS sample, where both the H I-to-stellar mass ratio and the total H I mass depend only weakly on the SFR, roughly at log 10 (M H I /M * ) ∼ −1.8 and log 10 (M H I /M e )∼8.7, respectively.At higher SFRs, the SC sample follows the M H I -SFR relation of the SDSS sample, but is biased toward relatively low log 10 (M H I /M * ).This should be produced by the relatively high stellar masses of the SC sample at which the average SFR is reduced by galaxies from the quiescent sequence (see Figure 3 below).When the MPA/JHU SFRs are used and at the low-SFR end, the SDSS sample presents a clumplike distribution in the M H I /M * versus SFR relation and a positive correlation between M H I and SFR, with a more limited SFR range and larger scatter in M H I /M * and M H I than when the GSWLC SFRs are adopted.The distributions of the RR, RN, and SC samples relative to the SDSS sample remain similar, however, regardless of the adopted SFR estimate.
One may worry that the RN sample is not substantially representative of the general population of red galaxies, as it is selected from the H I surveys but not directly from the parent SDSS sample.To address this possibility, we use the SDSS sample to construct another control sample (referred to as "CSC"), which consists of 423 galaxies and is matched with the RR sample in both M * and NUV − r.The distributions of the CSC sample are shown as filled gray contours in Figure 2. As expected, the CSC sample well represents the underlying red population of the SDSS sample.
The majority of the RR galaxies are located beyond the 95% contour of the SDSS sample, suggesting that red galaxies with large amounts of H I gas are rather rare.In other words, red/ quenched galaxies are mostly H I-poor, with only a small fraction having unusually large amounts of H I gas.In the next Section we will compare the RR and RN samples in detail, aiming to determine the origin of the high H I gas content of RR galaxies.

Optical Properties
In Figure 3 we examine the optical properties for the RR and RN samples, and compare with the SC sample and the SDSS volume-limited sample.We consider a variety of galaxy properties, each plotted in one panel as a function of log 10 M * .In each panel, the RR, RN, and SC galaxies are plotted as red squares, green triangles, and blue diamonds, respectively.The distribution of the SDSS sample is plotted as the background contours.The CSC sample is shown as the gray filled contours.A side panel is added to the right of each panel, showing the cumulative distribution of the corresponding property for the RR, RN, SC, and CSC samples, respectively.Overall, the well-known galaxy bimodality can be clearly seen from the SDSS sample in many properties, e.g., NUV − r, SFR, D 4000 , R 90 /R 50 , μ * , and the T-type.In all panels, the RN sample shows similar distributions to the CSC sample, demonstrating that the RN sample is a representative subsample of the red galaxy population of their mass.The RR galaxies are limited to relatively high masses, with M * 10 10 M e , a result that can be understood from their red colors.As can be seen from the color-mass diagram (the top-left panel), the general population of galaxies at NUV − r > 5 from the SDSS sample is predominantly massive.
We find that the RR and RN galaxies have similar distributions in all of the properties considered.Both types of galaxies fall in the sequences of red (by selection) and quiescent galaxies with relatively low SFRs and old stellar populations as indicated by their high D 4000 .The two types of galaxies are also similar in optical size (R 50 ) and structural properties (R 90 /R 50 and μ * ) at any given mass, with smaller sizes, more centrally concentrated light distributions, and higher stellar mass densities in comparison to the general galaxy population.In addition, they exhibit similarly high σ * , high B/T, and low T-type, all indicative of early-type morphologies.
In Figure 4 we further examine the merger probability (P merger ), the probabilities of hosting a galactic bar (P bar,GZ2 and P bar,N10 ), the (g − r) color indices measured at R 50 and 0.5R 50 , and the radial gradient of the color indices.For each property, we show both the differential and the cumulative distributions for the three samples: RR, RN, and SC.The RR and RN samples show very similar distributions for all properties considered except P bar,N10 , which appears to be smaller in the RR sample on average.The different distributions of P bar,GZ2 and P bar,N10 may be produced by the different definitions of the two parameters, which can also be seen from their distributions in the SC sample.Given this difference, both the large p-value of the Kolmogorov-Smirnov (K-S) test using P bar,GZ2 and the small p-value using P bar,N10 should be taken with caution.We note, however, that nearly all of the galaxies in both samples have P bar < 0.5, indicative of no/weak bars.This is a result that is true for both P bar,GZ2 and P bar,N10 .The majority of the RR and RN galaxies have low merger and bar probabilities, as well as negative color gradients, with P merger < 0.5, P bar,GZ2 < 0.2, P bar,N10 < 0.3, and Δ g−r < 0 for ∼80% of the sample galaxies.Nearly all of the RR and RN galaxies have (g − r) > 0.7 at both 0.5R e and R e , consistent with their globally red colors.When compared to the SC sample, both the RR and RN samples are similar in P merger and Δ g−r , smaller in P bar,GZ2 and P bar,N10 , and larger in (g − r) at both R 50 and 0.5R 50 .
Results as seen from Figures 3 and 4 strongly suggest that, when compared to the RN galaxies of similar stellar masses, the RR galaxies have nothing special in any of the properties we have examined, including both the properties related to stellar populations and those related to galaxy structure, morphology, and mergers.In fact, for each property we have performed a K-S test for the RR and RN samples.The resulting p-values are indicated above the panels in both figures, which confirm that the two samples are statistically drawn from the same parent sample.

Environments
The cold gas content of galaxies may be effectively reduced by environmental effects occurring in dense regions such as tidal stripping and ram-pressure stripping.These effects are expected to be more efficient for satellite galaxies than for central galaxies.To examine the central/satellite fractions of our galaxies and their environments at different scales, we use the SDSS galaxy group catalog of Yang et al. (2007) and the local environment density inferred from the density field reconstructed by Wang et al. (2016a).Additionally, we measure the projected two-point correlations and neighbor counts for our galaxy samples using the SDSS spectroscopic and photometric samples.First, we cross match the RR, RN, and SC samples with the group catalog constructed by Yang et al. (2007).A small number of our galaxies are not included in the group catalog, including 10, 16, and 3 galaxies from the RR, RN, and SC samples, respectively.For these galaxies, we use the NSA to examine their neighboring galaxies, and we find that most of them (10/10, 14/16, 3/3 for the RR, RN, and SC samples, respectively) are either the most massive or the only galaxy within a projected radius of 1 Mpc.We thus classify these "locally dominant galaxies" as centrals.The numbers and fractions of central/satellite galaxies in the three samples are listed in Table 2.We find that the majority of the RR galaxies are centrals in their host groups, with a very high central fraction of 89%, compared to 64% and 76% in the RN and SC samples, respectively.According to the group catalog, we find the median halo mass of the RR galaxies is M h ∼ 10 12 h −1 M e , with 91% of the sample (compared to 71% and 84% for the RN and SC samples) to be hosted by relatively low-mass halos with M h < 10 13 h −1 M e .For the RN galaxies, their host halos span a wide range of halo mass, with a median halo mass of M h ∼ 10 12.5 h −1 M e and a significant fraction (18%) above M h ∼ 10 14 h −1 M e .This result clearly shows that, different from the RN galaxies, the RR galaxies are mainly central (or isolated) galaxies in relatively low-mass halos, thus free of strong environmental effects and more capable of accreting and/or retaining their cold gas.
Next, we examine the 3D overdensity of the local environment of our galaxies, d r r º -1 ¯, where ρ and r āre the local and mean matter density, respectively.Taken from the "Exploring the Local Universe with the reConstructed Initial Density Field" project (Wang et al. 2016a), the matter density field can well reproduce the distribution of both galaxies and groups of galaxies in the local Universe as observed by SDSS.Overdensities can be estimated at different scales by smoothing the density field with Gaussian kernels of different sizes.In Figure 5 we compare the distribution of the local overdensity as estimated at scales of 2 Mpc (left panel) and 4 Mpc (right panel) for different samples.At both scales, we find a trend for the RR galaxies to be located in relatively underdense regions compared to the RN sample, although the  (Li et al. 2022).In this figure, some galaxies have NUV − r = 6.5; this is not the observed value but set artificially (see Section 2.4 for details).Solid/open markers represent H I detections/nondetections(upper limits).The red squares represent RR galaxies.The green downward triangles represent RN galaxies.The blue diamonds represent SC galaxies.The colored lines denote the normalized cumulative distribution of the y-axis galaxy properties (red for RR, green for RN, blue for SC, and black dashed for CSC).The gray filled contours represent the CSC sample.The p-value of the two-sample Kolmogorov-Smirnov (K-S) test conducted between the RR and the RN sample is indicated above each side panel.
p-value of the two-sample K-S test shows that this difference is not very significant.The difference is more pronounced at 4 Mpc than at 2 Mpc, as indicated by the p-values.
We estimate two more statistics to further probe the environment of our galaxies over a wider range of spatial scales.The first statistic is the projected cross-correlation function, w p (r p ), measured for each of our samples with respect to a reference sample of ∼5.3 × 10 5 galaxies from the SDSS spectroscopic survey.This statistic quantifies the clustering of galaxies over scales from a few tens of kiloparsecs up to a few tens of megaparsecs.The other statistic is the backgroundsubtracted neighbor count, N c ( < R p ), that is, the average number of neighboring galaxies within a projected distance R p around the galaxies in our samples, as estimated in the SDSS photometric sample down to the r-band limiting magnitude of = r 21 lim mag.The contribution by uncorrelated background galaxies is statistically estimated and subtracted.Detailed descriptions of the two statistics and tests/applications can be found in Li et al. (2006bLi et al. ( , 2008aLi et al. ( , 2008b) ) and Wang & Li (2019).Figure 6 shows our measurements of w p and N c for the different samples.First of all, for the CSC sample, both the w p (r p ) and N c measurements are in good agreement with those of the RN sample, demonstrating again that the RN sample is a random subset of the red population of their mass.
On scales larger than a few megaparsecs, the RN and SC samples present similar w p (r p ) amplitudes, implying similar dark matter masses of their host halos, while the RR sample shows weaker clustering (although with large errors).This result is consistent with both the smaller halo masses of the RR galaxies as found above from the SDSS group catalog and the lower local densities found at 4 Mpc from the reconstructed density field.On intermediate scales from a few megaparsecs down to ∼100 kpc, the RR and SC samples show similar clustering amplitudes, which are significantly lower than that of the RN sample.The stronger clustering of the RN sample than the SC sample can be attributed to the different NUV − r colors of the two samples.Previous studies of galaxy clustering have shown that, at scales below a few megaparsecs, red galaxies are more strongly clustered than blue galaxies of similar masses (e.g., Li et al. 2006a).Since the RR and RN samples are matched in both stellar mass and color, the weaker clustering of the RR sample at intermediate scales can be understood from the aforementioned higher fraction of central galaxies.As shown in Li et al. (2006b) and Wang & Li (2019), the dip in w p (r p ) at intermediate scales can well be reproduced by models in which a higher-than-average fraction of the galaxies are located at the centers of their host dark matter halos.
The N c measurements present similar behaviors to w p (r p ).This is expected, as the two statistics are closely related.Compared to w p (r p ), the values of N c are measured with smaller errors on the smallest scales (100 kpc) thanks to the much deeper photometric sample.We see a trend for a dip at R p ∼ 40h −1 kpc as well as an upturn at smaller scales in the average number of close companions around the RR galaxies.This result might indicate some role of tidal interactions or mergers for RR galaxies.Large samples and more studies are needed to test this hypothesis.

HI-to-optical Size Ratio
Previous studies have revealed a tight relation between the H I disk size and H I mass for nearby galaxies (e.g., Wang et al. 2016b, and reference therein).The unusually high H I mass fractions of the RR galaxies thus imply large H I disk sizes relative to galaxies with normal H I contents, if we assume that different types of galaxies follow the same H I size-mass relation.We estimate the radius of the H I disk (R HI ) for each From top left to bottom right, the quantities in each panel are probability of being a merger or projected pair (P merger ), the probability of having bar signatures (P bar,GZ2 and P bar,N10 ), g − r color gradient, g − r at 0.5R 50 , and g − r at R 50 .The colored solid lines show the normalized cumulative distribution of the histograms with the same color.The p-value of the two-sample K-S test conducted between the RR and the RN sample is indicated above each panel.galaxy in our samples using the H I disk size-mass relation from Wang et al. (2016b), where R HI is defined as the radius at which the surface H I mass density is equal to 1 M e pc −2 .The left-hand panel of Figure 7 shows the ratio of the estimated R HI to the optical radius R 90 as a function of stellar mass for the different samples of galaxies, with the side panel showing the histogram of R HI /R 90 of all of the galaxies in each sample.As expected, the RR galaxies have the largest H I-to-optical size ratios with an average of R HI /R 90 ∼ 4, compared to R HI /R 90 ∼ 1 for the RN sample and R HI /R 90 ∼ 1.7 for the SC sample.We see weak dependence on the stellar mass for all three samples, although there is considerable scatter at fixed mass.Note that the RN sample is dominated by H I upper limits, and so the estimated H I sizes are also upper limits.This means that the H I disk of red galaxies with normal H I content is comparable to or even smaller than their optical disk.One of our RR galaxies, UGC 1382 (ID = 0 in Table 1) with resolved H I observations available for measuring H I disk size (Hagen et al. 2016), is highlighted in Figure 7 as a red square with black borders.This galaxy is above the average H I sizemass relation from Wang et al. (2016b), but with a deviation that is only slightly larger than the 3σ scatter.Although with only one galaxy, this result supports the assumption that the galaxies from different samples follow the same H I size-mass relation.This assumption is actually consistent with previous studies of resolved H I observations of nearby galaxies.For instance, using H I maps of nearby galaxies obtained from the Westerbork Synthesis Radio Telescope, Wang et al. (2013) found that their "bluedisk" galaxies with unusually high H I mass fractions lie on the same H I mass versus H I size relation as the control galaxies with normal H I content, although the former have much larger H I-to-optical size ratios, with R HI extending to values as large as ∼100 kpc for M H I ∼2 × 10 10 M e (see Figure 6 of their paper).Therefore, applying the same size-mass relation to estimate H I radii for different samples appears to be a reasonable choice for our analysis here, before large samples of resolved H I observation become available in the future.

Uniqueness of the RR Galaxies
In this Section, we compare our RR sample with similar samples in previous studies.First, we consider the sample of ETGs from the ATLAS 3D H I survey by Serra et al. (2012).Similar to the RR galaxies studied here, the ETGs in ATLAS 3D also presented much larger H I disks than their optical sizes, and the H I-richest ETGs were also found in low-density environments (Serra et al. 2012).For quantitative comparison, we have attempted to match the ATLAS 3D sample with our SDSS galaxy catalogs.Due to the relatively low redshifts of the ATLAS 3D sample, none of the galaxies are included in GSWLC, and only a fraction have counterparts in the NSA and the MPA/JHU catalog.Therefore, for the ATLAS 3D galaxies, we opt for stellar masses derived by Cappellari et al. (2013) based on dynamical modeling, i.e., M JAM in that paper.By comparing M JAM with the NSA M * for those galaxies that have counterparts in NSA, we find M JAM is systematically larger by ∼0.3 dex than M * with no obvious dependence on M JAM .To take into account this difference, we use -M log 0.3 10 JAM as the stellar mass for the ATLAS 3D .The SFRs are estimated by Kokusho et al. (2017) from polycyclic aromatic hydrocarbons (PAH) luminosities, which were obtained from SED fitting to data from AKARI, the Widefield Infrared Survey Explorer (WISE), and the Two Micron All Sky Survey (2MASS).We have compared the SFRs for those with counterparts in the MPA/JHU catalog and found no systematic differences.In Figure 8 we compare the ATLAS 3D sample with the RR sample in terms of the distributions of M * , M H I /M * , NUV − r, and sSFR.We find that, when compared to the RR sample, the ETGs from ATLAS 3D mostly have similarly red colors and low sSFRs as expected, and significantly lower H I mass fractions.Only one galaxy from the ATLAS 3D meets our selection criteria of the RR sample, while most of them fall on or even below the contours of the general red population, indicating that the majority of red ETGs are H I normal or H I poor, as generally expected.The difference between the RR sample and the general ETG population can also be seen from the bottom-right panel in Figure 3, where the RR, RN, and CSC samples show very similar distributions in T-type, all dominated by early-types, but the RR sample is very different in H I mass fraction to the other two samples.The E+A galaxies previously studied are also different from our RR galaxies.In fact, Figure 3 of Zwaan et al. (2013) has shown that their E+A galaxies as well as the several from earlier studies are mostly "green" with 3 < NUV − r < 4 and H I-normal for their color, with log 10 (M H I /M * ) < − 1.
The samples of H I-rich galaxies with suppressed or low star formation (Lemonias et al. 2014;Parkash et al. 2019;Sharma et al. 2023) and massive red spirals (Guo et al. 2020;Wang et al. 2022a) are more similar to our RR sample.The selection criteria of these previous samples are mentioned in Section 1.The comparison of these samples with the RR sample is shown in Figure 8.As can be seen from the top-left panel, all of the samples in comparison are limited to relatively high stellar mass and high H I mass fraction with respect to the complete SDSS sample, although the FAST sample from Wang et al. (2022a) includes some galaxies with low H I mass fractions.Our RR sample is distinct from other samples by having the reddest colors (top-right panel) and the lowest sSFRs (bottomleft panel), both indicating the fully quenched status of star formation in these galaxies.In contrast, both the H I-rich galaxies with suppressed star formation from Lemonias et al. (2014) and the massive red spirals from Guo et al. (2020) and Wang et al. (2022a) are mostly "green" with intermediate colors, 3  NUV − r  5, and consistently they fall mainly in the transition region between the star-forming and quiescent sequences (bottom-left panel), indicating that these galaxies are not fully quenched.As pointed out in Zhou et al. (2021), red spirals selected by optical colors such as u − r are mostly green/blue in NUV − r, as the UV emission is more sensitive than optical colors to the ongoing weak star formation in galaxies.The fact that the sample of Parkash et al. (2019) is also limited to the transition region suggests that the criterion of mid-infrared color (W2 − W3 < 2) is unable to select fully quenched galaxies as well.Sharma et al. (2023) did not tabulate their galaxies, but it is natural to expect their sample to occupy similar regions in the sSFR versus M * diagram considering the same color selection as adopted in Parkash et al. (2019).We note that massive red spirals show quite different distributions in the sSFR versus M * diagram depending on the adopted SFR estimates.When the MPA/JHU catalog is used instead of GSWLC, these galaxies move from the transition region to the quiescent sequence, showing similarly low sSFRs to our RR galaxies.We argue that this discrepancy can largely be attributed to the larger uncertainties and biases in the SFRs from the MPA/JHU catalog than from the GSWLC at intermediate-to-low SFRs (see Li et al. 2023 for a recent comprehensive investigation and discussion).
We conclude that the RR sample selected in the M H I /M * versus NUV − r diagram represents a unique population of galaxies, which are fully quenched but contain unusually large amounts of H I gas.Previous samples selected by u − r (Guo et al. 2020;Wang et al. 2022a), or W2 − W3 (Parkash et al. 2019;Sharma et al. 2023), or by a fixed cut in sSFR (Lemonias et al. 2014) represent galaxies in the "green valley," i.e., in the transition phase to being fully quenched.Their H I contents can be partly (if not fully) explained by the blue/green NUV − r colors and the intermediate sSFRs.
Considering that the cessation of star formation in massive galaxies happens from the inside out (e.g., Li et al. 2015;Wang et al. 2018), it could be that the H I-rich galaxies in previous samples and the RR galaxies in our sample have the same origin and follow a common evolutionary path, but they are at different stages.In this case, the RR sample and the previous samples should share some common properties.In fact, like our RR galaxies, the H I-rich galaxies in Lutz et al. (2018Lutz et al. ( , 2020) ) were also found to have similar optical properties (e.g., the size-mass relation and radial distributions of gas metallicity) to that of H I-normal galaxies.In addition, the large H I-to-optical size ratios found for our RR galaxies were also found in Figure 6.Measurements of projected cross-correlation functions with respect to the SDSS reference sample (upper) and background-subtracted neighbor counts in the SDSS photometric sample down to the r-band magnitude limit of = r 21 lim (lower).Results for the different samples are plotted in different symbols/lines/colors, as indicated.Lemonias et al. (2014), with the help of resolved H I observations: most of the massive H I-rich galaxies with suppressed star formation in their sample had H I radii at least twice as large as their optical radii (see their Figure 7).The authors found all of their galaxies had substantially low H I surface densities, likely to be the main cause of their low sSFRs.Given the even larger H I-to-optical radius ratios of our RR galaxies, they are expected to have even lower H I surface densities.
On the other hand, however, our RR galaxies also show quite different properties in many other aspects from galaxies used in previous investigations.For instance, Lemonias et al. (2014) found some evidence that active galactic nuclei (AGNs) or bulges also contributed to the star formation suppression in their galaxies, giving support to AGN feedback and morphology quenching (Martig et al. 2009).Sharma et al. (2023) also found some evidence for AGN feedback and bar quenching for the H I-rich but low-SFR galaxies based on the IFS data from MaNGA.In addition, Parkash et al. (2019) found a high fraction (75%) of their galaxies to be LINERs or LIERs, and they speculated that the presence of H I is a precondition for LI(N)ER emission in galaxies.Sharma et al. (2023) found a similarly high fraction (77% ± 11%) of LIERs.However, they found their control sample of low-H I and low-SFR galaxies also has a high fraction of LIERs (60% ± 10%), and concluded that H I content is not required for the high LIER fraction to correlate with the low SFR.We find no evidence for either morphology quenching, given the similar distributions in B/T and σ * for the RR, RN, and CSC samples, as seen in Figure 3, or bar quenching, given the similar distributions in P bar , as seen in Figure 4. Selected from SDSS, most of our RR galaxies do not have IFS, but we have examined their single-fiber SDSS spectra.We find that the emission lines relevant for AGN identification are weak in most galaxies, indicative of no AGNs or LIERs in their central regions.This result is consistent with the global quiescence of the RR galaxies.
In Figure 9 we compare the projected 2PCF and local densities between our RR sample and the previous samples.The samples from Lemonias et al. (2014) and Parkash et al. (2019) are too small to obtain reliable measurements, while Sharma et al. (2023) did not tabulate their galaxies.In order to compare, we select a sample of "H I-rich low-SFR" galaxies from the SDSS to mimic their samples, by applying the following two criteria: log 10 M H I > 9.3 and - < <-11.5 log sSFR 10.4 10 ( ) .We further require this sample to have the same stellar mass distribution as our RR sample.The sample includes 2200 galaxies.The measurements for this sample are plotted as black symbols in the Figure, and those of a combined sample of massive red spirals with H I detections from ALFALFA (Guo et al. 2020) and FAST (Wang et al. 2022a) are plotted as cyan symbols.As can be seen, the "H I-rich low SFR" sample and the massive red spirals are very similar in both w p (r p ) and d + log 1 10 ( ), but when compared to them, the RR galaxies are less clustered at all scales except r p ∼ 0.2 − 0.3 h −1 Mpc and have lower densities on average.We have also examined the central and satellite fractions of the two samples based on the SDSS group catalog, and found a central fraction of 73% for the massive red spiral galaxies and 70% for the "H I-rich low SFR" sample.These results imply that low-density environments and high central fractions are necessary conditions to make the RR galaxies a unique population, which is unlikely to entirely follow the same evolutionary pathway as the H I-rich low-SFR samples previously studied.

Existence and Rareness of the RR Galaxies
In a recent review, based on the xGASS sample, Saintonge & Catinella (2022) concluded that there is no evidence for a significant population of passive galaxies with H I reservoirs  1), the only galaxy in the RR sample that has resolved H I observations available for measuring H I disk size (Hagen et al. 2016).
comparable to those of star-forming main-sequence (SFMS) galaxies.Our result is not contradictory to theirs.As can been seen from Figure 1, xGASS indeed contributes only a few galaxies to our RR sample.In fact, our RR sample is contributed mainly by ALFALFA, which is much larger in sample size, though less complete in H I mass fraction when compared to xGASS.This strongly suggests that the quenched but H I-rich galaxies are rather rare, and they can be identified as a significant population for statistical studies, as done in this work only if the parent sample is large enough.The rareness of the RR galaxies can be more clearly seen when compared to the SDSS volume-limited sample.As we have seen from Figure 2, the majority of our RR galaxies are found beyond the 95% contour of the SDSS sample.In addition, as can be seen in Figure 8, all of the previously studied samples including the H I-rich but low-SFR sample (Lemonias et al. 2014), the massive red spirals (Guo et al. 2020;Wang et al. 2022a), and the ATLAS 3D sample (Serra et al. 2012) contain a few galaxies satisfying our criteria of RR galaxies.This result is also consistent with both the existence and the rareness of the RR galaxies.
Using data from both xGASS and ALFALFA, however, Zhang et al. (2019) claimed that nearly all of the massive quiescent central disk galaxies are H I-rich, with M H I similar to that of star-forming galaxies.This result was questioned by Cortese et al. (2020) who found no passive disk galaxies in xGASS with H I reservoirs comparable to those typical of starforming galaxies.The authors suggested that the previous claim by Zhang et al. (2019) was due to the use of aperturecorrected SFR estimates from the MPA/JHU catalog, in which the global star formation of H I-rich galaxies with extended star-forming disks is significantly underestimated.Here we revisit this problem, taking advantage of the SDSS volumelimited sample for which we have estimated an H I mass for each galaxy.Adopting the same sample selection as in Zhang et al. (2019), we select a sample of 1609 massive (10 10.6 M e < M * < 10 11 M e ) central disk galaxies using the SDSS group catalog (Yang et al. 2007) and the SDSS morphology catalog from Domínguez Sánchez et al. (2018).We define central galaxies as the most massive galaxy of their group and select disk galaxies by requiring T-type > 0 and P disk > 0.5, where P disk is the probability for a galaxy to be a disk provided by the morphology catalog.The average and scatter of M H I /M * and M H of this sample as functions of NUV − r and SFR are plotted as the orange lines and shaded regions in Figure 2. We observe the massive central disk galaxies to roughly follow the SC and RN samples at the highand low-SFR ends in terms of M H I /M * , thus with much lower I mass fractions than the RR galaxies of similar SFRs. of their higher stellar masses, these galaxies have higher M H I than the RN and SC samples at a given NUV − r or SFR, but still have lower M H I than the RR galaxies.This result suggests both that the massive quiescent central disk galaxies investigated in Zhang et al. (2019) are H I-poor, as expected from their red colors and low SFRs, and that these galaxies are different from the RR galaxies in our sample, which are truly quenched, have unusually high H I mass fractions, but are rare.
Our result is consistent with that of Cortese et al. (2020), where both H I detections and nondetections from xGASS were used to quantify the H I reservoirs of central galaxies as a function of stellar mass, SFR, and B/T.Following the format of their Figures 1 and 2, but using our volume-limited SDSS sample instead of xGASS, we examine the M H I -SFR relationship for central galaxies in different intervals of M * and below different B/T upper limits, plotted as black contours in Figures 10 (with SFRs from GSWLC) and 11 (with SFRs from the MPA/JHU catalog).In each panel, we indicate the best-fitting SFMS and the 2σ lower limit of the SFMS, as well as the median H I mass of the SFMS galaxies (those falling within 1σ around the SFMS).The low-SFR but H I-rich galaxies are then defined to be those with SFR below the 2σ lower limit of the SFMS and above the median M H I , as indicated by the cyan area.We calculate the fraction of the SDSS galaxies falling in this area and indicate the result in each panel.As can be seen, the fractions are small in both figures and in all panels.In the case with the GSWLC SFRs, the fraction is at most (1.24 ± 0.17)% as found for all of the centrals in the most massive bin, and decreases as one goes to galaxies with lower M * and smaller B/T.Although the fractions become larger when the MPA/JHU SFRs are used, this type of galaxy is still a minor population with a fraction of (7.11 ± 1.03)% at most.For comparison, the RR, RN, and SC galaxies are plotted in both Figures.Most of the RR galaxies fall in the cyan areas, as expected.
We conclude that, although the quenched but H I-rich galaxies like those in our RR sample do exist in the nearby Universe, they contribute only a tiny fraction to the general population of their colors and SFRs.In other words, the majority of massive quiescent galaxies are H I-poor, as found in many previous studies (see Cortese et al. 2020;Saintonge & Catinella 2022, and references therein).

Origin of the RR Galaxies
The formation scenarios of the RR galaxies are not immediately clear based on our results.The high fraction of central galaxies in the RR sample as well as the relatively low density of the local environment are both helpful for the galaxies to sustain their cold gas against environmental effects such as ram-pressure and tidal stripping, which are stronger for satellite galaxies and in denser regions.However, the origin of the large amounts of H I gas in these galaxies is still unclear.As mentioned, the large H I-to-optical size ratios of the RR galaxies imply very extended H I disks with low surface densities of H I in the outskirt.It is thus natural to expect no or weak star formation even with a substantially large amount of total H I mass.This is consistent with the finding of Lemonias et al. (2014), where the resolved H observations revealed low H I densities in a sample of 20 H I-rich galaxies with suppressed SFRs.Resolved H I observations from The ATLAS 3D H I survey also show that about 20% of the nearby ETGs outside the Virgo Cluster are surrounded by large H I disks or rings of low column density (Serra et al. 2012(Serra et al. , 2014)), with signs of recent star formation similar to those in the outskirt of spiral galaxies (Yıldız et al. 2015(Yıldız et al. , 2017)).We have visually examined the optical image of our RR galaxies from the DESI Legacy Survey (see Figure 12 in the Appendix), finding about half of our galaxies to present faint and extended ring or spiral structures in the outskirt.This fraction could be even larger if deeper images are available.In addition, some of these galaxies show misalignment between the main stellar body and the extended structure, which was also observed for half of the ATLAS 3D ETGs that have large H I disks/rings (Serra et al. 2012(Serra et al. , 2014)).This result is consistent with the finding of the ATLAS 3D survey, supporting an external origin of the massive H I gas content, such as gas accretion and/or mergers.
The extended H I disks may be supported by high angular momentum, which is inherited from the high spin of host dark matter halos, as suggested by Lutz et al. (2018) for their H I- extreme galaxies.Mancera Piña et al. (2021a) and Mancera Piña et al. (2021b) recently measured the specific angular momentum for both stars and atomic gas in a sample of nearby disk galaxies, finding a strong correlation between the gas fraction and the specific angular momentum, with higher angular momenta and lower star formation efficiencies for galaxies of higher gas fractions.In particular, it was found that the H I-extreme galaxies from Lutz et al. (2018) tightly followed the same scaling correlations of angular momentum as galaxies of normal H I content.Although our RR galaxies differ from those previous samples in certain aspects, as discussed above, they could also have high angular momentum, due to the corresponding H I gas not being able to be efficiently converted into stars.In fact, the angular momentum of the circumgalactic medium has been found to either enhance or suppress star formation in simulated disk galaxies (Wang et al. 2022b;Lu et al. 2022).The ideal of "angular momentum quenching" has also been advocated for in empirical models of galaxy evolution (e.g., Obreschkow et al. 2016;Peng & Renzini 2020), although the model predictions are found to not fully match the current H I surveys, such as xGASS (Hardwick et al. 2022).The RR galaxies studied in our work represent a unique population useful to test such models.However, for this Figure 10.Distributions of central galaxies from the RR, RN, SC, and SDSS samples on the M H I vs. SFR diagram.Panels from left to right correspond to different stellar mass ranges, as indicated above the top panels, while panels from top to bottom correspond to different thresholds of the bulge-to-total ratio, as indicated on the left-hand side.In each panel, the red, green, and blue symbols represent the RR, RN, and SC galaxies, respectively, while the black contours enclose 20%, 45%, 70%, and 95% of the SDSS sample.The cyan, yellow, and black lines indicate the SFR on the best-fitting star-forming main sequence (SFMS), the 2σ lower limit of the SFMS, and the median H I mass of the main-sequence galaxies (those falling within 1σ around the SFMS), respectively.The fraction of H I-rich but low-SFR galaxies in the SDSS sample (i.e., those falling in the cyan region) and its Poisson error are indicated in each panel.
purpose, one would need larger samples with both resolved H I observations and angular momentum measurements.
We should point out that, when estimating the H I sizes for our galaxies, we have assumed the H I gas is distributed in a regular disk, as in most H I-bearing galaxies.In fact, the only galaxy in our sample with resolved H I observations, UGC 1382 (ID = 0 in Table 1), indeed shows a regularly rotating H I disk.However, irregular and even off-galaxy H I distributions have been observed for a handful of H I-bearing ETGs, e.g., the giant H I ring around NGC 3384 and M96 in the Leo group (Schneider et al. 1983;Schneider 1985Schneider , 1989;;Michel-Dansac et al. 2010), and the large H I ring around the quiescent galaxy AGC 203001 (Bait et al. 2020).As can be seen from Figure 8, both the Leo ring and AGC 203001 meet the selection criteria of our RR sample.Therefore, it is possible that some fraction of our RR galaxies also show similarly irregular H I distributions.Although its physical origin is still under debate, such H I-dominated rings are likely expelled gas from H I-rich ETGs caused by a collision with an intruder galaxy (e.g., Michel-Dansac et al. 2010;Bait et al. 2020) or tidal tripping by the group potential (e.g., Bekki et al. 2005;Serra et al. 2013;Corbelli et al. 2021), or falling back of the tidal tail of gas-rich major mergers (e.g., Morganti et al. 2003;Sameer et al. 2022).Considering the extreme paucity of such observational cases in the literature as well as the fact that the majority of our RR galaxies have no close companions, we would expect such an irregular or off-galaxy H I distribution to be not a common feature in the RR galaxies.Resolved H I observations for the RR galaxies would be necessary in order to test out this conjecture.

Summary
In this work we have identified and investigated a rare population of RR galaxies in the local Universe.We make use of three H I surveys including ALFALFA, xGASS, and H I-MaNGA, from which we select a sample of 47 RR galaxies on the diagram of log 10 (M H I /M * ) versus NUV − r.Our RR galaxies required to have NUV − r >5, and be outliers in the distribution of log 10 (M H I /M * ), located outside the 95% contour of a volume-limited sample of SDSS for which we have estimated an H I mass for each galaxy based on Bayesian inferences of its optical properties.The RR galaxies are relatively massive, with stellar mass M *  10 10 M e .For comparison, we have selected a control sample of "red and H I-normal" (RN) galaxies that are closely matched with the RR sample in M * and NUV − r, as well as a "stellar mass control" (SC) sample that is matched only in M * .We have demonstrated that the RN sample is representative to the massive red galaxy population (Figures 3 and 6).
We have examined a variety of optical properties, finding the RR and RN samples to be very similar in all properties, which are typical of massive red (quenched) galaxies in the SDSS volume-limited sample.We have also examined the environment of our samples, as quantified by three different statistics: overdensity of the local environment (1 + δ) as estimated at 2 and 4 Mpc, projected two-point cross-correlation function w p (r p ) with respect to a reference sample of half a million galaxies from the SDSS spectroscopic survey, and backgroundsubtracted neighbor counts N c ( < R p ) estimated with the SDSS photometric sample down to an r-band magnitude of = r 21 lim mag.We find the RR sample to be associated with lower-density environments at a scale of 4 Mpc, and have lower clustering amplitudes and smaller neighbor counts at intermediate scales from a few hundred kiloparsecs to a few megaparsecs, when compared to the RN sample.These results imply that the RR galaxies are preferentially located at the center of relatively low-mass dark matter halos.This is confirmed by the SDSS galaxy group catalog, which reveals that about 89% of our RR galaxies are indeed the central galaxy of their groups, compared to a central fraction of ∼60% for the RN sample.The median halo mass of the RR sample is found to be M h ∼ 10 12 h −1 M e , compared to M h ∼ 10 12.5 h −1 M e for the RN sample.Furthermore, we have estimated the H I disk size for each of our galaxies based on the tight H I size-mass relation from the literature.We find an average H I-to-optical size ratio of R HI /R 90 ∼ 4 for the RR sample, which is four times the average ratio found for the RN sample.We have compared our RR sample with similar samples investigated in previous studies.We find that the previous samples are dominated by star-forming galaxies with extremely high H I content, or quenched H I-poor galaxies, or transition galaxies with blue-to-green colors and intermediate SFRs, thus falling in between the SFMS and the quiescent sequence.Our RR galaxies are unique for their (NUV-to-optical) red colors, indicative of fully quenched star formation status and simultaneously high fractions of H I mass.We use the SDSS volume-limited sample with estimated H I masses to show that, although the RR galaxies form a unique population, they are rather rare, contributing only a tiny fraction of the quenched galaxy population.The majority of massive quenched galaxies in the low-redshift Universe are H I-poor.

Figure 1 .
Figure1.Sample selection.The gold solid dots represent the parent sample of all galaxies with H I observation from xGASS (upper-left), ALFALFA (upper-right), and H I-MaNGA (lower-left for detections and lower-right for upper limits), after removing potentially confused galaxies and galaxies with σ NUV−r > 0.3.Solid and open symbols represent H I detections and upper limits separately.Galaxies included in the final RR and RN sample are shown as red squares and green triangles, respectively.The color is set to NUV − r = 6.5 for a few galaxies that have large uncertainties in NUV (see the text in Section 2.4 for details).The gray contours represent the SDSS volume-limited sample, of which the H I mass of each galaxy is predicted by the estimator developed byLi et al. (2022).The outermost contour includes 95% of the total sample.

Figure 2 .
Figure 2. Diagrams of M H I /M * (the top row) and M H I (the bottom row), vs. NUV − r (the left column) and SFRs (with SFRs from GSWLC in the middle column, and SFRs from the MPA/JHU catalog in the right column).The RR, RN, and SC galaxies are shown as red squares, green triangles, and blue diamonds, with solid and open symbols for H I detections and nondetections (upper limits) separately.The color is set to NUV − r = 6.5 for a few galaxies that have large uncertainties in NUV (see the text in Section 2.4 for details).The gray contours represent the SDSS volume-limited sample, of which the H I mass of each galaxy is predicted by the estimator developed by Li et al. (2022).The outermost contour includes 95% of the total sample.The filled gray contours represent the CSC sample.The orange solid line and the shaded region show the median and the 1σ scatter of the H I-to-stellar mass ratio for massive disk central galaxies.

Figure 3 .
Figure 3. Galaxy properties as a function of stellar mass.The gray contours represent the SDSS volume-limited sample, of which the HI mass of each galaxy is predicted by the H I estimator(Li et al. 2022).In this figure, some galaxies have NUV − r = 6.5; this is not the observed value but set artificially (see Section 2.4 for details).Solid/open markers represent H I detections/nondetections(upper limits).The red squares represent RR galaxies.The green downward triangles represent RN galaxies.The blue diamonds represent SC galaxies.The colored lines denote the normalized cumulative distribution of the y-axis galaxy properties (red for RR, green for RN, blue for SC, and black dashed for CSC).The gray filled contours represent the CSC sample.The p-value of the two-sample Kolmogorov-Smirnov (K-S) test conducted between the RR and the RN sample is indicated above each side panel.

Figure 4 .
Figure4.Probability distribution of galaxy properties.From top left to bottom right, the quantities in each panel are probability of being a merger or projected pair (P merger ), the probability of having bar signatures (P bar,GZ2 and P bar,N10 ), g − r color gradient, g − r at 0.5R 50 , and g − r at R 50 .The colored solid lines show the normalized cumulative distribution of the histograms with the same color.The p-value of the two-sample K-S test conducted between the RR and the RN sample is indicated above each panel.

Figure 5 .
Figure 5. Histograms and cumulative distributions of d + log 1 10 ( ), where δ is the local overdensity estimated at scales of 2 Mpc (left) and 4 Mpc (right) by Wang et al. (2016a).Results for the RR, RN, SC, and CSC samples are plotted with different colors as indicated.For clarity, the histogram is not shown for the CSC sample.The p-values from K-S tests are indicated by comparing the RR sample with the RN and SC samples.

Figure 7 .
Figure7.The H I-to-optical radius ratio as a function of stellar mass.The red squares represent the RR sample.The green triangles represent the RN sample.The blue diamonds represent the SC sample.Solid and open symbols represent H I detections and nondetections (upper limits) separately.The dashed lines show the median value of the H I-to-optical radius ratio for the three samples, with the shaded region indicating the uncertainty due to the scatter of the H I size-mass relation.The histograms show the normalized distribution of the H I-to-optical radius ratio.The p-value of the two-sample K-S test conducted between the RR and the RN sample is indicated on top of the left panel.The red square is UGC 1382 (ID = 0 in Table1), the only galaxy in the RR sample that has resolved H I observations available for measuring H I disk size(Hagen et al. 2016).

Figure 8 .
Figure 8. Comparisons of the RR sample studied in this work with previous H I samples in the diagrams of M H I /M * vs. M * (top-left), M H I /M * vs. NUV − r (topright), and sSFR vs. M * (with SFRs from GSWLC in the bottom-left panel, and SFRs from the MPA/JHU catalog in the bottom-right panel).Solid and open symbols represent H I detections and nondetections (upper limits) separately.Red squares are the RR galaxies.Yellow dots are the H I-rich galaxies with suppressed star formation from Lemonias et al. (2014).Purple pluses and cyan dots are massive red spirals from Guo et al. (2020) and Wang et al. (2022a).Black solid dots with error bars are the median relation and scatter of the sample from Parkash et al. (2019).The yellow dashed lines indicate the sample selection in Sharma et al. (2023): M H I >10 9.3 M e (top-left) and sSFR<10 −10.4 yr −1 (bottom panels).Black solid and open stars represent the H I detections and nondetections (upper limits) of the ATLAS 3D sample, for which stellar masses are derived from M JAM given by Cappellari et al. (2013) based on dynamical modeling, and sSFRs are estimated from PAH luminosities by Kokusho et al. (2017) based on fits to data from AKARI, WISE, and 2MASS.The blue crosses and dots indicate the Leo ring and AGC203001 for which the H I and stellar masses and the colors are taken from Bait et al. (2020).The gray contours represent the SDSS volume-limited sample, of which the H I mass of each galaxy predicted by the estimator from Li et al. (2022).

Figure 9 .
Figure 9. Projected two-point correlation function w p (r p ) (left) and cumulative distributions of local density estimated at 4 Mpc (right), as measured for the RR and RN samples studied in this work, as well as the H I-detected massive red spirals from Guo et al. (2020) and Wang et al. (2022a) and a sample selected from SDSS to mimic the H I-rich low-SFR samples studied in previous studies (see the text for details).The value of "p" in the legend of the right panel is the p-value of the twosample K-S test conducted between the RR and the corresponding sample.

Figure 11 .
Figure 11.Same as Figure 10 except that the SFR estimates are take from the MPA/JHU catalog instead of GSWLC.

Figure 12 .
Figure12.The optical images of the RR galaxies.The galaxies are ordered by decreasing H I-to-stellar mass ratio ( M M log 10 HI * ).The galaxy ID from Table1and the value of M M log 10 HI * for each galaxy are indicated on top of the corresponding image.

Table 1
Basic Properties of the RR Sample Note.Source ID, R.A., decl., redshift, stellar mass from the NSA, NUV − r color index (those assigned a fixed NUV − r are shown as 6.50), Petrosian 50% and 90% radii in the r band, SFR taken from GSWLC and MPA/JHU (those not matched to these catalogs are denoted as "L"), H I mass, and source of H I mass (0:ALFALFA, 1:xGASS, 2:H I-MaNGA).

Table 2
Central/Satellite Numbers and Fractions