Molecular hydrogen in IllustrisTNG galaxies: carefully comparing signatures of environment with local CO & SFR data

The focus of this paper is on the results of the TNG100 simulation (Pillepich et al. 2018b; Nelson et al. 2018, 2019a; Marinacci et al. 2018; Naiman et al. 2018; Springel et al. 2018), part of the IllustrisTNG magnetohydrodynamic, cosmological simulation suite. The simulation was run with the arepo code (Springel 2010), and follows the standard $\Lambda$CDM cosmological model, with parameters based on the Planck Collaboration XIII (2016) results: $\Omega_{m}\!=\!0.3089$, $\Omega_{\Lambda}\!=\!0.6911$, $\Omega_{b}\!=\!0.0486$, $h\!=\!0.6774$, $\sigma_{8}\!=\!0.8159$, $n_{s}\!=\!0.9667$. The TNG model includes a range of sub-grid prescriptions that are designed to accommodate relevant astrophysical processes in the formation and evolution of galaxies: for example, gas cooling, star formation, massive black-hole growth, and feedback from both stars and active galactic nuclei (Weinberger et al. 2017; Pillepich et al. 2018a). These are based on the earlier Illustris simulation’s sub-grid models (see Vogelsberger et al. 2013, 2014a, 2014b; Genel et al. 2014; Torrey et al. 2014). TNG100 employs a periodic box of length $75\,h^{-1}\!\simeq\!110\,{\rm cMpc}$, containing $1820^{3}$ dark-matter particles of mass $7.5\!\times\!10^{6}\,{\rm M}_{\odot}$, and $1820^{3}$ initial gas cells of typical mass $\sim$$1.4\!\times\!10^{6}\,{\rm M}_{\odot}$. The minimum gravitational softening scale for gas — achieved only in galaxy centres — is 0.19 kpc.

Gravitationally bound structures in the simulations were found with the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). These ‘subhaloes’ exist within haloes that are built with a friends-of-friends (FoF) particle group finder. The most massive subhalo in a FoF group is dubbed the ‘central’, and the remainder are ‘satellites’. We follow the nomenclature of S19 when describing the integrated properties of galaxies in the simulation. That is, ‘inherent’ galaxy properties are defined as the subset of particles/cells in a subhalo that also lie within the spherical aperture given by the ‘baryonic mass profile’ criterion (‘BaryMP’ for short) of Stevens et al. (2014). ‘Mock’ properties depend on the dataset being compared to. In this paper, these relate to the xCOLD GASS survey and Sloan Digital Sky Survey (SDSS — York et al. 2000). We describe the mocking process for this in Section 2.4.

All gas cells in the simulation are post-processed in order to calculate their mass fractions in the form of atomic and molecular hydrogen. The method is described in S19, which builds from Lagos et al. (2015) and Diemer et al. (2018). In this paper, we follow the same three prescriptions used by S19, namely those based on the works of Gnedin & Kravtsov (2011, hereafter GK11), Krumholz (2013, hereafter K13), and Gnedin & Draine (2014, hereafter GD14). All three prescriptions rely principally on the local UV field strength (which we model for each galaxy in three dimensions — see Diemer et al. 2018 for details), local metallicity (a proxy for dust), and local surface density (translated from volumetric gas densities through the Jeans approximation) to calculate the H i/H₂ mass ratio of each gas cell.

For TNG100 results at $z\!=\!0$, the same galaxy sample analysed in S19 is analysed in this paper. That is, only subhaloes with $m_{*}\!\geq\!10^{9}\,{\rm M}_{\odot}$ and dark-matter fractions above 5% are included.

The xCOLD GASS^†††eXtended Carbon monOxide Legacy Database for the Galaxy evolution explorer Arecibo ‘Sloan digital sky survey’ Survey program comprises two observational surveys of the CO(1–0) line emission in galaxies with the IRAM^‡‡‡Institute for Radio Astronomy in the Millimetre range 30-m telescope. The surveys were designed to measure (or place upper limits on) the H₂ content of a representative sample of galaxies over $>$2 dex in stellar mass. The original COLD GASS survey (Saintonge et al. 2011) targeted 350 galaxies with $m_{*}\!\geq\!10^{10}\,{\rm M}_{\odot}$ in the redshift range $0.025\!\leq\!z\!\leq\!0.05$. These galaxies were a subset of those observed in 21-cm H i emission in the GASS survey (Catinella et al. 2010). Each galaxy was observed sufficiently long to either detect the CO(1–0) line or place an upper limit of $\sim$$0.015\,m_{*}$ on $m_{\rm H_{2}}$. The second survey — COLD GASS-low (Saintonge et al. 2017) — observed galaxies with $10^{9}\!<\!m_{*}/{\rm M}_{\odot}\!<\!10^{10}$ at $0.01\!\leq\!z\!\leq\!0.02$, with detections in CO for $m_{\rm H_{2}}\!\gtrsim\!0.025\,m_{*}$. Most galaxies in COLD GASS-low, but not all, also have H i data from GASS-low (Catinella et al. 2018). In total, xCOLD GASS comprises 532 galaxies, 477 of which are also in xGASS (GASS + GASS-low). As per Catinella et al. (2018), we refer to the subsample of galaxies in both surveys as ‘xGASS-CO’. All H₂ masses from xCOLD GASS shown in this work assume the CO-to-H₂ conversion factor of Accurso et al. (2017), which depends on both the metallicity and star formation rate (SFR) of a galaxy.

All xCOLD GASS data were sourced from the publicly available catalogue. We note that tags for galaxies being centrals and satellites are not included in this dataset, but they are for xGASS (see Section 2.3). Where we present centrals and satellites for xCOLD GASS, we therefore use the xGASS-CO subset. We also note that H₂ masses given in the public dataset and presented in Saintonge et al. (2017) are actually $1.36\,m_{\rm H_{2}}\!\simeq\!m_{\rm H_{2}}\!+\!m_{\rm He}$ (see section 1 of Accurso et al. 2017). In this work, $m_{\rm H_{2}}$ means the mass of molecular hydrogen only. We therefore subtract the helium contribution from all xCOLD GASS measurements.

Measured fluxes must undergo ‘beam corrections’ (often also referred to as ‘aperture corrections’ in the literature) to estimate the intrinsic flux of an object. The beam corrections assume that (i) the face-on, two-dimensional H₂ half-mass radius, $r_{\rm H_{2}}^{\rm half}$, and half-SFR radius, $r_{\rm SFR}^{\rm half}$, are equivalent for each galaxy (where the latter is measured independently); (ii) the H₂ surface density of every galaxy falls exponentially with radius; and (iii) the IRAM beam shape is a Gaussian. We can analytically write the beam correction factor as

where $i$ is the galaxy’s inclination, and ‘FWHM’ refers to the full width at half maximum of the beam — often simply referred to as the ‘beam size’ — where the beam size of IRAM is 22 arcsec at the frequency of the CO(1–0) observations. Further details regarding the beam correction method can be found in section 2.2 of Saintonge et al. (2012). Both the beam-corrected and uncorrected H₂ masses are used in this work. For a given galaxy, we assume the conversion factor from CO luminosity to $m_{\rm H_{2}}$ is the same for the corrected and uncorrected quantities.

Because xCOLD GASS is not complete in terms of stellar mass (even though it is representative), whenever we calculate statistics such as running percentiles or means, we assign weights to each galaxy to compensate for the relative galaxy count in adjacent 0.1-dex slices of stellar mass. The expected frequency is based on the ratio of galaxy counts in TNG100 for those same bins. This method follows S19, which is very similar to that used in Catinella et al. (2018).

We note that xCOLD GASS was not designed for environmental galaxy studies. In truth, neither the sample size nor detection threshold are ideal for this use. Nevertheless, there is a sufficiently representative range of environments covered by the two surveys (see the bottom panels of Fig. 1), and they presently remain unrivalled in their combined number counts and depth of CO surveys at $z\!\simeq\!0$. While we stress that we push these data to their usage limit, xCOLD GASS persists as the most relevant dataset in the literature for a point of comparison for this paper.

All stellar masses and SFRs from observational data used in this paper (i.e. for xCOLD GASS and the additional SDSS sample we describe below) originate from the MPA–JHU^§§§Max Planck institute for Astrophysics–Johns Hopkins University catalogue from SDSS Data Release 7 (Abazajian et al. 2009). These properties have been adjusted for a Chabrier (2003) stellar initial mass function (IMF) and Planck Collaboration XIII (2016) cosmology as necessary (to be consistent with the TNG model). Stellar masses are based on fits to each galaxy’s spectral energy distribution, as derived from SDSS’s broadband photometry (Salim et al. 2007; also see Kauffmann et al. 2003). SFRs rely on fibre spectroscopy, but have had a correction applied for the small fibre aperture (Salim et al. 2007). SFRs for actively star-forming galaxies use H$\alpha$ emission as their primary indicator, which canonically probe historical time-scales of order $10^{7}$ yr. For passive galaxies, estimated SFRs are related to the 4000-$\AA$ break (Brinchmann et al. 2004), which probes an average time-scale of order $10^{9}$ yr. Central/satellite assignments and host halo masses are based on the Yang et al. (2007) group finder. For xGASS, this includes tweaks outlined in Janowiecki et al. (2017). Otherwise, we use the ‘modelB’ catalogue.

In Section 4.4 and Appendix A.1, we use a subsample of galaxies from SDSS that is much larger than xCOLD GASS to look at environmental effects on galaxies’ specific star formation rates (sSFRs). This sample is volume-limited over the redshift range $0.02\!\leq\!z\!\leq\!0.05$, completely covers stellar masses from $10^{9}$ to $10^{11.5}$ M_⊙, and only includes galaxies with SFR_FLAG_tot = 0. This is taken from the sample used in Brown et al. (2017), which was also used in S19.

Throughout this paper, we adhere to the following conventions.

• •

  Whenever we refer to an ‘H₂ fraction’, we exclusively mean the ratio of H₂ mass to stellar mass of a galaxy (and likewise for ‘H i fractions’).

• •

  All cosmology-dependent quantities — both simulated and observational — assume (or have been adjusted to assume) $\Omega_{m}\!=\!0.3089$, $\Omega_{\Lambda}\!=\!0.6911$, and $h\!=\!0.6774$ (Planck Collaboration XIII 2016).

• •

  Likewise, all IMF-dependent quantities assume a Chabrier (2003) IMF. When converting from a Kroupa (2001) IMF, we assume $m_{\rm*,Chabrier}\!=\!\frac{0.61}{0.66}m_{\rm*,Kroupa}$ and ${\rm SFR}_{\rm Chabrier}\!=\!\frac{0.63}{0.67}\,{\rm SFR}_{\rm Kroupa}$ (see Madau & Dickinson 2014).

• •

  Zeroes are never discarded when calculating statistical quantities such as means, medians, and other percentiles.

• •

  Logarithms are always taken after the statistical quantity of interest is computed.

• •

  Upper-case $R$ refers to three-dimensional/spherical radii. Lower-case $r$ refers to two-dimensional/cylindrical radii.

• •

  Upper-case $M$ (or $\mathcal{M}$) refers to the mass of haloes. Lower-case $m$ refers to the mass of galaxies.

• •

  When binning data, we use bins of nominal (but not strict) width 0.2 dex. The absolute least number of galaxies we use in a bin is 15. No galaxies are discarded when binning. When necessary, bin widths automatically adapt to ensure both latter criteria are fulfilled.

• •

  There are no ‘helium contributions’ to H i or H₂ masses.

• •

  Whenever we refer to ‘halo mass’, we always mean that of the ‘host’ or ‘parent’ halo. We use these terms interchangeably.

In the top two panels of Fig. 3, we show the variation in the H₂ fraction of galaxies from both TNG100 and xCOLD GASS as a function of stellar mass, including the relation’s scatter (cf. fig. 3 of Diemer et al. 2019). The top panel makes no correction to the CO observations for the finite beam, and instead applies the same beam to the simulations to derive comparable ‘mock observed’ H₂ masses (as outlined in Section 2.4). Forward-modelled uncertainties on stellar masses are also included. By contrast, the second panel does employ the Saintonge et al. (2012) beam correction to the observations, while the inherent (gravitationally bound) properties of TNG galaxies are plotted. In both cases, non-detections in xCOLD GASS assume their upper-limit measurements for the percentiles shown with hexagons and error bars, and are instead set to zero for the medians shown by the plusses. Loosely, both datasets in both panels suggest that higher-mass galaxies have lower H₂ fractions, with significant scatter at all masses. It is important that we show both these panels, as it highlights the difficulty and care required to compare simulations with observations; the level of agreement is notably different between the two. At the low-mass end, the uncorrected and mock quantities seem better matched. But at the high-mass end, the choice of comparison either leads to the conclusion that TNG100 galaxies are systematically too gas-rich or contrarily that they are suggestively too gas-poor^******The detection fraction of xCOLD GASS is too low to claim this definitively. (modulo systematic uncertainties in xCOLD GASS masses, which are not shown).

Fig. 2 provides a deeper understanding of why there can be such significant differences between the mock and inherent H₂ properties of some TNG100 galaxies, even while other galaxies show negligible differences. The first two columns of images show two galaxies with comparable stellar masses and H₂ fractions, yet vastly different H₂ structure. If the first galaxy were observed face-on, its two-dimensional H₂ half-mass radius would be measured as $r_{\rm H_{2}}^{\rm half}\!=\!2.2$ kpc, well within the mock beam size. In fact, nearly all the H₂ in this galaxy fits within the mock beam. By contrast, the second galaxy has $r_{\rm H_{2}}^{\rm half}\!=\!7.3$ kpc, which exceeds the beam size. The outer half of its H₂ mass also lies in extended spiral arms, all of which is missed by the beam. The next three galaxies all have higher stellar masses, each closer to $10^{11}\,{\rm M}_{\odot}$. Again, these highlight a variety of situations where anywhere from 35% to 88% of the galaxies’ H₂ is missed by the mock beam. The most extreme of these three examples is actually representative of the average galaxy at this stellar mass, as discerned from the top-middle panel of Fig. 2.

All of this motivates the question: is one form of comparison (mock/uncorrected or inherent/corrected) more meaningful than the other? An argument can be made that, because the beam correction for xCOLD GASS requires a model for how H₂ in galaxies should be distributed (i.e. that it is distributed exponentially, with a half-mass radius equal to the radius enclosing half the galaxy’s star formation rate), the mock/uncorrected comparison is the most direct. Unfortunately, this implicitly relies on H₂ being distributed in TNG galaxies the same way it is in reality; otherwise the applied beam would exclude a different fraction of the true H₂ mass of the galaxy. As it happens, Diemer et al. (2019, see their fig. 7) have shown that H₂ in TNG galaxies can extend to higher radii than found in observations (cf. Bolatto et al. 2017). While the difference in extent is only of order tens of percent on average, there is a tail of galaxies that stretches out to factors of a few. We have since learned that higher-mass galaxies lie in this tail at greater frequency. Fig. 2 tells a consistent story, where several examples of galaxies with H₂ half-mass radii notably exceeding their stellar half-mass radii are shown. Based on this then, if the goal is to compare the integrated H₂ masses of galaxies, one can make a counter-argument that the inherent/corrected comparison might be more meaningful.

The truth is that neither form of comparison is perfect. But there is a test we can perform to help inform when one comparison might be more meaningful than the other. Namely, we can calculate an ‘effective beam correction factor’ for each TNG100 galaxy using Equation (1), and compare the distribution of these factors with (i) those of xCOLD GASS, and (ii) the ratio of inherent to mock H₂ masses of TNG100. Importantly here, Equation (1) is independent of any actual H₂ properties of the galaxies. We define $r_{\rm SFR}^{\rm half}$ for TNG galaxies using the instantaneous SFRs of their gas cells. If Equation (1) is an accurate method for beam correction, then $f^{\rm beam}_{\rm corr}$ should be almost the same as $m_{\rm H_{2}}^{\rm inherent}/m_{\rm H_{2}}^{\rm mock}$ for each TNG galaxy. To check this, we plot the running distributions of both $f^{\rm beam}_{\rm corr}$ and the ratio of it to $m_{\rm H_{2}}^{\rm inherent}/m_{\rm H_{2}}^{\rm mock}$ in Fig. 4 (both quantities are interpreted to be equal to 1 for galaxies with no star formation activity).

For $m_{*}\!\lesssim 10^{10.3}\,{\rm M}_{\odot}$, the bottom panel of Fig. 4 highlights that Equation (1) does a good job of predicting $m_{\rm H_{2}}^{\rm inherent}/m_{\rm H_{2}}^{\rm mock}$ for the median TNG galaxy. However, the 16th percentiles show that $f^{\rm beam}_{\rm corr}$ can underestimate this ratio by a factor of $\sim$2 at significant frequency. This notably depends on the H i/H₂ prescription; in the case of K13, H₂ does not extend as far out for each galaxy, and thus $f^{\rm beam}_{\rm corr}$ remains accurate within $\sim$30% (still as an underestimate) for most galaxies.^*†*†*†As a point of comparison, the uncertainty taken for beam corrections in Saintonge et al. (2017, section 2.5) is 15%. For the same mass range, the top panel of Fig. 4 shows that the running distributions of $f^{\rm beam}_{\rm corr}$ for TNG100 and xCOLD GASS are fairly similar. In principle, it does not matter if these are quantitatively identical or not; differences between distributions here are the result of differences in the distributions of $r_{\rm SFR}^{\rm half}$:^*‡*‡*‡Per Equation (1), the only other thing that could matter is differences in the distribution of $\sigma_{\rm beam}$ in physical units, but these are the same by construction. the point of the correction being that it accounts for this. This comparison simply tells us that the SFR sizes — and therefore the gas sizes — of TNG100 galaxies are physically reasonable for the majority with $m_{*}\!\lesssim 10^{10.3}\,{\rm M}_{\odot}$. Given this, and the unsurprising fact that $f^{\rm beam}_{\rm corr}$ is not identical to $m_{\rm H_{2}}^{\rm inherent}/m_{\rm H_{2}}^{\rm mock}$, forward-modelling the beam is, objectively, a better way to compare simulations with observations. In other words, the mock/uncorrected comparison should be more meaningful than the inherent/corrected comparison on the basis that it convolves a known gas distribution with a known beam to compare to a directly observed quantity, as opposed to assuming how gas is distributed in xCOLD GASS galaxies to effectively deconvolve the beam. The fact that we find $f^{\rm beam}_{\rm corr}$ to frequently underestimate the beam correction factor explains why TNG100 and xCOLD GASS are better aligned at low masses in the top panel of Fig. 3 than in the middle panel.

However, the situation is not so clear-cut at higher masses. The top panel of Fig. 4 shows wild differences in $f^{\rm beam}_{\rm corr}$ between TNG100 and xCOLD GASS for $m_{*}\!\gtrsim 10^{10.7}\,{\rm M}_{\odot}$. This is most likely symptomatic of gas being unrealistically distributed in the simulated galaxies. Moreover, as seen in the bottom panel of Fig. 3, $f^{\rm beam}_{\rm corr}$ and $m_{\rm H_{2}}^{\rm inherent}/m_{\rm H_{2}}^{\rm mock}$ start diverging at $m_{*}\!>\!10^{10.3}\,{\rm M}_{\odot}$. This implies that the radial gas profiles of these galaxies deviate significantly — i.e. more than one would expect — from a symmetric exponential. The major differences in the mock and inherent H₂ masses seen at these stellar masses are therefore somewhat artificial (or, in other words, not predictive of reality). In light of this, we suggest that the inherent/corrected comparison with xCOLD GASS is a safer choice than the mock/uncorrected comparison at high stellar masses, albeit still an imperfect one.

Given the uncertainties raised above, our philosophy is to present both the mock/uncorrected and inherent/corrected comparisons where relevant throughout this paper. Each requires a grain of salt in its interpretation, but both are needed to provide an as-complete-as-possible picture.

As shown explicitly in the bottom panel of Fig. 3, a significant number of the xCOLD GASS galaxies are not detected in CO; in fact, for the highest-mass bin, the majority are non-detections. This is also reflected in the shaded regions in the top two panels; these cover the area where there are more upper limits in the observational data than detections.^*§*§*§While the survey was designed with an H₂ fraction detection limit (mentioned in Section 2.2), in practice, some non-detections have higher upper limits on their H₂ fractions than other detections’ measurements at the same stellar mass. Because the detection fraction of xCOLD GASS galaxies is $<\!0.84$ in all mass bins, we cannot show the 16th percentile in Fig. 3 for xCOLD GASS H₂ fractions if non-detections are assumed to have no H₂ (those percentiles would all be zero, i.e. $-\infty$ in log space). In essence, this means the lower half of the scatter in the H₂ fraction–stellar mass relation is observationally unconstrained. We therefore cannot make meaningful statements about how representative the frequency of properly H₂-poor galaxies in TNG is of reality.

Nevertheless, we note that the 16th percentile of inherent H₂ fractions of TNG100 galaxies shows a dip over $10^{10.2}\!<\!m_{*}/{\rm M}_{\odot}\!<\!10^{11.3}$, reaching a trough outside the axes of Fig. 3. This lines up with the dip in H i fractions identified in S19. This feature effectively disappears in the mock properties, in part because the H₂ fractions of all the most massive galaxies drop after the application of the relatively small beam. As discussed in S19, this feature is most likely related to the onset of kinetic feedback from massive black holes / active galactic nuclei (AGN — for details on the implementation and result of said feedback, see Weinberger et al. 2017; Zinger et al. 2020, respectively).

In general, there is a slightly greater sensitivity to the choice of H i/H₂ prescription for TNG H₂ results than for H i (cf. Diemer et al. 2018, 2019; S19). This is perhaps unsurprising, as H i tends to dominate over H₂, meaning small fractional changes in H i between prescriptions are reflected by larger fractional changes in H₂ (as the total H i+H₂ is fixed). The source of this difference is most readily understood by looking at the $m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$ ratios of the simulated galaxies, which we assess next.

In Fig. 5, we show the H₂/H i mass ratios of all TNG100 and xGASS-CO galaxies as a function of stellar mass. As documented by Catinella et al. (2018), the H₂/H i mass ratios of xGASS-CO galaxies rise with stellar mass on average. But even over 2.5 dex in stellar mass, this rise is less than (or at least comparable to) the relation’s scatter. We again show two panels that compare the mock-observed simulation properties with uncorrected observations (top) and inherent simulation properties with corrected observations (bottom). Although, we note that no beam correction has been applied to the H i content of xGASS-CO galaxies in the bottom panel of Fig. 5 (it has only been applied for H₂, as given in the public xGASS and xCOLD GASS databases); see S19 for the potential significance of this (much larger) beam. We also need to be particularly careful about the treatment of non-detections; while some galaxies are detected in both H i and CO, some are only detected in one, some the other, and some neither. We still present running percentiles where non-detections (in both phases) are set to their upper limits. But it no longer makes sense to plot the equivalent when non-detections are set to zero; to do this for $m_{\rm H\,{\LARGE{\textsc{i}}}}$ would mean dividing by zero (for the same reason, only TNG100 galaxies with $m_{\rm H\,{\LARGE{\textsc{i}}}}\!>\!0$ are included in Fig. 5). Instead, we show a second set of points (exes) with dashed error bars, which exclusively consider galaxies detected in both phases. This again provides two imperfect means of comparing to simulation results. Somewhat reassuringly though, there is little difference between the two sets of observational points.

TNG100 mock results are broadly consistent with observations in the left half of the top panel of Fig. 5, but are again hampered by the overextension of H₂ at high stellar masses in the simulation. Here, the H₂/H i mass ratios show a systematic drop from $m_{*}\!\gtrsim\!10^{10.3}\,{\rm M}_{\odot}$, before rising again at $m_{*}\!\gtrsim\!10^{11}\,{\rm M}_{\odot}$. While the depth of this feature is tied to our mock procedure, comparison with the bottom panel of Fig. 5 shows that the feature still exists for the inherent galaxy properties, just more weakly. As such, there must be an underlying physical cause for this feature in the simulation. This suggests that gas from the galaxies’ centres — i.e. where the majority of H₂ but only a fraction of H i is (see e.g. Leroy et al. 2008; Obreschkow et al. 2009; Stevens et al. 2019b) — is being removed or depleted. Based on differences in mock/inherent H i masses alone, especially when contrasting satellites and centrals, S19 hypothesized that AGN feedback was likely responsible for displacing gas from galaxies’ centres in TNG100 at the mass scale in question. Having imaged several galaxies, we can confirm there are clear holes and asymmetries in the gas structure of many galaxies at these masses, as would be expected from ejective feedback at each galaxy’s centre (for a related analysis of galactic outflows in the TNG model, see Nelson et al. 2019b). The fourth galaxy in Fig. 2 is a prime example of this. Most of the H₂ on one side of the galaxy has been blown out, with the other side clearly disrupted, as seen by its lack of alignment with the stellar disc. This physical effect is then compounded when mock measurements are made, as it is then only the central region of the galaxy that would contribute meaningfully to the H₂ mass, but that gas has been displaced.

While we cannot claim that TNG100 is a 100% faithful representation of reality when it comes to galaxies’ H₂ content — indeed, we have identified some issues — Figs 3 & 5 at least give us confidence that the simulation is ‘in the right ball park’, so to speak (especially for $10^{9}\!\leq\!m_{*}/{\rm M}_{\odot}\!\lesssim\!10^{10.2}$). Similar statements could be made about the EAGLE^*¶*¶*¶Evolution and Assembly of GaLaxies and their Environments simulations (cf. fig. 5 of Lagos et al. 2015; also see Davé et al. 2020). It will likely be years before next-generation simulations and CO surveys will allow us to make comparative statements that are more grand and robust than this. In the interim, the simulations that we have on hand still provide worthwhile synthetic laboratories to study the effects of environment on H₂, bearing in mind the caveats that carry through from the above.

As outlined in Section 1, one of the main aims of this paper is to investigate how centrals and satellites differ in their H₂ content at fixed stellar mass. In Fig. 6, we show precisely this for xGASS-CO (i.e. the xCOLD GASS galaxies that have central/satellite status determined as part of xGASS, per Section 2.2) and TNG100 at $z\!=\!0$. Similar to Fig. 3, we compare the mock/uncorrected values in the top panels, and the inherent/corrected values in the middle panels. Because we are now subsampling the observational data, we are dealing with weaker statistics. This, combined with there being visually more information in the figure, means we have shown the running medians and 84th percentiles in separate columns. We exclude a column of panels for the 16th percentile here primarily (but not exclusively — see below) due to the low detection fraction in the observations. Importantly, we can trust that the distribution of parent halo masses of satellites in xGASS-CO is very similar to that of TNG100 (Fig. 1), meaning their gas properties at fixed stellar mass should be directly comparable.

Looking at the observational data alone in any of the main four panels of Fig. 6, there is little evidence to support the idea that satellites are significantly depleted of their H₂ at a particular mass scale. But some degree of systematic difference does appear to be present, as 5–6 of the seven stellar-mass bins show satellites to have lower median and 84th-percentile H₂ fractions than centrals. While the bootstrapping errors on these percentiles imply this is not statistically significant for any singular bin, the likelihood of randomly observing this result for all bins simultaneously is low (we elaborate below). Nevertheless, the results of Fig. 6 contrast with the directly comparable plot for H i in xGASS (fig. 5 of S19, ), where there is a clearer distinction between the distribution of H i fractions of satellites and centrals for all fixed $m_{*}\!\lesssim\!10^{10.75}\,{\rm M}_{\odot}$ (see Appendix A.2). We have tested that those H i results remain when we select xGASS-CO galaxies exclusively, implying that the relative lack of observed difference in H₂ for centrals and satellites cannot be fully explained away by the subsample being small and/or potentially unrepresentative.

By contrast, satellites in TNG100 show a clear depletion in their H₂ fractions relative to centrals of the same stellar mass. The inherent results of the simulation find the median satellite to be $\sim$0.6 dex (a factor of $\sim$4) poorer in H₂ than the median central at fixed $m_{*}\!<\!10^{10.3}\,{\rm M}_{\odot}$. From a theoretical, qualitative stand-point, this is not an unreasonable result. One can imagine that satellites’ molecular gas would continue to be consumed in star formation (and accompanied feedback), but would be unable to efficiently replenish, due to the denial of cosmological gas accretion (similar to the picture outlined by Stevens & Brown 2017). If the depletion of H₂ were to be approximately linear with time, we should see a greater difference between centrals and satellites in log space for low percentiles, and likewise a smaller difference at high percentiles. Indeed, there is only a 0.2–0.3 dex difference (a factor of $<$2) in the inherent 84th percentile of centrals and satellites in TNG100 (lower-right panel of Fig. 6). At the same time, for those satellites that experience a pericentric passage close to the halo centre, ram pressure from the relative motion of the intrahalo medium would be so strong that it could entirely remove the interstellar medium (i.e. including H₂). Another reason we opted to exclude panels for the 16th percentile in Fig. 6 is that it equals zero ($-\infty$ in log) for TNG100 satellites $\forall m_{*}\!\geq\!10^{9}\,{\rm M}_{\odot}$.

It is starkly obvious that the separation between centrals and satellites diminishes significantly when we look at the mock-observed simulation outputs. The separation between running medians has been reduced from $\sim$0.6 to $\sim$0.2 dex. From testing (not shown here), we learned that while the application of the mock beam (Section 2.4.1) has an impact on decreasing this separation, the most important factor is the cross-contamination of centrals and satellites (Section 2.4.3). That is, if we apply the contamination to the inherent properties, we find a separation $\sim$0.3 dex. Similarly, had we excluded the cross-contamination from our mock procedure (but kept the other aspects the same), we would have found a separation of $\sim$0.4 dex. Satellite–central cross-contamination is always present in the xGASS-CO data, regardless of whether beam corrections are considered or not. There is no way to correct for this contamination without applying a completely different group-finding algorithm to SDSS; even then, a ‘perfect’ group finder does not exist to our knowledge, nor might it ever. The best option (arguably, the only option) is to forward-model the contamination, like we do with the mock results of TNG100. As such, the mock/uncorrected comparison is always more meaningful than the inherent/corrected comparison when investigating the relative influence of environment, despite the concerns about the absolute mock H₂ properties at high stellar masses raised in Section 3.1.

A natural question to ask is: do the xGASS-CO data favourably support our prediction that the survey should see a 0.2-dex systematic difference between the H₂ fractions of centrals and satellites? Close inspection of the top panels of Fig. 6 suggest that the data are consistent with this prediction. But the uncertainties in each stellar-mass bin for xGASS-CO make it difficult to read how quantitatively precise this apparent consistency is. By combining the data from all bins, however, we can make a more rigorous comparison. To this end, we calculate the difference in H₂ fraction for each galaxy from the median of centrals in each bin, and thus define

In Fig. 7, we show the cumulative distributions of this property for both centrals and satellites from each of TNG100 and xGASS-CO. The top panel shows the predicted distributions from TNG100 to lie almost on top of the uncertainty range of xGASS-CO. A Kolmogorov–Smirnov (KS) test of the observational data is enough to rule out the possibility that satellites and centrals follow the same distribution in $\Delta\log_{10}\!\left(m_{\rm H_{2}}/m_{*}\right)$.^*∥*∥*∥The two-sided, two-sample KS statistic for the two xGASS-CO dot-dashed curves in the top panel of Fig. 7 is 0.31, with a $p$-value of order $10^{-3}$. Not only does this show that the survey is consistent with the prediction of TNG100 regarding the systematic difference between centrals and satellites, but this also means TNG100 predicts very similar scatter in H₂ fractions to what xGASS-CO observes for both satellites and centrals. The bottom panel of Fig. 7 highlights how different the ‘true’ distributions of H₂ fractions are for TNG100 satellites. We include this panel for completeness, but remind the reader that the top panel is the definitive comparison to be made with observations.

As has already been discussed here and in many other works, H₂ in satellite galaxies is expected to be stripped by ram pressure and/or tides systematically later after infall than H i. Based on this idea, one might expect that the H₂/H i mass ratios of satellites should be higher than centrals of the same stellar mass on average (although, this implicitly neglects other environmental and secular effects that might alter the H i and H₂ reservoirs of satellites differently). We therefore compare the H₂/H i mass ratios of TNG100 centrals and satellites with xGASS-CO data in Fig. 8. At first glance, the median $m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$ points appear to be higher for satellites than centrals in xGASS-CO, albeit with potential overlap in some bins. To confirm this, we can follow the same procedure as we did for $m_{\rm H_{2}}/m_{*}$, i.e. we can define a quantity $\Delta\log_{10}\!\left(m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\right)$ that is analogous to Equation (3), swapping all instances of ‘$m_{\rm H_{2}}/m_{*}$’ with ‘$m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$’. While we have omitted a figure analogous to Fig. 7 for $\Delta\log_{10}\!\left(m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\right)$ for brevity, we can again rule out the possibility that centrals and satellites have consistent distributions of $\Delta\log_{10}\!\left(m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\right)$ through a KS test. This holds regardless of whether we include non-detections at their upper limit or exclude them entirely.^*********The two-sided, two-sample KS statistic of $\Delta\log_{10}\!\left(m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\right)$ for xGASS-CO centrals versus satellites when we include (exclude) non-detections in CO and H i is 0.31 (0.17), with a corresponding $p$-value of order $10^{-6}$ ($10^{-2}$). These values are for the uncorrected H₂ masses, but including beam corrections does not lead to significant difference. The data therefore support the theory that satellites’ H₂/H i mass ratios should be elevated relative to centrals.

Contrary to expectation, the mock-observed H₂/H i mass ratios of TNG100 satellites are lower than centrals on average, as seen in the top panel of Fig. 8. This is, however, merely a reflection of the mocking procedure. The applied beam for H i is typically greater in size than the extent of the simulated galaxies’ H i, while the opposite is true for H₂. For satellites in particular, this artificially drives $m_{\rm H\,{\LARGE{\textsc{i}}}}$ up (see S19), and thus drives the $m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$ ratio down. This is clarified by the bottom panel of Fig. 8, where the inherent $m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$ ratios of TNG100 galaxies are shown. For fixed H i/H₂ prescription, satellites in fact have systematically higher inherent $m_{\rm H_{2}}$/$m_{\rm H\,{\LARGE{\textsc{i}}}}$ ratios than centrals of the same mass. This difference is notably greater at $m_{*}\!\lesssim\!10^{10}\,{\rm M}_{\odot}$ than at higher stellar masses, going from $\sim$0.2 dex to effectively zero eventually.

Let us now assess the influence of environment on H₂ in a more quantitative sense. Several metrics exist for quantifying galaxy environment. We use host halo mass as our metric, because (i) this is consistent with our previous work with TNG100 (Stevens et al. 2019a, b) and (ii) observational trends are evidently stronger with this metric than other popular alternatives, such as those that use the distance to the $N$th nearest neighbour (Brown et al. 2017).

In Fig. 9, we break TNG100 satellites into four bins of parent halo mass (in line with S19), plotting the mean H₂ fraction as a function of stellar mass in the top panel. For a given H i/H₂ prescription, there is a clear decline in $m_{\rm H_{2}}$ in higher-mass haloes (denser environments) at fixed stellar mass; satellites hosted in haloes with $M_{\rm 200c}\!\geq\!10^{14}\,{\rm M}_{\odot}$ have a factor of $\sim$10 less H₂ than those in haloes with $M_{\rm 200c}\!<\!10^{12}\,{\rm M}_{\odot}$. This result is very similar to that found for H i in S19 (see their fig. 6), implying that environmental effects in TNG100 have a clear and comparable signature on both atomic and molecular gas.

Despite the signature of environment appearing similar for both H₂ and H i in TNG100, we know from Section 4.1 that the H₂/H i mass ratios of satellite galaxies in TNG100 are raised relative to centrals of the same stellar mass. For this to be consistent, we should expect the average H₂/H i mass ratios of TNG100 satellites to have a residual dependence on parent halo mass as well. Indeed, this dependence is present, as shown in the middle panel of Fig. 9. The quantity on the $y$-axis of this panel has a somewhat subtle but deliberate definition. For satellites in each halo- and stellar-mass bin, we measure the mean H₂ and mean H i mass, take the ratio of those means, log it, and from this subtract the same log ratio for all satellites at the same stellar mass. That is,

This definition removes systematic differences between the H i/H₂ prescriptions (which can already be seen in the top panel) and allows for all galaxies with gas masses of zero to be naturally included [the latter would not be the case were we to plot $\langle m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\rangle$ or ${\rm Median}\!\left(m_{\rm H_{2}}/m_{\rm H\,{\LARGE{\textsc{i}}}}\right)$ instead, for example]. Evidently, parent halo mass only plays a noteworthy role in this quantity for $m_{*}\!<\!10^{10.5}\,{\rm M}_{\odot}$, and its role is rather noisy for $m_{*}\!\gtrsim\!10^{9.7}\,{\rm M}_{\odot}$. The former resembles the fact that high-mass satellites are generally less prone to stripping, as they can often be comparable in size to their central. Even for $m_{*}\!<\!10^{9.5}\,{\rm M}_{\odot}$, $\Delta\log_{10}\!\left(\langle m_{\rm H_{2}}\rangle/\langle m_{\rm H\,{\LARGE{\textsc{i}}}}\rangle\right)$ of satellites in the densest environments differs from those in the lowest-mass haloes by only $\sim$0.2 dex.

What really drives the separation (or lack thereof) in the lines of different halo mass in the top two panels of Fig. 9 is the fraction of galaxies that lack gas entirely. This fraction is plotted in the bottom panel. High-density environments in TNG100 are eventually able to fully remove the gas reservoirs of satellites. An increase in the ‘gas-free’ fraction has a clear signature in mean H₂ mass, while making it difficult for there to be significant differences between H i and H₂. We should be clear that stripping is not actually binary in the simulation. Rather, once the externally forced removal of gas reduces a galaxy to small number of cells, resolution effects take over, dooming the galaxy to lose its remaining few cells. Up until this point, gas should be stripped smoothly.

Previous studies of galaxies’ cold-gas content across environments have typically included much larger samples of observed galaxies than we have presented with xGASS-CO (e.g. Stark et al. 2016; Brown et al. 2017; Stevens & Brown 2017); the important difference is they only had H i data on those galaxies. With only 109 xGASS-CO satellites (including non-detections in CO), it is difficult to split the sample into bins of parent halo mass without losing statistical significance. We therefore leave the results of Fig. 9 as predictions for future observations to test.

A primary utility of cosmological simulations is that galaxies can be tracked through time: something that is observationally impossible for real galaxies. With TNG100, we can directly see how the H i and H₂ reservoirs of galaxies change after infall (i.e. once they have become satellites). To this end, we select a sample of TNG100 satellite galaxies at $z\!=\!0$ with $m_{*}\!<\!10^{10.2}\,{\rm M}_{\odot}$ (avoiding AGN complications) that reside within $1.5\,R_{\rm 200c}$ of their parent halo, have been a satellite for $<\!3$ Gyr, and have an H i mass that is between 0.4 and 2.0 dex below the median for galaxies with the same stellar half-mass radius. This last criterion means there is a nominal expectation that some stripping should have taken place (as the galaxies are ‘H i deficient’ — e.g. Boselli et al. 2014), which we can test, while also ensuring they still have some gas left. The requirement to have only been a satellite for 3 Gyr minimizes any potential changes to the galaxies’ gas properties associated with morphological transformation (see the related discussions of Cortese et al. 2019; Joshi et al. 2020). We track those galaxies back through the baryonic SubLink merger trees (Rodriguez-Gomez et al. 2015) for TNG100, calculating their H i and H₂ content at each snapshot.^{*††*††}*††We note that not all snapshots have the fields saved that we require as input for our gas-phase decomposition methods. We have therefore approximated the values for these missing properties, as described in Appendix B, which we apply to all snapshots for Fig. 10 (exclusively) to ensure we can fairly calculate changes to the satellites’ gas content over time. The uncertainty introduced by these approximations is entirely negligible.

In the top panel of Fig. 10, we show the ratio of the sample satellites’ H i mass at infall to that at redshift zero, versus the same ratio for H₂. Compared in the bottom panel are three example evolutionary tracks of galaxies from their time of infall. We define the ‘time of infall’ here as that of the first snapshot at which the galaxy was both a satellite according to subfind and within $1.5\,R_{\rm 200c}$ of its host halo. We chose $1.5\,R_{\rm 200c}$ rather than $1.0\,R_{\rm 200c}$ because environmental effects are not confined to $R_{\rm 200c}$; in fact, we could have selected an even larger radius, as environmental effects can extend out to $5\,R_{\rm 200c}$ in clusters (see Bahé et al. 2013; Behroozi et al. 2014; Ayromlou et al. 2019). Galaxies that have lost H i since infall will move to the right, while those that have lost H₂ will move upward.

Fig. 10 makes it clear that there is inhomogeneous evolution in satellites’ H i and H₂ content after infall. A relatively simple example of post-infall evolution is seen by the galaxy represented by hexagons in the bottom panel of Fig. 10; this galaxy maintains its gas for several snapshots before gradually losing H i and H₂ at an increasing rate, moving towards the top-right of the panel. Visual inspection of this galaxy (not shown here) clarifies that it is experiencing ram-pressure stripping. The other two galaxies in the panel do not evolve so simply, reminding us that there are many physical processes that take place after infall, not just environmental effects. The galaxy represented by the diamonds was in the process of experiencing a minor merger during infall. This merger brought gas into the galaxy, causing both $m_{\rm H\,{\LARGE{\textsc{i}}}}$ and $m_{\rm H_{2}}$ to increase, shifting the galaxy to the bottom-left of the panel. After this, its gas started to be stripped, causing a reversal in its parameter-space trajectory. The galaxy represented by the squares began to be stripped after infall, but twice collided with intervening gas in the intrahalo medium (which likely would have otherwise been accreted by the central galaxy of the halo). In the first instance, the galaxy was unable to hold the gas for longer than a single snapshot interval. In the second instance, where the galaxy was closer to pericentre, it successfully accreted the gas it collided with, rejuvenating the galaxy for the better part of a dynamical time, before being stripped once again.

Despite the diversity of post-infall changes to satellites’ gas, an underlying trend presents itself. This is readily identified by the running median in the top panel of Fig. 10, which is shallower than the 1:1 line, highlighting that losses in H i are generally more pronounced than in H₂. In S19, we argued that H₂ in TNG100 satellites might be stripped too efficiently, based on their SFRs being too strongly coupled to their H i content. Here, we revise that interpretation somewhat; the preference for stripping H i appears to be present, but it may well be weaker than in the real Universe.

It is important to emphasize that H₂ need not be lost from satellites specifically by direct stripping. The three examples we described above were selected to be diverse rather than representative. For other galaxies, it could feasibly be that H i is lost from ram pressure, but H₂ is consumed in star formation and simply unable to replenish (e.g. Stevens & Brown 2017). However, the facts that (i) stripping manifests as the complete removal of satellites’ gas in the densest environments in many cases (bottom panel of Fig. 9), (ii) depletion of H i in TNG100 satellites has a parallel effect on SFR (S19), and (iii) hydrodynamical forces in the simulation have no explicit knowledge of whether gas is predominantly atomic or molecular (as the phase decomposition is calculated in post-processing) imply that stripping affects all gas phases in the simulation. Regardless of the physical mechanism(s) at play, Fig. 10 definitively shows us that H i is more readily depleted than H₂ (on average) once galaxies become satellites in TNG100.

A wealth of empirical evidence points to a tight connection existing between the amount of molecular gas in a galaxy and its star formation rate, applying not only to integrated properties (captured by a common H₂ depletion time — e.g. Saintonge et al. 2017; Tacconi et al. 2018), but also to localised surface densities (Bigiel et al. 2008; Leroy et al. 2008). After all, both molecule formation and star formation should occur in relatively dense, cold gas clouds. Observational measurements of galaxies’ star formation rates are orders of magnitude more numerous than those of galaxies’ H₂ content in the local Universe. Comparing the star formation activity of TNG galaxies against those numerous data provides a less direct, but more statistically meaningful test of the simulated galaxies’ molecular-gas content and its relation to galaxy environment.

Because the H i/H₂ modelling was done in post-processing, the star formation law applied in TNG and the H₂ fractions we derive for the simulated galaxies do not have a precise one-to-one connection. Unsurprisingly however, one does still find that TNG galaxies’ star formation rates and H₂ masses are closely correlated (Diemer et al. 2019), as is the case in other cosmological, hydrodynamic simulations (Lagos et al. 2015; Rodríguez Montero et al. 2019; Davé et al. 2020). Comparison with xCOLD GASS shows that the SFR–$m_{\rm H_{2}}$ relationship in TNG is consistent with the real, local Universe (fig. 4 of Diemer et al. 2019). We have further found that there is no environmental dependence on this relationship in TNG100 (not shown here), which is important for the below.

To provide context for the results that follow, we present and discuss the $m_{*}$–sSFR plane of TNG100 galaxies in Appendix A.1. We also refer the reader to the closely related work by Donnari et al. (2020, 2021), where the quenched fractions of TNG galaxies are investigated in detail, including environmental influences.

In the two large panels of Fig. 11, we show running means separately for satellites’ and centrals’ sSFRs as a function of stellar mass, for SDSS and both the mock (left panels, as described in Section 2.4.2) and inherent (right-hand panels) properties of TNG100. The left panels remain the definitive comparison with SDSS though. It is clear that satellites have systematically lower sSFRs than centrals of the same stellar mass (or, rather, a greater fraction of satellites are passive), both in the observed and simulated universes. Despite TNG100 generally having a smaller quenched population at most stellar masses (as can be inferred from Fig. 12) — which leads to differences in the shape of the mean $m_{*}$–sSFR curves between the TNG100 mock predictions and SDSS — both datasets show similar relative behaviour between centrals and satellites. In particular, we can quantify the typical central–satellite separation in $\langle{\rm sSFR}\rangle$ across the full stellar-mass range as 0.19 and 0.15 dex for SDSS and TNG100, respectively.

The key result of this subsection can be found in the bottom-left panel of Fig. 11, where we split satellites by environment (quantified by their parent halo mass), and show the difference between the mean sSFR of those satellites and that of all satellites at the same stellar mass:

At fixed stellar mass, the mean sSFR of satellites has a visible dependence on parent halo mass both in SDSS and TNG100 (for a directly related assessment of SDSS, see Wetzel, Tinker & Conroy 2012). If we are to trust that sSFR is a fair proxy for H₂ fraction, then Fig. 11 provides evidence — in addition to and independent from Fig. 6 — that environment-driven losses of H₂ in TNG are in quantitative agreement with the real Universe. Comparison of the two bottom panels of Fig. 11 highlights how, between satellites hosted in haloes with $M_{\rm 200c}\!<\!10^{12}\,{\rm M}_{\odot}$ and those with $M_{\rm 200c}\!\geq\!10^{14}\,{\rm M}_{\odot}$, a ‘true’ difference in $\langle{\rm sSFR}\rangle$ of almost an order of magnitude can be reduced to a factor of $\sim$2.5 ($\sim$0.4 dex) through the interpretation of and uncertainties associated with observational data. Comparison of the bottom-right panel of Fig. 11 with the top panel of Fig. 9 also highlights that sSFR and H₂ are indeed very similarly affected by environment in TNG100.

Simulation outcomes are always at the mercy of resolution. Of particular relevance here is the effect this can have on gas stripping in satellites, both from tides and ram pressure. For example, a finite resolution limits the accuracy of hydrodynamical-force calculations between the interstellar medium of a satellite and the intrahalo medium it moves through. While idealized simulations with arepo suggest that the moving-mesh hydrodynamics scheme generally captures effects of ram pressure and induced Kelvin–Helmholtz instabilities more accurately than traditional SPH^{*‡‡*‡‡}*‡‡Smoothed-Particle Hydrodynamics codes (Springel 2010; Heß & Springel 2012; Sijacki et al. 2012), the median neutral gas cell number (neutral-gas mass divided by the baryonic mass resolution of $1.4\!\times\!10^{6}\,{\rm M}_{\odot}$) for TNG100 galaxies used in our study is only $\sim$1200. The gas in many of the TNG100 galaxies is therefore more poorly resolved than the regimes of these tests. The high fraction of gas-free satellites in TNG100 discussed in Section 4.2 could potentially be influenced by this at some level.

Furthermore, the potential-well depth of galaxies can be affected by a host of numerical effects, including gravitational softening, the unequal mass partitioning of dark matter and baryons, the imposed density threshold for star formation, and the implementation of feedback (see the recent results of Dutton et al. 2019; Ludlow et al. 2019, 2020). This can artificially alter how susceptible a galaxy is to being stripped. For example, we raised in Section 2.4 that the inner 2 kpc of TNG100 galaxies are effectively ‘unresolved’. Yet we showed a galaxy in Fig. 2 with an H₂ half-mass radius of 2.2 kpc, scarcely larger than the ‘resolution limit’. Fig. 1 also shows that mock-observed H₂ masses can be sensitive to as little as the inner few kpc for a fraction of low-mass galaxies in TNG100. These examples are by no means representative of all galaxies, but they anecdotally highlight that we are pushing TNG100 close to the limit of its capabilities. In these situations, it also becomes especially difficult to capture environment-driven gas compression, where ram pressure might not only strip gas, but also lead to the conversion of H i into H₂ (e.g. Henderson & Bekki 2016), which empirical evidence shows can happen (e.g. Lee et al. 2017a).

To elaborate on the non-trivial topic of resolution-based force limitations in the context of TNG would be worthy of a paper by itself, if not several. Where me may glean some insight regarding the impact of these limitations in future is through applying a similar analysis in this work (and S19) to the higher-resolution TNG50 simulation (Nelson et al. 2019b; Pillepich et al. 2019).

Our post-processing models for performing the H i/H₂ breakdown in TNG rely on common assumptions. Namely, (i) they all assume some degree of equilibrium; (ii) they rely on the same principle ingredients: local density, temperature, metallicity, and UV intensity; (iii) they are inherently two-dimensional models that have been translated to three dimensions using the imprecise Jeans approximation; and (iv) they are instantaneous in their application, meaning it is possible for discontinuities in the H₂/H i mass ratio of a given galaxy to occur from snapshot to snapshot. The lattermost of these could artificially amplify the ‘random walk’ nature of the tracks shown in the bottom panel of Fig. 10. Many of the limitations and assumptions of our H i/H₂ modelling are discussed in greater detail in Diemer et al. (2018).

It should not be forgotten that observational data also carry potentially important limitations too. Most notably, uncertainties in observed H₂ masses are dominated by systematics associated with the CO-to-H₂ conversion factor. While the conversion factor adopted for xCOLD GASS is physically motivated and empirically calibrated (Accurso et al. 2017), it principally depends on metallicity, and secondarily depends on sSFR relative to the main sequence (which traces the UV field strength). Ideally, the conversion would be done on local scales, but it is instead based on the galaxies’ ‘integrated’ properties, where ‘integrated’ really means ‘inferred from the 3-arcsec SDSS fibre’. The CO-to-H₂ conversion therefore implicitly assumes that the gas-phase metallicity and UV field are uniform in a given galaxy, which is almost certainly false. As discussed in Section 2.4.2, the SDSS fibre size is far smaller than the beam size used to measure CO, and smaller still relative to the total extent of the xCOLD GASS galaxies. Moreover, the diagnostics used to infer metallicity (oxygen abundance) from the SDSS emission line ratios carry significant systematic uncertainty themselves (see the review by Kewley, Nicholls & Sutherland 2019). Further systematics lie in the assumption that the gas-phase metallicity traces the dust-to-gas ratio (incidentally, this also applies to the H i/H₂ prescriptions used for the simulation).

We note that we neglected to forward-model any H₂ mass uncertainties in our simulation predictions. Had we done this, the scatter in H₂ fractions in the ‘mock’ results would have increased to reflect the typical statistical error in the xCOLD GASS $m_{\rm H_{2}}$ measurements. Given that systematics are even more important, it matters less that a predicted line goes through the data, and more that it is close and parallel.

In the spirit of forward-modelling, models that can be overlaid on simulations to infer CO emission provide a means to avoid converting CO flux from observations to H₂ mass entirely. While this requires thoughtful consideration of radiative transfer and chemistry to model in detail (see e.g. Glover et al. 2010), even relatively straightforward methods to forward-model CO emission from cosmological-simulation galaxies can lead to different conclusions than H₂-based comparisons with observations (Davé et al. 2020). Some level of systematic uncertainty will likely still be present, but as with the mock-observation strategy employed in this paper, such a method is at least self-consistent and circumvents assumptions regarding structure and gradients in the observed galaxies.