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

It has been shown in previous publications that the TNG100 simulation quantitatively reproduces the observed reduction in each of the total atomic and total molecular hydrogen gas for galaxies within massive halos, i.e., dense environments. In this Letter, we study how well TNG50 reproduces the resolved effects of a Virgo-like cluster environment on the gas surface densities of satellite galaxies with m * > 109 M ⊙ and star formation rate > 0.05 M ⊙ yr−1. We select galaxies in the simulation that are analogous to those in the HERACLES and VERTICO surveys and mock-observe them to the common specifications of the data. Although TNG50 does not quantitatively match the observed gas surface densities in the centers of galaxies, the simulation does qualitatively reproduce the trends of gas truncation and central density suppression seen in VERTICO in both H i and H2. This result promises that modern cosmological hydrodynamic simulations can be used to reliably model the post-infall histories of cluster satellite galaxies.


INTRODUCTION
In recent years, the realism of cosmological hydrodynamic simulations has grown to the point where the empirical effects of galaxy environment on the global gas properties of galaxies are demonstrably reproducible at low redshift (e.g.Stevens et al. 2019bStevens et al. , 2021;;Yun et al. 2019).Indeed, this is also true for some semi-analytic models of galaxy formation (e.g.Stevens & Brown 2017;Xie et al. 2020).These tests of forefront models in the literature have become possible thanks to statistically significant and representative observational surveys that have traced the emission from both the atomic and molecular gas in galaxies that cover the spectrum of environments from field isolation to massive galaxy clusters at low redshift (e.g.Catinella et al. 2018).
While the broad effects of environment on galaxies' gas is now reasonably well understood (see the review by Cortese et al. 2021), surveys have started directing their attention to the gas properties on local scales inside galaxy discs.The 'Virgo Environment Traced In CO' survey (VERTICO, Brown et al. 2021) has observed the molecular hydrogen gas content [H 2 , traced via the CO(2-1) line] at sub-kpc resolution across 51 late-type galaxies in the Virgo cluster, all of them with existing multi-wavelength maps of atomic hydrogen gas (H i), stellar emission, and star formation activity.Recently, Watts et al. (2023) used this dataset to show that, as galaxies are processed by the cluster, environmental mechanisms drive a continuous decrease in both H i and H 2 local surface densities (Σ H i and Σ H2 , respectively) at fixed local stellar surface density (Σ * ) with respect to the field.From a theoretical perspective, it is clearly an important next step to establish if the latest simulations are capable of reproducing (i) the scaling relations of Σ H i and Σ H2 values with Σ * observed in nearby field and cluster populations, and (ii) the systematic influence of environment upon those relationships.
Taking that step in this Letter, we investigate the resolved gas scaling relations of individual galaxies from the TNG50 simulation (Nelson et al. 2019a;Pillepich et al. 2019), mock-observing them to be directly comparable to VERTICO and the 'Heterodyne Receiver Array CO Line Extragalactic Survey' (HERACLES; Leroy et al. 2009).

Observations
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 SFR.The H iobservations for VERTICO and HERACLES galaxies are drawn from the VLA Imaging of Virgo in Atomic gas survey (VIVA; Chung et al. 2009) and The HI Nearby Galaxy Survey (THINGS; 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 (z0MGS; 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) andJiménez-Donaire et al. (2023), respectively.VERTICO is typically 2-3 times more sensitive than HERACLES in both H i and H 2 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 9 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 resolutionmatched 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 −1 .
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, Jiménez-Donaire et al. (2023) for the star formation rates, and 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 2 .

Simulations
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 (AGN), and more, all within a ΛCDM cosmological, magnetohydrodynamic framework with the movingmesh 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 4 M ⊙ .Galaxy subhalos in TNG are identified with subfind (Springel et al. 2001;Dolag et al. 2009).
The H i and H 2 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 2 properties of TNG100 galaxies and their trends with environment have been explored in depth (Diemer et al. 2019;Stevens et al. 2019bStevens et al. , 2021)), but the same is not true in the literature 1 for TNG50 (but see Boselli et al. 2023).In brief, after mock-observing the simulated galaxies to match survey specifications, Stevens et al. (2019bStevens et al. ( , 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 2 fractions of the simulations in Fig. 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 9 M ⊙ and SFR ≥ 0.05 M ⊙ 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 (cf. the criteria in Nelson et al. 2019b).This totals 2479 galaxies from TNG50.Other than a small systematic increase in H 2 fraction for TNG50, the two simulations are well aligned.There is generally good agreement between global H i and H 2 gas fractions in TNG100 and low-redshift observations (e.g., ALFALFA, xGASS, and xCOLD GASS; Diemer et al. 2019;Stevens et al. 2019a), 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 10 M ⊙ .Due to the small number statistics in this regime, we caution against over interpretation, yet it is possible that this is due to a systematic difference 1 Figures published in Diemer et al. (2019) with TNG100 and TNG300 have been reproduced with TNG50 at http://www.benediktdiemer.com/data/hi-h2-in-illustris/.
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.We create mock surface density maps of H i, H 2 , and stellar mass for each TNG galaxy in our sample.These maps are matched to the specifications of the HERA-CLES 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 sub-samples 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 −2 per galaxy.The 'cluster' sample, comparable to VERTICO, contains only satellites (galaxies that are not centrals) in the simulation's two clusters 2 with M 200c ≳ 10 14 M ⊙ .The most massive cluster, with M 200c = 1.8×10 14 M ⊙ , is very similar in mass to Virgo (1.4-4.2×10 14M ⊙ ; 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 13 M ⊙ ), but has a different dynamical state, and has already been justified as a Virgo analogue by Joshi et al. (2021).The cluster sample totals 41 galaxies with an average of 1725 pixels with Σ * > 1 M ⊙ pc −2 per galaxy.The higher pixel count per galaxy in the cluster sample reflects its higher average stellar mass (cf. the distribution of TNG50 points in Fig. 1).

RESULTS
We explore how the resolved gas-stellar surface density scaling relations of galaxies are affected by environment, which we do in two steps.First, we compare sample-averaged scaling relations between the field and Figure 1.The H i (top panel) and H2 (bottom panel) fractions of TNG50 and TNG100 galaxies as a function of stellar mass at z = 0.Only galaxies with m * ≥ 10 9 M⊙ and SFR ≥ 0.05 M⊙ yr −1 included (without any environmental sub-sampling).Hex bins show the number density of TNG50 galaxies.Lines are the running median (thick) and 16th and 84th percentiles (thin) for TNG50 (solid) and TNG100 (dashed).Points with approximate errors compare the VERTICO and HERACLES galaxies that we use in this paper (a subset of the full surveys; cf.fig. 1 of Zabel et al. 2022).We show these observational data for reference, but we do not necessarily expect the simulation medians to align closely with them (but they should be closer to HERACLES than VERTICO).Individual points from our TNG50 cluster and field samples, shown, can be respectively compared to VERTICO and HERACLES.
cluster samples in TNG50 to those of HERACLES and VERTICO.Second, we present the resolved relations of individual cluster galaxies, assessing how the relations change with global H i deficiencies, and comparing analogous galaxies between TNG50 and VERTICO.

Sample averages
Our first aim is to see if the average behavior of kpcscale gas in TNG50 reflects that of reality in both the field and cluster environments.In Fig. 2, we present the pixel-based relations of Σ H i and Σ H2 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 HER-ACLES and VERTICO.To demonstrate that the preference of HERACLES galaxies to have m * ≳ 10 10 M ⊙ (seen in Fig. 1) does not affect our comparison, we add results to Fig. 2 for a TNG50 field sub-sample where we have excluded galaxies with m * ≤ 10 10 M ⊙ .We exclude the resolved molecular Kennicutt-Schmidt relation (Σ H2 -Σ SFR ) from Fig. 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 (Σ * -Σ H2 , addressed further in Brown et al. 2023).See Motwani et al. (2022) for further analysis on kpc-scale Σ SFR in TNG50.
The left column of panels in Fig. 2 accounts for all pixels with Σ * detections, regardless of whether they were detected in H i or H 2 .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 Σ H2 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 Σ H2 .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 Σ H2 become increasingly divergent between the TNG50 cluster and field samples at Σ * ≲ 30 M ⊙ pc −2 .This average behavior is typical of disc truncation, which one would Shaded regions cover the 16th to 84th percentiles (not shown for the high-mass field sub-sample).The left column accounts for all pixels in both the observations (provided they were detected in stellar emission) and simulations, setting non-detections in either Σ H i or ΣH 2 in the observations to zero, which are accounted for in the percentiles.The right column removes any non-detections in the observations by cutting out any pixels that would fall below the axes as plotted.The lower boundary of each axis represents the 1st percentile of all gas-detected pixels (irrespective of whether the pixel is a detection in Σ * ) across both observational surveys.This boundary also represents the cut in gas surface density applied to TNG50 in the right-hand panels, as to emulate a detection threshold.The top two panels show the one-dimensional histograms of Σ * for pixels in each sample, normalised by the number of galaxies in that sample.The y-axis in the top-right panel is stretched by a factor of two for clarity.
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) describe 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. 2019bStevens et al. , 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 Fig. 2.This implies that galaxy centers are generally under-dense 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 Fig. 2. The low gas densities, combined with this systematic underestimation of Σ * in the centers of galaxies, is 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 Fig. 2. We apply a lower limit for Σ H i and Σ H2 equal to the minimum value in the respective axes in the figure to emulate a detection threshold for TNG.While the median Σ * -Σ H2 relations for Σ H2 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 kpc-scale gas is absent in the simulation, as we explore in the next subsection.

Individual cluster galaxies
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 Figure 3. Individual H sequences for the nine most H i-deficient galaxies among VERTICO analogues in TNG50.Points are pixels from the TNG50 galaxies, with thick solid lines the running median of those points.The thin, solid, red line that repeats in each panel is the median for the TNG50 field sample.Each dashed line is the running median for the VERTICO galaxy that the TNG50 galaxy is matched to, based on its stellar mass and distance from the star-forming main sequence.NGC4533 lacks any detected resolved H i. The lower bound of the y-axis in each panel is ∼0.5 dex lower than what is detected in VERTICO galaxies.
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 Σ * -Σ H2 sequences.In a similar fashion to Watts et al. (2023), we rank order our TNG50 cluster sample by H i deficiency (H idef for short), and show their individual Σ * -Σ H i and Σ * -Σ H2 sequences in Figs 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 Fig. 1).3Each pixel for each TNG50 galaxy in Figs 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 analogue 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: For TNG galaxies, we iteratively perform a least-squares straight-line fit to log 10 (m * )-log 10 (sSFR), removing outliers of > 3σ each time to ensure we do not accidentally include quenched galaxies, until the fit converges.The resulting fit is 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 analogues, of which the nine with the highest H i-def are shown in Figs 3 and 4. Fig. 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 discs 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 Fig. 2 because 28 of the 41 galaxies in the TNG50 cluster sample are either H i-normal, i.e. they have H idef < 0.3, and/or have no (unique) analogue in VER-TICO.
Fig. 4 shows the Σ * -Σ H2 sequences for the same galaxies as in Fig. 3.We see truncation occurring at the same Σ * in both Σ H i and Σ H2 , and we see a suppression in central Σ H2 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,4 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.
4. SUMMARY TNG50 quantitatively reproduces the Σ * -Σ H2 relation found in observations in both cluster and field samples.The Σ * -Σ H i relation, on the other hand, is found to be endemically gas-poor at fixed Σ * with respect to the observations.In addition, we find that the kpc-scale effects of a Virgo-like environment on satellite galaxies' H i and H 2 gas content is qualitatively recovered by TNG50 at z = 0. Gas discs are not only truncated more the more they have been affected by their environment, but their central gas densities are also relatively suppressed.However, this effect is likely quantitatively weaker than in reality, because the central gas surface densities of TNG50 galaxies, particularly in H i, are systematically low relative to observations, irrespective of environment.
With this baseline performance of TNG50 established, it opens the door to using the TNG simulations to show how the resolved gas scaling relations of galaxies change after infall into a cluster or otherwise dense environment.Such an experiment is crucial to reinforce the theoretical interpretation of the empirical results of VER-TICO, which must rely on conjecture in describing cluster galaxies' histories.

Figure 2 .
Figure2.Resolved scaling relations for TNG50 galaxies per their field and cluster samples (and a high-mass field sub-sample), compared respectively with HERACLES and VERTICO.Lines are running medians in 0.2-dex bins of Σ * .Shaded regions cover the 16th to 84th percentiles (not shown for the high-mass field sub-sample).The left column accounts for all pixels in both the observations (provided they were detected in stellar emission) and simulations, setting non-detections in either Σ H i or ΣH 2 in the observations to zero, which are accounted for in the percentiles.The right column removes any non-detections in the observations by cutting out any pixels that would fall below the axes as plotted.The lower boundary of each axis represents the 1st percentile of all gas-detected pixels (irrespective of whether the pixel is a detection in Σ * ) across both observational surveys.This boundary also represents the cut in gas surface density applied to TNG50 in the right-hand panels, as to emulate a detection threshold.The top two panels show the one-dimensional histograms of Σ * for pixels in each sample, normalised by the number of galaxies in that sample.The y-axis in the top-right panel is stretched by a factor of two for clarity.

Figure 4 .
Figure 4. Same as Fig. 3 but for the resolved molecular-gas main sequence.