HI-MaNGA: Tracing the physics of the neutral and ionized ISM with the second data release
David V. Stark, Karen L. Masters, Vladimir Avila-Reese, Rogemar Riffel, Rogerio Riffel, Nicholas Fraser Boardman, Zheng Zheng, Anne-Marie Weijmans, Sean Dillon, Catherine Fielder, Daniel Finnegan, Patricia Fofie, Julian Goddy, Emily Harrington, Zachary Pace, Wiphu Rujopakarn, Nattida Samanso, Shoaib Shamsi, Anubhav Sharma, Elizabeth Warrick, Catherine Witherspoon, Nathan Wolthuis
aa r X i v : . [ a s t r o - ph . GA ] J a n MNRAS , 1–23 (2020) Preprint 1 February 2021 Compiled using MNRAS L A TEX style file v3.0
HI-MaNGA: Tracing the physics of the neutral and ionized ISM with thesecond data release
David V. Stark, ★ Karen L. Masters, Vladimir Avila-Reese, Rogemar Riffel, , Rogerio Riffel, , Nicholas Fraser Boardman, Zheng Zheng, Anne-Marie Weijmans, Sean Dillon, Catherine Fielder, Daniel Finnegan, Patricia Fofie, Julian Goddy, Emily Harrington, Zachary Pace, Wiphu Rujopakarn, , , Nattida Samanso, Shoaib Shamsi, Anubhav Sharma, Elizabeth Warrick, Catherine Witherspoon, Nathan Wolthuis, Department of Physics and Astronomy, Haverford College, 370 Lancaster Avenue, Haverford, PA 19041, USA Instituto de Astronomía, Universidad Nacional Autónoma de México, A.P. 70-264, 04510, Mexico, D.F., México Laboratório Interinstitucional de e-Astronomia, 77 Rua General José Cristino, Rio de Janeiro, 20921-400, Brasil Departamento de Física, CCNE, Universidade Federal de Santa Maria, 97105-900, Santa Maria, RS, Brazil Universidade Federal do Rio Grande do Sul, IF, CP 15051, Porto Alegre 91501-970, RS, Brazil Department of Physics and Astronomy, University of Utah, 115 S. 1400 E., Salt Lake City, UT 84112, USA National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, China School of Physics and Astronomy, University of St Andrews, North Haugh, St. Andrews KY16 9SS, UK Department of Physics and Astronomy, University of Texas at San Antonio, One UTSA Circle, San Antonio, TX 78249, USA PITT PACC, Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA Department of Physics, Siena College, 515 Loudon Road, Loudonville, New York 12211, USA Williams College Department of Physics, 33 Lab Campus Drive, Williamstown, MA 01267, USA Department of Physics, Bryn Mawr College, 101 N Merion Ave, Bryn Mawr, Pennsylvania 19010, USA Department of Astronomy, University of Wisconsin-Madison, 475N. Charter St., Madison WI 53703, USA Department of Physics, Faculty of Science, Chulalongkorn University, 254 Phayathai Road, Pathumwan, Bangkok 10330, Thailand National Astronomical Research Institute of Thailand, Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand Kavli IPMU, Todai Institutes for Advanced Study (WPI),The University of Tokyo, Kashiwa, Japan 277- 8583
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present the second data release for the HI-MaNGA programme of H i follow-up observations for the SDSS-IV MaNGAsurvey. This release contains H i measurements from 2135 Green Bank Telescope (GBT) observations and an updated crossmatchof the SDSS-IV MaNGA sample with the ALFALFA survey. We combine these data with MaNGA spectroscopic measurementsto examine relationships between H i -to-stellar mass ratio ( M H i / M ∗ ) and average ISM/star formation properties probed by opticalemission lines. M H i / M ∗ is very weakly correlated with the equivalent width of H 𝛼 , implying a loose connection between theinstantaneous star formation rate and the HI reservoir, although the link between M H i / M ∗ and star formation strengthens whenaveraged even over only moderate timescales ( ∼
30 Myrs). Galaxies with elevated H i depletion times have enhanced [O i ]/H 𝛼 and depressed 𝐻𝛼 surface brightnesses, consistent with more H i residing in a diffuse and/or shock heated phase which is lesscapable of condensing into molecular clouds. Of all optical lines, M H i / M ∗ correlates most strongly with oxygen equivalentwidth, EW(O), which is likely a result of the existing correlation between M H i / M ∗ and gas-phase metallicity. Residuals in the M H i / M ∗ − EW(O) relation are again correlated with [O i ]/H 𝛼 and H 𝛼 surface brightness, suggesting they are also driven byvariations in the fraction of diffuse and/or shock-heated gas. We recover the strong anti-correlation between M H i / M ∗ and gas-phase metallicity seen in previous studies. We also find a positive relationship between M H i / M ∗ and [O i ]6302/H 𝛼 , independentof metallicity, which may suggest that higher fractions of diffuse and/or shock heated gas are more prevalent in gas-rich galaxies. Key words: galaxies: ISM – radio lines: galaxies – catalogs – surveys
Neutral atomic hydrogen (H i ) is typically the dominant phase ofcold gas found within galaxies at the present epoch and plays ★ E-mail: [email protected] an important role in the evolution of galaxies. Although stars arenot expected to form out of H i directly, it is still the large-scalegas reservoir out of which the dense, cold, star-forming molecularclouds condense (e.g. Elmegreen 1993). Generally, galaxies withhigher H i -to-stellar mass ratios ( M H i / M ∗ ) have higher specificstar formation rates, sSFR (Doyle & Drinkwater 2006; Wang et al. © D. V. Stark et al. M H i / M ∗ has a particularly strong rela-tionship with galaxy color, which traces the relative growth ratesof galaxies over ∼
100 Myr to Gyr timescales (Kannappan et al.2013), such that galaxies with bluer colors/higher long-term growthrate tend to have large H i reservoirs, while galaxies with red col-ors/lower long-term growth rates tend to have little or no detectableH i (Kannappan 2004; Catinella et al. 2013; Kannappan et al. 2013;Eckert et al. 2015; Catinella et al. 2018). The meaning of the H i -color relationship has been called into question by Jaskot et al. (2015)however, who argue it is primarily driven by internal dust extinction,not recent star formation.The ability of H i reservoirs to contribute to star formation isstrongly linked to the local conditions within the interstellar medium(ISM), namely those which promote the formation of molecularclouds (e.g. Krumholz et al. 2009a,b). The typical conditions of theISM likely vary between low stellar mass (typically H i -rich) galax-ies and high stellar mass (typically H i -poor) galaxies, which canin turn impact the ability of their H i reservoirs to contribute tostar formation. For instance, Dalcanton et al. (2004) find that theISM is more confined to dense thin disks in high-mass galaxies,while the ISM is more diffuse with a higher vertical scale length,in lower mass galaxies. This result further implies lower turbulentvelocities and more disk instabilities in higher mass galaxies. Con-sistently, several studies argue low-mass galaxies are generally sta-ble against fragmentation (Meurer et al. 1996; van Zee et al. 1997;Hunter et al. 1998). Furthermore, spectroscopic observations havealso shown that gas-phase metallicity and H i content are inverselyrelated (Lequeux et al. 1979; Skillman et al. 1996; Peeples et al.2008; Robertson et al. 2012; Moran et al. 2012; Hughes et al. 2013;Bothwell et al. 2013; Brown et al. 2018; Zu 2020). The combineddynamical stability and lower metal content of low mass gas richgalaxies is expected to make it harder for them to meet the condi-tions necessary for molecular cloud formation and subsequent starformation.Optical spectroscopy provides a powerful means of probing howthe conditions of the ISM vary between gas-rich and gas-poor galax-ies, not only by revealing metal abundances, but also ionization frac-tions, properties of ionizing sources, gas densities, and gas temper-atures. With the exception of metallicity, many of these propertieshave not been explored in great detail, and prior studies which relategas content to optical spectroscopic properties have been limited in anumber of ways. Traditionally, spatially resolved spectroscopy mea-sured over large fractions of galaxy disks (i.e., IFU spectroscopy)has been difficult to obtain, and such data have been limited tosmall samples. Sample sizes can be dramatically increased by usinglarge fibre-fed optical spectroscopic surveys (e.g., SDSS, York et al.2000), but these observations only capture spectra in the nuclei ofnearby galaxies. The current generation of large IFU surveys likeMaNGA (Bundy et al. 2015), SAMI (Croom et al. 2012), and CAL-IFA (Sánchez et al. 2012) address these limitations by providingspatially resolved spectroscopy out to large radii for thousands ofgalaxies.Complementing the optical IFU observations of ∼ i follow-up program collecting single-dish 21cm data for all MaNGA galaxies lacking coverage from otherH i surveys. In this paper, we present the second data release for HI-MaNGA which adds 1846 new observations (a factor of 5.5 timeslarger than the previous data release) and an updated crossmatchagainst the ALFALFA survey, yielding a sample of 3698 uniquegalaxies with global H i and optical IFU spectroscopic data. We usethis powerful data set to investigate how the average ISM properties vary as a function of gas-richness, providing additional importantclues as to what affects the ability of different types of galaxies toconvert their ISM into stars.In Section 2, we present our H i observations and updated catalog, adiscussion of how H i upper limits are incorporated into our analysis,and a description of our analysis sample and derived data products.In Section 3, we present scaling relations between M H i / M ∗ and ISMproperties derived from MaNGA optical spectroscopy. In Section 4,we interpret our findings and discuss potential systematic errors. Ourconclusions are summarized in Section 5.Throughout this work, we assume a cosmology with 𝐻 =
70 km s − Mpc − , Ω 𝑚 = . , and Ω Λ = . . HI-MaNGA (Masters et al. 2019) is an H i follow-up program for theSDSS-IV Mapping Nearby Galaxies at Apache Point Observatory(MaNGA) survey (Gunn et al. 2006; Smee et al. 2013; Bundy et al.2015; Blanton et al. 2017; Yan et al. 2016a). In this paper, we releasenew Robert C. Byrd Green Bank Telescope (GBT) observations of1846 galaxies that were taken over 1079 hours between 2017 and2018 (all under proposal code AGBT17A-012). Details about targetselection, observing strategy, and data reduction can be found in thedata release 1 (DR1) paper (Masters et al. 2019). To briefly sum-marize, all MaNGA galaxies lacking H i data are targeted, with nopre-selection based on other properties such as stellar mass, color, ormorphology, although there is a 𝑧 < . restriction which introducesa bias against the highest stellar mass galaxies due to the stellar mass-redshift relation built into the MaNGA sample definition (Wake et al.2017). All targets are observed using standard position-switching tothe same depth ( ∼ ∼
15 minutes on-source). All scans areaveraged, any strong RFI is removed, and the spectra are boxcar andhanning smoothed to a final velocity resolution of ∼
10 km s − . Apolynomial is fit to remove any leftover baseline variations, afterwhich source parameters are measured (see Section 2.4).There are few modifications to the derived data products in thisrelease compared to DR1. First, we have incorporated cosmologicalcorrections (factors of + 𝑧 ) to relevant quantities. Second, all fluxeshave been multiplied by a factor of 1.2 to correct for an underestima-tion of the default flux calibration scale built into GBTIDL , the datareduction software used to process our data (see Goddy et al. 2020 forfurther details). Third, linewidths have been corrected for instrumen-tal broadening. Finally, the catalog contains additional informationto aid users, including additional general information for galaxies,an assessment of source confusion, and flags to highlight potentiallyunreliable observations. Further details are given in Section 2.4. A significant fraction of MaNGA ( ∼ , which includes all MaNGA galaxies). We assignthe closest optical match within one half-power beam width to theH i source. Cases with multiple optical matches, including any outto 1.5 times the half-power beam width, are flagged as confused. http://gbtidl.nrao.edu/ MNRAS , 1–23 (2020)
I-MaNGA Further details on the criteria used for this matching, including howmatching is done in velocity space, are given in Section 2.1. Forany MaNGA galaxies which fall inside the ALFALFA footprint butdo not have matches within the a.100 catalog, we extract its spec-trum from the ALFALFA data cubes using a ′ × ′ window at eachgalaxy’s location. A first order baseline is subtracted, and the rmsnoise is measured enabling an estimate of an H i flux upper limit.ALFALFA is an untargeted H i survey, and as such requires a rel-atively high signal-to-noise ratio ( 𝑆 / 𝑁 ) for detections in order toconsider them reliable enough to incorporate into the final a.100 cat-alog. In contrast, a pointed survey (like our GBT observations) canaccept lower 𝑆 / 𝑁 on detections due to the prior knowledge of wherethe galaxies are located both on the sky and in redshift space. Thea.100 catalog contains code 1 sources with 𝑆 / 𝑁 > . , and code2 sources with lower 𝑆 / 𝑁 if they coincide with an existing opticalcounterpart with known redshift. Although there is no minimum 𝑆 / 𝑁 required for the code 2 sources, their frequency declines rapidly be-low 𝑆 / 𝑁 ∼ . . This 𝑆 / 𝑁 corresponds to an integrated flux approx-imately matching the “detection limit" of the survey (Haynes et al.2011). While non-detections in our GBT observations are used toderive 𝜎 upper limits, using 𝜎 upper limits for objects lackingmatches in the a.100 catalog would be inappropriate because there isthe possibility that they are weak sources with < 𝑆 / 𝑁 < . thatwere missed. Therefore, all ALFALFA non-detections are treated as . 𝜎 upper limits in order to more accurately match the completenesslevel of the a.100 catalog. The GBT and Arecibo beams are 9 ′ and 3.5 ′ , respectively, at 21cm.Thus, there is a strong potential for H i source confusion caused bymultiple galaxies within the beam. We now flag potential cases ofconfusion for all H i detections by identifying all known galaxiesfrom the NSA within 1.5 times the half-power beam width (HPBW)from the beam center and at similar velocity as our primary target.Specifically, assuming the primary galaxy (that being targeted) hasa linewidth 𝑊 and central velocity 𝑉 𝐻 𝐼 equivalent to what wasmeasured, and the secondary galaxies (those in close proximity tothe primary) have linewidths of 𝑊 =
200 km s − and central veloc-ities equivalent to their known optical redshifts, we look for overlapbetween their H i profiles in velocity space. Any cases with veloc-ity overlap are flagged as potentially confused (see the confused column in the final catalog (Section 2.4). For this analysis, we as-sess velocity overlap using 𝑊 +
20 km s − to account for the factthat profiles typically extend slightly beyond 𝑊 . This adjustmentis more akin to using a 𝑊 measurement (Kannappan et al. 2002),but we use the approximation instead of actual measurements of 𝑊 because the measurements may be unreliable for a significant frac-tion of our data set with low 𝑆 / 𝑁 . We also enforce a minimum 𝑊 of 20 km s − to account for turbulent motions.The above analysis attempts to identify all cases of confusion, butlikely over-estimates the true rate of confusion because some sec-ondary galaxies may be inherently gas-poor relative to the primarytarget and/or far enough away from the beam center that they do notcontribute significantly to the total measured H i flux. To refine ourestimate of source confusion, we perform an additional analysis tak-ing into account each galaxy’s likely H i mass and position relativeto the beam center. For each potentially confused observation (iden-tified as described above), we estimate the likelihood distributions of 𝑀 𝐻 𝐼 for the primary and secondary galaxies using their color andsurface brightness as described in Section 2.7.2 and Appendix B.We estimate the relative flux likelihood distributions of each galaxy by scaling the H i mass likelihood distributions based on the relativebeam power at their angular distance from the beam center, where thisscale factor is estimating assuming a Gaussian beam with FWHMappropriate for Arecibo of GBT observations. We then calculate thelikelihood distribution of the ratio, 𝑅 , of the flux, 𝐹 HI , from from allsecondary galaxies to that of all galaxies: 𝑅 = Í 𝐹 HI , secondary 𝐹 HI , pr imary + Í 𝐹 HI , secondary (1)We then estimate 𝑃 𝑅> . , the probability that 𝑅 is greater than 0.2,i.e., the probability that the companions are contributing at least 20%of the total measured flux (a value of 20% is chosen as it is comparableto typical flux calibration uncertainties). 𝑃 𝑅> . serves as means ofincorporating galaxies back into an analysis which likely have reliableflux measurements even in the presence of nearby companions. Forexample, only considering galaxies with 𝑃 𝑅> . > . as confusedlowers the total confusion rates of GBT and ALFALFA data by12% and 22%, respectively. We incorporate 𝑃 𝑅> . into our catalog(Section 2.4) so that any user can make their own judgment as towhat value is appropriate for their analysis. The HI-MaNGA DR2 catalog can be found at https://greenbankobservatory.org/science/gbt-surveys/hi-manga/ .This release includes the galaxies from HI-MaNGA DR1(Masters et al. 2019), and all data products here supersede thosefrom the earlier catalogue. The data will also be released as part ofa Value Added Catalogue with SDSS-IV DR17, which will enableeasy access alongside all MaNGA data through the
Marvin pythonpackage (Cherinka et al. 2019).The DR2 catalog contains H i information for MaNGA galaxiesin the MaNGA Product Launch (MPL) 8 sample, although does notprovide H i data for all 𝑆 / 𝑁 , whereas ALFALFA data have a minimum 𝑆 / 𝑁 of 4.5. Thereare 149 galaxies with both GBT and ALFALFA data (these GBTdata were mostly taken prior to the final release of the ALFALFAcatalog when the final footprint was unknown), and we include bothobservations for these galaxies in this catalogue. See Section 2.5for a description of the criteria used in this work to select whichobservation is incorporated into our final sample.A complete description of the quantities in the final catalog is asfollows: MaNGA-ID, Plate-IFU:
Internal MaNGA-ID and Plate-IFUand designations. 𝜶 , 𝜹 : Galaxy Right Ascension and Declination, J2000 (takenfrom manga DRPALL catalog ). log 𝑴 ∗ : Base-10 log of the stellar mass. Taken from the DRPALLcatalog, which is turn is taken from the NSA catalog. MNRAS , 1–23 (2020)
D. V. Stark et al. sin 𝒊 : An estimate of the sin of the inclination using the axial ratio, 𝑏 / 𝑎 :sin 𝑖 = s − ( 𝑏 / 𝑎 ) − 𝑞 (2)where 𝑞 = . 𝑏 / 𝑎 is takenfrom the DRPALL catalog. 𝑽 opt : Systemic velocity derived from the optical spectrum of thistarget (taken from MaNGA DRPALL catalog). Session:
Session IDs during which observations were con-ducted. For GBT observations, the format is [programID]_[session_number] (e.g.,
AGBT17A_012_01 ). For alldata from the ALFALFA survey, the session is listed as
ALFALFA . 𝒕 𝒊𝒏𝒕 : Total on-source integration time in seconds. 𝒓𝒎𝒔 : rms noise in mJy measured in a signal free part of the H i spectrum. log 𝑴 HI , lim , For non-detections, the base-10 log of the H i mass upper limit (in 𝑀 ⊙ ). For GBT data, this is calculated as 𝑀 𝐻 𝐼 𝑀 ⊙ = × . × ( + 𝑧 ) (cid:18) 𝐷 Mpc (cid:19) (cid:18) 𝑟𝑚𝑠𝐽𝑦 (cid:19) s(cid:18) 𝑊 km s − (cid:19) (cid:18) 𝑑𝑉 km s − (cid:19) (3)where 𝑊 is the assumed linewidth, 𝐷 is the luminosity distance to thegalaxy, and 𝑑𝑉 is the velocity resolution. For ALFALFA data, the up-per limits are estimated by reordering Equation (4) from Haynes et al.(2018) to determine the flux of a 𝑆 / 𝑁 = . 𝑆 / 𝑁 for an detection with an existing optical red-shift in the a.100 catalog), then converting to HI mass: 𝑀 𝐻 𝐼 𝑀 ⊙ = . × . × ( + 𝑧 ) (cid:18) 𝐷 Mpc (cid:19) (cid:18) 𝑊 km s − (cid:19) (cid:18) 𝑟𝑚𝑠 Jy (cid:19) 𝑤 − / 𝑠𝑚𝑜 (4)where 𝑤 𝑠𝑚𝑜 = 𝑊 /
20. We assume 𝑊 =
200 km s − for all upperlimits. 𝒇 peak : The peak flux density of the detected H i emission line. 𝑺 / 𝑵 : The peak signal-to-noise ratio defined as 𝑓 peak / 𝑟𝑚𝑠 . 𝑭 HI For detections, the integrated flux of the H i line in Jy km s − . log 𝑴 𝑯𝑰 : For detections, base-10 log of the H i mass (in solarmasses) defined as 𝑀 𝐻 𝐼 𝑀 ⊙ = . × ( + 𝑧 ) (cid:18) 𝐷 Mpc (cid:19) (cid:18) 𝐹 𝐻 𝐼
Jy km s − (cid:19) , (5) 𝑾 𝑴 : Width of the H i line measured at 50% of the mean fluxdensity of the two peaks (or the overall peak flux density in thecase of a single-peaked profile). The locations of the 50% peak fluxdensity levels on either side of the H i profile are estimated usinglinear interpolation. 𝑾 𝑷 : Width of the H i line measured at 50% of the overall peakflux density of the H i line. The locations of the 50% peak fluxdensity levels on either side of the H i profile are estimated usinglinear interpolation. 𝑾 𝑷 : Width of the H i line measured at 20% of the overall peakflux density of the H i line. The locations of the 20% peak fluxdensity levels on either side of the H i profile are estimated usinglinear interpolation. 𝑾 𝑷 : Width of the H i line measured at 50% of the peak fluxdensity on either side of the profile (or the overall peak flux densityin the case of a single-peaked profile). The location of the 50% peakflux density level is estimated using linear interpolation. 𝑾 𝑭 : Width of the H i line measured at 50% of the peak fluxdensity on either side of the profile. The 50% flux density level isestimated using linear fits ( 𝑓 = 𝑎 + 𝑏𝑣 ) to each side of the profilebetween 15-85% of the peak flux (minus the rms noise) on eitherside of the profile. 𝚫 𝒔 : Correction for linewidth overestimation due to instrumentalbroadening following Springob et al. (2005), defined as Δ 𝑠 = Δ 𝑣𝜆 + 𝑧 (6)This factor is already subtracted from the linewidths tabulated here,but is provided here in case the user wants to remove it and applyalternative corrections to the GBT data. We do not provide thisvalue for ALFALFA data since Haynes et al. (2018) only provideslinewidths that are already corrected. 𝑽 HI : Centroid of the H i line detection, calculated by taking themidpoint of the velocities at the 50% peak flux density level oneither side of the H i profile determined when calculating 𝑊 𝐹 . 𝑷 𝒓 : The peak H i flux at the high velocity peak (identical to 𝑃 𝑙 for a single-peaked profile). 𝑷 𝒍 : The peak H i flux in the low velocity peak (identical to 𝑃 𝑟 fora single-peaked profile). 𝒂 𝒓 : 𝑦 -intercept fit parameter for high velocity side of the H i profile, used to calculate 𝑊 𝐹 . 𝒃 𝒓 : slope fit parameter for high velocity side of the H i profile,used to calculate 𝑊 𝐹 . 𝒂 𝒍 : 𝑦 -intercept fit parameter for low velocity side of the H i profile, used to calculate 𝑊 𝐹 . 𝒃 𝒍 : slope fit parameter for low velocity side of the H i profile,used to calculate 𝑊 𝐹 . 𝝈 𝑽 Hi Error on 𝑣 H i , defined as 𝜎 𝑉 𝐻𝐼 = s(cid:18) 𝑟𝑚𝑠𝑏 𝑙 (cid:19) + (cid:18) 𝑟𝑚𝑠𝑏 𝑟 (cid:19) (7) confusion flag: Flag to indicate possible confusion of the H i profile (0 = not confused, 1 = possibly confused). 𝑷 𝑹> . : Probability that more than 20% of the total measuredflux comes from galaxies other than the primary. See Section 2.3 fordetails. MNRAS , 1–23 (2020)
I-MaNGA z l o g M H I ( M ⊙ ⊙ GBT upper limitALFALFA upper limitGBT detectionALFALFA detection
Figure 1.
Derived H i mass versus redshift for the HI-MaNGA DR2 catalog.Upper limits from ALFALFA are 4.5 𝜎 while GBT upper limits are 3 𝜎 . Themore conservative upper limits for ALFALFA data are due to the imposedminimum 𝑆 / 𝑁 in the a.100 catalog of detections. OFF detection flag:
Flag to indicate the presence of an H i source in the OFF observation which has a similar velocity as theredshift of our target. The presence of H i emission in the OFF beammay lead to an underestimation of the flux from our primary target. Baseline structure flag : Flag indicating the presence of strongbaseline variations which limit our ability to measure H i sourceproperties. Such structure can hamper fitting a smooth baseline tothe data, and may have structure on scales similar to typical galaxylinewidths, making differentiating baseline structure from true signalchallenging.Any missing measurements (e.g., linewidths for HI non-detections)are represented by -999 in the catalog. All linewidths include a red-shift stretching correction factor of ( + 𝑧 ) − , as well as an instrumen-tal broadening correction following Springob et al. (2005). We havenot applied corrections for turbulent motion or galaxy inclination dueto their potentially large uncertainties, but for convenience we haveincluded inclination estimates in our catalog. Multiple linewidth esti-mates are provided to enable comparison with other studies, althoughwe consider 𝑊 𝐹 the most reliable because it is the least sensitiveto low 𝑆 / 𝑁 data at the edges of H i profiles (Giovanelli et al. 1997;Springob et al. 2005). However, all linewidths can become unreliableat 𝑆 / 𝑁 . i profiles (Stark et al. 2013).Figure 1 shows the derived H i mass as a function of redshift for ourcatalog, including both new GBT data and the ALFALFA crossmatch.Fig 2 shows the sky coverage of HI-MaNGA observations.We also make available fits and csv files containing the spectrumof each galaxy. These files contain the following columns: 𝒗 𝑯𝑰 : Barycentric recession velocity for each channel (opticalconvention). 𝑭 𝑯𝑰 : Channel flux density after baseline subtraction. 𝑩 𝑯𝑰 : Channel flux density prior to baseline subtraction; providedso a user can subtract their own baseline if desired. Our parent sample is from the internal MaNGA Product Launch 8(MPL-8), which is composed of 6583 unique galaxies with a roughlyuniform stellar mass distribution from log 𝑀 ∗ / 𝑀 ⊙ = . − . i data. TheMaNGA survey is designed such that more massive galaxies arechosen to be at larger distances, ensuring they can have spectrosopiccoverage out to the same relative radius as the rest of the sample. Thisdistance-mass dependence, combined with the fact that HI-MaNGAonly observes out to 𝑧 ∼ .
05, means the data set with H i is somewhatbiased against the most massive galaxies (log 𝑀 ∗ / 𝑀 ⊙ >
11; seeFig. 1 of Masters et al. 2019).For galaxies which have data from both GBT and ALFALFA, weuse the following criteria to determine which observation to use. Ifboth observations are detections, and neither are confused, we usethe higher 𝑆 / 𝑁 observation. If one of the observations is confused,we use the non-confused observation. A case where both observa-tions are confused is irrelevant because these are not included inany analysis. If both observations are nondetections, we use the ob-servation with the lowest rms noise level. If one observation is adetection and the other a non-detection, we use the detection as longas it is not flagged as confused. Otherwise, the upper limit is used.Any galaxies with observational artifacts which may impact the H i mass measurement, such as strong baseline oscillations or “negative"detections caused by a galaxy in the OFF beam at similar redshift asthe primary target, are rejected from our analysis.We use the Primary+ MaNGA subsample for our analysis, includ-ing the Primary sample with a flat stellar mass selection, as well as the“color-enhanced" sample designed to more evenly incorporate less-populated regions in color-stellar mass parameter space (Wake et al.2017). The IFU bundles measure spectra out to ∼ . 𝑅 𝑒 in this sub-sample, where 𝑅 𝑒 is the 𝑟 -band elliptical Petrosian half-light radius.We do not include the Secondary sample, with flat stellar mass se-lection and IFU coverage out to ∼ . 𝑅 𝑒 , as it lies at larger distanceand has a substantially higher H i non-detection rate.Our analysis focuses only on star forming galaxies, specifically anygalaxy which falls in the star-forming region of both the [O iii ]/H 𝛽 vs.[N ii ]/H 𝛼 and [O iii ]/H 𝛽 vs. [S ii ]/H 𝛼 Baldwin, Phillips & Telervich(BPT) diagrams of Kewley et al. (2006b), where the line ratios areintegrated out to 𝑅 𝑒 . We further require galaxies to have integratedH 𝛼 equivalent widths, EW(H 𝛼 ), measured over the same area, of > < EW ( 𝐻𝛼 ) <
6, making them “greenvalley" galaxies possibly transitioning between the star forming andpassive populations (Sánchez et al. 2014; Cano-Díaz et al. 2019). Wedo not include this subsample in any statistical analysis althoughthese data are shown in figures throughout this paper to illustratehow relationships between H i content and optical ISM diagnosticsmay evolve as galaxies migrate off the star-forming sequence. Toavoid active galactic nuclei (AGN), including those embedded in star-forming disks, we also require that the line ratios measured withinthe central 2.5 ′′ fall within the star-formation region of the BPTdiagrams. Future work will explore M H i / M ∗ in AGN and passivegalaxies.After these various selections, we are are left with a sample of 874galaxies, including 614 detections (591 of which have 𝑆 / 𝑁 >
5) and258 non-detections. These numbers do not include the 21 green valleygalaxies described above. Figure 3 illustrates the parent MaNGA
MNRAS , 1–23 (2020)
D. V. Stark et al. α [deg] −10010203040506070 δ [ d e g ] DR1 GBT dataDR2 GBT dataALFALFA data
ALFALFAApertif
Figure 2.
Distribution of observed MaNGA plates with follow-up H i observations. The footprints of ALFALFA and the expected Apertif Medium Deep surveysare overlaid. log M * /M ⊙ −15−14−13−12−11−10−9−8 s S F R [ y r − ] HI-MaNGA sampleAnalysis sample (detections)Analysis sample (non-detections)
Figure 3.
The HI-MaNGA sample and the analysis sample for this paper insSFR vs. log 𝑀 ∗ parameter space. Black contours show the distribution forthe parent MaNGA sample. SFRs taken from the Pipe3d analysis of MaNGAdata cubes (see Section 2.6. sample, the full HI-MaNGA catalog sample, and the sample used forour analysis in sSFR vs. 𝑀 ∗ space. Details of MaNGA instrumentation, observing strategy, data reduc-tion pipeline (DRP), data analysis pipeline (DAP), and data prod- ucts can be found in Drory et al. (2015), Law et al. (2015),Law et al.(2016), Yan et al. (2016b), Westfall et al. (2019), Aguado et al.(2019), and Belfiore et al. (2019). Although our parent sample istaken from MPL-8, we use more up-to-date data products from MPL-9. All measurements are derived from maps of each spectral quantityusing the
SPX binning scheme, where each 0 . ′′ × . ′′ spaxel isanalyzed independently. Below we describe each of the parametersmeasured and their use in characterizing the ISM. Unless otherwisenoted, spectral properties are all measured within 𝑅 𝑒 in order toprovide consistent measurements across our full sample. In caseswhere we measure line ratios, we measure the ratio spaxel-by-spaxeland then take the median within our region of interest. When mea-suring line ratios, we require all relevant lines to have 𝑆 / 𝑁 > kcorrect with
𝐺 𝐴𝐿𝐸 𝑋 +SDSSphotometry (Blanton & Roweis 2007).In part of our analysis, we use parameters measured from ananalysis of MaNGA data using
Pipe3d . We refer the reader toSánchez et al. 2016a and Sánchez et al. 2016b for details. Specif-ically, we make use of two estimates of star formation rate, onemeasured from extinction-corrected H 𝛼 flux and the other from stel-lar population synthesis modeling. A constant factor of − .
24 dex issubtracted from all
Pipe3d
SFRs to convert from a Salpeter IFM toa Chabrier initial mass function (IMF) used by the MaNGA DAP.
Equivalent width (EW) measures emission line flux normalized bythe strength of stellar continuum. When examining how M
H i / M ∗ scales with the strength of individual emission lines, it is appropriateto normalize emission line fluxes by the stellar continuum level in away that is analogous to how M H i / M ∗ normalizes the H i mass bythe stellar mass.The DAP measures equivalent widths for each bin, which in thecase of the SPX binning scheme, is each 0.5 ′′ × ′′ spaxel. The con- MNRAS , 1–23 (2020)
I-MaNGA tinuum flux level used in the calculation for each emission line andspaxel is recorded. Therefore, the integrated equivalent width is cal-culated simply by dividing the summed emission line fluxes by thesummed continuum flux densities within 𝑅 𝑒 . Our analysis considersall measured lines from 3727Å ([O ii ]) to 9548Å (P 𝜖 ). The metal abundance (metallicity) has a profound impact on the ISM,dramatically increasing the cooling rate and increasing the ability ofclouds to self-shield from ionizing radiation. There are many cali-brations between gas-phase metallicity and various strong emissionline ratios. Unfortunately, there are also significant disagreementsin the absolute metallicities provided by different calibrations, al-though relative metallicity variations tend to be more consistent(Kewley & Ellison 2008). For this study, we primarily care about rel-ative trends, and we use the N2O2 calibrator from Kewley & Dopita(2002) which predicts the metallicity (12 + log O/H) as a functionof the [N ii ]6585Å/[O ii ]3727Å ratio. The strong advantage of thiscalibration is its insensitivity to other local ISM properties, notablythe ionization parameter, 𝑞 (see Section 2.6.3), and contaminationby emission from diffuse interstellar gas, or DIG (Zhang et al. 2017).As will be discussed in Section 4.4.3, variations in the scaling ofM H i / M ∗ against EW may be associated with the fraction of emis-sion from DIG, necessitating the need for a metallicity calibrationthat is insensitive to exactly where the line emission is originating. Adisadvantage of the N2O2 method is it can only be used for metallic-ities above ∼ . 𝑍 / 𝑍 ⊙ . However, we expect all metallicities in oursample fall above this threshold.Since [O ii ] 3727Å and [N ii ] 6585Å lie at significantly differentwavelengths, we apply internal extinction corrections before measur-ing their ratio. We use the Balmer decrement to estimate 𝐴 𝑣 assumingan intrinsic H 𝛼 /H 𝛽 ratio of 2.86 (Osterbrock & Ferland 2006) andthe extinction curve of O’Donnell (1994). We require H 𝛼 and H 𝛽 tobe detected to 𝑆 / 𝑁 > 𝑞 The ionization parameter is the maximum velocity of an ionizationfront driven by a local radiation field, typically defined as, 𝑞 = 𝑆 𝐻 / 𝑛 𝐻 (8)where 𝑆 𝐻 is the flux of ionizing photons per unit area and 𝑛 is thehydrogen number density, including all ionized and neutral hydrogen(Kewley & Dopita 2002). This parameter serves as a useful indicatorof the ionization state of the ISM. We estimate 𝑞 as a function of[O iii ] 5008Å/[O ii ] 3727,3729Å and 12 + log O / H using the cali-bration of Kewley & Dopita (2002). We apply extinction correctionsto the [O iii ] and [O ii ] lines in the same manner as described inSection 2.6.2. Electron density, 𝑛 𝑒 , affects the collision rates of particles in theISM, which can in turn affect the strengths of the forbidden linesused in our analysis. The ratio [S ii ]6718Å/[S ii ]6732Å provides anestimate of electron density due to the fact that these lines are emittedby two different energy levels with very similar excitation energy,making their relative strengths depend primarily on the collisionstrength. The ratio is sensitive to variations in 𝑛 𝑒 between ∼ ∼ cm − corresponding to line ratios between ∼ . ∼ . II ]/H 𝛼 , [S II ]/H 𝛼 , [O I ]/H 𝛼 The [N ii ]6586Å, [S ii ]6718,6732Å, and [O i ]6302Å lines are lowionization lines which are emitted in partially ionized regions ofthe ISM. Around star-forming regions, these lines are emitted fromthe transition zone between the ionized and neutral gas at theedges of star-forming clouds. The strengths of these lines rela-tive to 𝐻𝛼 are thought to be significantly enhanced in the pres-ence of a harder ionizing field (relative to HII regions) and shocks,and they often coincide with regions ionized by AGN or massiveevolved stars (Veilleux & Osterbrock 1987; Stasińska et al. 2008;Cid Fernandes et al. 2011; Belfiore et al. 2016). These line ratios arealso strengthened in DIG by factors of several relative to HII re-gions (Reynolds 1985b,a; Reynolds et al. 1998; Haffner et al. 1999;Hoopes & Walterbos 2003; Madsen et al. 2006; Voges & Walterbos2006; Oey et al. 2007; Zhang et al. 2017). 𝛼 surface brightness, 𝜇 𝐻 𝛼
While H 𝛼 emission from HII regions is a known tracer of star for-mation (Kennicutt & Evans 2012), when measured over large scalesit can be contaminated by emission from DIG residing beyond HIIregions (Oey et al. 2007). The effective MaNGA PSF ( ∼ ′′ ) cor-responds to spatial scales of ∼ 𝛼 emission aris-ing from HII regions versus DIG in our data. However, H 𝛼 surfacebrightness, 𝜇 H 𝛼 can be used as a an indicator of the contribution ofDIG to the overall H 𝛼 emission, as supported by observations whichshow increasing low-ionization line strength with decreasing 𝜇 H 𝛼 (Oey et al. 2007; Zhang et al. 2017). We measure the median valueof 𝜇 𝐻 𝛼 within 𝑅 𝑒 to characterize the average impact of DIG in arelative sense throughout our sample. It is not corrected for internalreddening. Even the deepest H i surveys do not detect all optically identifiedtargets, but the H i mass upper limits derived from non-detectionsstill hold valuable information. Ignoring non-detections when ex-amining relationships between H i mass and other galaxy propertiescan be misleading, causing artificially strong or weak correlations,depending on the exact sampling function (e.g. Calette et al. 2018).HI-MaNGA detects approximately 55% of all observed targetsand 70% in our analysis sample. In the analyses presented in this pa-per, we employ two distinct methods to study population trends in thepresence of non-detections: (1) survival analysis with the generalizedKendall’s 𝜏 and the Akritas-Theil-Sen Estimator to conduct correla-tion tests and perform linear fits in the presence of non-detections,and (2) the Photometric Gas Fractions technique to replace upperlimits with H i mass estimates based on other galaxy properties. 𝜏 andAkritas-Theil-Sen Estimator Making inferences from data in the presence of non-detections (or“left censored" data) has been a pursuit of the specific branch ofstatistics known as “survival analysis", and there exist a numberof studies where these methods have been applied to astronomi-cal data (Feigelson & Nelson 1985; Schmitt 1985; Isobe et al. 1986;Akritas et al. 1995; Akritas & Siebert 1996; Feigelson & Babu 2012;Calette et al. 2018; Yesuf & Ho 2019; Rodríguez-Puebla et al. 2020).
MNRAS , 1–23 (2020)
D. V. Stark et al.
Our goal is to test the strengths and statistical significances of cor-relations between M
H i / M ∗ and other parameters, while also deter-mining the best linear relations, all with a significant fraction ofH i non-detections. To achieve this goal, we follow the methodologyoutlined by Akritas et al. (1995). First, to test correlation strengthsand significances, we employ a generalized version of Kendall’s 𝜏 correlation test which incorporates upper limits. Building directlyon the generalized Kendall’s 𝜏 is the Theil-Sen estimator to deter-mine a best-fit line (as Akritas et al. 1995 developed the generalizedTheil-Sen estimator used here, we refer to the method as the Akritas-Theil-Sen estimator, or ATS). This method determines the slope, 𝑎 ,between the response variable, 𝑦 , and the covariate, 𝑥 , by finding thevalue of 𝑎 where 𝜏 calculated between the residuals, 𝑦 − 𝑎𝑥 , and 𝑥 isapproximately zero. The 𝑦 -intercept is then found by estimating themedian of 𝑦 − 𝑎𝑥 . To find the median when the residuals contain upperlimits, we use the Kaplan-Meier estimator (Kaplan & Meier 1958) toestimate their survival function, 𝑆 ( 𝑥 ) . The survival function providesthe likelihood that a distribution has a value above 𝑥 , and relates tothe cumulative distribution function, 𝐹 ( 𝑥 ) , by 𝑆 ( 𝑥 ) = − 𝐹 ( 𝑥 ) . Theintercept is thus where the Kaplan-Meier estimator equals 0.5. TheKaplan-Meier estimator of the fit residuals also allows us to estimatethe 1- 𝜎 scatter around the mean relation, which is done by findingwhere it equals 0.16 and 0.84In practice we use the R packages cenken and survfit to conductthese analysis. To quantify uncertainties on the estimated slopes andintercepts, we use bootstrapping following the recommendation ofWilcox (2010). We resample 𝑁 pairs of data points from the 𝑁 original pairs (preserving whether each pair was an H i detectionor upper limit) and re-estimate the fit parameters. This process isrepeated 600 times and the middle 68% of the returned distributionsare used to define the confidence intervals on our fit parameters.The ATS estimator as applied here assumes all relationships arelinear. This assumption generally appears to be valid, although theremay be some exceptions, notably analyses involving ionization pa-rameter, 𝑞 . Kendall’s 𝜏 assumes a monotonic relationship betweenvariables, not necessarily a linear one, so should be a generally reli-able indicator of correlation strength.To test that the ATS estimator gives reliable and consistent linearfits in the presence of upper limits, we have conducted a series offits between M H i / M ∗ and 𝑔 − 𝑟 color for a mock data set where weincrease the fraction of censored data in a manner consistent withobserving H i in MaNGA galaxies to progressively shallower depth.We find the linear fit parameters returned by the ATS estimatorare very robust to variations in the fraction of censored data andconsistently agree well with the “true" linear fit to the uncensoreddata, whether it is determined using the Theil-Sen estimator or amore traditional least-squares approach, as long as the censoringfraction does not exceed ∼ i detections. Further details are provided in Section 3.The tests described here are by no means an exhaustive analysis ofthe limitations of the ATS estimator, but do assess its applicabilityto our own data set, and we find that it should perform reliably for our purposes. The ATS estimator is the primary tool used to conductstatistical analysis in the presence of H i non-detections throughoutthis work. An alternative approach to handing H i upper limits is to replace themwith estimates of H i content using other galaxy properties. In Sec-tion 4.1 we examine how our results may change using this approachas opposed to using the ATS estimator. We follow a method very simi-lar to that described in Eckert et al. (2015) to assign gas masses to ournon-detections. Using data from the RESOLVE survey , Eckert et al.(2015) fit a model to the 2D probability distribution of galaxies inthe M H i / M ∗ vs. modified color plane ( 𝑃 ( 𝑀 Hi / 𝑀 ∗ | color ) , wheremodified color refers to a linear combination of color with anothergalaxy parameter. The model is composed of two major compo-nents: H i detections which follow a linear relationship with mod-ified color with some scatter, and non-detections which begin toappear at larger modified colors and cluster around M H i / M ∗ = . 𝑃 ( 𝑀 Hi / 𝑀 ∗ | color ) , upper limits can be assigned new estimatesof M H i / M ∗ as long as their modified color is measured. However,MaNGA photometry comes from the NSA, and there are likely sys-tematic differences between NSA and RESOLVE photometry, so wedo not use the same modified color and fitted model parameters fromEckert et al. (2015) to predict M H i / M ∗ for non-detections. Instead,we crossmatch the RESOLVE survey with the NSA, then rerun theanalysis of (Eckert et al. 2015) to determine 𝑃 ( 𝑀 H i / 𝑀 ∗ | color ) . Thefull details of our analysis are given in Appendix B. We now examine how M
H i / M ∗ is correlated with different param-eters derived from MaNGA optical spectroscopy. In Section 3.1we characterize scaling relationships between M H i / M ∗ and emis-sion line equivalent widths, followed by testing correlations betweenM H i / M ∗ and other ISM diagnostics in Section 3.2. In Section 3.3,we examine whether the scatter in the M H i / M ∗ -EW relations can beattributed to variations in their typical ISM properties. H i / M ∗ vs. EW relations We begin by simply testing the correlation strengths between all mea-sured optical emission lines and M
H i / M ∗ . Instead of using emissionline luminosities, we use emission line equivalent widths (EWs),which normalize the emission line strength by the stellar continuumstrength in a manner analogous to the gas-to-stellar mass ratio. Ta-ble 1 summarizes the correlations between each emission line EWand M H i / M ∗ . The table includes the generalized Kendall’s 𝜏 corre-lation coefficient and the corresponding p-value, the linear fit coef-ficients and their uncertainties determined using the ATS estimator,and the scatter above and below the best-fit line ( 𝜎 ℎ , 𝜎 𝑙 ).As discussed in Section 2.7.1, we limit our statistical analysis toavoid regions of parameter space dominated by upper limits. Theseregions are identified by first binning along the 𝑥 -axis variable using7 bins centered on the median 𝑥 value and spanning the inner 98%of all data. Any bin with a non-detection rate of more than 50% isexcluded from the ATS estimator. We also only include bins withat least 10 data points to avoid regions of parameter space which https://resolve.astro.unc.edu MNRAS , 1–23 (2020)
I-MaNGA log EW([OII] 3729Å) (Å) −1.5.1.0.0.50.00.5 l o g M H I / M * max medianmin medianupper limitdetectionsupper limitsgreen valley ATS fit −0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50 log EW([OIII] 5008Å) (Å) −1.5−1.0−0.50.00.5 l o g M H I / M * −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 log EW([OI] 6302Å) (Å) −1.5−1.0−0.50.00.5 l o g M H I / M * log EW(Hα 6564Å) (Å) −1.5−1.0−0.50.00.5 l o g M H I / M * Figure 4. 𝐸𝑊 vs. M H i / M ∗ for a subset of optical emission lines analyzed in this work. The three oxygen lines ([O i ],[O ii ],[O iii ]) show the strongest correlationsamong all emission lines. 𝐻 𝛼 (bottom-right) is significantly weaker. Black lines show the linear fit derived from the Akritas-Theil-Sen estimator which accountsfor upper limits (see Section 2.7.1). Large green filled and black unfilled points represent binned medians under two different assumptions about the upper limits:either the true M H i / M ∗ is equal to the measured upper limit, or it is zero. Inverted triangles represent bins which have at least 50% non-detections, such thatthere minimum binned median is zero. Thick orange points represent galaxies in the green valley. are poorly sampled. In the vast majority of cases, excluding thesedata have little impact on our results, but we nonetheless adopt thisapproach out of caution. This fitting approach is adopted throughoutour paper.All lines, with the exception of [N i ]5201Å have statistically sig-nificant correlations with M H i / M ∗ , although not all correlations areparticularly strong. In all cases, the scatter around the best-fit lineis typically asymmetric, skewing towards larger values, by as muchas a factor of ∼ H i ; we have rerun our analysis settingthe minimum required 𝑆 / 𝑁 at 3, 5, and 10, and the scatter decreasesby ∼
5% at most. Alternatively, some of the scatter may be driven bythe significantly different spatial scales over which optical and 21cmline emission is measured in our sample. In Section 4.3, we explorethis issue further.Notably, M
H i / M ∗ is very weakly correlated ( 𝜏 = .
13) withEW(H 𝛼 ), despite the fact that both 𝐸𝑊 ( 𝐻𝛼 ) and M H i / M ∗ are knownto relate to sSFR (this is further discussed in Section 4.4.1). In con-trast, the equivalent widths of the Oxygen forbidden lines ([O i ],[O ii ], and [O iii ]) all show the strongest correlation coefficients, with 𝜏 = . − . H i / M ∗ -EW(O)relation in Section 4.4.2. Figure 4 illustrates the observed scaling re-lations between M H i / M ∗ and equivalent widths of 𝐻𝛼 , [O i ], [O ii ],and [O iii ], with binned medians and the best fit line using the ATSestimator overlaid.M H i / M ∗ also shows a strong correlation ( 𝜏 = .
35) with[Ne iii ]3869/3968Å equivalent width. The [Ne iii ] and [O iii ] linesare already known to strongly correlate due to their similar ionizationstructure and constant Ne/O abundances (Pérez-Montero et al. 2007).Thus, since EW([O iii ] correlates well with M
H i / M ∗ , EW([Ne iii ])is almost guaranteed to as well.As discussed in Section 2.5, there are a small number of galaxiesthat have star formation as the primary ionization source, but theyhave 3 < 𝐸𝑊 ( 𝐻𝛼 ) < MNRAS , 1–23 (2020) D. V. Stark et al.
Table 1. M H i / M ∗ correlation statisticscovariate 𝜏 p-value ∗ slope intercept 𝜎 † ℎ 𝜎 † 𝑙 Equivalent Widths [ O ii ] < . × − ± ± [ O ii ] < . × − ± ± < . × − ± ± < . × − ± ± 𝜃 < . × − ± ± 𝜂 < . × − ± ± [ Ne iii ] < . × − ± ± [ He i ] < . × − ± ± 𝜁 < . × − ± ± [ Ne iii ] < . × − ± ± 𝜖 < . × − ± ± 𝛿 < . × − ± ± 𝛾 < . × − ± ± [ He ii ] . × − ± ± 𝛽 . × − ± ± [ O iii ] < . × − ± ± [ O iii ] < . × − ± ± [ N i ] . × − ± ± [ N i ] . × − ± ± [ He i ] < . × − ± ± [ O i ] < . × − ± ± [ O i ] < . × − ± ± [ N ii ] . × − -0.54 ± ± 𝛼 . × − ± ± [ N ii ] . × − -0.54 ± ± [ S ii ] < . × − ± ± [ S ii ] < . × − ± ± [ He i ] < . × − ± ± [ Ar iii ] < . × − ± ± [ Ar iii ] < . × − ± ± 𝜂 . × − ± ± [ S iii ] < . × − ± ± 𝜁 < . × − ± ± [ S iii ] < . × − ± ± 𝜖 . × − ± ± [ N ii ]/ 𝐻 𝛼 -0.43 < . × − -1.80 ± ± [ S ii ]/ 𝐻 𝛼 < . × − ± ± [ O i ]/ 𝐻 𝛼 < . × − ± ± [ O iii ]/ 𝐻 𝛽 < . × − ± ± [ S ii ] /[ S ii ] . × − ± ± 𝑞 -0.04 . × − -0.32 ± ± 𝜇 𝐻 𝛼 -0.12 . × − -0.22 ± ± + log 𝑂 / 𝐻 (N2O2) -0.46 < . × − -2.70 ± ± ∗ Values of < . × − represent the limits of numerical precision in the R code used for thisstatistical test. † These refer to the scatter above and below the fitted line. low equivalent widths (as one would expect) but they generally fallalong the fitted relation. This behavior will change somewhat inSection 3.2. H i / M ∗ and ISM diagnostics Next we examine how M
H i / M ∗ correlates with the other ISM di-agnostics discussed in Section 2.6. The parameters describing eachcorrelation are provided in Table 1 and all the relations are plotted inFigure 5.Again, statistically significant correlations are found with everysingle parameter, even though several correlation coefficients arequite small. M H i / M ∗ has only a very mild dependence on [S ii ]/H 𝛼 , [S ii ]6718Å/[S ii ]6732Å, 𝑞 , and 𝜇 𝐻 𝛼 . Our analysis also highlightshow the range of [S ii ]6718Å/[S ii ]6732Å is very narrow, implyingroughly constant typical electron densities across all galaxies in oursample. Similarly, our sample falls within a narrow range of 𝑞 , exceptfor a small subset of galaxies (typically at M H i / M ∗ <
1) which skewtowards higher values.Much stronger correlations are found with [N ii ]/H 𝛼 , [O i ]/H 𝛼 ,[O iii ]/H 𝛽 , 12+log O/H. The strong relationship between M H i / M ∗ and 12+log O/H is consistent with previous studies, and can be un-derstood as relating to the correlation between M H i / M ∗ and stellarmass, the stellar mass-metallicity relation (MZR), and anticorre-lation of the residuals in the MZR with H i mass (Lequeux et al.1979; Skillman et al. 1996; Peeples et al. 2008; Robertson et al. MNRAS , 1–23 (2020)
I-MaNGA −1.0 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 log [NII]/Hα −1.5−1.0−0.50.00.5 l o g M H I / M * −0.60 −0.55 −0.50 −0.45 −0.40 −0.35 −0.30 −0.25 log [SII]/Hα −1.5−1.0−0.50.00.5 l o g M H I / M * −1.8 −1.6 −1.4 −1.2 −1.0 log [OI]/Hα −1.5−1.0−0.50.00.5 l o g M H I / M * −0.6 −0.4 −0.2 0.0 0.2 log [OIII]/Hβ −1.5−1.0−0.50.00.5 l o g M H I / M * [SII]6718Å/[SII]6732Å −1.5−1.0−0.50.00.5 l o g M H I / M * log q (cm s −1 ) −1.5−1.0−0.50.00.5 l o g M H I / M *
12 + log O/H (N2O2) −1.5−1.0−0.50.00.5 l o g M H I / M * log μ Hα (erg s −1 kpc −2 ) −1.5−1.0−0.50.00.5 l o g M H I / M * Figure 5. M H i / M ∗ vs. various ISM diagnostics based on optical emission lines ratios. Symbols are the same as in Figure 4.MNRAS , 1–23 (2020) D. V. Stark et al. i mass,or M H i / M ∗ is beyond the scope of this work.Both [N ii ]/H 𝛼 and [O iii ]/H 𝛽 are strongly correlated with metal-licity, so [N ii ]/H 𝛼 , [O iii ]/H 𝛽 , and 12+log(O/H) are all effectivelyillustrating the same M H i / M ∗ -metallicity relation. The metallicitydependence of [O i ]/H 𝛼 is significantly weaker, and disappears be-low 12 + log 𝑂 / 𝐻 ∼ .
05 in our data. We rerun the correlation testsfor [O i ]/H 𝛼 , limiting to 12 + log 𝑂 / 𝐻 < .
05, and find that theM
H i / M ∗ vs. [O i ]/H 𝛼 correlation persists at a statistically signifi-cant level (but with slightly lower correlation coefficient of 0.23 andp-value 9 . × − ).In contrast to Figure 4, the green valley points often skew be-low the best-fit line to star-forming galaxies. This behavior providessome idea of how the observe relations may vary as galaxies migrateoff the star-forming sequence. However, examining the green valleypopulation in further detail is beyond the scope of this work. H i / M ∗ − EW(O) residuals
Given that the emissivity of optical emission lines are typically de-pendent on the local conditions of the ISM, we explore whether thescatter in the M
H i / M ∗ vs. EW relations can be linked to variationsin the average properties of the ISM. For this analysis, we focus ourattention on the residuals, 𝛿 M H i / M ∗ , of the strongest correlationsusing the three strong oxygen lines. We test for correlations between 𝛿 M H i / M ∗ and each of the ISM properties described in Section 2.6.Table 2 provides Kendall’s 𝜏 and p-value for correlations between 𝛿 M H i / M ∗ and the various ISM diagnostics, while Figure 6 plotsthese residuals for the M H i / M ∗ -EW([O iii ]) relation. Correlation co-efficients are generally low but statistically significant. The strongestcorrelations ( 𝜏 ∼ . − .
3) are consistently found with 𝜇 𝐻 𝛼 and[O i ]/H 𝛼 , regardless of which M H i / M ∗ -EW(O) relation is consid-ered. In some select cases, other ISM diagnostics show similarlystrong correlations with 𝛿 M H i / M ∗ (e.g., the residuals when usingthe [O iii ] line are strongly correlated with [S ii ]/H 𝛼 ) likely owingto the different conditions under which these lines originate. In con-trast, correlations of 𝛿 M H i / M ∗ with 𝑞 , [S ii ]6732Å/[S ii ]6718Å, and[O iii ]/H 𝛽 are consistently the weakest, with | 𝜏 | . . In Section 3, we showed that out of all observed optical emissionlines, M
H i / M ∗ in star-forming galaxies correlates best with the EWof [O ii ], [O iii ], and [O i ], all to a similar degree. Correlation strengthwith most other lines, most notably H 𝛼 , is significantly weaker. Wealso recover a strong anti-correlation between M H i / M ∗ and gas-phase metallicity, as well as a a positive correlation between M H i / M ∗ and [O i ]/H 𝛼 (independent of metallicity). Furthermore, we foundthat the residuals of these scaling relations, 𝛿 M H i / M ∗ , correlate bestwith the mean 𝐻𝛼 surface brightness and [O i ]/H 𝛼 .In the following section we investigate and interpret these resultsin further detail. In Section 4.1, we compare our results to anotheranalysis approach where H i upper limits are replaced with estimatesbased on photometric parameters. In Section 4.2, we discuss whetheroptical spectroscopic data may be a useful means of indirectly es-timating M H i / M ∗ , similar to how broadband colors are used. InSection 4.3, we explore how the aperture over which emission lineproperties are measured affects the resulting relations, specificallywhether tighter scaling relations can be obtained with larger aper-tures. Lastly, in Section 4.4, we discuss what the observed trends may be telling us about the properties of the ISM in galaxies as a functionof their gas-richness. As described in Section 2.7.2, another approach to incorporatingupper limits into our analysis is to replace them with estimatesof M
H i / M ∗ based on existing scaling relations with already mea-sured photometric properties. Using the methods described in Sec-tions 2.7.2 and Appendix B, we replace H i upper limits with newestimates and compare the derived correlation properties to thosedetermined using the ATS estimator and our original upper limits.Figure 7 shows the M H i / M ∗ -EW(O) relations using the replacedupper limits. The black lines in Figure 7 represent the new linear fits,while the red lines represent the original linear fits from Section 3and Figure 4. Although using photometric gas fractions removesupper limits and the need for survival analysis, we still use the ATSestimator for the linear fits (with all data now treated as detections) tobe as consistent as possible with our earlier analysis. The two fits arenot drastically different. The values of 𝜏 when using the photometric-gas fraction approach are slightly higher (0.40 −− 𝜎 ℎ , 𝜎 𝑙 ), which are still asymmetric. The locusof points around M H i / M ∗ ∼ − . H i / M ∗ assumed by the photometric gas fractionmethod used here. If we rerun the ATS fit treating these values asupper limits, we obtain very similar results. We avoid consideringeither the ATS or photometric gas fraction approach to upper limitsto be “best", but our analysis illustrates that both approaches giveconsistent results. H I / M ∗ A number of studies have identified tight correlations betweenM
H i / M ∗ and broadband photometric properties, most notably color,but even tighter scaling relations have been formulated by incorpo-rating additional parameters as well (e.g., surface brightness, axialratio; see Section 1). Here, we ask the question of whether similarlytight scaling relations can be formulated using spectroscopic infor-mation. One such attempt was made by Brinchmann et al. (2013),who was able to estimate local gas surface densities to within a fac-tor of two using optical spectroscopy. However, they were unable toextend their approach to provide reliable estimates of integrated gasmass; their spectroscopic data was limited to the central 3 ′′ of eachgalaxy, making the difference in spatial scales between the opticaland 21cm data even more pronounced than in our study. The authorsexplored using aperture corrections, but found the systematic errorswere too large to reliably estimate gas masses.The tightest photometric gas fraction estimators have scatters of ∼ .
3, or a factor of two (e.g. Catinella et al. 2013; Eckert et al.2015). A clear disadvantage of the scaling relations in this work isthe asymmetric scatter. Although the upper scatter approaches 0.3dex in some cases, the corresponding lower scatter typically > . H i / M ∗ and optical emission line prop-erties. For this purpose we combine EW(O) and 𝜇 𝐻 𝛼 , which con-sistently correlates most strongly with the residuals of the M
H i / M ∗ - MNRAS , 1–23 (2020)
I-MaNGA −1.0 −0.8 −0.6 −0.4 −0.2 log [NII]/Hα −1.0−0.50.00.5 l o g M H I / M * s l o g E W ([ O III ] Å ) r e s i d u a l s −0.6 −0.5 −0.4 −0.3 −0.2 log [SII]/Hα −1.0−0.50.00.5 l o g M H I / M * v s l o g E W ([ O III ] Å ) r e s i d a l s −2.0 −1.8 −1.6 −1.4 −1.2 −1.0 log [OI]/Hα l o g M H I / M * v s l o g E W ([ O III ] Å ) r e s i d u a l s −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 log [OIII]/Hβ l o g M H I / M * v s l o g E W ([ O III ] Å ) r e s i d u a l s q (cm/s) −1.00−0.75−0.50−0.250.000.250.500.75 l o g M H I / M * v s l o g E W ([ O III ] Å ) e s i d u a l s [SII]6718Å/[SII]6732Å −1.0−0.50.00.5 l o g M H I / M * v l o g E W ([ O III ] Å ) r e i d u a l
12 + log O/H (N2O2) −1.0−0.50.00.5 l o g M H I / M * v l o g E W ([ O III ] Å ) r e i d u a l log μ Hα (e g s −1 kpc −2 ) −1.0−0.50.00.5 l o g M H I / M * v s l o g E W ([ O III ] Å ) e s i d u a l s Figure 6.
Correlations between the residuals of the M H i / M ∗ -EW([O iii ]) relation vs. various ISM diagnostics. Symbols are the same as in Figure 4.MNRAS , 1–23 (2020) D. V. Stark et al.
Table 2.
Correlation test 𝜏 and (p-value) for M H i / M ∗ vs. EW residuals [ O ii ] [ O iii ] [ O i ] [ N ii ]/ 𝐻 𝛼 -0.14( . × − ) -0.13( . × − ) -0.19( . × − ) [ S ii ]/ 𝐻 𝛼 . × − ) 0.22( < . × − ) 0.17( . × − ) [ O i ]/ 𝐻 𝛼 < . × − ) 0.26( < . × − ) 0.20( . × − ) [ O iii ]/ 𝐻 𝛽 . × − ) 0.07( . × − ) 0.16( . × − ) [ S ii ] /[ S ii ] -0.06( . × − ) -0.08( . × − ) -0.08( . × − ) 𝑞 -0.07( . × − ) -0.14( . × − ) -0.06( . × − ) + log 𝑂 / 𝐻 (N2O2) -0.14( . × − ) -0.14( . × − ) -0.20( . × − ) 𝜇 𝐻 𝛼 -0.23( < . × − ) -0.28( < . × − ) -0.25( < . × − ) log EW [OII] 3729Å −1.5−1.0−0.50.00.5 o g M H I / M * detectionsreplaced upper imitsnew ATS fitorigina ATS fit −0.5 0.0 0.5 1.0 1.5 2.0 og EW [OIII] 5008Å −0.8 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 og EW [OI] 6302Å Figure 7. M H i / M ∗ vs. EW(O) correlations where upper limits have been replaced by new estimates of M H i / M ∗ using the photometric gas fraction technique(orange points). The black line shows the ATS fit to these data (where all points are treated as detections), while the red line shows the ATS fit to the data withthe original upper limits. EW(O) relations. The final relations, shown in Figure 8, are:log 𝑀 H i / 𝑀 ∗ = ( . ± . ) log 𝐸𝑊 ([ O ii ] )+ (− . ± . ) log 𝜇 𝐻 𝛼 + ( . ± . ) , (9)log 𝑀 H i / 𝑀 ∗ = ( . ± . ) log 𝐸𝑊 ([ O iii ] )+ (− . ± . ) log 𝜇 𝐻 𝛼 + ( . ± . ) , (10)log 𝑀 H i / 𝑀 ∗ = ( . ± . ) log 𝐸𝑊 ([ O i ] )+ (− . ± . ) log 𝜇 𝐻 𝛼 + ( . ± . ) . (11)The upper/lower scatter ( 𝜎 ℎ / 𝜎 𝑙 ) are 0.31/0.42, 0.29/0.40, and0.31/0.41 for equations 9, 10, and 11, respectively. The correlationcoefficients are 0.4 − H i / M ∗ -[N ii ]/H 𝛼 and M H i / M ∗ -(12 + log 𝑂 / 𝐻 ) relations have similar values of 𝜏 and comparablescatter using only one parameter, the relations using EW(O) and 𝜇 𝐻 𝛼 have some notable advantages: they require no extinction cor-rections, the lines are less likely to be blended, and they are notdependent on specific photoionization models. Nonetheless, for thepurposes of using optical properties to predict the H i content ofgalaxies, existing relations in the literature using broadband colorare likely more reliable at this time. Improvements may be possible,however. For example, deeper H i data may allow us to identify re-gions of parameter space where the scatter is more symmetric, andusing optical emission line measurements over larger areas (or ap-plying aperture corrections) is likely to tighten these relations (seeSection 4.3). Throughout our analysis, we have been using H i measurements fromsingle-dish observations, where the beams (FWHM=3.5 ′ –9.1 ′ ) aresignificantly larger than the individual galaxies and likely contain all of their H i emission. In contrast, our emission line measurementsare conducted out to 𝑅 𝑒 . Although the MaNGA IFU coverage ofthe Primary+ sample nominally extends to 1.5 𝑅 𝑒 , there are slightvariations in coverage, so 𝑅 𝑒 was chosen as a scale over which mea-surements can be conducted consistently for all galaxies. However,H i disks are frequently observed to extend well beyond this radius(Broeils & Rhee 1997). It is possible that the difference in spatialscales being probed is contributing to at least some of the scatter inthe relationships presented in Section 3, unless all H i disks are self-similar such that they always contains the same fraction of H i inside 𝑅 𝑒 , but this is unlikely. Although Broeils & van Woerden (1994) andSwaters et al. (2002) do find similarities in radial H i distributions atfixed morphology, there is significant scatter around the mean distri-butions, with an additional dependence on environment also possible(Cayatte et al. 1994).To explore this issue further, we compare the M H i / M ∗ vs.EW([O i ]) correlation using emission line measurements conductedout to 𝑅 𝑒 / 𝑅 𝑒 , and 2 𝑅 𝑒 . The line measurements out to 2 𝑅 𝑒 arenot done uniformly; we simply take any data in the IFU which lieswithin 2 𝑅 𝑒 . Typically, the data at larger radius falls along the minoraxis and is visible within the IFU due to the a galaxy’s inclination.Although these measurements are not uniform, if a larger apertureimproves the M H i / M ∗ scaling relations, we should still see somesign of it, although the overall improvement may be underestimated.[O i ] is used here as an example; the results of this analysis also applyto the other lines. MNRAS , 1–23 (2020)
I-MaNGA −1.50−1.25−1.00−0.75−0.50−0.25 0.00 0.25 Hα + 12.73 −1.5−1.0−0.50.00.5 l o g M H I / M * −1.25 −1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 Hα + 15.16 −2.0 −1.5 −1.0 −0.5 0.0 0.5 Hα + 16.63 Figure 8.
Final M H i / M ∗ scaling relations combining EW(O) with 𝜇 𝐻 𝛼 . Symbols are the same as in Figure 4.
Table 3.
Fit parameters for different IFU radial coverage 𝑅 𝑒 / 𝑅 𝑒 𝑅 𝑒 M H i / M ∗ vs. EW([O i ])6302Å 𝜏 ± ± ± ± ± ± 𝜎 ℎ 𝜎 𝑙 𝛿 M H i / M ∗ vs. 𝜇 𝐻 𝛼 𝜏 -0.26 -0.25 -0.25slope -0.39 ± ± ± ± ± ± 𝜎 ℎ 𝜎 𝑙 M H i / M ∗ vs. EW([O i ]6302Å) + 𝜇 𝐻 𝛼 𝜎 ℎ 𝜎 𝑙 Table 3 gives the linear fit parameters to the M
H i / M ∗ -EW([O i ])and the 𝛿 M H i / M ∗ - 𝜇 𝐻 𝛼 relations, as well as the final upper andlower scatter in the combined scaling relation. Increasing the emis-sion line measurement aperture does appear to mildly strengthenthe M
H i / M ∗ -EW([O i ]) correlation strength, although it levels offonce an aperture of 𝑅 𝑒 is reached. The derived slope and interceptalso change significantly. The strength of the correlation between 𝛿 M H i / M ∗ and 𝜇 𝐻 𝛼 is essentially independent of IFU radial cov-erage, and while the derived slopes/intercepts change, they are inagreement within their error bars. Notably, increasing IFU radialcoverage decreases both the size of the scatter and the asymmetry inthe relations shown in Table 3, implying the different spatial scalesover which 21cm and optical emission is measured is contributing tothe scatter.As a consistency check, we conduct the same analysis as abovebut using the MaNGA Secondary sample (with radial coverage outto ∼ . 𝑅 𝑒 ) where spectroscopic measurements to 2 𝑅 𝑒 can includeall data along the major axis, but with the caveat that this subsamplesuffers from a smaller sample size and a larger rate of non-detectionsdue to its higher redshift distribution. We find similar qualitativeresults as described above, although all fit parameters suffer fromlarger errors.From this analysis, we conclude that the correlations identified inSection 3 persist regardless of the area over which the emission linesare measured, but measuring emission line properties over larger areas does appear to increase the correlation strength and decreasethe scatter. The improvement when using an aperture of 2 𝑅 𝑒 mayactually be underestimated due to our non-uniform measurementsout to this radius. Given the tightening of the final relation whenlarger apertures are used, it may be valuable to explore in futurework whether using even larger-area IFU bundles or applying carefulaperture corrections can yield even stronger relations. Our results provide a means of understanding the typical ISM condi-tions in star-forming galaxies as a function of their gas-richness. Weagain stress that we are discussing average
ISM properties, and thelocal conditions within galaxies almost certainly vary significantlyaround the average behavior. Furthermore, as was discussed in Sec-tion 4.3, the optical emission and H i emission are being measuredover significantly different areas. While the optical emission is mea-sured within 𝑟 -band effective radius ( 𝑅 𝑒 ), Broeils & Rhee (1997)find that the typical H i effective radius roughly equals the 𝐵 -band 25mag arcsec − isophotal radius, and thus the majority of H i emissioncomes from beyond 𝑅 𝑒 . Thus, we are not measuring the conditionsof all the H i in our sample. Rather, we are examining the ISM con-ditions in the star-forming disks of galaxies, and assessing how theserelate to H i -richness.In agreement with previous work, we find a clear anticorre-lation between gas content and metallicity (Lequeux et al. 1979;Skillman et al. 1996; Peeples et al. 2008; Robertson et al. 2012;Hughes et al. 2013; Brown et al. 2018). This result is implied boththrough our estimate of metallicity using the N2O2 predictor aswell as [N ii ]/H 𝛼 and [O iii ]/H 𝛽 which are very sensitive to metal-licity. The gas-phase metallicity shows the strongest variation withM H i / M ∗ of all properties estimated from emission lines. In contrast,[S ii ]6718Å/[S ii ]6732Å varies only slightly, with most of the datalying around ∼ ii ]6718Å/[S ii ]6732Å is minimally depen-dent on electron density at these values (Osterbrock & Ferland 2006),but most galaxies are consistent with having average electron den-sities .
100 cm − . Systematic variations in electron densities mayhave impacted other estimates of ISM properties which make as-sumptions about electron density, but thankfully this is not the case.Ionization parameter, 𝑞 , similarly tends to cluster around a narrowrange of values ( ∼ . cm s − ). The few galaxies which have 𝑞 athigher levels notably tend to only exist at M H i / M ∗ < H i / M ∗ and [O i ]/H 𝛼 , even MNRAS , 1–23 (2020) D. V. Stark et al. after removing the metallicity dependence of [O i ]/H 𝛼 . [O i ]/H 𝛼 issensitive to ionization fraction and gas temperature (Reynolds et al.1998). Its enhancement in gas-rich galaxies may indicate an on-average harder ionizing spectrum incident on the emitting gas,or more excitation from shocks (Dopita et al. 2000; Kewley et al.2006a). Increased shock excitation would be consistent withthe higher turbulent velocities expected in lower mass galaxies(Dalcanton et al. 2004). This trend may also be consistent with alarger contribution of optical line emission from DIG, especiallygiven the additional (albeit very weak) correlation between M H i / M ∗ and 𝜇 𝐻 𝛼 . HI / M ∗ correlate weakly with H 𝛼 emission? M H i / M ∗ is only weakly correlated with 𝐸𝑊 ( 𝐻𝛼 ) , a finding thatis consistent with previous work by Jaskot et al. (2015). This re-sult implies a relatively loose connection between the instantaneousspecific star formation rate (H 𝛼 traces star formation on ∼ i reservoir. Theweak correlation between M H i / M ∗ and EW(H 𝛼 ) is in stark contrastto the significantly tighter relationship between M H i / M ∗ and color,even though both color and EW(H 𝛼 ) can be thought of as tracers ofsSFR. We consider a number of possible explanations for the largescatter in the M H i / M ∗ -EW(H 𝛼 ) relation, especially in contrast tothe M H i / M ∗ -color relation: (1) scatter in the relationship betweenEW(H 𝛼 ) and sSFR, (2) 𝐸𝑊 ( 𝐻𝛼 ) and color measured over differentapertures, (3) dust extinction artificially tightening the M H i / M ∗ -color relation, and (4) the different timescales over which 𝐸𝑊 ( 𝐻𝛼 ) and color trace sSFR.We have been assuming that EW(H 𝛼 ) has a 1:1 relationship withsSFR. EW(H 𝛼 ) is the H 𝛼 line flux normalized by the stellar contin-uum flux density, while sSFR is star formation rate normalized bystellar mass. The mapping between stellar continuum flux densityand stellar mass may depend on the details of the stellar population(age, metallicity, IMF). Similarly, the H 𝛼 luminosity to SFR con-version makes a number of assumptions about the stellar IMF, gasmetallicity, gas temperature, and gas density (Kennicutt 1998), andif any of these assumptions are not true, then they might contributeto a more scattered M H i / M ∗ -EW(H 𝛼 ) relation. A full analysis ofthe possible error in the H 𝛼 luminosity to SFR conversion is beyondthe scope of this work, but we do analyze the sSFR-M H i / M ∗ corre-lation using extinction-corrected H 𝛼 -based SFR from Pipe3D andstellar masses from the NSA, and we find the correlation coefficientis almost exactly the same as when using EW(H 𝛼 ).The limited radial coverage of our spectroscopic measurementsmay drive at least some scatter in the M H i / M ∗ -EW(H 𝛼 ) relation.Our analysis in Section 4.3 implies that increasing the area overwhich optical spectroscopic properties are measured should tightenrelations with M H i / M ∗ . However, increasing the radial coverageof the spectrosopic measurements is unlikely to account for all thescatter; Jaskot et al. (2015) find a similarly weak M H i / M ∗ -EW(H 𝛼 )relation using narrow band H 𝛼 imaging where emission can be mea-sured globally. Nonetheless, we examine the change in correlationstrength when it is measured out to 1.5 𝑅 𝑒 , i.e., the radial extent ofthe MaNGA Primary+ IFU bundles. 𝜏 increases only slightly from0.13 to 0.16, but we cannot obtain a truly global EW(H 𝛼 ) measure-ment due to the limited IFU sizes. Therefore, we reverse the problem,and examine whether the M H i / M ∗ -color relation becomes a muchweaker correlation if the colors are measured within 𝑅 𝑒 , using the 𝑔 − 𝑟 color derived from images from the NSA (Figure 9). The origi-nal relationship using global 𝑔 − 𝑟 has 𝜏 = .
46 while the version with 𝑔 − 𝑟 measured within 𝑅 𝑒 has 𝜏 = .
41. We conclude that the limited radial coverage of our spectroscopic measurements can explain somebut not all of the weakness of the M
H i / M ∗ -EW(H 𝛼 ) relationship.Jaskot et al. (2015) argue that dust extinction artificially tightensthe M H i / M ∗ -color relationship, and when applying internal extinc-tion corrections, they obtain a correlation with M H i / M ∗ that hasscatter similar to that of the M H i / M ∗ -EW(H 𝛼 ) relation. However,we argue that Jaskot et al. (2015) may have over-corrected the broad-band colors for internal dust extinction. Specifically, they estimateinternal extinction from the WISE mid-infrared (MIR) dust emission,which will be weighted towards the densest, and thus more extincted,regions of the ISM. However, broadband colors can reflect long-termgrowth rates of galaxies on timescales of up to Gyrs (Kannappan et al.2013), so the stellar light is coming from stars which are very likelyno longer embedded in their birth clouds and not as extincted as theemission coming from more active star forming regions.An alternative approach to determine internal extinction correc-tions is to conduct a statistical analysis of galaxy colors as a functionof axial ratio (a proxy for inclination) and absolute magnitude (orstellar mass), under the assumption that the reddening will be morepronounced for edge-on galaxies and more massive galaxies withhigher metallicities (Hernández-Toledo et al. 2008; Masters et al.2010). We apply the statistical internal extinction corrections ofHernández-Toledo et al. (2008) to the 𝑔 − 𝑟 colors of our sample. Thecorrelation strength is only marginally weakened with the addition ofdust corrections, with 𝜏 changing from -0.47 to -0.42, and the weakercorrelation is still much stronger than the M H i / M ∗ -EW(H 𝛼 ) relation(Figure 9). As a final test, we assess the strength of the M H i / M ∗ -EW(H 𝛼 ) relation with both the dust corrections and colors measuredout to 𝑅 𝑒 , but find the combination still does not significantly affectthe correlation strength (Figure 9).Having largely ruled out more mundane explanations of theweak relationship between M H i / M ∗ and EW(H 𝛼 ), we now explorewhether there is a more physical explanation for the loose associationbetween H 𝛼 emission and H i content. A weak relationship betweenthese quantities may not actually be surprising when considering pre-vious work on the star formation law ( Σ SFR vs. Σ gas , where Σ refersto surface density) in nearby galaxies, where the results clearly showactive star formation rate is most closely correlated with the amountof molecular hydrogen (H ), not the amount of H i (Bigiel et al.2008). This result holds even when the two phases of hydrogen arecompared at fixed surface density (Schruba et al. 2011). If the globalH /H i were constant, or smoothly varying as a function of M H i / M ∗ ,then the M H i / M ∗ -EW(H 𝛼 ) relation might appear tighter. How-ever, Catinella et al. (2018) and Calette et al. (2018) illustrate thatH /H i scaling relations have very large scatter, and H /H i are moreclosely correlated with properties like stellar surface mass densityand bulge strength. Other studies have also highlighted how H /H i can change significantly due to interactions (Kenney & Young 1989;Braine & Combes 1993; Lisenfeld et al. 2011; Stark et al. 2013) andbar inflows (Sakamoto et al. 1999; Sheth et al. 2005). On small scaleswithin galaxies, molecular cloud formation is likely dependent ona number of local conditions, including the overall gas density,details of the ionizing radiation field, and gas-phase metallicity(Krumholz et al. 2009a). In summary, although H i “fuels" star for-mation by providing a gas reservoir from which molecular clouds canform, the detailed and complicated physics involved in converting H i into H on short timescales lead to a very messy M H i / M ∗ -EW(H 𝛼 )relation.Instead, we suggest the difference in timescales probed by EW(H 𝛼 )and color explains their different correlation strengths with M H i / M ∗ .EW(H 𝛼 ), being sensitive on 5 Myr timescales, most strongly corre-lates with the properties of the gas where those stars are actively MNRAS , 1–23 (2020)
I-MaNGA forming, namely the overall mass density of H . Meanwhile, globalcolors typically trace star formation over much longer timescales. Forexample, Kannappan et al. (2013) show that U-NIR color (referringto a near ultraviolet minus a near infrared color, e.g., 𝑢 − 𝐽 ) is agood proxy for the relative increase in stellar mass over the last Gyr .Therefore, although stars do not form in H i , it is the H i reservoir(and its likely regular replenishment) which sustains star formationby providing the material out of which molecular clouds and starsform, thus driving the M H i / M ∗ -color relationship.We do not actually need to average over ∼ Gyr timescales to see atightening in the relationship between M
H i / M ∗ and sSFR. Figure 10plots M H i / M ∗ vs. sSFR , where sSFR is the sSFR overthe last 32 Myr derived from Pipe3d stellar population synthesismodeling of the MaNGA datacubes. SFR is calculated by es-timating the SFR history from fitting a composite stellar populationto each spaxel spectrum, taking into account the contribution of dustextinction. The obtained SFR history is then averaged within a pe-riod of 32 Myr for each spaxel, then coadded across the full FOV inorder to obtain SFR . Although sSFR does not trace starformation over the extremely long timescales often probed by color,it still averages over timescales several times that of the H 𝛼 emis-sion. The resulting correlation is notably stronger ( 𝜏 = .
35) thanthe M
H i / M ∗ -EW(H 𝛼 ) relation. However, we must also consider thepossibility that at least some of the scatter in the M H i / M ∗ -EW(H 𝛼 )relation is due to the aforementioned potential systematic uncertain-ties in the H 𝛼 − SFR calibration. Discrepancies between the assumedand true conditions of star-forming regions across a wide range ofgalaxies may contribute to a more scattered relation of M
H i / M ∗ with H 𝛼 emission, while we recover cleaner relation with sSFR where the SFR is based on stellar population modeling. HI / M ∗ more strongly correlate with oxygenemission? M H i / M ∗ is most strongly correlated with EW(O), and the correla-tion strengths are approximately equal regardless of which availableoxygen line is considered. A correlation between H i and [O i ] isexpected to some extent, as oxygen has an ionization energy verysimilar to hydrogen (both approximately 13.6 eV), so they are ex-pected to co-exist within the ISM. Indeed [O i ] is typically observedin partially neutral regions of the ISM like the outer edges of neb-ula (Reynolds et al. 1998; Kewley et al. 2019). However, the actualemissivity of the [O i ] line is closely tied to the excitation mecha-nism, not just its presence in neutral gas. Furthermore, the [O iii ] and[O ii ] lines are found throughout the inner regions of star-formingnebulae (Kewley et al. 2019) where hydrogen is mostly ionized, yetthey show similarly strong correlations with neutral hydrogen.The link between M H i / M ∗ and EW(O) can more simply be under-stood as a result of the relation between M H i / M ∗ and 12 + log 𝑂 / 𝐻 .Simply put, at high enough metallicities, oxygen line emission will inversely correlate with oxygen abundance as a result of the highermetallicity increasing the cooling efficiency of the gas and lower-ing its overall temperature (this trend reverses at low metallicity,but we are well above that regime; McGaugh 1991). EW(O) de-pends on the total line flux and the stellar continuum flux density, The growth rate used in Kannappan et al. (2013) is defined as new starsformed in the last Gyr divided by the pre-existing stellar mass, so it is notexactly sSFR like used throughout our work where the denominator is totalstellar mass including recently formed stars. and as M
H i / M ∗ increases, two changes occur. The metallicity de-creases, thereby increasing the relative oxygen emission from gas inthe appropriate conditions, and the stellar mass/surface density willon-average decrease(Catinella et al. 2013), thus lowering the stellarcontinuum flux density and further boosting the EW. If we assumestar forming regions are the primary source of optical line emission,the oxygen line luminosity can be expressed as 𝐿 O ∝ (cid:18) 𝑂𝐻 (cid:19) 𝑛 𝑀 𝑔 (cid:18) 𝑀 SF 𝑀 𝑔 (cid:19) (12)where 𝑂𝐻 is the oxygen abundance relative to hydrogen, 𝑛 is somevalue less than zero to reflect the anti-correlation between oxygenabundance and relative line strength, 𝑀 SF is the gas mass in starforming regions, and 𝑀 𝑔 is the total gas mass. Since 𝐸𝑊 is theemission line flux normalized by the stellar continuum flux density,we can rewrite this as 𝐸𝑊 ( ) ∝ (cid:18) 𝑂𝐻 (cid:19) 𝑛 (cid:18) 𝑀 SF 𝑀 ∗ (cid:19) ∝ (cid:18) 𝑂𝐻 (cid:19) 𝑛 𝑠𝑆𝐹 𝑅 (13)where we have assumed 𝑀 SF is proportional to SFR. The assumptionthat oxygen line emission is only coming from star-forming regionsmay be an over-simplification, as will be discussed in the followingsection. Optical line emission will not only arise from star forming regions,but also from DIG (also known as the warm ionized medium, orWIM). The contribution of DIG to the overall emission line lu-minosity of a galaxy can be substantial. For instance, Oey et al.(2007) find on average ∼
60% of H 𝛼 emission arises from DIG,with the fraction being highest for low surface brightness galax-ies (they do not find any significant trend with H i fraction, butalso do not incorporate H i non-detections into their analysis). Metalline emission has also been observed from DIG, and in some casesis enhanced relative to Balmer hydrogen emission (e.g., [O ii ]/H 𝛽 ,[O i ]/H 𝛼 , [N ii ]/H 𝛼 , [S ii ]/H 𝛼 ; Reynolds 1985b,a; Reynolds et al.1998; Haffner et al. 1999; Hoopes & Walterbos 2003; Madsen et al.2006; Voges & Walterbos 2006; Zhang et al. 2017), which may im-ply an even larger fraction of the integrated metal line emission arisesfrom the DIG.To account for DIG, Eq 13 can be modified to read: 𝐸𝑊 ( 𝑂 ) = (cid:18) 𝑂𝐻 (cid:19) 𝑛 (cid:20) 𝐴 (cid:18) 𝑀 SF 𝑀 ∗ (cid:19) + 𝐵 (cid:18) 𝑀 DIG 𝑀 ∗ (cid:19) (cid:21) (14)where 𝐴 and 𝐵 are unknown constants, and M DIG is the amount ofemitting gas associated with DIG. We make no attempt to constrainthe unknowns in this equation here, and they probably vary withgalaxy properties (mass, surface brightness, etc.), but it highlightsour key takeaways: (1) the oxygen emission is correlated with gas-phase metallicity (2) both star-forming clouds and DIG contribute tothe overall emission. The relative contributions of HII regions andDIG to the M
H i / M ∗ -EW(O) relation may not be fixed, and likely varysmoothly on average along the trend for it to remain well-behaved(such smooth variations should just be absorbed into the linear fitcoefficients). Evidence of a growing contribution from DIG as H i fraction increases is seen in the correlation between M H i / M ∗ and[O i ]/H 𝛼 (and perhaps to a weaker extent, the anti-correlation with 𝜇 𝐻 𝛼 ). Additionally, in Section 3.3, we found that the scatter in theM
H i / M ∗ vs. EW(O) relations correlates with [O i ]/H 𝛼 and 𝜇 𝐻 𝛼 ,which suggests the scatter may be driven by variations in the DIGfraction. However, we caution that similar results could be explained
MNRAS , 1–23 (2020) D. V. Stark et al. g-r −1.5−1.0−0.50.00.5 l o g M H I / M * uncorrectedτ = -0.44 g-r r < r e , uncorrectedτ = -0.38 −0.2 0.0 0.2 0.4 0.6 0.8 g-r −1.5−1.0−0.50.00.5 l o g M H I / M * dust_correctedτ = -0.41 −0.2 0.0 0.2 0.4 0.6 0.8 g-r r < r e , dust-correctedτ = -0.40 Figure 9. M H i / M ∗ vs. 𝑔 − 𝑟 before and after correcting for internal extinction and/or limiting the color measurements to within 𝑟 𝑒 . Neither correction, nor thecombination of the two, lead to a relationship as weak as the M H i / M ∗ -EW(H 𝛼 ) relation. −11.25−11.00−10.75−10.50−10.25−10.00 −9.75 −9.50 −9.25 log sSFR [Gyr −1 ] −1.5−1.0−0.50.00.5 l o g M H I / M * Figure 10. M H i / M ∗ vs. sSFR measured over the last 32 Myr. Symbols are thesame as in Figure 4. This correlation is notably stronger than the relationshipbetween M H i / M ∗ and EW(H 𝛼 ) which traces SFR over ∼ M H i / M ∗ is a strong influence on long-term averagedstar formation. by increased shock heating in more gas-rich galaxies. DIG line ratiosmay themselves be driven by shocks, but this is a subject of debate(Haffner et al. 2009).If our observations can be explained by more DIG in gas-rich galaxies, one might also wonder if this implies the presence ofmore diffuse H i as well. Analyses of the neutral content of theDIG implies that it is indeed mostly ionized, with neutral fractions <
10% (Reynolds et al. 1998; Hausen et al. 2002). However, kine-matic studies find that DIG is almost always found to coincide in bothspace and velocity with H i , specifically the warm neutral medium(WNM; the warmer, diffuse, more extended neutral hydrogen) asopposed to the cold neutral medium (CNM; the cooler, denser hy-drogen clouds) (Spitzer & Fitzpatrick 1993; Reynolds et al. 1995;Hartmann & Burton 1997; Howk et al. 2003; Haffner et al. 2009). Therefore, it is possible that galaxies with more DIG may also have alarger fraction of their HI residing in the WNM.
A further implicationof this relationship would be that the most gas-rich galaxies (whichshow evidence of a higher DIG fraction) would host a larger WNMfraction. Similarly, the scatter in the M
H i / M ∗ -EW(O) relation, whichwe attributed to variations in the fraction of DIG emission, may alsoreflect variations in the amount of H i in the CNM vs. WNM phases.Whether an H i reservoir has a large diffuse component or is subjectto more shock heating, such a change in the ISM may impact globalstar formation efficiencies. The H i will not be able to condense intomolecular clouds nearly as easily as it would otherwise, increasing theoverall H i depletion time, 𝑡 dep , defined as 𝑀 H i / 𝑆𝐹 𝑅 . We explore thispossibility by plotting 𝑀 H i vs. SFR in Figure 11 and comparingto lines of constant 𝑡 dep . We specifically use SFR in order tosmooth over stochastic variations in the SFR and look at the long-term averaged behavior. Points in Fig 11 are color-coded by [O i ]/H 𝛼 (left) and 𝜇 𝐻 𝛼 (right), and we restrict our analysis to the subset ofgalaxies with 12 + log 𝑂 / 𝐻 < .
05 where there is no metallicity
MNRAS , 1–23 (2020)
I-MaNGA dependence on [O i ]/H 𝛼 . Longer H i depletion times at fixed SFRcoincide with enhanced [O i ]/H 𝛼 and depressed 𝜇 H 𝛼 , supporting theidea that longer H i depletion times are caused by a larger fraction ofH i residing in the WNM and/or heated by shocks. Using the second data release of the HI-MaNGA survey, an H i follow-up program for the SDSS-IV MaNGA survey, we conductedan analysis of the scaling relations between galaxy H i -to-stellar massratio and average ISM properties derived from optical spectroscopyfor star forming galaxies. Our key results are as follows: • We recover the strong anti-correlation between M
H i / M ∗ andgas-phase metallicity seen in previous studies. We also find a mildcorrelation between M H i / M ∗ and [O i ]/H 𝛼 , which suggests the pres-ence of a harder ionizing field, more shock excitation, and/or a largerfraction of diffuse ionized gas (DIG) in gas-rich galaxies. A largerDIG fraction may in turn imply a larger fraction of H i residing in themore diffuse warm neutral medium (WNM). • M H i / M ∗ is weakly correlated with EW( 𝐻𝛼 ) and other opticalhydrogen lines. This weak connection implies the global H i reser-voir does not have a strong impact on the star formation on short-timescales. However, the stronger link between H i and SFR whenmeasured over longer timescales implies the H i reservoir still playsan important role in sustaining the average SFRs of galaxies. • Of all optical emission lines, M
H i / M ∗ correlates most stronglywith Oxygen lines ([O i ], [O ii ], and [O iii ]). This result is likely drivenby the existing anti-correlation between gas fraction and metallicity. • The residuals in the M
H i / M ∗ -EW(O) relations are most stronglycorrelated with [O i ]/H 𝛼 and mean H 𝛼 surface brightness, 𝐻𝛼 . Thisresult suggests the scatter is driven by variations in the amount of gasassociated with the warm neutral medium/DIG. • Galaxies with longer than average depletion times also haveelevated [O i ]/H 𝛼 and depressed 𝜇 𝐻 𝛼 , consistent with long depletiontimes occurring when a significant fraction of H i is in a diffuse phaseand/or subject to more shock heating, making it more difficult tocondense into molecular clouds. • We make a first attempt to calibrate multi-parameter M
H i / M ∗ scaling relations using optical spectroscopic information, but findthat the asymmetric scatter makes them less preferable than existingscaling relations using broadband photometric properties. However,the relations in this work appear to strengthen, with notably lessasymmetric scatter, if emission line properties are calculated usinglarger apertures. More data with larger area coverage – or carefulaperture corrections applied to our current MaNGA sample – mayprove valuable when revisiting these scaling relations.Our results highlight how ISM properties vary from typically low-mass gas-rich galaxies, to typically more massive gas-poor galaxies.In many ways, gas-rich dwarf galaxies represent aspects of the pro-genitors of the Milky Way and other more massive galaxies. Under-standing the properties of their ISM, and how they might impact starformation, is a key ingredient to our general understanding of galaxyformation. Analyses of the relative fractions of diffuse and denseH i may be particularly useful for simulations which aim to properlyrecreate the true breakdown of the ISM into its different phases.In a practical sense, this work also demonstrates ways to con-duct statistical analysis in the presence of H i upper limits. Boththe Akritas-Theil-Sen estimator and substituting upper limits withphotometric gas fractions are useful approaches and give generallyconsistent results. A third approach, which we have not demonstrated here, is stacking of H i data, with the caveat that stacking removesinformation about the distribution of data around mean trends. Re-gardless, H i non-detections contain valuable information, it is crucialto incorporate them into any analysis of gas content in order to avoidbiased results. ACKNOWLEDGEMENTS
DATA AVAILABILITY
The HI-MaNGA catalog used in this study can be found at https://greenbankobservatory.org/science/gbt-surveys/hi-manga/ .The MaNGA MPL-9 data products can be gener-ated by the public using the raw data (available at with DRP v2.7.1, DAP v2.4.3, and Pipe3D v2.7.1.
REFERENCES
Aguado D. S., et al., 2019, ApJS, 240, 23Akritas M. G., Siebert J., 1996, MNRAS, 278, 919Akritas M. G., Murphy S. A., LaValley M. P., 1995, J. Amer. Stat. Assn, 90,170 MNRAS , 1–23 (2020) D. V. Stark et al. −1.5 −1.0 −0.5 0.0 0.5 1.0 log SFR [M ⊙ yr −1 ] l o g M H I [ M ⊙ ] . G y r . G y r G y r G y r G y r −1.50−1.45−1.40−1.35−1.30−1.25−1.20−1.15 log [OI]⊙Hα −1.5 −1.0 −0.5 0.0 0.5 1.0 log SFR [M ⊙ yr −1 ] l o g M H I [ M ⊙ ] . G y r . G y r G y r G y r G y r μ Hα [erg s −1 kpμ −⊙ ] Figure 11. 𝑀 H i versus SFR for galaxies with + log 𝑂 / 𝐻 < . . Longer depletion times coincide with higher [O i ]/H 𝛼 and depressed 𝜇 𝐻 𝛼 , consistent witha larger fraction of DIG/WNM and/or gas heated by shocks.Belfiore F., et al., 2016, MNRAS, 461, 3111Belfiore F., et al., 2019, AJ, 158, 160Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B.,Thornley M. D., 2008, AJ, 136, 2846Blanton M. R., Roweis S., 2007, AJ, 133, 734Blanton M. R., et al., 2017, AJ, 154, 28Bothwell M. S., Maiolino R., Kennicutt R., Cresci G., Mannucci F., MarconiA., Cicone C., 2013, MNRAS, 433, 1425Braine J., Combes F., 1993, A&A, 269, 7Brinchmann J., Charlot S., Kauffmann G., Heckman T., White S. D. M.,Tremonti C., 2013, MNRAS, 432, 2112Broeils A. H., Rhee M. H., 1997, A&A, 324, 877Broeils A. H., van Woerden H., 1994, A&AS, 107, 129Brown T., Cortese L., Catinella B., Kilborn V., 2018, MNRAS, 473, 1868Bundy K., et al., 2015, ApJ, 798, 7Calette A. R., Avila-Reese V., Rodríguez-Puebla A., Hernández-Toledo H.,Papastergis E., 2018, Rev. Mex. Astron. Astrofis., 54, 443Cano-Díaz M., Ávila-Reese V., Sánchez S. F., Hernández-Toledo H. M.,Rodríguez-Puebla A., Boquien M., Ibarra-Medel H., 2019, MNRAS,488, 3929Catinella B., et al., 2013, MNRAS, 436, 34Catinella B., et al., 2018, MNRAS, 476, 875Cayatte V., Kotanyi C., Balkowski C., van Gorkom J. H., 1994, AJ, 107, 1003Cherinka B., et al., 2019, AJ, 158, 74Cid Fernandes R., Stasińska G., Mateus A., Vale Asari N., 2011, MNRAS,413, 1687Croom S. M., et al., 2012, MNRAS, 421, 872Dalcanton J. J., Yoachim P., Bernstein R. A., 2004, ApJ, 608, 189Dopita M. A., Kewley L. J., Heisler C. A., Sutherland R. S., 2000, ApJ,542, 224Doyle M. T., Drinkwater M. J., 2006, MNRAS, 372, 977Drory N., et al., 2015, AJ, 149, 77Eckert K. D., Kannappan S. J., Stark D. V., Moffett A. J., Norris M. A., SnyderE. M., Hoversten E. A., 2015, ApJ, 810, 166Elmegreen B. G., 1993, ApJ, 411, 170Feigelson E. D., Babu G. J., 2012, Modern Statistical Methods for AstronomyFeigelson E. D., Nelson P. I., 1985, ApJ, 293, 192Giovanelli R., Haynes M. P., Herter T., Vogt N. P., Wegner G., Salzer J. J., daCosta L. N., Freudling W., 1997, AJ, 113, 22Goddy J., Stark D. V., Masters K. L., 2020,Research Notes of the American Astronomical Society, 4, 3Gunn J. E., et al., 2006, AJ, 131, 2332Haffner L. M., Reynolds R. J., Tufte S. L., 1999, ApJ, 523, 223Haffner L. M., et al., 2009, Reviews of Modern Physics, 81, 969Hartmann D., Burton W. B., 1997, Atlas of Galactic Neutral Hydrogen Hausen N. R., Reynolds R. J., Haffner L. M., 2002, AJ, 124, 3336Haynes M. P., et al., 2011, AJ, 142, 170Haynes M. P., et al., 2018, ApJ, 861, 49Hernández-Toledo H. M., Vázquez-Mata J. A., Martínez-Vázquez L. A., AvilaReese V., Méndez-Hernández H., Ortega-Esbrí S., Núñez J. P. M., 2008,AJ, 136, 2115Hoopes C. G., Walterbos R. A. M., 2003, ApJ, 586, 902Howk J. C., Sembach K. R., Savage B. D., 2003, ApJ, 586, 249Huang S., Haynes M. P., Giovanelli R., Brinchmann J., 2012, ApJ, 756, 113Hughes T. M., Cortese L., Boselli A., Gavazzi G., Davies J. I., 2013, A&A,550, A115Hunter D. A., Elmegreen B. G., Baker A. L., 1998, ApJ, 493, 595Isobe T., Feigelson E. D., Nelson P. I., 1986, ApJ, 306, 490Jaskot A. E., Oey M. S., Salzer J. J., Van Sistine A., Bell E. F., Haynes M. P.,2015, ApJ, 808, 66Kannappan S. J., 2004, ApJ, 611, L89Kannappan S. J., Fabricant D. G., Franx M., 2002, AJ, 123, 2358Kannappan S. J., et al., 2013, ApJKaplan E. L., Meier P., 1958, Journal of the American Statistical Association,53, 457Kenney J. D. P., Young J. S., 1989, ApJ, 344, 171Kennicutt Robert C. J., 1998, ARA&A, 36, 189Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531Kewley L. J., Dopita M. A., 2002, ApJS, 142, 35Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183Kewley L. J., Geller M. J., Barton E. J., 2006a, AJ, 131, 2004Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006b, MNRAS,372, 961Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511Krumholz M. R., McKee C. F., Tumlinson J., 2009a, ApJ, 693, 216Krumholz M. R., McKee C. F., Tumlinson J., 2009b, ApJ, 699, 850Law D. R., et al., 2015, AJ, 150, 19Law D. R., et al., 2016, AJ, 152, 83Lequeux J., Peimbert M., Rayo J. F., Serrano A., Torres-Peimbert S., 1979,A&A, 500, 145Leroy A. K., et al., 2012, AJ, 144, 3Lisenfeld U., et al., 2011, A&A, 534, A102Madsen G. J., Reynolds R. J., Haffner L. M., 2006, ApJ, 652, 401Masters K. L., et al., 2010, MNRAS, 405, 783Masters K. L., et al., 2019, MNRAS, 488, 3396McGaugh S. S., 1991, ApJ, 380, 140Meurer G. R., Carignan C., Beaulieu S. F., Freeman K. C., 1996, AJ, 111, 1551Moran S. M., et al., 2012, ApJ, 745, 66O’Donnell J. E., 1994, ApJ, 422, 158Oey M. S., et al., 2007, ApJ, 661, 801MNRAS , 1–23 (2020)
I-MaNGA Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae andactive galactic nucleiPeeples M. S., Pogge R. W., Stanek K. Z., 2008, ApJ, 685, 904Pérez-Montero E., Hägele G. F., Contini T., Díaz Á. I., 2007, MNRAS,381, 125Reynolds R. J., 1985a, ApJ, 294, 256Reynolds R. J., 1985b, ApJ, 298, L27Reynolds R. J., Tufte S. L., Kung D. T., McCullough P. R., Heiles C., 1995,ApJ, 448, 715Reynolds R. J., Hausen N. R., Tufte S. L., Haffner L. M., 1998, ApJ, 494, L99Robertson P., Shields G. A., Blanc G. A., 2012, ApJ, 748, 48Rodríguez-Puebla A., Calette A. R., Avila-Reese V., Rodriguez-Gomez V.,Huertas-Company M., 2020, Publ. Astron. Soc. Australia, 37, e024Sakamoto K., Okumura S. K., Ishizuki S., Scoville N. Z., 1999, ApJ, 525, 691Sánchez S. F., et al., 2012, A&A, 538, A8Sánchez S. F., et al., 2014, A&A, 563, A49Sánchez S. F., et al., 2016a, Rev. Mex. Astron. Astrofis., 52, 21Sánchez S. F., et al., 2016b, Rev. Mex. Astron. Astrofis., 52, 171Schmitt J. H. M. M., 1985, ApJ, 293, 178Schruba A., et al., 2011, AJ, 142, 37Sheth K., Vogel S. N., Regan M. W., Thornley M. D., Teuben P. J., 2005, ApJ,632, 217Skillman E. D., Kennicutt Robert C. J., Shields G. A., Zaritsky D., 1996, ApJ,462, 147Smee S. A., et al., 2013, AJ, 146, 32Spitzer Lyman J., Fitzpatrick E. L., 1993, ApJ, 409, 299Springob C. M., Haynes M. P., Giovanelli R., Kent B. R., 2005, ApJS,160, 149Stark D. V., Kannappan S. J., Wei L. H., Baker A. J., Leroy A. K., EckertK. D., Vogel S. N., 2013, ApJ, 769, 82Stark D. V., et al., 2016, ApJ, 832, 126Stasińska G., et al., 2008, MNRAS, 391, L29Swaters R. A., van Albada T. S., van der Hulst J. M., Sancisi R., 2002, A&A,390, 829Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295Voges E. S., Walterbos R. A. M., 2006, ApJ, 644, L29Wake D. A., et al., 2017, preprint, ( arXiv:1707.02989 )Wang J., et al., 2011, MNRAS, 412, 1081Westfall K. B., et al., 2019, AJ, 158, 231Wilcox R., 2010, Fundamentals of Modern Statistical Methods: Sub-stantially Improving Power and Accuracy. Springer New York, https://books.google.com/books?id=uUNGzhdxk0kC
Yan R., et al., 2016a, AJ, 151, 8Yan R., et al., 2016b, AJ, 152, 197Yesuf H. M., Ho L. C., 2019, ApJ, 884, 177York D. G., et al., 2000, AJ, 120, 1579Zhang K., et al., 2017, MNRAS, 466, 3217Zu Y., 2020, MNRAS, 496, 111van Zee L., Haynes M. P., Salzer J. J., Broeils A. H., 1997, AJ, 113, 1618
APPENDIX A: PERFORMANCE OF THEAKRITAS-THEIL-SEN ESTIMATOR
To test the performance of the ATS estimator, we examine its con-sistency and accuracy when applied to a mock data set where thefraction of censored data is progressively increased. For this analysiswe use data from the xGASS survey (Catinella et al. 2018), chosenbecause it has deep H i observations and a roughly uniform stel-lar mass distribution similar to our MaNGA data set. We start withonly the xGASS detections, so that we have a sample where true H i masses are known for all galaxies. Fig. A1 shows the M H i / M ∗ vs. 𝑔 − 𝑟 relation for this sample, as well as the ordinary unweightedleast-squares (OLS) fit and the ATS fit. These two fitting algorithmsare in good agreement within their errors. g − r −2.0−1.5−1.0−0.50.00.51.0 l o g M H I / M * OLS fitAk itas-Theil-Sen (ATS) fit
Figure A1.
The M H i / M ∗ vs. 𝑔 − 𝑟 relation using only H i detections fromthe xGASS survey (Catinella et al. 2018). The red and blue lines representthe OLS and ATS fits to the data. The M
H i / M ∗ vs. 𝑔 − 𝑟 relation is a ideal test of the ATS esti-mator for this work because when we simulate different observingdepths and introduce mock upper limits, the distribution of upperlimits is similar to the distribution seen in the relations presentedin Section 3. Specifically, the range of M H i / M ∗ for detections andnon-detections overlaps significantly at fixed 𝑥 -axis value. To createmock non-detections in this data set in a manner consistent with thenon-detections in the HI-MaNGA survey, we assign each xGASSgalaxy a new distance based on its 𝑖 -band absolute magnitude, 𝑀 𝑖 ,mimicking the relationship between 𝑀 𝑖 and redshift which is explic-itly built into the MaNGA Primary sample design (Wake et al. 2017).For a given 𝑀 𝑖 , a galaxy in MaNGA will only fall within a set rangeof redshifts, and we assign each galaxy a random value within thatrange. To ensure consistency with MaNGA data, we use the ellipicalPetrosian magnitudes from the NSA catalog. Each galaxy’s observedH i flux is recalculated using its known H i mass and this new dis-tance. We then assume all galaxies are observed down to a fixed 𝑟𝑚𝑠 noise level, which we use to calculate the integrated flux 𝑆 / 𝑁 foreach galaxy. Anything with a mock 𝑆 / 𝑁 < 𝜎 upper limit assuming a linewidth of 200km s − . We start with an rms of 1.5 mJy (the nominal HI-MaNGAsurvey depth) and progressively increase it up to 10 mJy, each timerunning an OLS fit on the mock detections and an ATS fit on all thedata, including the upper limits.Figure A2 shows the results of this analysis. Unsurprisingly, theOLS fit on just the detections becomes increasingly biased as themock survey depth decreases. Meanwhile, the ATS estimator fitchanges very little, until the fraction of censored data exceeds ∼ ∼ MNRAS000
H i / M ∗ vs. 𝑔 − 𝑟 relation is a ideal test of the ATS esti-mator for this work because when we simulate different observingdepths and introduce mock upper limits, the distribution of upperlimits is similar to the distribution seen in the relations presentedin Section 3. Specifically, the range of M H i / M ∗ for detections andnon-detections overlaps significantly at fixed 𝑥 -axis value. To createmock non-detections in this data set in a manner consistent with thenon-detections in the HI-MaNGA survey, we assign each xGASSgalaxy a new distance based on its 𝑖 -band absolute magnitude, 𝑀 𝑖 ,mimicking the relationship between 𝑀 𝑖 and redshift which is explic-itly built into the MaNGA Primary sample design (Wake et al. 2017).For a given 𝑀 𝑖 , a galaxy in MaNGA will only fall within a set rangeof redshifts, and we assign each galaxy a random value within thatrange. To ensure consistency with MaNGA data, we use the ellipicalPetrosian magnitudes from the NSA catalog. Each galaxy’s observedH i flux is recalculated using its known H i mass and this new dis-tance. We then assume all galaxies are observed down to a fixed 𝑟𝑚𝑠 noise level, which we use to calculate the integrated flux 𝑆 / 𝑁 foreach galaxy. Anything with a mock 𝑆 / 𝑁 < 𝜎 upper limit assuming a linewidth of 200km s − . We start with an rms of 1.5 mJy (the nominal HI-MaNGAsurvey depth) and progressively increase it up to 10 mJy, each timerunning an OLS fit on the mock detections and an ATS fit on all thedata, including the upper limits.Figure A2 shows the results of this analysis. Unsurprisingly, theOLS fit on just the detections becomes increasingly biased as themock survey depth decreases. Meanwhile, the ATS estimator fitchanges very little, until the fraction of censored data exceeds ∼ ∼ MNRAS000 , 1–23 (2020) D. V. Stark et al. g − r l o g M H I / M * −2.0−1.5−1.0−0.50.00.51.0 Censored fraction: 0.07 Censored fraction: 0.20
Censored fraction: 0.32
Censored fraction: 0.59 original OLS fitOLS fit to detectionsoriginal ATS fitATS fitdetectionsupper li its
Figure A2. M H i / M ∗ vs. 𝑔 − 𝑟 color relation for xGASS data where galaxies are redistributed in redshift space to match the MaNGA survey design, then mockobserved to progressively lower depths. The thick lines represent the original fits to the uncensored data from Figure A1. The Akritas-Theil-Sen line remainsextremely stable until the censoring fraction exceeds ∼ APPENDIX B: PHOTOMETRIC GAS FRACTIONCALIBRATION
To determine 𝑃 ( M H i / M ∗ | color ) , the 2D probability distribution ofgalaxies as a function of M H i / M ∗ and modified color, we followthe methodology of Eckert et al. (2015) using the RESOLVE surveyH i catalog (Stark et al. 2016). RESOLVE is the ideal data set todetermine this probability distribution as it is a highly complete,closed volume survey with H i data that is uniformly deep (as afraction of stellar mass). We select only data from the RESOLVE-A volume above the baryonic mass completeness limit, enforce alldetections to have 𝑆 / 𝑁 >
5, and reject any confused targets whosesystematic error due to confusion is >
25% of the integrated flux.The RESOLVE catalog is crossmatched with the NSA using a searchradius of 5 ′′ in order for us to use photometry and stellar massesconsistent with our MaNGA sample.We first determine which color provides the tightest correlationwith M H i / M ∗ when using NSA data, which we find to be 𝑢 − 𝑖 . Wenext determine the optimum third parameter, testing both axial ratio( 𝑏 / 𝑎 ) and 𝑟 -band surface brightness calculated within 𝑅 𝑒 ( 𝜇 𝑟 ). Thefits are limited to the H i detections between 0 . < 𝑢 − 𝑖 < 𝑏 / 𝑎 to be the optimal third parameter, wefind it has very little correlation with the residuals in the M H i / M ∗ - ( 𝑢 − 𝑖 ) correlation. A more significant correlation is found with 𝜇 𝑟 .Based on the fit between M H i / M ∗ , 𝑢 − 𝑖 , and 𝜇 𝑟 , our final modifiedcolor is defined as 0 . ( 𝑢 − 𝑖 ) − . 𝜇 𝑟 +
4. The additive factor of
Table B1.
Best fit values to the 𝑃 ( 𝑀 H i | color ) model 𝐴 𝐴 𝐴 𝐴 -1.10 𝐴 𝐴 𝐴 𝐴 𝐴 >
0, which is neededfor the full fit to 𝑃 ( M H i / M ∗ | color ) . 𝑃 ( M H i / M ∗ | color ) is fit using Eqs. 1–4 from Eckert et al. (2015).To briefly summarize, the model assumes there are two distinctpopulations: (1) H i -rich galaxies (typically detections) which fol-low a linear relation with modified color and have Gaussian scatterwhich broadens at redder colors, and (2) gas-poor (typically non-detections) which are assumed to all fall around M H i / M ∗ ∼ .
05. Weuse M
H i / M ∗ and color bin sizes of 0.2 and weight the data by 1 / 𝑁 ,where 𝑁 is the number of galaxies per bin. The final fit parametersare given in Table B1.For each non-detection in our MaNGA sample, we determine itsM H i / M ∗ probability distribution based on its modified color and themodel fit parameters given in Table B1. We set all probabilities tozero above the measured upper limit. The new estimate of M H i / M ∗ for this galaxy is determined by randomly drawing from its M H i / M ∗ probability distribution. MNRAS , 1–23 (2020)
I-MaNGA This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000