Dominant Role of Coplanar Inflows in Driving Disk Evolution Revealed by Gas-Phase Metallicity Gradients

Our parent sample is derived from the Sloan Digital Sky Survey IV (SDSS IV; Blanton et al., 2017) Mapping Nearby Galaxies at Apache Point Observatory (MaNGA; Bundy et al., 2015) survey. MaNGA provides two-dimensional spectra for approximately 10,000 galaxies with redshifts in the range $0.01<z<0.15$, utilizing the spectrographs on the 2.5-meter Sloan Telescope at Apache Point Observatory (Gunn et al., 2006; Smee et al., 2013). The survey covers wavelengths from 3600 to 10300 Å with a spectral resolution of $R\sim 2000$. Spatial coverage typically extends beyond 1.5 $R_{\mathrm{e}}$ for individual galaxies, with a spatial resolution of 1–2 kpc. Detailed descriptions of the MaNGA instrument design, testing, and assembly are provided in Drory et al. (2015).

We adopt the emission line flux and stellar and gas kinematics measurements from the MaNGA Data Analysis Pipeline (DAP) (Westfall et al., 2019; Belfiore et al., 2019b). The DAP employs adaptive Voronoi spatial binning (Cappellari & Copin, 2003) to achieve a target signal-to-noise ratio (S/N) of approximately 10 in the stellar continuum. After measuring stellar kinematics, the DAP fits Gaussian profiles to derive emission line fluxes for each spaxel using the full spectral fitting code pPXF (Cappellari & Emsellem, 2004; Cappellari, 2017).

We cross-match the MaNGA galaxy sample with the GALEX-SDSS-WISE Legacy Catalog (GSWLC; Salim et al., 2016, 2018), which provides galaxy stellar masses derived from spectral energy distribution (SED) fitting of ultraviolet, optical, and infrared photometry. Our selection criteria include galaxies with stellar mass and star formation rate (SFR) data available in the GSWLC-M2 catalog. We adopt the SERSIC_TH50 as effective radii ($R_{\mathrm{e}}$) from the NASA-Sloan Atlas (NSA; Blanton et al., 2011) and convert it to physical scale assuming a flat cold dark matter cosmology model with $\Omega_{m},\Omega_{\Lambda},h=0.27,0.73,0.7$. We select galaxies with axis ratios ($b/a$) larger than 0.4, where $b/a$ are taken from SERSIC_BA in NSA catalog. Stellar masses are rescaled to a Kroupa (2001) initial mass function (IMF). We focus on star-forming galaxies, selecting those with $\Delta$MS $>-1$ dex. Here, $\Delta$MS represents the vertical offset from the star-forming main sequence (SFMS):

which is defined in Ma et al. (2024) for MaNGA sample based on an iterative algorithm (Woo et al., 2013; Donnari et al., 2019; Wang & Lilly, 2023). Specifically, the algorithm begins with an initial guess of a linear function and iteratively selects galaxies within 1 dex of the SFMS. It then refits the slope and intercept until convergence is achieved. This selection results in a final sample of 3,808 galaxies.

The gas velocity dispersion ($\sigma_{\mathrm{gas}}$) is also taken from the MaNGA DAP product. For each galaxy, we select spaxels with good data quality (indicated by QUALDATA in MaNGA DAP) and valid H$\alpha$ fluxes within galactocentric radii of 0.5–1.5 $R_{\mathrm{e}}$. After correcting for instrumental dispersion, we use the median H$\alpha$ velocity dispersion in annuli 0.5–1.5 $R_{\mathrm{e}}$ for each galaxy, as an indicator of gas turbulence.

For each galaxy, we identify star-forming and composite spaxels within galactocentric radii of 0.5–1.5 $R_{\mathrm{e}}$ using the BPT-NII diagnostic diagram (Baldwin et al., 1981; Kewley et al., 2001; Kauffmann et al., 2003). Selected spaxels meet the criteria S/N $>3$ for the fluxes of H$\alpha$, H$\beta$, [OIII]$\lambda$5007, [NII]$\lambda$6583, and [OII]$\lambda$3737, 3729. Emission line fluxes are corrected for interstellar reddening using the Cardelli et al. (1989) reddening law with $R_{V}=3.1$. Extinction is calculated by comparing the measured H$\alpha$/H$\beta$ flux ratio to the theoretical value of 2.86 (Osterbrock & Ferland, 2006), assuming case B recombination. Spaxels classified as Seyfert or LINERs are excluded. For star-forming and composite spaxels, we determine the H$\alpha$-flux-weighted gas-phase metallicity (12+log(O/H)) using the N2S2H$\alpha$ calibrator proposed by Dopita et al. (2016), which is insensitive to the reddening and ionization parameter via the [NII]/H$\alpha$ term (e.g., Zhang et al., 2017; Hwang et al., 2019). In addition, Easeman et al. (2024) proposed that the N2S2H$\alpha$ is preferred when studying the distribution of metals within galaxies because of the near-linear relation with $T_{\mathrm{e}}$-based measurement. Alternative calibrations, such as S-calibration (Scal, Pilyugin & Grebel, 2016) and N2O2 (Dopita et al., 2013; Zhang et al., 2017) also confirm the trends in this work.

We focus on disk-dominated regions. For each galaxy, the median gas-phase metallicity is calculated by considering all spaxels within each annulus, ranging from 0.5 $R_{\mathrm{e}}$ to 1.5 $R_{\mathrm{e}}$ with each annulus having a width of 0.2 $R_{\mathrm{e}}$, requiring at least 70% of spaxels to be available within each annulus. Gas-phase metallicity radial gradients, expressed in dex/$R_{\mathrm{e}}$ ($\nabla\log(\mathrm{O/H})$), are derived using linear fits to the radial profiles in this range. This approach minimizes the impact of smearing effects and avoids significant deviations from linear fits (Sánchez et al., 2014; Sánchez-Menguiano et al., 2016; Belfiore et al., 2017; Boardman et al., 2022; Barrera-Ballesteros et al., 2023). Metallicity diagnostics are typically calibrated based on HII region models. However, one may be concerned that the spaxel sample selection criteria (e.g., BPT-NII) could result in non-negligible contamination from diffuse ionized gas (DIG) when measuring chemical abundances from emission line ratios (e.g., Belfiore et al., 2017; Vale Asari et al., 2019). To address this issue, we have tested the requirement that spaxels must have a minimum equivalent width of H$\alpha$ (EW(H$\alpha$)) of 10 Å (e.g., Vale Asari et al., 2019) to derive the metallicity profiles and gradients of each galaxy. Combined with our sample selection criteria, this EW(H$\alpha$) cut removes approximately 24% of the used spaxels. However, since we calculate the metallicity gradient using the median metallicity value of each annulus and subsequently apply linear fitting, this cut has a negligible impact on our final valid galaxy sample (it only reduces the sample by $\sim$8.1%, from 3,808 to 3,499 valid galaxies). Furthermore, the impact of this cut on our results is minimal.

We utilize two-dimensional stellar mass maps from Lu et al. (2023), derived using pPXF. Stellar population synthesis was performed using the FSPS model (Conroy et al., 2009; Conroy & Gunn, 2010), encompassing 43 age and 9 metallicity grids. The datacubes are Voronoi binned to achieve S/N $\sim 30$ prior to spectral fitting. For a given Voronoi bin, spaxels share the same stellar mass.

Stellar mass surface density ($\Sigma_{*}$) is calculated for each spaxel by dividing the stellar mass by the physical area of the pixel, corrected for inclination effect. To avoid biases from incomplete coverage, the median $\Sigma_{*}$ is calculated by considering all spaxels within each annulus, ranging from 0.5 $R_{\mathrm{e}}$ to 1.5 $R_{\mathrm{e}}$ with each annulus having a width of 0.2 $R_{\mathrm{e}}$. While the annuli at least 70% spaxels available are finally adopted to build the radial profile. The radial gradients of $\Sigma_{*}$, in units of dex/$R_{\mathrm{e}}$ ($\nabla\log\Sigma_{*}$), are determined by linear fits to the radial profiles in the same range.

To explore the potential impact of gas coplanar inflow on metallicity distribution, we examine the gas-phase metallicity gradient $\nabla\log(\mathrm{O/H})$ as a function of stellar mass and local gas-phase metallicity at the effective radius, in the top panels of Figure 2. Overall, the local gas-phase metallicity at the effective radius increases with stellar mass, reflecting a local manifestation of the well-established mass-metallicity relation (MZR, e.g., Tremonti et al., 2004; Mannucci et al., 2010; Steidel et al., 2014; Wuyts et al., 2014; Sánchez et al., 2019; Ma et al., 2024). At a given stellar mass, galaxies with lower gas-phase metallicity tend to exhibit more negative or steeper $\nabla\log(\mathrm{O/H})$. This trend aligns with the predictions of the modified accretion disk model proposed by Wang & Lilly (2022a), wherein inflowing gas is progressively enriched by in situ star formation as it slowly spirals inward toward the disk center. Consequently, higher metallicity in the disk reduces the difference between the metallicities of the outer and inner regions, leading to a relatively flatter $\nabla\log(\mathrm{O/H})$. Jara-Ferreira et al. (2024) found that low-mass galaxies with strong positive metallicity gradients are less enriched than the median mass-metallicity relation, while those with strong negative gradients are more enriched. This connection between metallicity and metallicity gradient is observed in both the MaNGA sample and the EAGLE simulations, and is consistent with our results. Furthermore, this mass-metallicity-gradient relation also follows the mass-size-metallicity relationship (e.g., Ellison et al., 2008; D’Eugenio et al., 2018) combined with the mass-size-gradient trend reported in Boardman et al. (2021), where at a given stellar mass, larger galaxies on average possess lower gas-phase metallicities and display steeper metallicity gradient than smaller galaxies. The relation between stellar mass, gas-phase metallicity, and metallicity gradient for MaNGA data has been also presented by Franchetto et al. (2021). They found that, at a given stellar mass, metallicity profiles tend to flatten as metallicity increases, which is consistent with the trends shown in Figure 2. Additionally, they observed that the increasing trend of metallicity is also found to be similar to the decreasing trend of gas fraction on the stellar mass–metallicity gradient plane.

To account for the mass dependence of the metallicity gradient, we calculate the metallicity excess ($\Delta\log(\mathrm{O/H})$) for each galaxy, defined as the deviation of the metallicity from the median value for galaxies of the same stellar mass (represented by the black curve and dots in the top panels of Figure 2). By removing the mass dependence, $\Delta\log(\mathrm{O/H})$ provides a more direct measure of the relative enrichment of the interstellar medium (ISM) at the effective radius of each galaxy. Although its value can be influenced by calibration methods and measurement indicators, the metallicity excess reflects the balance of metal production and transport within star-forming galaxy disks. As such, it is a fundamental parameter for understanding galactic chemical evolution on the disk.

To intuitively illustrate the relationship between metallicity excess and the radial metallicity profiles of galaxies, we present the median metallicity radial profiles for the $R_{\mathrm{e}}$ range in the bottom panels of Figure 2. These profiles are shown across different stellar mass bins (individual panels) and metallicity excess bins (distinguished by colors). From left to right, as stellar mass increases, the overall normalization of the metallicity profiles rises systematically. Within each stellar mass bin, the profiles become progressively flatter as the metallicity excess increases (from blue to red curves). This trend demonstrates the strong correlation between local metallicity excess and the metallicity gradient in the disk.

Figure 3 examines the metallicity gradients as functions of both $\Delta\log(\mathrm{O/H})$ and the gradient of $\Sigma_{*}$. The results show that $\nabla\log(\mathrm{O/H})$ is only weakly correlated with $\nabla\log\Sigma_{*}$, exhibiting a Pearson correlation coefficient of $\sim$0.058 (with $p$-value of 0.0011), compared to its much stronger correlation with $\Delta\log(\mathrm{O/H})$. This distinction is further supported by the angles of partial correlation (e.g., Bait et al., 2017; Bluck et al., 2019; Baker et al., 2023), as indicated by the black arrows. Furthermore, the weak correlation between $\Delta\log(\mathrm{O/H})$ and $\nabla\log\Sigma_{*}$ shows that the two quantities are not degenerate. Therefore, the metallicity gradient is not closely linked to the stellar mass surface density gradient, effectively ruling out the dominant role of Discrete Evolution Mode in driving disk formation.

The gas disk can be highly turbulent due to stellar winds and supernova explosions, which contribute to radial gas mixing and reduce metallicity gradients. Processes such as turbulence, stellar feedback, accretion, mergers, and satellite interactions mix metals in the ISM, generally flattening metallicity gradients. This effect is particularly pronounced at early times when turbulence is higher and the ISM is geometrically thicker (Bird et al., 2012, 2021; McCluskey et al., 2024; Graf et al., 2024). Additional detailed mechanisms that influence metallicity gradients include galactic winds, which eject metals out of galaxies (e.g., Hayward & Hopkins, 2017; Pandya et al., 2021), and radial advection, which transports pristine or metal-rich gas across the disk in bulk flows (e.g., Sharda et al., 2021). Together, these processes highlight the dynamic and interconnected nature of gas mixing and its role in shaping metallicity profiles.

Motivated by these processes, we define $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$ as a proxy for local gas turbulence within the galaxy disk. This quantity encapsulates the strength of gas turbulence at galactic scale for individual galaxies, which serves as an imprint of various physical processes. Gas turbulence tends to flatten metallicity gradients in star-forming galaxies through various physical processes that drive gas mixing and enrichment. However, distinguishing between these processes is challenging without numerous physical assumptions. In Figure 4, we plot the metallicity gradient as a function of metallicity excess and $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$. We find that $\nabla\log(\mathrm{O/H})$ correlates strongly with both $\Delta\log(\mathrm{O/H})$ and $\log(\sigma_{\mathrm{gas}}/R_{\mathrm{e}})$, consistent with the diffusion effect described by Wang & Lilly (2022a), where larger turbulent velocities lead to more flattened metallicity profiles. These trends are clearly depicted through the median metallicity profiles across different ranges of $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$ and $\Delta\log(\mathrm{O/H})$, as presented in the bottom panel of Figure 4. We further quantify the combined effect of $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$ and $\Delta\log(\mathrm{O/H})$ in shaping metallicity gradients in Appendix A and consider alternative combinations of physical properties in predicting metallicity gradient in Appendix B. The residual scatters in Figure 7 within Appendix B further demonstrate that $\Delta\log(\mathrm{O/H})$ is strongly related to the metallicity gradient, with $\sigma_{\mathrm{gas}}$ and $R_{\mathrm{e}}$ being comparably informative on the metallicity gradient.

These findings suggest that the gas-phase metallicity gradient results from competing effects of coplanar gas inflow and local gas turbulence on the disk. These two physical quantities co-evolve through the interplay of processes, such as gas inflow, star formation, metal enrichment, gas outflows, and gas mixing. Consequently, there is also a positive correlation between metallicity excess and $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$, as shown in Figure 4. This relationship underscores the complex interactions between turbulence and chemical enrichment in galaxy disks.

Previous studies have inferred the roles of various physical processes in shaping metallicity gradients by examining correlations with simple, reliable galaxy properties. Observational evidence indicates that metallicity gradients correlate with a range of galaxy characteristics, including stellar mass (e.g., Belfiore et al., 2017; Mingozzi et al., 2020; Schaefer et al., 2020; Franchetto et al., 2021), galaxy size (e.g., Carton et al., 2018; Boardman et al., 2021, 2022, 2023; Lin et al., 2024), gas fraction (e.g., Franchetto et al., 2021), bulge-to-total (B/T) ratio (e.g., Moran et al., 2012), and environmental factors (e.g., Lian et al., 2019; Franchetto et al., 2021). Identifying a concise set of physical quantities that encapsulate or reflect the impact of these processes is critical to understanding the factors driving metallicity gradients.

To investigate whether metallicity excess and local gas turbulence sufficiently describe metallicity gradients, we employ a Random Forest model. This machine-learning approach allows us to assess the contributions of different galaxy properties to the metallicity gradient. In this context, relative feature importance quantifies the contribution of input features to the model’s predictions. Using the scikit-learn implementation, feature importance is ranked based on the Gini importance (Pedregosa et al., 2011). Figure 5 indicates that metallicity excess and $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$ dominate the feature importance rankings. These two quantities significantly outperform other commonly used galaxy properties, such as stellar mass, effective radius, star formation rate, and effective radius, in predicting metallicity gradients. While the $\nabla\log\Sigma_{*}$ plays a relatively minor role, consistent with the trends in Figure 3. We have also checked to perform an alternative Random Forest model with independent parameters ($\sigma_{\mathrm{gas}}$, $M_{*}$, $R_{\mathrm{e}}$, $\Delta\log\mathrm{(O/H)}$, SFR, redshift, $\nabla\log\Sigma_{*}$, and a random variable). The analysis consistently identifies $\Delta\log\mathrm{(O/H)}$, $\sigma_{\mathrm{gas}}$, and $R_{\mathrm{e}}$ as the top three most important features, which aligns with the results in Figure 5. Additionally, we have tested alternative metallicity indicators and calibrations to calculate both metallicity gradients and metallicity values, finding consistent results (see Appendix C).

This consistency supports the robustness of our findings that metallicity excess and $\sigma_{\mathrm{gas}}/R_{\mathrm{e}}$ as predictive features reinforce the role of coplanar gas inflow in disk formation and growth. These results also support the hypothesis that coplanar gas inflow plays a dominant role in driving the disk evolution of star-forming galaxies.