THE FUELING DIAGRAM: LINKING GALAXY MOLECULAR-TO-ATOMIC GAS RATIOS TO INTERACTIONS AND ACCRETION

Our primary sample comes from the NFGS (Jansen et al. 2000a, 2000b; see also Jansen & Kannappan 2001), a set of ∼200 galaxies selected to span a broad range of B-band luminosities and morphologies. The data products of the original survey include UBR surface photometry and optical spectroscopy (Jansen et al. 2000a, 2000b; Kannappan & Fabricant 2001; Kannappan et al. 2002, 2009). The sample also includes extensive supporting data relevant to this study. Roughly 90% of galaxies have Sloan Digital Sky Survey (SDSS) DR8 optical imaging (Aihara et al. 2011). All have near-infrared (NIR) data from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) while 55% have deeper Spitzer Infrared Array Camera (IRAC; Fazio et al. 2004) imaging, mainly from Sheth et al. (2010), Moffett et al. (2012) and Kannappan et al. (2013, hereafter K13). In addition, all galaxies have single dish 21 cm data (Wei et al. 2010a, K13).

For this study, we analyze the 35 out of 39 NFGS galaxies with CO data that pass the usability cuts applied to our final sample (see Section 2.3.7). This subset of the NFGS has a stellar mass range of 10^8.8–10^10.5 M_☉, morphologies ranging from E to Sbc, and diverse states of star formation (e.g., starbursting, post-starburst, quiescent). Most of the CO data for this sample came from new observations, and unlike many previous investigations of molecular gas in galaxies, our NFGS CO observations were not in general designed to emphasize CO-bright galaxies, but were instead designed to reach strong, scientifically useful upper limits in the case of CO non-detections.

To strengthen our results, we supplement our sample with galaxies from the literature with available CO, H i, and multi-band imaging data.

Our literature sample includes galaxies from three large surveys: the Spitzer Near Infrared Galaxy Survey (SINGS; Kennicutt et al. 2003), ATLAS-3D (Young et al. 2011), and COLD GASS (Saintonge et al. 2011). These surveys are dominated by high mass galaxies, so to supplement the low mass end of our data set, we add galaxies from Barone et al. (2000), Garland et al. (2005), Leroy et al. (2005),⁶ Taylor et al. (1998), and Kannappan et al. (2009). Most of these references are themselves the sources of the CO data, although some CO data for galaxies in SINGS and Kannappan et al. (2009) come from Leroy et al. (2009), Albrecht et al. (2007), Sage et al. (2007), Zhu et al. (1999), or our own observations (Section 2.2). H i data often come from the same source as the CO data, or else from alternate sources in the literature or HyperLeda (Huchtmeier et al. 1995; Smoker et al. 2000; Salzer et al. 2002; Paturel et al. 2003; Garland et al. 2004; Springob et al. 2005; Meurer et al. 2006; Walter et al. 2008; Catinella et al. 2010; Haynes et al. 2011; Catinella et al. 2012; Serra et al. 2012). All optical data come from the SDSS DR8, except for a subset of the SINGS sample outside the SDSS footprint that has BVRI photometry (Muñoz-Mateos et al. 2009). NIR imaging data is available for all galaxies from 2MASS.

Combined, these data sources bring our full sample (only considering galaxies that have all necessary data) to 627 galaxies. However, our total decreases to 323 galaxies after we institute a number of usability cuts (see Section 2.3.7).

Optical/NIR photometry is needed to estimate stellar masses and track recent central star formation enhancements. For our analysis, we do not use products of the SDSS pipeline since it is prone to shredding extended sources and does not make use of the most recent background subtraction algorithm (Blanton et al. 2011). We instead recalculate total magnitudes using our own custom pipeline, described in greater detail by K. D. Eckert et al. (in preparation). After the downloaded images are co-added, bright stars and interloping galaxies are masked. The masking process is automated with the aid of SExtractor (Bertin & Arnouts 1996), but each mask is inspected by eye and adjusted when there is clear over- or under-masking using the automatic routine. The ELLIPSE task in IRAF is used to extract surface brightness profiles of constant center, position angle (P.A.), and ellipticity, and to sum up the flux within each isophote. For the NFGS sample, we adopt the same P.A.s and ellipticities used for the UBR photometry (Jansen et al. 2000b). For our literature sample, we adopt the method of K. D. Eckert et al. (in preparation), who use ELLIPSE to determine the best P.A. and ellipticity from the low surface brightness outer disk of each galaxy using the co-added gri images. Models of each galaxy are created with the resulting surface brightness profiles, and are used to fill in masked regions to correct the total flux.

Total magnitudes are calculated two ways. First, we adopt a curve of growth technique very similar to the one outlined in Muñoz-Mateos et al. (2009), where the outer disk values of the enclosed magnitude and its radial gradient are fit with a linear function. The y-intercept of this line (i.e., where the enclosed magnitude is no longer increasing) is the total magnitude. Total magnitudes are also calculated by fitting the outer disk of the surface brightness profile with an exponential function, similar to the method in Jansen et al. (2000b), except that outlier points are rejected from the fit. The total flux is summed up to the last isophote used in the fit, after which the fit itself is used to estimate the remaining outer flux.

To estimate systematic uncertainties in our total magnitudes, we perform each of these total magnitude extrapolations using slightly differently defined fit ranges (between 1 and 8 times the sky noise, between 3 and 10 times the sky noise, within a 1 mag arcsec⁻² range ending at 1 or 3 times the sky noise, and finally using the last 5 data points above 1.5 times the sky noise) and then average the results (ignoring >6σ outliers) to obtain our final magnitude. We take the difference of the maximum and minimum magnitude estimates divided by two (also ignoring outliers) as our official systematic uncertainty. This uncertainty is added in quadrature with the Poisson statistical uncertainty.

The ellipses used for total magnitude calculation best match the shape of the far outer disk or halo, and are therefore not ideal for calculating galaxy inclinations. To estimate photometric inclinations, ELLIPSE is run a second time where ellipticity is allowed to vary while P.A. is kept fixed to the previously used value. The ellipticity at a surface brightness of 22.5 mag arcsec⁻² in the co-added gri images usually reliably traces the higher surface brightness inner disk and provides accurate inclinations, which are estimated using

$$\begin{equation} \cos {i}=\sqrt{\frac{(b/a)^2-q^2}{1-q^2}}. \end{equation} \tag{ 3 }$$

Here, q is the intrinsic disk thickness (assumed to be q = 0.2), and b/a is the minor-to-major axis ratio derived from the ellipticity. For the NFGS sample, we use the same inclinations as K13 for consistency, which are based on the same equation, but not using the b/a measured from SDSS photometry.

JHK magnitudes are also recalculated using the 2MASS imaging data. Here, we redo the background subtraction for the 2MASS images using a method similar to that used in the original 2MASS pipeline (Jarrett et al. 2000), where the sky was fit by 3rd order polynomials. However, we fit the background of each relevant frame only in the region local to the galaxy of interest (within ∼5 × R₂₅), with the galaxy itself and any stars or background galaxies masked. We impose the parameters from our first set of optical surface brightness profiles (center, P.A., and ellipticity) to extract NIR surface brightness profiles.

The two methods for calculating total optical magnitudes described above are again used to calculate the total NIR magnitudes. However, tests comparing the output of these two methods against magnitudes calculated from much deeper Spitzer IRAC 3.6 μm imaging reveal that the exponential fit method performs systematically worse than the curve of growth method when applied to the relatively shallow 2MASS data. The fits also perform best when using a fit range between 3 and 10 times the sky noise. We thus solely use the curve-of-growth derived total magnitudes determined using this fit range as our final estimates, although we still rely on the difference between the curve-of-growth and exponential fit derived magnitudes to estimate the systematic error for most galaxies. For any galaxy for which our different magnitude estimation techniques yield very different results (disagreement of more than 0.5 mag), we resort to taking an aperture magnitude of the galaxy, and then infer the total magnitude by multiplying by the total-to-aperture flux ratio of a higher S/N passband (another NIR band like J if possible, otherwise the i band).

All optical and NIR magnitudes are corrected for foreground extinction using the dust maps of Schlegel et al. (1998) and the extinction curve of O'Donnell (1994).

Stellar masses are calculated with an improved version of the method described by Kannappan et al. (2009). Taking as inputs a combination of UBR, ugriz, JHK, and 3.6 μm photometry, along with global optical spectra (or any subset of these inputs that are available; see K13 and references therein), the stellar mass estimation code fits mixed young + old stellar populations built from pairs of simple stellar population (SSP) models from Bruzual & Charlot (2003) assuming a Salpeter initial mass function (IMF). Output stellar masses are scaled by 0.7 to match a "diet" Salpeter IMF containing fewer low mass stars (Bell et al. 2003). The only significant change compared to Kannappan et al. (2009) is the addition of a very young SSP with an age of 5 Myr to the suite of models. Note also that model SSP pairs can include a "middle-aged" young SSP, as long as its age is younger than the old SSP, which was also true for the Kannappan et al. (2009) model grid. In addition, the UBR zero points have been adjusted for consistency (see K13). The final stellar mass estimate is defined by the median and 68% confidence interval of the likelihood weighted mass distribution over the full model grid.

Stellar masses for our literature sample are calculated using only SDSS and 2MASS photometry. To determine whether the lack of spectroscopy leads to any systematic differences in our stellar mass estimates, we compare the high quality stellar masses from the NFGS against a second set calculated using only our custom SDSS and 2MASS photometry. As shown in Figure 1, the two methods of estimation are in good agreement. There is no statistical significance in the weak linear trend between the two mass estimates as a function of stellar mass, and the 1σ scatter between the two estimates is 0.04 dex, much less than the typical uncertainty of ∼0.2 dex for stellar masses.

Zoom In Zoom Out Reset image size

To track enhancements in recent central star formation, we rely on a simple measure of the color gradient referred to as blue-centeredness (ΔC), defined as the outer disk color from the half-light radius (r₅₀) to the 75% light radius (r₇₅) minus the color from the center to r₅₀ (Jansen et al. 2000b).¹⁰ For future reference, general discussion of blue-centeredness will use the notation, ΔC, while specific discussion involving a particular color will note that color explicitly, e.g., Δ(g − r). The radii in this definition reliably separate bulge/disk colors without being sensitive to variations in the bulge-to-disk ratio (Kannappan et al. 2004). The half-light radius is a natural separator of inner/outer galaxy growth; for example, Peletier & Balcells (1996) note shifts in colors and ages at approximately this radius. We also stress that this simple measure appears robust to variations in the dividing radii. For example, using the 40% and 90% enclosed light radii returns approximately the same results.

Before blue-centeredness is used to measure enhancements in recent central star formation, we account for other galaxy properties that may affect the color gradients. Kannappan et al. (2004) find that blue-centeredness correlates with galaxy luminosity, likely due to the fact that both metallicity and stellar population age can influence galaxy color gradients, reddening galaxy centers with increasing luminosity. These authors remove the luminosity trend by subtracting off the fitted relation between Δ(B − R) and luminosity. We follow the same approach, performing an ordinary least-squares fit (Isobe et al. 1990) of ΔC against stellar mass (Figure 2). We derive the following relations:

$$\begin{eqnarray} &&\Delta (u-r)=\left(-0.114\pm 0.025\right)\log {M_*}+\left(0.983\pm 0.249\right) \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\Delta (u-g)=\left(-0.064\pm 0.017\right)\log {M_*}+\left(0.556\pm 0.167\right) \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\Delta (g-r)=\left(-0.049\pm 0.009\right)\log {M_*}+\left(0.417\pm 0.095\right). \end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

The use of an ordinary least-squares fit of y versus x is crucial, since in order to remove any trend with stellar mass, we must minimize the scatter in ΔC relative to stellar mass. Following Kannappan et al. (2004), we calibrate this mass correction using the NFGS sample (not just our subset with CO data) since it is our most representative data set and has the most robust stellar mass estimates. The fit is limited to star-forming disk galaxies, identified as galaxies with morphology S0 or later and either detected Hα emission extending beyond the nucleus or U − K < 4, which helps exclude quenched galaxies (Kannappan & Wei 2008; K13). We reject galaxies with known strong active galactic nuclei (AGNs) and restrict the fit to M_*  >  10^8.5 M_☉, avoiding the low M_* tail of the NFGS where all blue-centeredness states may not be evenly represented. In total, 135 NFGS galaxies were used to derive these relations. The residuals of this fit give mass-corrected blue-centeredness (denoted by ΔC^m). In simple terms, ΔC^m is a measure of the color difference between the inner and outer regions of a galaxy, relative to the typical color difference for galaxies at a given stellar mass, where higher values imply bluer centers. More physically, this parameter tracks recent central star formation relative to the typical star-forming galaxy at a given stellar mass.

We find that galaxies with and without strong peculiarities (i.e., signs of recent interactions) as defined by Kannappan et al. (2004) have different underlying distributions of Δ(g − r)^m at a 99% confidence level, so our definition of mass-corrected blue-centeredness successfully recreates the correlation found by these authors using luminosity-corrected blue-centeredness (their Figure 8). However, we stress that while our mass correction is physically motivated, it is modest, and the distributions of uncorrected Δ(g − r) for peculiar/non-peculiar galaxies differ at a confidence level comparable to that for distributions using Δ(g − r)^m. Our general results presented in Section 3 are still found without the correction in place, as discussed in Section 3.2.2.

Within this paper, the total error budget on ΔC^m includes contributions from the measurement error of ΔC, as well as additional uncertainties due to the error in r₅₀ and r₇₅, the error in the stellar mass, and the uncertainty in the fitted relation between ΔC and stellar mass. In the majority of cases, the Poisson error is dominant.

We note that optical color gradients have advantages over direct star formation rate tracers (e.g., FUV + 24 μm) for the purposes of our study. A practical advantage is that the data needed to compute color gradients are readily available for most galaxies. In addition, the longer timescales over which optical indicators remain sensitive to enhanced star formation are more useful for studying the extended evolution of galaxies. They enable us to analyze the stellar populations beyond the timescale of the star formation and gas consumption itself, allowing us to note the longer term consequences of recent star formation episodes. We defer a detailed analysis of the timescales of these optical indicators to future work, but Figure 11 of Kannappan et al. (2004) shows a simple model wherein the blue-centeredness fades on the order of 0.5–2 Gyr after the central star formation episode ends. It should be noted that this model represents a single case, and the evolution of blue-centeredness may vary depending on the duration and size of the burst as well as the composition of the preexisting stellar populations in the bulge and disk.

For consistency and since the vast majority of our galaxies have SDSS data, we always quote ΔC^m using SDSS-equivalent colors. A handful of galaxies from our NFGS sample do not have SDSS data, but do have UBR photometry. We derive the following conversions between the NFGS UBR and ugriz systems using colors measured within the B-band 25 mag arcsec⁻² isophote for 183 NFGS galaxies:

$$\begin{eqnarray} &&g-r=0.62(B-R)-0.18;\sigma =0.04 \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &&u-g=1.16(U-B)+1.08; \sigma =0.16 \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&u-r=0.91(U-R)+0.61; \sigma =0.21. \end{eqnarray} \tag{ 9 }$$

These relations are useful only for the NFGS due to the possibly different UBR zero points used in the NFGS relative to other photometric studies (see K13). The relations have been corrected for foreground extinction, although as we are measuring color differences, both foreground extinction corrections and k-corrections cancel. Where these conversions are used, their uncertainties are incorporated into the final error on ΔC^m.

Since we make use of optical color gradients, internal dust extinction is a concern. A systematic effect from dust may manifest itself as a dependence on inclination, since higher inclination galaxies show more dust along the line of sight. To determine whether internal extinction is systematically biasing our results, we run Spearman rank tests on all star-forming disk galaxies in the NFGS with measured inclinations >15° after excluding AGNs. We find no significant correlation between inclination and ΔC^m. Shown in Figure 3, the same rank test for galaxies with stellar masses above ∼10^9.7 M_☉—roughly equivalent to the 120 km s⁻¹ rotation velocity threshold above which dust lanes become more prominent in galaxies (Dalcanton et al. 2004; see also K13)—implies somewhat significant correlations (3%, 5%, and 4% chance of being random for Δ(u − r)^m, Δ(u − g)^m, and Δ(g − r)^m, respectively). Though weak correlations may exist, these trends are much smaller than the scatter (Figure 3). Thus we conclude that systematic internal extinction effects should have a minimal effect on our results.

Zoom In Zoom Out Reset image size

We should note that while we find only small systematic effects due to inclination, dust within galaxies will inevitably alter our color gradients if it is present. There are galaxies in our sample for which visual inspection has revealed significant dust features. These galaxies and their impact on our results are discussed in Sections 3 and 4.

To ensure that our data are as uniform as possible, we recalculate all gas masses using the same formulae. The H i mass is calculated as (Haynes & Giovanelli 1984):

$$\begin{equation} M_{\rm H\,{\scriptsize {I}}}=2.36\times 10^5D^2F_{21}. \end{equation} \tag{ 10 }$$

Here, F₂₁ is the measured 21 cm flux in Jy km s⁻¹ and D is the distance to the galaxy in Mpc. The molecular gas mass is estimated as (Sanders et al. 1991):

$$\begin{equation} M_{\rm H_2}= 1.18\times 10^4 \left(\frac{X_{\rm CO}}{3\times 10^{20}} \right)D^2 F_{\rm CO}. \end{equation} \tag{ 11 }$$

Here, F_CO is the measured CO flux in Jy km s⁻¹. We assume a constant X_CO of 2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ (Strong & Mattox 1996; Dame et al. 2001); see Section 3.1.3 for an analysis of how assuming constant X_CO may affect our results. In further analysis, we multiply all gas masses by a factor of 1.4 to account for helium.

The majority of our CO flux measurements come from single dish telescopes with single pointings at galaxy centers. Unlike 21 cm beams, which are typically several times the size of our galaxies, the CO beams are often comparable to or smaller than our galaxies. However, since CO scale lengths are typically smaller than optical scale lengths, this size mismatch is not always a major issue when trying to estimate total H₂ mass. By modeling the beam pattern of the telescope and assuming a CO flux distribution, one can estimate the additional flux not detected by a single pointing. We apply beam corrections to our single-dish CO data following the same method as Lisenfeld et al. (2011). The aperture correction is defined as the ratio of the total-to-observed CO intensity,

$$\begin{equation} f = I_{\rm CO,total}/I_{\rm CO,observed}. \end{equation} \tag{ 12 }$$

We assume that the CO distribution I_CO follows an exponential disk,

$$\begin{equation} I_{\rm CO}(r)=I_0e^{-r/h_{\rm CO}}, \end{equation} \tag{ 13 }$$

where h_CO is the scale length of the CO emission. The total CO intensity is simply the integral of this CO distribution to infinity:

$$\begin{equation} I_{\rm CO,total}=\int _0^{\infty }{2\pi I_{0}re^{-r/h_{\rm CO}}dr} = 2\pi I_0 h_{\rm CO}^2. \end{equation} \tag{ 14 }$$

The observed CO intensity is the integral of the CO distribution convolved with the beam pattern. Assuming a Gaussian beam with HPBW, Θ, this is written:

$$\begin{eqnarray*} I_{\rm CO,observed}&=& 4I_{\rm CO}\int ^{\infty }_0{\int _0^{\infty }{\exp {\left(-\frac{\sqrt{x^2+y^2}}{h_{\rm CO}}\right)}}}\nonumber \\ && \times \exp \left(-\ln {2}\left[\left(\frac{2x}{\Theta }\right)^2 \nonumber +\left(\frac{2y\cos {i}}{\Theta }\right)^2\right]\right)dxdy. \end{eqnarray*}$$

A major uncertainty in this calculation is the value of h_CO. Studies of star-forming late-type galaxies have shown h_CO ∼ (0.2 ± 0.05) R₂₅, where R₂₅ is the B-band 25 mag arcsec⁻² isophotal radius (Young et al. 1995; Leroy et al. 2008; Schruba et al. 2011). We assume this relation holds for all late type galaxies in our sample. Whether the same is true for E/S0 galaxies is uncertain, as their CO scale lengths have not been as thoroughly studied in the literature. We examine the CO profiles of Wei (2010; see also Wei et al. 2010b), who mapped the CO(1–0) distribution in a sample of E/S0 galaxies using the CARMA array. We find an average value of h_CO = 0.1 R₂₅ with a standard deviation of 0.05 R₂₅. We adopt this scale length for all E/S0 galaxies in our sample. It should be noted that these radial profiles are primarily for blue-sequence E/S0s, not the more common red-sequence E/S0s, simply because CO is more often detected in the former. The two red-sequence E/S0 galaxies in the Wei (2010) sample have CO scale lengths close to 0.05 R₂₅, but two data points are not enough to warrant separate definitions of h_CO for red- and blue-sequence E/S0s.

We assume uncertainties in h_CO of 0.05 R₂₅ for both spiral and E/S0 galaxies. The additional uncertainty that propagates into the beam-corrected H₂ mass depends on the size of the galaxy relative to the beam. Within our final sample (see Section 2.3.7), additional uncertainties are on the order of 5%, though some are as large as 25%. We also note that for galaxies in our sample for which we have both single dish CO fluxes and resolved CO maps from Wei (2010), we adopt the CO fluxes from the single dish observations and beam-correct these fluxes using the measured scale-lengths from the CO maps.

After identifying from our NFGS+literature compilation all galaxies with the necessary optical/NIR imaging as well as 21 cm and CO(1–0) flux measurements (627 galaxies), we institute a number of usability criteria that galaxies must pass before they are included in our final sample. These are: (1) SDSS r-band half-light radii larger than 5'' to ensure blue-centeredness calculations are not compromised by variations in the point-spread function between the different SDSS passbands; (2) minimal beam corrections to detected CO fluxes, so that the corrected CO fluxes are no larger than 1.5 times the measured fluxes (equivalent to a change in ${\rm H_2/{\rm H\,{\scriptsize {I}}}}<0.2$ dex), ensuring the estimated H₂ masses are minimally dependent on the assumed model of the CO distribution; (3) strong upper limits on total H₂ masses, i.e., where CO detections are missing but H i detections are not, H₂/H i must be <0.05. We do not remove any galaxies with H i upper limits from our sample. With these cuts, our NFGS+literature sample totals 323 galaxies. Included in this tally are 11 galaxies we judged to be highly peculiar/interacting whose derived P.A.s and ellipticities carry little meaning, and whose blue-centeredness calculations are therefore untrustworthy; we consider their properties in Section 3.1.2. We also include an additional 122 "quenched" galaxies (two of which are also counted as highly peculiar), which pass all usability criteria except for having upper limits in both H i and H₂ mass; these are discussed in Section 3.2.3. For our NFGS sample, upper limits from the IRAM 30 m telescope are preferred over the ARO 12 m telescope since they are always stronger. Weak upper limits led to the removal of 4 out of 39 NFGS galaxies with CO data from our final sample, though 3 of the 4 were upper limits from prior CARMA observations (Wei et al. 2010b). Where we have both IRAM 30 m telescope and ARO 12 m telescope detections for our NFGS sample, we prefer the IRAM 30 m data as long as the beam correction is less than 1.5. Otherwise, we use the ARO 12 m data.

It is important to note that while our sample is diverse, it is not statistically representative of the galaxy population since our combined data set is subject to whatever biases exist in past CO studies. Nonetheless, our usability criteria, which enforce minimum and maximum apparent radii for many galaxies in our sample, do not appear to bias us toward galaxies of a certain physical size thanks to the wide variety of distances within our sample. The distributions of physical sizes of our galaxies before and after we institute our usability criteria are not significantly different, as confirmed by a K-S test. Our final sample includes some galaxies taken from the NFGS, which were in fact originally chosen to be representative of the galaxy population (see Section 2.1.1). We will use this subset to investigate our results in the context of a broadly representative data set.

Our final sample, including the quenched and peculiar galaxies, is shown on the u − r color versus stellar mass plane in Figure 4. This sample includes some galaxies with M_* < 10^8.5 M_☉, although the blue-centeredness mass-correction was calibrated with only galaxies above this mass. These lower mass galaxies do not show unusual behavior within our results, and are therefore kept in our final sample (see Section 3.2.2 for further discussion).

Zoom In Zoom Out Reset image size

Our final compiled data set is made available in machine readable format in the online edition of our paper. A summary of the data is given in Table 3.

The fueling diagram is shown in Figure 5, plotted using mass-corrected blue-centeredness based on u − r, u − g, and g − r colors, all of which recreate the same basic structure. The regions of parameter space roughly defining the left, right, and bottom branches are overlaid on the Δ(g − r)^m data in Figure 6.

Zoom In Zoom Out Reset image size

The majority of our sample falls on the left branch, which shows a positive correlation between H₂/H i and mass-corrected blue-centeredness. There is considerable scatter which appears to increase above H₂/H i ∼ 0.5, predominantly biasing galaxies toward lower ΔC^m. We attribute at least some of this scatter to centrally concentrated dust. Visual inspection has shown several of the galaxies that scatter to the left of the main trend have distinct dust features (e.g., M82, which lies at Δ(g − r)^m = −0.14 and H₂/H i = 3.6). Galaxies in the "dusty zone" do not have preferentially high inclination (Section 2.3.4). Rather, their centrally concentrated dust may reflect a shared evolutionary state (see Section 4.1.1). The right branch of the fueling diagram, which begins to appear at Δ(g − r)^m ∼ 0.05–0.1, is far less populated but still well defined. It shows a relationship opposite to the left branch, where mass-corrected blue-centeredness increases while H₂/H i decreases. The bottom branch, which encompasses galaxies below H₂/H i ∼ 0.06, shows no clear correlation between H₂/H i and mass-corrected blue-centeredness.

Among the galaxies in our sample with gas detections and reliable ΔC^m measurements, ∼75% are on the left branch, ∼5% are on the right branch, and ∼20% are on the bottom branch. However, these percentages are very crude since the branches connect with each other and their definitions are not exact. We also stress that our sample is not statistically representative of the galaxy population, largely because of the lack of CO measurements for low-mass galaxies. Thus, the fractions of galaxies falling on the left, right, and bottom branches in Figure 5 cannot be used to infer the true frequency with which galaxies fall on these parts of the fueling diagram. If instead we limit ourselves to the more representative NFGS subsample, then the fractions of galaxies on the left, right, and bottom branches are roughly 65%, 10%, and 25%. But again, this subsample is only approximately representative.

Between these three branches exists a region where no data points lie. This hole is clearest in the version of the fueling diagram using Δ(g − r)^m because the right/bottom branches are most evenly populated. For this reason, we default to using the Δ(g − r)^m version of the plot to show further results.

The structure seen in the fueling diagram—particularly the unoccupied region between the three branches—does not appear to be artificially created by our sample restrictions. As described in Section 2.1, our full NFGS+literature sample was unrestricted in terms of stellar mass, gas content, and star formation properties. Subsequent restrictions applied were mostly related to the usability of the data. The only selection criterion we can reasonably relax is our limit on the level of permitted CO beam-correction used to estimate flux missed by the telescope beam. As discussed in Section 2.3.6, we enforce the corrected-to-measured CO flux to be less than 1.5. After relaxing the corrected-to-measured CO flux to be less than 2, 5, and 10, the three branches of the fueling diagram remain well defined (although with increased scatter) and no new galaxies appear to invalidate the existence of the empty region between the branches.

It is interesting to examine whether peculiar or actively interacting galaxies could potentially fill the empty region. These galaxies lack reliable ΔC measurements, so we do not plot them within the fueling diagram. However, we mark their measured H₂/H i to the left of the y-axis in Figure 5. Their H₂/H i values do not exclusively cluster in the range spanned by the hole, although some have the proper values to fall within it. We conclude that if galaxies like those represented in our sample ever fall in the empty region, they must do so only briefly during a phase of rapid evolution. The specific reasons why galaxies rarely settle in this region of parameter space are not immediately apparent. Detailed modeling of galaxies may be necessary to explain this phenomenon and will be the focus of future work.

The structure of galaxies in the fueling diagram appears robust against possible variations of X_CO due to its dependence on metallicity. To test this, we use O/H to estimate X_CO for each galaxy, employing the calibration from Obreschkow & Rawlings (2009). Estimates of O/H from optical line ratios are only available for our NFGS (Kewley et al. 2005) and SINGS (Moustakas et al. 2010) samples and the studies of Barone et al. (2000) and Taylor et al. (1998). To be able to carry out this analysis for our full sample, we also estimate X_CO for each galaxy from B-band luminosity, again using the calibration of Obreschkow & Rawlings (2009), who find it to be the next most reliable estimator of X_CO after metallicity. The values of 12+log O/H in our sample (where known) range from 7.9–9.2, yielding X_CO estimates between 0.4 and 10 times the Milky Way value. The calibration using B-band luminosity yields a similar range of X_CO.

Figure 7 shows the effect that variable X_CO due to metallicity has on the fueling diagram. Even with the variation in X_CO, the right/bottom branches remain distinct from the left branch and the hole remains intact, although some of the upper limits on the bottom branch show significant increases in H₂/H i due to their estimates of X_CO being several times the Milky Way value. These galaxies are mainly dwarfs with stellar masses ≲ a few × 10⁸ M_☉, where this sort of deviation might be expected. In general, however, the structure of the fueling diagram remains unchanged.

Zoom In Zoom Out Reset image size

The above analysis ignores other physics that may alter X_CO. Most relevant to our study is the possible decrease in X_CO in very high gas surface density regimes where the interstellar medium turns almost entirely molecular. Such a situation could occur as the result of inflows that drive large amounts of gas to the centers of galaxies, and has been observed in many systems from ULIRGs to ordinary spirals (Downes et al. 1993; Regan 2000). The gas in these extremely high surface density regions is expected to have increased temperature relative to molecular clouds in less dense environments (Narayanan et al. 2012), which would be coupled with an increase in the F_CO(2–1)/F_CO(1–0) ratio. Where available, we compare F_CO(2–1)/F_CO(1–0) (both beam corrected as in Section 2.3.6) with mass-corrected blue-centeredness and H₂/H i, but we find no correlations. Any evidence of increased central gas temperature is likely being washed out in our global CO measurements. Even if some of the increase in H₂/H i with ΔC^m along the left branch is the result of overestimated X_CO, this is still consistent with our physical interpretation that the left branch of the fueling diagram is the result of galaxy interactions and inflows (see Section 4.1.1).

For the remainder of this paper, we use molecular gas masses estimated as described in Section 2.3.5, assuming constant X_CO.

Figure 8(a) displays galaxy morphologies within the fueling diagram. All galaxies are classified by eye using SDSS g-band images, but we check our classifications against previously published types when available. Galaxies are separated into two categories: E/S0s (including S0a) and spirals (Sa–Sd). Using the distinction between S0a and Sa as the separation between early- and late-type galaxies may be somewhat sensitive to classification error, but this division is useful since the presence or absence of extended spiral arms represents a basic transition in structure likely strongly linked to star formation history.

We find a clear bimodality in the distribution of E/S0 and spiral morphologies. Spiral galaxies almost exclusively fall on the left branch, and in fact, a Spearman rank test on all spiral galaxies with H₂/H i > 0.06 (to avoid confusion with the bottom branch) confirms a correlation between H₂/H i and ΔC^m, with ∼5σ confidence (for u − r, u − g, and g − r, the probabilities of the distributions being random are 5 × 10⁻⁸, 8 × 10⁻⁷, and 1 × 10⁻⁷, respectively). E/S0s show a much more diverse distribution throughout the fueling diagram: some occupy the left branch with spirals, but also and more noticeably the rest almost completely define the right branch and much of the bottom branch. Several of the E/S0s on the right/bottom branches are also centrally concentrated blue compact dwarf galaxies (BCDs), which are lumped with traditional E/S0s in our simplified morphological classification scheme. Toward the left end of the bottom branch, there is a strong shift in morphology from E/S0 to spiral.

Barred galaxies are also shown in Figure 8(a). Except for galaxies that have existing bar classifications from the NFGS or Nair & Abraham (2010), bar classifications were done by eye using SDSS i-band imaging, since bars are best identified in bands that trace the stellar light (Eskridge et al. 2000). Bars are noticeably absent from all galaxies above Δ(g − r)^m ∼ 0.15, which includes most of the right branch, as well as sections of the left and bottom branches. Within the left branch alone, a Spearman Rank test does not suggest a smooth correlation of bar fraction with ΔC^m or with H₂/H i. Although our full sample is not statistically representative of the galaxy population, we speculate that the same processes that produce extreme blue-centeredness may destroy bars, while milder processes associated with evolution along the left branch do not. If we limit our examination of bars to only the more representative NFGS subsample, we still see no significant trend between ΔC^m and bar fraction within the left branch. Bars are important galactic structures to put into the context of our study, since like galaxy interactions, they are thought to enable inflows of gas to the centers of galaxies. We discuss bars and interpret their role in Section 4.3.

Figure 8(b) displays the distribution of stellar masses within the fueling diagram. Instead of examining the continuous distribution of stellar masses, we divide the data into three characteristic mass regimes: (1) M_* > M_b = 3 × 10¹⁰ M_☉, where M_b is the bimodality mass, a stellar mass scale above which the population of galaxies goes from being typically star-forming with disk-like morphology to typically non-star-forming with spheroidal morphology (Kauffmann et al. 2003), (2) M_* < M_t = 5 × 10⁹ M_☉, where M_t is the gas-richness threshold mass, below which there is a significant increase in gas-dominated galaxies (Kannappan et al. 2009; K13), and (3) M_t < M_* < M_b, where bulged spirals with intermediate mass content are the norm (K13). Galaxies with M_* < 10^8.5 M_☉ are also denoted with an extra "+" symbol in this figure. These galaxies technically fall outside the stellar mass range used to calibrate the blue-centeredness mass-correction (Section 2.3.3), and Kannappan et al. (2004) suggest their possibly irregular patterns of star-formation may make color gradients hard to interpret in relation to higher mass galaxies. However, we include them in our plots because they show no unusual distribution relative to the rest of the sample below M_t. Among these galaxies, the most extreme blue-centeredness is found for NGC 3738, with Δ(g − r)^m = 0.19. This galaxy is a centrally concentrated BCD, making it completely consistent with its neighboring galaxies in the fueling diagram.

The different stellar mass regimes display patterns within the fueling diagram. While galaxies in all three regimes span the full range of H₂/H i, we find almost no galaxies above M_b on the right/bottom branches. Instead, most of the galaxies on these branches fall below M_t. There is also a tendency for galaxies below M_t to cluster in the lower-left corner of the fueling diagram (many having H₂/H i upper limits), but elsewhere along the rising branch they appear to spread roughly evenly, as do galaxies in the higher stellar mass regimes. We note that with our low number statistics, the apparent tendency of galaxies on the upper right branch to have stellar masses between M_t and M_b is not significant and should not be over-interpreted given that the range between M_t and M_b is very narrow, and stellar mass estimation involves typical errors of ∼0.2 dex.

In Section 2.3.3, we described our motivation and methods for subtracting the underlying dependence of blue-centeredness on stellar mass, yielding mass-corrected blue-centeredness, or ΔC^m. To determine how dependent our results are on this mass correction, we plot the fueling diagram using the uncorrected Δ(g − r), rather than Δ(g − r)^m, in Figure 9. The general structure of the fueling diagram remains intact, specifically the presence of three branches with a hole between them. The distribution of different mass regimes illustrates the tendency of high-mass galaxies to have more red-centered color gradients, which originally motivated our mass correction.

Zoom In Zoom Out Reset image size

Figure 8(c) shows the distribution of red and blue-sequence E/S0s and spirals in the fueling diagram, which are classified based on their positions within the u − r versus stellar mass plane shown in Figure 4. The left branch is composed of a mix of red and blue sequence galaxies. Intriguingly, the right and bottom branches are almost completely dominated by blue sequence galaxies, despite the common assumption that E/S0 galaxies always fall on the red sequence. High-mass, red-sequence E/S0 galaxies in our sample are often quenched, i.e., have no detected CO or H i emission and total gas-to-stellar mass ratio upper limits less than 0.04 M_☉. These systems are shown below the x-axis in Figures 5, 8, 9, and 10.¹¹

Zoom In Zoom Out Reset image size

Figure 8(d) shows the distribution of the total gas-to-stellar mass ratio (total gas = H i+H₂ with a 1.4× mass correction for He) within the fueling diagram. Broadly speaking, H i-to-stellar mass ratios and total gas-to-stellar mass ratios on the left branch are both anti-correlated with H₂/H i. However, there is significant diversity below H₂/H i ∼ 0.5, where galaxies with both very high and very low gas fractions fall. Above H₂/H i ∼ 0.5, the gas fractions on the left branch are consistently low.

Most of the galaxies comprising the right and bottom branches have substantial gas content, between 10% and 100% of their stellar mass. Along the bottom branch, both H i-to-stellar and total gas-to-stellar mass ratios appear anti-correlated with mass-corrected blue-centeredness, that is, gas content is higher on the left, more red-centered side. Using a Spearman rank test on all unquenched bottom branch galaxies (defined as H₂/H i < 0.06), we find the negative correlation of the total gas-to-stellar mass ratio and Δ(g − r)^m has a 1.6% chance of being random. This probability drops to 0.2% when restricting the test to M_* > 10^8.5 M_☉ (i.e., the mass above which the blue-centeredness mass correction was calibrated). A few red-sequence galaxies with very low gas fractions also lie on the bottom branch, which are among the minority of massive (M_*  >  M_t) E/S0s whose gas data are not upper limits.

While our sample lacks uniform/complete data for AGN classification, we display known AGNs in Figure 10 and briefly discuss them for two reasons. First, it is important to note that AGNs themselves are not the cause of strongly blue-centered color gradients, as their light contribution is typically too small to have any significant effect on mass-corrected blue-centeredness, which uses colors measured over large regions of galaxies. Second, we note that AGNs are predominantly seen among the high mass galaxies in our sample (see also Figure 8(b)), with 25 out of 27 AGNs hosted by galaxies with M_* > M_t. Their presence, particularly among the high mass elliptical galaxies which show AGNs in the quenched regime but also among galaxies with detected gas, affects our physical interpretation of the fueling diagram, and is discussed in Sections 4.1.1 and 4.1.3.

To summarize, the main properties of each branch of the fueling diagram are as follows.

The left branch is where spiral galaxies are primarily found. Many galaxies here are barred, but without any discernible pattern with respect to ΔC^m or H₂/H i. Galaxies on the left branch cover the full range of stellar masses, along with a wide range of H i- and total gas-to-stellar mass ratios, which broadly behave inversely to H₂/H i.

The right and bottom branches are almost entirely populated by unbarred E/S0 galaxies, although this type of distribution largely transitions into spirals on the left side of the bottom branch, which is also where we begin to see barred galaxies again. Stellar masses fall predominantly below the bimodality mass, and most fall below the threshold mass. Gas fractions are moderate to high, between 10% and 100% of the stellar mass, and show a negative correlation with mass-corrected blue-centeredness along the bottom branch.

We also note that the right/bottom branches are dominated by blue-sequence E/S0 galaxies, i.e., galaxies with spheroidal morphologies that fall on the blue sequence in color versus stellar mass space, which are thought to represent a transitional phase (Kannappan et al. 2009). This fact, as well as the properties of galaxies within the different branches of the fueling diagram, may reveal important clues as to how these branches relate to different evolutionary states, as discussed in Section 4.

The very existence of the left branch provides support for the idea that local galaxy interactions play a key role in replenishment of molecular gas: since mass-corrected blue-centeredness is a signpost of recent galaxy interactions (Kannappan et al. 2004), and it is correlated with H₂/H i, galaxy interactions appear to be linked to H₂/H i in a systematic way. The analysis of Δ(u − g)^m* and Δ(g − r)^m* in Section 3.3 implies that galaxies on the left branch may evolve along it in both directions: up as a result of each galaxy encounter, boosting the galaxy to higher H₂/H i and ΔC^m (with a possible contribution due to overluminous CO), then down as the molecular reservoir is consumed and the young central stellar population fades. Total gas-to-stellar mass ratios support this scenario. We expect H i and total gas-to-stellar mass ratios to decrease at higher ΔC^m since the gas is being converted into H₂ and then stars, but we also expect a mix of gas fractions at low ΔC^m since galaxies settle here before or after each burst event while accreting fresh disk gas. The bulk of the evolution along the left branch must be driven by minor rather than major mergers or interactions, since most galaxies appear to retain their spiral morphologies. An alternate possibility is that the relation along the left branch reflects bar-driven inflows. However, bars do not show any statistically significant correlation with ΔC^m along the left branch. The possible role of bars is explored further in Section 4.3.

Major mergers in the high stellar mass and/or low gas fraction regime may help to explain the E/S0s at high H₂/H i on the left branch (in contrast, gas-rich, low mass, but comparable mass ratio mergers appear to follow a different path; see Sections 4.1.2–4.1.3). Alternatively, these galaxies could be formed by repeated minor mergers rather than a single event (Bournaud et al. 2007), and they could therefore be galaxies that have made several oscillations along the left branch. Either way, many of the E/S0s at high H₂/H i on the left branch are likely moving into the quenched regime, up and off the plot, as their remaining gas reservoirs have become almost entirely molecular and may soon be completely depleted. In addition, some of the high stellar mass E/S0s at the peak of the rising branch may have previously been quenched but recently experienced small gas accretion events, possibly associated with satellite accretion (e.g., Martini et al. 2013). Some of these E/S0s host AGNs, as do many of the high mass E/S0s in the quenched regime. A gas accretion event could provide fuel for AGN activity coupled with a central molecular gas concentration, resulting in a high H₂/H i ratio with minimal change to ΔC^m (since most of the star formation would be occurring very close to the AGN). Galaxies undergoing this process would jump vertically in the fueling diagram, between the quenched regime and the peak of the left branch. This path could place them in the hole of the fueling diagram, although they would likely move through it relatively quickly.

As previously noted in Section 3.1.1, there is substantial scatter above the left branch toward more red-centered color gradients, likely caused by internal extinction. While dust effects may be altering the measured color gradients, they usefully highlight galaxies in early stages of star-forming events when the young stars are still heavily embedded in dust clouds. We suspect galaxies in this "dusty zone" represent a stage very soon after mergers or interactions that are mild enough not to have driven the galaxy into the peculiar morphology region off the plot, shown on the left hand side of the fueling diagram. As the star formation progresses, the dust is likely to clear, allowing these galaxies to develop the bluer-centered color gradients expected for their H₂/H i ratios.

We suggest galaxies on the right branch may be the result of gas-rich mergers, specifically between galaxies of comparable stellar mass. Along with their high values of mass-corrected blue-centeredness (suggestive of a strong central starburst), the spheroidal morphologies of galaxies on the right branch are consistent with recently violent histories. Most of these galaxies are classified as blue-sequence E/S0s, and several are also classified as BCDs (e.g., Haro 2, NGC 7077, UM 465). The existence of blue-sequence E/S0s and BCDs in the same regime of the fueling diagram argues in favor of them experiencing similar evolutionary processes. Previous observational and theoretical studies of BCDs (e.g., Östlin et al. 2001; Pustilnik et al. 2001; Bekki 2008) also support merger driven evolution.

While most of these galaxies do not show obvious outward signs of a recent strong interaction in their optical images, such as irregular structure or tidal features (the lack of these features is actually built into our analysis since we do not plot highly peculiar/interacting systems), smooth optical morphology is not inconsistent with a merger having recently occurred. Merger simulations find the strongest morphological disturbances before the peak of induced star formation, which in turn typically occurs before galaxies land on the right branch, and the complete coalescence of the two merging galaxy nuclei commonly happens several hundred Myr before the main starburst event (Lotz et al. 2008, 2010). The diverse arrow directions in Figure 11 on the right branch suggest these galaxies are still actively forming stars well after the merger remnant has settled. Signatures of recent mergers may be more obvious in observations of gas morphology and kinematics. H i maps can be extremely useful since they trace extended structure and can retain signatures of interactions as long as 1 Gyr after the events (Holwerda et al. 2011). For example, the high resolution H i map of Haro 2 (Bravo-Alfaro et al. 2004) shows the H i kinematic and optical major axes to be almost perpendicular, consistent with a merger or recent accretion event.

Blue-sequence E/S0s are known to emerge primarily below M_b, and become abundant below M_t (Kannappan et al. 2009). As seen in Figure 8, the blue E/S0s on the right branch are consistent with this pattern. However, the existence of the right branch cannot be driven by stellar mass alone. Low stellar mass galaxies are found throughout the fueling diagram, and if ΔC^m were dictated solely by stellar mass (i.e., equal size bursts occurring in higher/low mass galaxies yielding lower/higher ΔC^m), then the hole seen in the fueling diagram should not exist. A merger origin for the right branch is more consistent with such a large gap between the left and right branches. Furthermore, gas richness likely produces distinct evolutionary tracks for galaxies: gas-rich mergers drive galaxies along the right branch, while gas-poor mergers drive galaxies into the quenched regime, up and off the plot, as discussed in Section 4.1.1. The association of the gas-rich merger track with galaxies below M_b and especially M_t is a simple consequence of increasing gas richness below those scales (Kannappan et al. 2009; K13).

The bottom branch appears to be part of the same evolutionary sequence as the right branch but at a later stage. Galaxies here show many of the same properties as those on the right branch in terms of stellar mass, gas richness, and prominence of blue-sequence E/S0s. However, they show depressed H₂/H i ratios while uniformly evolving leftward in the fueling diagram (Section 3.3). These observations are all consistent with a scenario where these are post-starburst galaxies with depleted central gas and fading young central stellar populations.