VERTICO and IllustrisTNG: The Spatially Resolved Effects of Environment on Galactic Gas

The observational samples of field and cluster galaxies used in this paper are drawn from HERACLES (Leroy et al. 2009) and VERTICO (Brown et al. 2021), respectively. The analysis sample used in this paper is a superset of that used in Watts et al. (2023), and there is a detailed discussion of the observational biases present in comparing these data in that work. Briefly, the two data sets contain star-forming late-type galaxies that are well matched in global stellar mass and specific star formation rate (SFR). The H iobservations for VERTICO and HERACLES galaxies are drawn from the Very Large Array  Imaging of Virgo in Atomic gas survey (VIVA; Chung et al. 2009) and The H i Nearby Galaxy Survey (Walter et al. 2008), respectively. All molecular-gas information is derived from the public data cubes using the methodology described in Brown et al. (2021). Similarly, global stellar-mass and SFR estimates are drawn from the z = 0 Multiwavelength Galaxy Synthesis database (Leroy et al. 2019), while resolved stellar and SFR surface densities are derived from identical multiwavelength data for both surveys following the methods outlined in Villanueva et al. (2022) and Jiménez-Donaire et al. (2023), respectively. VERTICO is typically 2–3 times more sensitive than HERACLES in both H i and H₂ densities, which ensures that any observed differences with environment are not driven by sensitivity differences between the surveys.

For this work, we select we only select galaxies with stellar mass m_* > 10⁹ M_⊙ and inclinations less than 70°, ensuring the galaxies' surface density maps can be reliably deprojected with a simple cosine-of-inclination correction factor. These selections yield a final resolution-matched sample consisting of 10 field galaxies from HERACLES and 33 cluster galaxies from VERTICO. These galaxies all have star formation rates SFR ≳ 0.05 M_⊙ yr⁻¹.

Very briefly, each galaxy has molecular gas, stellar mass, and SFR surface density maps that have been smoothed to a spatial resolution of ∼1.2 kpc with a pixel size of ∼650 pc to approximately Nyquist-sample the smoothing kernel [or CO(2–1) resolving beam]. This approximately matches the resolution of the H i data from VIVA (including an update where the D-configuration data were removed, leaving only the higher-resolution C-configuration data). The derived data products for VERTICO and HERACLES are produced using a near-identical procedure to that described in Brown et al. (2021) for the molecular-gas surface densities, in Jiménez-Donaire et al. (2023) for the star formation rates, and in Watts et al. (2023) for the stellar surface densities. The only difference between the data presented in previous VERTICO papers and here is that we have reduced the molecular-gas surface densities by a factor of 1.36 to remove the "helium contribution," as we are specifically interested in H₂.

The IllustrisTNG model of galaxy formation (Weinberger et al. 2017; Pillepich et al. 2018) comprises key descriptions of astrophysical processes, including gas cooling, star and black hole formation, stellar evolution and feedback, feedback from active galactic nuclei (AGNs), and more, all within a ΛCDM cosmological, magnetohydrodynamic framework with the moving-mesh arepo code (Springel 2010). Simulations of various box sizes and resolutions have been run with the IllustrisTNG model, all of which carry identical parameters for both the sub-resolution physical prescriptions and cosmology, with the latter based on Planck Collaboration et al. (2016). We use the TNG50 simulation (Nelson et al. 2019a; Pillepich et al. 2019) in this paper, which has a periodic box length of ∼50 cMpc and baryonic mass resolution of 8.5 × 10⁴ M_⊙. Galaxy subhalos in TNG are identified with subfind (Springel et al. 2001; Dolag et al. 2009).

The H i and H₂ properties of gas cells in TNG were calculated in post-processing (Diemer et al. 2018; Stevens et al. 2019b). In this paper, we use the decomposition based on Gnedin & Draine (2014) as described in Stevens et al. (2019b).

The integrated H i and H₂ properties of TNG100 galaxies and their trends with environment have been explored in depth (Diemer et al. 2019; Stevens et al. 2019b, 2021), but the same is not true in the literature ¹⁷ for TNG50 (but see Boselli et al. 2023). In brief, after mock-observing the simulated galaxies to match survey specifications, Stevens et al. (2019b, 2021) show that TNG100 quantitatively reproduces gas-fraction trends seen with ALFALFA, xGASS, and xCOLD GASS (i.e., from the results of Brown et al. 2017; Saintonge et al. 2017; Catinella et al. 2018). To show that the gas fractions of TNG50 and TNG100 are consistent, we compare the H i and H₂ fractions of the simulations in Figure 1. To approximately match the galaxies in the observational samples described in Section 2.1, we exclusively consider TNG galaxies at z = 0 with m_* ≥ 10⁹ M_⊙ and $\mathrm{SFR}\geqslant 0.05\,{{M}_{\odot }\,{yr}}^{-1}$ (based on measurements internal to the "BaryMP" radius of Stevens et al. 2014, also referred to as the "inherent" properties in Stevens et al. 2019b) in this figure and throughout this Letter. We also exclude any galaxy with a dark-matter fraction below 5% to conservatively remove non-cosmological objects (see the criteria in Nelson et al. 2019b). This totals 2479 galaxies from TNG50. Other than a small systematic increase in H₂ fraction for TNG50, the two simulations are well aligned. There is generally good agreement between global H i and H₂ gas fractions in TNG100 and low-redshift observations (e.g., ALFALFA, xGASS, and xCOLD GASS; Stevens et al. 2019a; Diemer et al. 2019), especially after mock-observing TNG galaxies to survey specifications. For the purposes of this paper, we infer by extension that TNG50 sufficiently agrees with observations. We note that in Figure 1 the TNG50 cluster sample has higher gas fractions than their VERTICO counterparts at fixed stellar mass, particularly where m_* ≤ 10¹⁰ M_⊙. Due to the small number statistics in this regime, we caution against overinterpretation, yet it is possible that this is due to a systematic difference between the simulated and observed samples. For example, environmental processes in TNG50 may not be regulating gas content in lower stellar mass simulated galaxies to the same extent as for the observed galaxies.

Zoom In Zoom Out Reset image size

We create mock surface density maps of H i, H₂, and stellar mass for each TNG galaxy in our sample. These maps are matched to the specifications of the HERACLES and VERTICO data by convolving each map at the simulation's native resolution with a Gaussian kernel with a full width at half maximum of 1.2 kpc, and then resampling with square pixels of length 0.65 kpc. This physical scale matches the resolution of the observational data. Each map is made face-on using the angular-momentum vector of all neutral gas (any gas that is not ionized) inside the "BaryMP" radius (Stevens et al. 2014), under the assumption that the observations have been correctly deprojected.

We separate TNG50 galaxies into two subsamples after applying the above cuts. The "field" sample, intended to be comparable to HERACLES, is selected to contain only central galaxies (those in the most massive subfind subhalo) in haloes with virial mass M_200c ≤ 10^12.2 M_⊙. The field sample totals 1739 galaxies with an average number of 1514 pixels with Σ_* > 1 M_⊙ pc⁻² per galaxy. The "cluster" sample, comparable to VERTICO, contains only satellites (galaxies that are not centrals) in the simulation's two clusters ¹⁸ with M_200c ≳ 10¹⁴ M_⊙. The most massive cluster, with M_200c = 1.8 × 10¹⁴ M_⊙, is very similar in mass to Virgo (1.4–4.2 × 10¹⁴ M_⊙; see Table 1 of Boselli et al. 2018 and references therein). The other TNG50 cluster is about half this mass (M_200c = 9.4 × 10¹³ M_⊙) but has a different dynamical state and has already been justified as a Virgo analog by Joshi et al. (2021). The cluster sample totals 41 galaxies with an average of 1725 pixels with Σ_* > 1 M_⊙ pc⁻² per galaxy. The higher pixel count per galaxy in the cluster sample reflects its higher average stellar mass (see the distribution of TNG50 points in Figure 1).

Our first aim is to see if the average behavior of kiloparsec-scale gas in TNG50 reflects that of reality in both the field and cluster environments. In Figure 2, we present the pixel-based relations of Σ_{H I} and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ as a function of Σ_*. The pixels for all galaxies in the TNG50 cluster sample are grouped together, as are those in the field sample. The same grouping strategy applies for HERACLES and VERTICO. To demonstrate that the preference of HERACLES galaxies to have m_* ≳ 10¹⁰ M_⊙ (seen in Figure 1) does not affect our comparison, we add results to Figure 2 for a TNG50 field subsample where we have excluded galaxies with m_* ≤ 10¹⁰ M_⊙. We exclude the resolved molecular Kennicutt–Schmidt relation (${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$–Σ_SFR) from Figure 2 as we find that the relation is independent of environment to first order in TNG50. Indeed, Jiménez-Donaire et al. (2023) have also shown that this is true in observations. Because of this, we also choose not to include the resolved star-forming main sequence (Σ_*–Σ_SFR) as it presents similar information to the resolved molecular-gas main sequence in TNG50, as it does in observations (Σ_*–${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$; addressed further in Brown et al. 2023). See Motwani et al. (2022) for further analysis on kiloparsec-scale Σ_SFR in TNG50.

Zoom In Zoom Out Reset image size

The left column of panels in Figure 2 accounts for all pixels with Σ_* detections, regardless of whether they were detected in H i or H₂. For the purposes of calculating percentiles, gas non-detections are assumed to have zero mass. The sharp downturn in the VERTICO and HERACLES medians from right to left is simply representative of the threshold surface densities needed for gas to be detected in those surveys.

From these panels, we can immediately identify that TNG50 quantitatively reproduces the correlation between ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ and Σ_*, while the Σ_{H I}–Σ_* relation in the simulation is systematically offset to lower Σ_{H I} values. VERTICO cluster galaxies are systematically suppressed in Σ_{H I} by a factor of ∼2 at fixed Σ_* relative to HERACLES but not significantly for ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$. This is in line with the findings of Watts et al. (2023) even though we are including lower-mass galaxies than in that study.

By contrast to the observations, the Σ_{H I} medians are not parallel for the field and cluster samples in TNG50, instead converging at high Σ_*. Both Σ_{H I} and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ become increasingly divergent between the TNG50 cluster and field samples at Σ_* ≲ 30 M_⊙ pc⁻². This average behavior is typical of disk truncation, which one would expect from processes like ram pressure stripping that affect the cluster galaxies.

TNG50 galaxies in field and cluster samples are systematically low in Σ_{H I} relative to observations and exhibit little to no environmental dependence at high Σ_* (the centers of galaxies). This central deficit in H i appears to be widespread in TNG50, a finding described in Gebek et al. (2023), and also seen in TNG100 by Diemer et al. (2019) and in Figures A1 and A2 of Stevens et al. (2019a). A similar outcome has also been found for the EAGLE simulation (Bahé et al. 2016). Discussion in Section 4.5.2 of Diemer et al. (2019) and Section 4.3 of Gebek et al. (2023) describes how (some) central H i deficits in TNG are likely the result of high ionized fractions in the interstellar medium and AGN feedback removing gas from galaxy centers (also see Stevens et al. 2019b, 2021). We also note that among all the post-processing methods used to decompose neutral gas in TNG into its atomic and molecular components, all have molecular fractions that asymptote to 1 at high densities, evidently leading to a relatively low saturation density for H i (see Diemer et al. 2018).

There are also fewer pixels per galaxy at high Σ_* in TNG50 than observations, as seen in the top panels of Figure 2. This implies that galaxy centers are generally underdense in the simulation. Numerical heating of stellar particles plays a role here in artificially decreasing central stellar densities (see Ludlow et al. 2021). This same effect has a minimal influence on gas, suggesting pixels in TNG50 galaxies are leftward of where they should be in Figure 2. The low gas densities, combined with this systematic underestimation of Σ_* in the centers of galaxies, are likely also the reason that we do not see an environmental influence on Σ_{H I} at high values of Σ_* compared to observations.

It is important to assess what the role of pixels undetected in gas have on (i) the average behavior of the observations and (ii) the comparison with TNG described above. To this end, we compare gas-detected pixels in both the observations and simulations in the right-hand panels of Figure 2. We apply a lower limit for Σ_{H I} and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ equal to the minimum value in the respective axes in the figure to emulate a detection threshold for TNG. While the median Σ_*–${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ relations for ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ detections in TNG50 and observations are in reasonable agreement, the systematic deficit in TNG's Σ_{H I} relative to observations is even more stark than before. The latter suggests that TNG galaxies are already suppressed in Σ_{H I} at high Σ_* before falling into the cluster. There is minimal suppression in their central Σ_{H I} thereafter, while the observations instead show a clear suppression. However, because the TNG galaxies are selected to be star forming, and only 41 of those meet our cluster criteria, this does not automatically mean the correct qualitative effect of environment on kiloparsec-scale gas is absent in the simulation, as we explore in Section 3.2.

The ensemble relations above demonstrate the systematic underestimation of Σ_{H I} at fixed stellar density in TNG50 with respect to the observations. But ensemble relations are merely the superposition of individual relations. Thus, to establish if TNG50 can reproduce the resolved effects of environmental mechanisms on individual galaxies' gas content reported by Watts et al. (2023), we need to explore their individual Σ_*–Σ_{H I} and Σ_*–${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ sequences.

In a similar fashion to Watts et al. (2023), we rank-order our TNG50 cluster sample by H i deficiency (H i-def for short) and show their individual Σ_*–Σ_{H I} and Σ_*–${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ sequences in Figures 3 and 4, respectively. H i deficiency is often used as a proxy for how processed a galaxy has been by its environment, provided H i-def >0.3. Here, H i-def is defined as the logarithmic deviation in expected H i mass at fixed stellar mass, relative to the median of TNG50 field galaxies (shown in Figure 1). ¹⁹ Each pixel for each TNG50 galaxy in Figures 3 and 4 is shown along with their running median. The running median for the whole TNG50 field sample (a "control") is shown along with one analog VERTICO galaxy, matched according to its stellar mass and distance from the (integrated) star-forming main sequence. For the real galaxies, we use the main sequence of Leroy et al. (2019), which we write in terms of specific star formation rate:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{\mathrm{sSFR}}_{\mathrm{MS}}^{{\rm{L}}19}}{{\mathrm{yr}}^{-1}}\right)=-0.346\,{\mathrm{log}}_{10}\left(\displaystyle \frac{{m}_{* }}{{M}_{\odot }}\right)-6.652.\end{eqnarray} \tag{ 1 }$$

For TNG galaxies, we iteratively perform a least-squares straight-line fit to ${\mathrm{log}}_{10}({m}_{* })$–${\mathrm{log}}_{10}(\mathrm{sSFR})$, removing outliers of >3σ each time to ensure we do not accidentally include quenched galaxies, until the fit converges. The resulting fit is

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{\mathrm{sSFR}}_{\mathrm{MS}}^{\mathrm{TNG}}}{{\mathrm{yr}}^{-1}}\right)=-0.159\,{\mathrm{log}}_{10}\left(\displaystyle \frac{{m}_{* }}{{M}_{\odot }}\right)-8.427.\end{eqnarray} \tag{ 2 }$$

The vector length between the TNG50 and VERTICO matched pairs is always shorter than 0.3 dex in this parameter space. Although H i deficiency and distance from the main sequence are physically correlated, we note that this method does not mean that the H i deficiencies of the paired VERTICO and TNG50 galaxies are the same. There are 13 galaxies in the TNG50 cluster sample that both have H i-def >0.3 dex and have VERTICO analogs, of which the nine with the highest H i-def are shown in Figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 demonstrates that when examining the most environmentally affected TNG galaxies, the simulation does reproduce the qualitative results of VERTICO per Watts et al. (2023), i.e., that (i) gas disks are truncated to higher Σ_* for higher H i-def and (ii) central Σ_{H I} decreases for high H i-def. This did not appear to be the case in Figure 2 because 28 of the 41 galaxies in the TNG50 cluster sample are either H i-normal, i.e., they have H i-def < 0.3, and/or have no (unique) analog in VERTICO.

Figure 4 shows the Σ_*–${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ sequences for the same galaxies as in Figure 3. We see truncation occurring at the same Σ_* in both Σ_{H I} and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$, and we see a suppression in central ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}$ for the most H i-deficient galaxies. These results are again qualitatively in line with that reported for VERTICO in Watts et al. (2023).

The fact that the majority of galaxies in the TNG50 cluster sample are unaffected by their environment (i.e., 28 of the 41 galaxies are H i-normal, a much larger fraction than in the observations) is because many TNG galaxies that are strongly affected by a cluster environment belong to a subhalo that is devoid of gas cells entirely, ²⁰ especially at low stellar masses (see Section 5.2 of Diemer et al. 2019; Figure 9 of Stevens et al. 2021). Naturally, such galaxies, where the environmental influence on gas content is greatest, cannot be included in this work. While the improved resolution of TNG50 relative to earlier simulations in the suite certainly combats this issue, numerical simulations are always limited by discretization. Only with enough dense gas elements are the hydrodynamical forces of environment reliable in the simulation. In essence, this means there is a narrow window of opportunity to catch simulated galaxies experiencing the onset of environmental effects. No analogous limitation exists for observed galaxies. This might explain why the TNG50 cluster sample is small in number, despite having two clusters of comparable mass to Virgo. In future work, this issue can be mitigated by using more snapshots to catch the moment of interest for each infalling satellite galaxy.