The relationship between cluster environment and molecular gas content of star-forming galaxies in the eagle simulation Open Access

ABSTRACT

We employ the eagle hydrodynamical simulation to uncover the relationship between cluster environment and H₂ content of star-forming galaxies at redshifts spanning 0 ≤ z ≤ 1. To do so, we divide the star-forming sample into those that are bound to clusters and those that are not. We find that, at any given redshift, the galaxies in clusters generally have less H₂ than their non-cluster counterparts with the same stellar mass (corresponding to an offset of ≲0.5 dex), but this offset varies with stellar mass and is virtually absent at M_⋆ ≲ 10^9.3 M_⊙. The H₂ deficit in star-forming cluster galaxies can be traced back to a decline in their H₂ content that commenced after first infall into a cluster, which occurred later than a typical cluster galaxy. Evolution of the full cluster population after infall is generally consistent with ‘slow-then-rapid’ quenching, but galaxies with M_⋆ ≲ 10^9.5 M_⊙ exhibit rapid quenching. Unlike most cluster galaxies, star-forming ones were not pre-processed in groups prior to being accreted by clusters. For both of these cluster samples, the star formation efficiency remained oblivious to the infall. We track the particles associated with star-forming cluster galaxies and attribute the drop in H₂ mass after infall to poor replenishment, depletion due to star formation, and stripping of H₂ in cluster environments. These results provide predictions for future surveys, along with support and theoretical insights for existing molecular gas observations that suggest there is less H₂ in cluster galaxies.

1 INTRODUCTION

Star formation proceeds in cold and dense molecular clouds (McKee & Ostriker 2007; Kennicutt & Evans 2012), where high densities are conducive for processes that protect the gas against ultraviolet (UV) photodissociation (Glover & Clark 2012), such as dust shielding and self-shielding within molecular clouds. Observations show that the star formation surface density is well correlated with the molecular gas surface density within galaxies in the local Universe, both on global (Kennicutt 1998) and local scales (Bigiel et al. 2008). This implies that the evolution of galaxies is regulated by the efficiency with which they convert molecular gas into stars (known as the ‘star formation efficiency’). As a result, one of the primary objectives of extragalactic astrophysics has been to understand the factors that govern the star formation efficiency of galaxies. The two key parameters that are relevant in this regard are the star formation rate (SFR) and molecular gas content of galaxies.

Since molecular hydrogen (H₂) is the most abundant molecule in the interstellar medium (ISM), understanding its evolution is crucial for developing a comprehensive picture of the formation and evolution of galaxies. H₂, however, is not a good observational tracer of molecular gas because it lacks a permanent dipole moment, which renders its dipolar rotational transitions forbidden and precludes observing it in emission. Fortunately, molecular gas clouds also contain carbon and oxygen, which combine to form the CO molecule under the conditions that are prevalent in molecular clouds (Van Dishoeck & Black 1988). This molecule possesses a weak permanent dipole moment that allows for its excitation to higher rotational states, causing transitions like the J = 1 → 0 emission which – owing to its high atmospheric transparency and intensity – is often used as an indirect tracer of H₂.

In practice, the H₂ fraction of a cloud is gauged by multiplying the CO luminosity by a CO-to-H₂ conversion factor (α_CO), which depends on ionisation conditions in the cloud and its environment (Feldmann, Gnedin & Kravtsov 2012; Narayanan et al. 2012; Clark & Glover 2015; Gong, Ostriker & Kim 2018; Hu et al. 2022). Initial conversion factors adopted in the literature were either calibrated to the Milky Way or other galaxies in the local Universe (see the review by Bolatto, Wolfire & Leroy 2013). Some studies have even extended this to z ≈ 3 (e.g. Genzel et al. 2015). However, the factors that govern α_CO are difficult to accurately account for, and the derived H₂ content of observed galaxies can carry significant intrinsic uncertainty. This is an important hindrance in accurate estimation of the H₂ content of metal-poor systems, where much of the carbon is ionised due to significant photodissociation of CO molecules by far-UV photons (Bolatto et al. 2013). CO can also be destroyed through collisions with |$\rm {He}\, \small {\rm II}$| present in cosmic rays (Bisbas, Papadopoulos & Viti 2015). For systems at high redshifts (z > 0.5), there is an additional uncertainty introduced due to the fact that only higher CO rotational transitions are currently observable, which have to be corrected for gas excitation states to obtain the equivalent J = 1 → 0 flux. The resulting uncertainties can be reduced to some extent by using complementary observations of [|$\rm {C}\, \small {\rm I}$|]¹(e.g. Valentino et al. 2018; Heintz & Watson 2020), |$\rm {C}\, \small {\rm II}$| (e.g. Zanella et al. 2018), dust extinction in optical/UV regime (e.g. Concas & Popesso 2019; Yesuf & Ho 2019; Piotrowska et al. 2020), and dust emission (e.g. Scoville et al. 2016), although the latter may not be suitable for high redshifts due to the decline in dust content with metallicity.

It is well established that galaxies in dense environments, like clusters, have different evolutionary tracks compared to their non-cluster counterparts. This occurs because cluster galaxies are subject to mechanisms such as the stripping of gas through ram pressure exerted by the intracluster medium (Gunn & Gott III 1972; Quilis, Moore & Bower 2000), thermal evaporation of gas due to interaction with hot intergalactic gas (Cowie & Songaila 1977), ‘strangulation’ – i.e. the cessation of cold gas replenishment due to the removal of hot circumgalactic gas (Larson, Tinsley & Caldwell 1980; Kawata & Mulchaey 2008), viscous stripping of cold gas due to the flow of hot gas past the galaxy (Nulsen 1982), and tidal effects due to interactions with other galaxies in the cluster or the cluster potential (Moore et al. 1996; Moore, Lake & Katz 1998; Smith et al. 2015).

The effect of these processes on the neutral atomic hydrogen (⁠|$\mathrm{H}\, \small {{\rm I}}$|⁠) of galaxies is well documented (e.g. Yun, Ho & Lo 1994; Verdes-Montenegro et al. 2001; Kenney, Van Gorkom & Vollmer 2004; Cortese et al. 2010; Yoon et al. 2017; Džudžar et al. 2019; Reynolds et al. 2021; López-Gutiérrez et al. 2022). Cluster galaxies, in particular, have been shown to have lower |$\mathrm{H}\, \small {{\rm I}}$| content (Giovanelli & Haynes 1985; Solanes et al. 2001; Cortese et al. 2011; Brown et al. 2017), severely truncated |$\mathrm{H}\, \small {{\rm I}}$| discs (Warmels 1988; Bravo-Alfaro et al. 2000), and disturbed |$\mathrm{H}\, \small {{\rm I}}$| morphologies (Chung et al. 2009; Fumagalli et al. 2014). However, there is no overall consensus regarding the impact of cluster environment on the molecular gas content of galaxies. This is despite a plethora of centimetre and (sub)millimetre observations of galaxies carried out after the advent of telescopes like the Atacama Large Millimeter/Submillimeter Array (ALMA; Wootten & Thompson 2009), Jansky Very Large Array (VLA; Perley et al. 2011), and the IRAM NOrthern Expanded Millimeter Array [NOEMA; an extension of the original Plateau de Bure interferometer (Guilloteau et al. 1992)].

Many early as well as recent studies have claimed that galaxies in clusters contain similar (e.g. Casoli et al. 1991; Boselli, Casoli & Lequeux 1995; Rudnick et al. 2017; Darvish et al. 2018; Wu et al. 2018; Lee et al. 2021), or even more molecular gas (Noble et al. 2017; Hayashi et al. 2018; Tadaki et al. 2019) than those not in clusters. On the other hand, there is also strong support for a deficit of molecular gas in such systems (e.g. Vollmer et al. 2008; Jablonka et al. 2013; Wang et al. 2018; Spérone-Longin et al. 2021; Alberts et al. 2022). However, the interpretation of these results is not so straightforward. Galaxies detected in molecular gas observations must possess a decent amount of H₂ in the first place, and therefore, tend to be star forming (hereafter SF), while the most affected systems are deficient in H₂, form a larger fraction of the cluster population, and evade detection. As such, statistically significant samples of SF systems in clusters spanning a wide range in cluster mass are required for conclusive results. Observational studies focussed on clusters, though, generally suffer from poor statistics and inhomogeneous sampling, which can introduce unintended biases and lead to conflicting results.

It is true that some studies focussing on nearby clusters have shown direct evidence of ram-pressure stripping (Vollmer et al. 2008, 2009; Zabel et al. 2019; Moretti et al. 2020), but they disagree on its impact on the molecular gas content. This is partially because, before cluster galaxies begin to lose molecular gas, they may gain it due to, for instance, momentary compression of galactic gas by the ram pressure (Fujita & Nagashima 1999; Kronberger et al. 2008; Kapferer et al. 2009; Bekki 2014; Henderson & Bekki 2016; Mok et al. 2016; Steinhauser, Schindler & Springel 2016; Lee et al. 2017; Ramos-Martínez, Gómez & Pérez-Villegas 2018; Vulcani et al. 2018; Safarzadeh & Loeb 2019; Troncoso-Iribarren et al. 2020; Roberts et al. 2021; Lee et al. 2022). On long timescales, however, the molecular gas content is expected to reduce. Hence, different results regarding the H₂ content of such galaxies can be obtained depending on the evolutionary stage that is observed. Moreover, since H₂ is predominantly concentrated in the inner regions of a galaxy (e.g. see Sun et al. 2020), where the gravitational potential is deepest, it is rare to find cases where there is ongoing stripping of H₂, and big samples are required in order to derive robust conclusions.

Furthermore, it is important to consider that galaxy populations in clusters differ between redshifts. Butcher & Oemler Jr (1978), and many others thereafter (Balogh et al. 1999; Poggianti et al. 2006; Saintonge, Tran & Holden 2008; Li, Yee & Ellingson 2009; Raichoor & Andreon 2012; Brodwin et al. 2013; Hennig et al. 2017; Wagner et al. 2017; Jian et al. 2018; Pintos-Castro et al. 2019), showed that clusters at high redshift have a higher fraction of blue SF galaxies compared to similar-mass clusters at z = 0. The morphology–density relation in clusters also gets shallower at high z (e.g. Postman et al. 2005), and galaxies in cluster cores sometimes show higher SFRs than the field (e.g. Tran et al. 2010). Taken together, these trends are attributed to the fact that clusters at high redshifts are newly formed and expected to go through an initial merger-driven phase, accompanied by the accretion of ample cool gas to support star formation. These factors, combined with poor statistics, complicate any investigation of the impact of environment on galaxies in high-z clusters.

Cosmological hydrodynamical simulations can help circumvent the issues present in observational studies by providing sufficient statistics, complete samples of galaxies, and a diverse range of environments. More importantly, such simulations enable a proper investigation of temporal variation in galaxy properties, and its relation to external processes. For example, Stevens et al. (2021) used the IllustrisTNG simulations (Marinacci et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018) to extensively explore the impact of environment on the H₂ content of z = 0 galaxies. They found that the H₂ content of satellites is, on average, lower than centrals by ≈0.6 dex. They also showed that the mean H₂ content of satellites in haloes with masses above |$10^{14}\, {\rm M_\odot }$| is lower by ≈1 dex relative to satellites in haloes with masses less than |$10^{12}\, {\rm M_\odot }$|⁠. This was attributed to mass loss post infall, which occurs for both |$\mathrm{H}\, \small {{\rm I}}$| and H₂, but at a faster rate for the former.

In this paper, we use the 100 Mpc ‘REFERENCE’ box of the eagle suite of hydrodynamical simulations (Crain et al. 2015; Schaye et al. 2015) to investigate whether SF galaxies in clusters differ in their H₂ content in comparison to other SF galaxies, and what is the underlying physics governing this difference, or lack thereof. This work complements Stevens et al. (2021) and extends it out to z = 1, roughly corresponding to half the age of the Universe. The paper is structured as follows. We describe the simulation and the methodology in Section 2, where Section 2.1 delineates details about the simulation data used in this work, and Section 2.2 explains the galaxy selection. In Section 3, we compare the H₂ content of SF galaxies in clusters with those outside clusters, comment on the factors that modulate the differences, and discuss the implications for galaxy quenching in clusters. We then extend the analysis on SF cluster galaxies in Section 4, where we demonstrate the influence from various modes of molecular gas transfer on the H₂ content with respect to their infall, and also assess the contribution from feedback. We summarise our main findings and some important caveats of the study in Section 5. Additionally, for the interested reader, we describe our approach for modelling the H₂ in Appendix  A, and the tests performed for identifying the valid prescriptions for eagle in Appendix  B.

2 SIMULATION AND METHODS

2.1 The eagle simulation

Evolution and Assembly of GaLaxies and their Environments, or eagle² (Crain et al. 2015; Schaye et al. 2015), is a suite of hydrodynamical simulations that were evolved with smoothed-particle hydrodynamics using a modified version of the gadget-3 code (a successor to the gadget-2 code described in Springel 2005). Here, we briefly describe the salient aspects of the simulation pertinent to the analysis in this paper. For a more generalised description of eagle, the interested readers should refer to Schaye et al. (2015), Crain et al. (2015), and the references therein.

In this work, we use the ‘REFERENCE’ eagle simulation corresponding to the box with a comoving side length |$L=100\, {\rm Mpc}$| and N = 2 × 1504³ particles (including both dark matter and baryons). The simulation is based on a flat Lambda cold dark matter (ΛCDM) cosmology with cosmological parameters from the Planck Collaboration XVI (2014) results. Each dark matter and (primordial)³ gas particle has mass |$m_\text{dm} = 9.70 \times 10^6\, {\rm M}_\odot$| and |$m_\text{g} = 1.81\times 10^6\, {\rm M}_\odot$|⁠, respectively. eagle lacks the resolution required for modelling fine-scale processes, which are rather implemented using the subgrid physics described in Crain et al. (2015). This, however, does not include the physics of the cold ISM. We compute the H₂ content for gas particles in post-processing, i.e. after the simulation was run. To this end, we first derive the neutral hydrogen content (Appendix A1) using the prescription of Rahmati et al. (2013), and then subsequently partition it into atomic and molecular components (Appendix A2). For the latter, we test five different prescriptions that have been widely used in the literature, namely Blitz & Rosolowsky (2006), Leroy et al. (2008), Krumholz (2013), Gnedin & Kravtsov (2011), and Gnedin & Draine (2014); hereafter BR06, L08, K13, GK11, and GD14, respectively. The description and implementation is provided in Appendix A2. We compare the outputs from these prescriptions against observational data and infer that only K13 and GD14 are sufficiently reliable for eagle (see Appendix  B). The results in this paper are based on GD14.⁴

Dark matter haloes in the simulation are first identified using a friends-of-friends (FoF) algorithm (Davis et al. 1985), and the star, BH and gas particles are considered to be part of the FoF halo containing the nearest dark matter particle, if the particle belongs to one. Each FoF halo along with its baryons is a ‘group,’ and the self-bound substructures (or subhaloes) within each group are identified using subfind (Springel et al. 2001; Dolag et al. 2009). We consider the virial mass of a group as M₂₀₀, which is the total mass enclosed within the radius where the mean density is 200 times the critical density of the universe. We define clusters as the groups with |$M_{200}\gt 10^{13.8}\, {\rm M}_\odot$| (see Section 3). The most massive subhalo in each group is referred to as the ‘central’ subhalo (which, if massive enough, may host the central galaxy), and the remaining subhaloes which orbit the central are ‘satellite’ subhaloes and may host the satellite galaxies. For each subhalo – and the corresponding galaxy – the centre is considered to be the position of the particle with the minimum gravitational potential according to subfind. The centre of the central subhalo/galaxy is considered to be the centre of the group. The merger tree for each halo constitutes the complete set of branches of its surviving subhaloes (Qu et al. 2017), and we use them to construct histories for several properties of dark matter haloes and the galaxies therein.

2.2 Selecting star-forming galaxies

We intend to adopt a selection strategy that results in populations that are (or close to) comparable to the heterogeneous observed samples across redshifts. As stated earlier, observed galaxies with molecular gas detections, especially at high redshifts (z > 0.5), tend to be SF or star bursting. This is unsurprising, since the systems that emit enough flux to allow for their detection by current telescopes are likely to be rich in molecular gas, and thus also have high SFRs. Considering this, we focus our analysis only on SF galaxies in eagle.

The star formation activity of a galaxy is usually characterised based on its position relative to the star formation main sequence, which is the average SFR–M_⋆ or sSFR–M_⋆ relation (sSFR is the specific star-formation rate or SFR/M_⋆). We use the sSFR–M_⋆ plane obtained for eagle to select SF galaxies at four different redshifts: z = 0, 0.3, 0.7 and 1. Each panel in Fig. 1 is the sSFR vs M_⋆ scatter plot for eagle and observed galaxies at the redshift ±0.1 of the redshift written in the bottom-left corner are shown. The latter only includes the galaxies with molecular gas measurements: Geach et al. (2011), Bauermeister et al. (2013), Tacconi et al. (2013), Scoville et al. (2014), Decarli et al. (2016), Grossi et al. (2016), Johnson et al. (2016), Saintonge et al. (2017), Wu et al. (2018), Aravena et al. (2019), Betti et al. (2019), Cairns et al. (2019), Castignani et al. (2019), Freundlich et al. (2019), Moretti et al. (2020), Castignani et al. (2020a, b, c), Belli et al. (2021), Dunne et al. (2021), Freundlich et al. (2021), Bezanson et al. (2022), Castignani et al. (2022), and Morokuma-Matsui et al. (2022). We omit eagle galaxies with |$M_\star \lt 10^9\, {\rm M}_\odot$|⁠, since the observed masses, particularly at high redshifts, tend to lie above this threshold. This criterion also happens to avoid galaxies that are not reasonably resolved in eagle (Schaye et al. 2015). The solid curve shows the main sequence (medians) at each redshift obtained for eagle galaxies, and the dashed curve shows the 16th percentiles of the sSFR for 12 bins of stellar mass between M_⋆/M_⊙ = [10⁹, 10^12.5] (the galaxies with |${\rm sSFR}=0\, {\rm yr}^{-1}$| are excluded while calculating the percentiles). The curves are only shown for the bins with >10 galaxies. The observed galaxies with molecular gas detections and upper limits are shown as filled and empty symbols, respectively.

Recall that each eagle galaxy is part of a gravitationally bound substructure identified using subfind, and this is in turn located within a group which may also contain other galaxies. For example, in a system with a central galaxy surrounded by multiple satellite galaxies, each galaxy is surrounded by its own DM subhalo, and the whole system is inside a group with virial mass of M₂₀₀. We define ‘clusters’ as those eagle groups that have mass |$M_{200}\gt 10^{13.8}\, {\rm M}_\odot$|⁠, and the rest as ‘non-clusters’. This threshold is somewhat arbitrary, as the transition mass is not a well-defined quantity. We adopt this value simply because it (approximately) delineates the lower bound of (proto-)cluster masses in molecular gas observations (e.g. Wagg et al. 2012; Tadaki et al. 2014; Casey 2016; Noble et al. 2017; Castignani et al. 2018; Coogan et al. 2018; Strazzullo et al. 2018; Wang et al. 2018; Zabel et al. 2019; Alberts et al. 2022; Castignani et al. 2022). The clusters thus selected have masses M₂₀₀/M_⊙ ≲ 10^14.6, 10^14.4, 10^14.3, 10^14.1 at z = 0, 0.3, 0.7, 1, respectively. Fig. 1 shows the eagle galaxies residing in clusters as black points, and the others as grey points.

The results show that, at each redshift, detections in the observational data sets predominantly lie on or above the dashed curve. Therefore, for each of these four redshifts, we select galaxies from eagle that have sSFR above the 16th percentile for their stellar mass.⁵ The breakdown of galaxies between SF and non-SF (including |${\rm sSFR}=0\, {\rm yr}^{-1}$|⁠) for |$M_\star \gt 10^{9}\, {\rm M}_\odot$| is shown in the first column of Table 1. Since |${\rm sSFR}=0\, {\rm yr}^{-1}$| galaxies do not constitute a substantial fraction, SF galaxies account for ≳65 per cent of the sample at each of the four redshifts. The second column of Table 1 shows the number of SF galaxies that reside in clusters, and those that do not. Only ≲3 per cent of all the SF galaxies in our sample are located in clusters. The last row shows that the number of SF cluster galaxies is <100 for z = 1, and we find it to reduce further for higher redshifts. This is the reason for limiting our analysis to z ≤ 1.

We realise that our adopted criteria are far from perfect if we are to match the observations exactly. Even if we were to attempt that, it would not be accurate because the observed samples are heterogeneous and a single selection function cannot mimic the full data set. In addition, the maximum sSFR of eagle galaxies is almost an order of magnitude below some observations (Fig. 1), which is – in part – due to the limited volume of eagle, which does not sample enough rare, star-bursting galaxies. However, we would like to emphasise that an exact match with observations is not our intention here. Instead, the aim is to exclude quiescent galaxies while also maintaining a sample that is fairly representative.⁶

Also, while there is no explicitly enforced limit on the gas content, the selection of SF galaxies inherently means that they must possess ample cold gas to sustain their star formation. Note that imposing a strict limit would not be appropriate if we want our selection to be guided by observations, as detections and upper limits overlap in the |$M_\mathrm{H_2}$|–M_⋆ space; this can be seen, for example, in Fig. A1 for xCOLD GASS. Moreover, the strategies for computing upper limits differ between surveys.

We note that subfind occasionally identifies dense star-clusters within galaxies as subhaloes (McAlpine et al. 2016). These structures have little stellar mass (⁠|$M_\star \lt 10^8\, {\rm M}_\odot$|⁠), but anomalously high metallicity or BH mass. Our sample is devoid of such spurious subhaloes due to the adopted stellar mass cut.

3 MOLECULAR GAS CONTENT OF SF GALAXIES IN CLUSTERS AND NON-CLUSTERS

In this section, we address the main aim of this study: the impact of cluster environment on the molecular gas content of SF galaxies, and its dependence on redshift. For this, at each redshift, we compare the SF eagle galaxies that are located within clusters against those that are not. Each of these two samples contains both central and satellite galaxies; the exact numbers are given in Table 1. The number of clusters containing these galaxies is equal to the number of centrals, except for z = 0 where the central in one of the clusters is not SF – that is, there are actually 14 clusters with SF galaxies at this redshift. Satellites are generally considered to be more susceptible to environment than centrals, but we do not show results for the two classes separately here for the sake of consistency with observational studies, which cannot reliably distinguish between them in most cases. Notwithstanding, we do mention the satellite fractions if required for proper interpretation of the results. Overall, our selection causes the non-cluster sample to be predominantly composed of centrals (≈70 per cent), and the cluster sample to be dominated by satellites (≈95 per cent), for all redshifts. Note that each cluster galaxy below |$M_\star \approx 10^{11.5}\, {\rm M}_\odot$| is a satellite, which is expected for haloes at this mass scale.

In Fig. 2, we show the |$M_{\mathrm{H_2}}$|–M_⋆ relations for the two samples at each of the four redshifts, along with observations. The median relations for eagle (based on GD14) are shown as the outlined solid orange curve for SF galaxies in clusters (⁠|$M_{200}\gt 10^{13.8}\, {\rm M}_\odot$|⁠) and dashed blue curve for non-clusters (⁠|$M_{200}\lt 10^{13.8}\, {\rm M}_\odot$|⁠), using 0.3 dex wide bins that lie within M_⋆/M_⊙ = [10⁹, 10^12.3]. The solid grey curve shows the medians computed if we consider all the cluster galaxies. The shaded regions with the same colour coding show 1σ bootstrapping errors on the medians, i.e. they span the 16th and 84th percentiles of the medians of the bootstrap samples. The dotted magenta curve shows the median based on satellites in the non-cluster sample.⁷ The filled and open squares show the detections and upper limits, respectively, for galaxies in (proto-)clusters taken from Geach et al. (2011), Corbelli et al. (2012), Cybulski et al. (2016), Grossi et al. (2016), Johnson et al. (2016), Saintonge et al. (2017), Wu et al. (2018), Betti et al. (2019), Cairns et al. (2019), Castignani et al. (2019), Zabel et al. (2019), Moretti et al. (2020), Castignani et al. (2020a, b, c), Dunne et al. (2021), Castignani et al. (2022), Morokuma-Matsui et al. (2022), and Zabel et al. (2022). For systems that are redundant between the data sets, the most recent values are shown. The scaling relation for the observed star formation main sequence by Tacconi et al. (2018) is displayed as the solid green line, and is based on a representative sample of SF galaxies between z = 0 and  4.4. This relation is applicable to molecular gas masses based on both CO and dust emission.

The z = 0 results for eagle indicate that, except for the lowest and highest stellar masses, SF galaxies in clusters typically have a deficit of H₂ compared to their non-cluster counterparts. The difference – at ≈0.5 dex – is most pronounced for |$M_\star \approx 10^{10.5}\, {\rm M}_\odot$|⁠. The trend also exists at z > 0, albeit with offsets that get smaller with increasing redshift at a given stellar mass. We find that, at each redshift, there are certain stellar mass bins where differences between the medians are significant even after accounting for their bootstrapping errors. Hence, these results establish that (simulated) SF galaxies in clusters have lower H₂ content than those not in clusters, provided they lie in a certain mass range.

The situation, though, is not so clear for observed galaxies, and the data shown in Fig. 2 help demonstrate why that is the case. In observations, the impact of the cluster environment is usually inferred by comparing the molecular gas content of cluster galaxies to that of galaxies on the main sequence at the cluster’s redshift. If we assume that the main-sequence relation from Tacconi et al. (2018) represents non-cluster galaxies (which is a reasonable assumption because they dominate the statistics), we can infer an enhancement or deficit of molecular gas in clusters depending on the data set to which we choose to compare. For example, at z = 0, all of the Moretti et al. (2020) detections lie above the observed main sequence, whereas most of the data from Cairns et al. (2019) lie below it. To add to this, an accurate determination of the difference in the H₂ content using observational data is limited by the presence of non-detections, which do not contaminate the simulated data.

Fig. 2 also shows that SF cluster galaxies are not representative of the whole cluster population. For any stellar mass bin with reliable statistics (>10 galaxies), an average cluster galaxy contains less H₂ than a typical SF one, and this difference increases at lower masses, indicating that SF systems sample a smaller proportion of the population. This is somewhat expected, given that susceptibility to environment increases as a galaxy’s mass decreases, which then reduces the likelihood of preserving the cold gas reservoir (e.g. Marasco et al. 2016). We predict that, for the cluster masses explored here, observations carried out using a radio/IR telescope with sufficient sensitivity to detect gas-poor systems would produce a |$M_\mathrm{H_2}$|–M_⋆ relation that is (at least qualitatively) similar to the solid grey curve.

Note that the difference in H₂ content of the cluster and non-cluster SF samples could be a representation of the difference in their satellite fraction, which is evidently greater for the former (Table 1). We investigate this by plotting the median H₂ content of satellites in the non-cluster sample in Fig. 2 (dotted curve). The offsets still persist relative to the full cluster sample, and at similar magnitudes as for the full non-cluster sample. Therefore, the offsets between the cluster and non-cluster sample are also indicative of differences between the satellites in these two kinds of groups, and not simply driven by the difference in satellite fraction.

The main-sequence relation itself is broadly consistent with eagle at z = 0 except at |$M_\star \gt 10^{11}\, {\rm M}_\odot$| (dashed blue and violet curves), but its normalisation gets progressively higher than that of eagle as redshift increases. Although the exact reason is uncertain, there are a number of factors that can potentially contribute to this. One possibility is that the detections are biased towards galaxies that are more gas-rich as redshift increases (Malmquist bias). Another is that the uncertainties regarding the CO-to-H₂ conversion factor can, in principle, lead to a systematic effect on the inferred masses. It may also reflect the limitations of the post-processing schemes adopted to compute H₂ masses for eagle galaxies at higher redshifts. Nevertheless, our focus is on understanding the offset in molecular gas content in different environments, which is more likely to be a function of the physics in the simulation rather than any post-processing scheme. Thus, the results in Fig. 2 are useful.

We note that eagle clusters contain more hot gas compared to observed clusters, due to AGN feedback not being efficient enough to eject the gas out of their host haloes (Barnes et al. 2017). Although this does lead to more massive centrals (see Bahé et al. 2017), our cluster sample is dominated by satellites instead. Barnes et al. (2017) showed that the feedback in centrals is only strong enough to transport the gas out to 0.5 ≲ r/R₅₀₀ ≲ 1, which means elevated density of intracluster gas, and therefore, stronger ram pressure acting on satellites at such radii. However, since the stellar mass function of cluster satellites in eagle is consistent with observations (Bahé et al. 2017), this does not seem to have a significant impact on the global properties of such systems. The exact physics involved remains elusive, and a proper investigation would require a detailed analysis, which is beyond the scope of this paper.

3.1 Potential factors modulating the offset in |$M_\mathrm{H_2}$|

Though the offset in |$M_\mathrm{H_2}$| is visible at all redshifts, its magnitude at a given stellar mass increases towards lower redshifts. The trend is particularly striking for the sample with all the cluster galaxies (solid grey curve in Fig. 2). One potential reason for this could be that galaxies residing in clusters at lower redshifts have been bound to a cluster, and thus influenced by the cluster environment, for longer durations. We now investigate the validity of this scenario.

Each galaxy in eagle has formed through mergers of smaller subhaloes, referred to as the galaxy’s ‘progenitors’. The progenitors at a given snapshot along the most massive branch of its merger tree is called the ‘main progenitor’ (see Qu et al. 2017). We study the acquisition of our galaxies into their groups and its impact on the H₂ content by tracking their main progenitors, noting their masses and that of their host group. For each cluster galaxy, we interpolate the group mass history to determine the look-back time (t_c) when |$M_{200}=10^{13.8}\, {\rm M}_\odot$|⁠, i.e. the first instant it was bound to a cluster.

Fig. 3 shows the |$M_\mathrm{H_2}$| histories for the cluster and non-cluster samples, where each panel corresponds to the redshift written in the bottom right corner. t is the look-back time corresponding to the progenitor’s redshift, and t₀ is the look-back time at the sample redshift; cosmic time moves forwards from left to right on the x-axis. The solid and dashed curves in each panel show the median |$M_\mathrm{H_2}$| history for the cluster and non-cluster samples, respectively. The corresponding shaded regions enclose the 20–80th percentiles. The dotted curve is the median for satellite galaxies in the non-cluster sample. The downward-pointing arrow is the median t₀ − t_c, corresponding to first infall into a cluster, and the horizontal error bar shows its 20–80th percentile range. Given that the |$M_\mathrm{H_2}$| offset varies with M_⋆, it is important to control for stellar mass for proper correspondence with the results in Fig. 2. Here, we show the results for the galaxies with 10^10.5 < M_⋆/M_⊙ < 10^10.8 just as an example, but have verified that similar results apply to other stellar mass bins that show noticeable H₂ deficits.

These results show that, at some point in the past, both the samples possessed similar quantities of H₂, which provides an important insight: the galaxies outside clusters did not start off with more H₂ content to begin with. Later on, the cluster galaxies incurred a steep decline in their |$M_\mathrm{H_2}$| relative to the non-cluster sample, and this evolution eventually led to the offset seen at t = t₀. A similar history is seen for the full cluster sample, but note the steeper decline. The non-cluster sample with only satellites also shows an average decline in H₂ mass relative to the full non-cluster sample, but – as expected – not as steep as the cluster samples.

Interestingly, the instant when the SF cluster sample deviated seems to have roughly coincided with the time when a cluster first accreted these galaxies. However, the wide error-bars on the arrows mean that the curves are averaging over a wide range of infall times, and are therefore not ideal for inferring the typical behaviour of a galaxy with respect to its infall into a cluster. Hence, we estimate the global significance of this trend by examining the quantitative relationship between the offset in |$M_\mathrm{H_2}$| and t_c − t₀, or the time elapsed since infall. The scatter plot is shown in Fig. 4, where each data point shows the offset between the median logarithmic H₂ masses (underpinning the dashed blue and solid orange curves in Fig. 2), and the median t_c − t₀, for a particular stellar mass bin and redshift, depicted using a specific colour and symbol, respectively. The labels on the colour bar are the bin edges used to compute the medians, and are the same as those used for Fig. 2. We only show the results for the bins that contain >10 galaxies, which effectively excludes all the bins beyond |$M_\star =10^{11.1}\, {\rm M}_\odot$|⁠. The error-bars are the bootstrapping uncertainties.⁸

The figure shows that the H₂ offset is indeed positively correlated with the time elapsed since infall for cluster galaxies; the Spearman rank coefficient and p-value for the whole sample is ≈0.57 and ≈0.002, respectively. The scatter at a fixed t_c − t₀ is likely due to some additional factors. For example, a higher group mass and proximity to group’s centre generally leads to a lower gas accretion rate (e.g. van de Voort et al. 2017). For a fixed mass, group haloes at different redshifts have varied crossing times, which means that galaxies accreted at those redshifts have traversed different number of orbits for the same time elapsed since infall. These haloes also have different intrahalo gas densities at their virial radius, which translates to different strengths of ram pressure for the same galaxy velocity. We assess the importance of potential modulating factors through a principal component analysis (PCA), which is a statistical technique that derives orthogonal vectors (called principal components) along which the variance is maximised in a multidimensional parameter space. This is utilised in galaxy evolution studies to find useful trends that may otherwise remain undiscovered due to noise in the data (e.g. Costagliola et al. 2011; Bothwell et al. 2016; Lagos et al. 2016; Cochrane & Best 2018). For our purposes, we conduct the analysis⁹ on the data set containing the current group mass, gas density at R₂₀₀ at infall (taken as that of dark matter; ρ_{200, inf}), crossing time-scale of group at infall (t_{cross, inf}), ratio of the galaxy’s mass to the group’s mass at infall (M_gal/M₂₀₀)_inf, t_c − t₀, and the H₂ offset. Fig. 4 can be envisaged as a projection of this hyperspace on the plane defined by the H₂ offset and t_c − t₀.

All the parameters are standardised before performing the PCA to ensure that the significance of a quantity is agnostic to its dynamic range. The results from this analysis are presented in Table 2, where the first column states the principal component’s rank, the second column shows the percentage of the total variance of the data explained by the component, and the rest of the columns show the coefficients for the various parameters in the data set. The coefficient’s absolute value indicates the level of importance of the corresponding parameter for that principal component. Most of the data can be explained by the first two components, accounting for ≈88 per cent of the variance. We focus on these two for the following.

The first component represents ≈61 per cent of the variance, and the coefficients show that most of it is contributed by M₂₀₀, ρ_{200, inf}, t_{cross, inf} and t_c − t₀. The similar coefficient magnitudes suggest significance of parameters other than t_c − t₀, and the dominance of four parameters implies that the component is useful for understanding trends mainly amongst them. The t_c − t₀ is positively correlated to M₂₀₀ and t_{cross, inf} but negatively correlated to ρ_{200, inf}. The latter hints towards the importance of ram pressure, as a galaxy loses cold gas quicker and sustains its star formation for shorter durations under stronger pressures. ρ_{200, inf} and t_{cross, inf} have the same strengths, partially because both quantities are just functions of the critical density of the universe.

The second component explains ≈27 per cent of the variance and is primarily composed of (M_gal/M₂₀₀)_inf and the H₂ offset, with coefficients almost twice those for other parameters. The coefficients suggest that these quantities are positively correlated for some systems, which is contrary to the expectation if gas stripping dominates the offset, as its efficacy is anti-correlated to this ratio (e.g. Tormen, Moscardini & Yoshida 2004; Bahé & McCarthy 2015; Yun et al. 2019). This can nevertheless be caused by feedback events through the proverbial ‘mass quenching,’ which is more effective at higher masses due to hotter circumgalactic gas (e.g. Gabor & Davé 2015). We will investigate this further in later sections. Overall, these results showcase the parameters that are contributing to the scatter in Fig. 4.

The correlation between H₂ offset and the time since infall is also apparent in Fig. 4 for points corresponding to a given redshift (same symbol) or stellar mass (same colour). For the latter, note that t_c − t₀ always decreases at higher z regardless of how the H₂ offset behaves for that stellar mass, because of less cosmic time elapsed since the cluster’s formation due to hierarchical assembly of haloes,¹⁰ and shorter quenching time-scales for cluster galaxies at higher redshifts (e.g. Baxter et al. 2022). In addition, we find that these galaxies were always centrals before infall (t = t_c), meaning they were not pre-processed in groups. Overall, these results support the notion that the trends in Fig. 2 can be attributed to the impact of the cluster environment on galaxies.

It is worth noting that the H₂ offset is non-existent at |$M_\star \lesssim 10^{9.3}\, {\rm M}_\odot$|⁠, which may seem peculiar at first glance, as one might expect such low-mass galaxies to be most susceptible to environmental processes due to their shallow gravitational potential wells. Even the H₂ histories of these galaxies are very similar to their counterparts outside clusters, as if their environment had no influence on them. We know that such low-mass satellites in z = 0 eagle clusters generally get quenched by ≲1 Gyr after infall (Wright et al. 2019). In fact, low-mass SF galaxies are rare in clusters, constituting ≈10 per cent of all the low-mass cluster satellites at z = 0 (in agreement with Shao et al. 2018), and ≈35 per cent at z = 1. Moreover, the cluster galaxies at this M_⋆ tend to be devoid of H₂, regardless of the redshift (see Fig. 2). This implies that such galaxies can be SF while also being in a cluster only if they were accreted recently, as they would not be able to preserve enough cold gas otherwise. As it turns out, all the SF cluster galaxies with |$M_\star \lt 10^{9.3}\, {\rm M}_\odot$| are satellites, and most of these have only been bound to a cluster for ≲300 Myr.

Similarly, Fig. 2 shows that the H₂ content of SF galaxies with |$M_\star \gtrsim 10^{11.2}\, {\rm M}_\odot$| is oblivious to the cluster environment. However, this cannot be simply attributed to a short time span post infall in all cases. For example, the z = 0 cluster galaxies at |$M_\star \approx 10^{11.2}\, {\rm M}_\odot$| are satellites with a median t_c − t₀ of ≈2 Gyr, which is similar to the quenching time-scale for such systems (Wright et al. 2019). Though, considering that the efficacy of gas stripping decreases with the satellite-to-central mass ratio (e.g. Tormen et al. 2004; Bahé & McCarthy 2015; Yun et al. 2019), it is plausible for a few massive galaxies in clusters to protect their gas reservoirs against such mechanisms for durations longer than usual. At |$M_\star \gtrsim 10^{11.4}\, {\rm M}_\odot$|⁠, the SF cluster galaxies are mostly centrals at all redshifts, and are therefore relatively immune to environmental stripping. We, however, would also like to qualify that the data suffers from poor statistics in this range.

3.2 Implications for galaxy quenching in clusters

While it is well accepted that infall into a cluster exacerbates quenching, there has been an ongoing debate on the details of how SFR declines post infall. Deriving an accurate model for this is important for a comprehensive understanding of galaxy evolution, because the rate of decline of a galaxy’s SFR provides crucial insights for the dominant quenching mechanism in clusters; for instance, ram-pressure stripping causes a more abrupt decline than reducing gas accretion.

In a seminal study, Wetzel et al. (2013) showed that observed properties of cluster galaxies (like colours and SFRs) cannot be explained by a model that only assumes rapid or slow quenching. In an attempt to resolve this, they proposed the ‘delayed-then-rapid’ quenching scenario, wherein a satellite’s SFR remains unaffected by the environment for a few Gyr after infall, and then incurs a rapid drop. This idea has amassed significant support over the last decade (e.g. McGee, Bower & Balogh 2014; Mok et al. 2014; Muzzin et al. 2014; Fossati et al. 2017; Foltz et al. 2018; Moutard et al. 2018; Rhee et al. 2020; Oman et al. 2021; Sampaio et al. 2022). Haines et al. (2015) purported a similar ‘slow-then-rapid’ model – involving an initial phase of gradually reducing SFR – as a more accurate model, because of its utility in explaining the observations of galaxies transitioning between SF and passive. This model has recently been gaining advocacy (Maier et al. 2019; Roberts et al. 2019; Kipper et al. 2021), and is generally understood physically as the ‘slow’ phase caused by strangulation, and the ‘rapid’ phase corresponding to the galaxy reaching inner regions of its group halo (within 0.5R₂₀₀) and experiencing strong cold gas stripping due to ram pressure and tidal forces. Futhermore, it has been suggested that the time spent in the ‘rapid’ phase reduces with stellar mass (e.g. Rhee et al. 2020), which is usually interpreted as feedback events being more important for massive galaxies. Therefore, it is warranted to investigate the evolution of our cluster samples in the broader context of galaxy quenching.

In Fig. 5, we show the SFR and |$M_\mathrm{H_2}$| normalized by their values at infall for all the cluster galaxies (grey) and the SF subsample (orange) with 10^10.5 < M_⋆/M_⊙ < 10^10.8. The solid grey curves show that the SFR for a typical cluster galaxy started declining ≳2 Gyr prior to infall, and the decline got stronger from infall onwards. This indicates that the cluster environment expedites quenching of galaxies at all redshifts, confirming the findings from previous studies based on both observations (e.g. Chung et al. 2011; Haines et al. 2015; Pintos-Castro et al. 2019; Stevens et al. 2021), and simulations (e.g. Bahé & McCarthy 2015; Pallero et al. 2019; Donnari et al. 2021; Stevens et al. 2021). The top panels show that the slow decline persisted for about few Gyr post infall, followed by a rapid drop in SFR, demonstrating that these galaxies underwent ‘slow-then-rapid’ quenching.

Based on similar curves for different mass bins, we find that the average active time, t_act¹¹ – i.e. the duration between t = t_c and the moment of sharpest drop in the SFR (apparent in Fig. 5) – is not constant for all stellar masses and redshifts. At a fixed redshift, the median t_act increases monotonically with mass until |$M_\star \approx 10^{9.9}\, {\rm M}_\odot$| and remains nearly constant for higher masses. The constant is lower for higher redshifts, ranging from the median t_act ≈ 3 Gyr for z = 0 to t_act ≈ 1.35 Gyr for z = 0.7, deviating by ≈0.25 Gyr between masses. This is, in fact, consistent with the evolution of t_cross that reduces from ≈2.9 Gyr at z = 0 to ≈1.6 Gyr at z = 0.7, in agreement with observations (Mok et al. 2014; Fossati et al. 2017; Lemaux et al. 2019; Rhee et al. 2020). We would like to point out, however, that ‘slow-then-rapid’ quenching is not an accurate description of the evolutionary histories of galaxies with |$M_\star \lesssim 10^{9.5}\, {\rm M}_\odot$|⁠, where the ‘slow’ phase is practically absent with t_act ≲ 1 Gyr.

These results favour a quenching mechanism that is largely independent of redshift and mass. Since it takes some time after infall for a galaxy to reach the dense regions of a cluster, it is expected to be influenced more by starvation/strangulation until the first pericentric passage, and then undergo ram-pressure and tidal stripping. Overall, this is expected to quench the galaxy over a duration similar to the crossing timescale, which is evident in the variation of t_act with redshift for |$M_\star \gtrsim 10^{9.9}\, {\rm M}_\odot$|⁠. Apart from indicating mass-invariance, the constant t_act in this regime for a fixed redshift rules out AGN feedback as the primary quenching process for cluster satellites, especially at |$M_\star \gtrsim 10^{10.5}\, {\rm M}_\odot$|⁠, where one might expect it to be relevant (Mitchell, Schaye & Bower 2020). At low masses (like |$M_\star \lesssim 10^{9.5}\, {\rm M}_\odot$|⁠), a mass-dependence emerges likely because gravitational potentials become considerably shallower, which is conducive for efficient ram-pressure stripping soon after infall, and expected to quench within ≈1 Gyr (Roediger & Hensler 2005; Steinhauser et al. 2016).

Note that the SF subsample has not evolved like a typical cluster galaxy. As shown by the solid orange curves in Fig. 5, their SFRs did not decline until infall, and the rate of decline was always shallower than the average cluster galaxy. Furthermore, we find some SF cluster galaxies that have not quenched even after being bound to clusters for durations longer than the typical active time. By selection, SF galaxies have been more robust against external quenching mechanisms than any typical galaxy that is exposed to cluster environment. We examine the group membership of our cluster galaxies across time, and find that most of them were centrals ≳ 2 Gyr before infall, whether we consider the full cluster sample or the SF one. However, while the SF systems remained centrals until infall into a cluster, an average cluster galaxy was a satellite in a lower mass group for about 2 Gyr before being accreted by a cluster. This means that, unlike the SF subsample, most cluster galaxies have been pre-processed in such groups, which explains the earlier commencement of their molecular gas loss and SFR’s decline.

In addition, the dashed curves in Fig. 5 show that the |$M_\mathrm{H_2}$|’s evolution always traces that of SFR, regardless of stellar mass or redshift, which means that the impact of molecular gas loss on the SFR is nearly instantaneous. Another way to interpret this is that the cluster environment has negligible impact on the star formation efficiency of galaxies, at least for the temporal resolution of the simulation. We find this to hold true across stellar masses.

4 UNCOVERING THE PHYSICAL ORIGIN OF THE DECLINE IN THE H₂ CONTENT OF SF CLUSTER GALAXIES

In this section, we investigate how various processes influenced the H₂ content of SF cluster galaxies before and after infall. Through this analysis, we aim to highlight the mechanisms that have a stronger bearing on the evolution of H₂ in such galaxies.

4.1 Contributions from phenomenological modes

In eagle, when a gas particle gets converted into a star particle, all of its mass gets transferred to the star particle, not just the mass that is contained in H₂ (which is post-processed). This means that, even if star formation is the only process affecting the ISM, a reduction in |$M_\mathrm{H_2}$| is not equal to the increase in M_⋆. Fortunately, since every particle has a unique ID in eagle, it is possible to identify the gas particles that were converted into stellar particles by comparing their IDs between successive snapshots. Likewise, we can also identify the gas particles that were removed from the galaxy, those that were accreted by the galaxy, and the ones that remained in the galaxy but incurred a change in their H₂ mass due to, for example, variation in the local ionisation field, mass-loss from stars, or radiative cooling/heating. Thus, particle-tracking is a viable technique for understanding how distinct modes govern gas flows to-, from- and within galaxies in eagle.

Note that the mass added via gas accretion cannot be estimated by simply identifying the gas particles that did not exist in the subhalo in the previous snapshot, because this does not account for changes incurred due to transport of gas from the circumgalactic medium (CGM) to the ISM. Thus, we need to distinguish between the two for the analysis that follows. Since there is no physical boundary that separates them, this can be carried out in a variety of ways. The interstellar gas is generally denser, cooler, star-forming, and located closer to the halo centre. Consequently, methodologies adopted in the literature usually involve enforcing a cut on star formation rate (e.g. Wijers, Schaye & Oppenheimer 2020) and/or radial distance (e.g. van de Voort et al. 2017; DeFelippis et al. 2020; Fielding et al. 2020) coupled with that on density (e.g. Anglés-Alcázar et al. 2017; Correa et al. 2018; Hafen et al. 2019; Mitchell et al. 2020), temperature (e.g. Correa et al. 2018; Mitchell et al. 2020; Wright et al. 2022) or rotational support (e.g. Mitchell et al. 2018; Garratt-Smithson et al. 2021). We follow the approach delineated in Wright et al. (2022), where we first determine the radial extent of the galaxy, R_bary: the smallest radius¹² at which the cumulative baryonic-mass profile based on the stars and all the gas that is cool (⁠|$T\lt 10^{4.3}\, {\rm K}$|⁠) or star-forming (SFR > 0) in the subhalo resembles an isothermal gas halo (i.e. the density follows ρ_bary∝r⁻²; see Stevens et al. 2014). The ISM is then defined as the cool or star-forming gas that lies within R_bary, and the rest of the subhalo gas is the CGM. Our visual inspection of the resultant gas distributions suggests that this method is indeed effective.

Once we have categorised the particles into the two classes, we track the stellar, gas and BH particles of the main progenitor and, for each pair of successive snapshots, determine the amount of H₂ that is

• lost due to star formation,

• lost when particles move from the ISM to the CGM,

• lost due to particles escaping the system (or getting unbound),¹³

• lost due to a decrease in the H₂ content of particles in the ISM,

• added due to introduction of fresh gas in the system,

• added when particles in the CGM get accreted by the galaxy, and

• added due to increase in the H₂ content of particles in the ISM.

For (i), we identify the progenitor’s gas particles that exist as star particles in the descendant, and determine their cumulative H₂ mass. For (ii), we select the gas particles in the progenitor’s ISM that are in the descendant’s CGM, identify the ones whose H₂ mass (⁠|$m_{\rm H_2}$|⁠) had reduced, and add up their deficits. For (iii), we identify the progenitor’s gas particles that are not present in the descendant as either gas, stellar or BH particles, and compute their cumulative H₂ mass. For (iv), we select the gas particles that are common in the ISMs of the progenitor and the descendant, identify the ones that lost their H₂, and calculate the cumulative deficit. For (v), we identify the descendant’s gas particles that are not present in the progenitor, and add up their H₂ mass. For (vi), we identify the CGM particles in the progenitor that are present in the descendant’s ISM, select the ones whose |$m_{\rm H_2}$| had increased, and compute their cumulative mass-credit. Lastly, for (vii), we select the ISM particles that are common between the progenitor and the descendant, identify those that had an increase in |$m_{\rm H_2}$|⁠, and add up their credits.

One may notice that we have not alluded to the role of BHs above. Given that BHs accrete material that lies in their proximity, some H₂ can be lost to this mechanism. Likewise, some mass may be lost due to conversion of gas particles to BHs. It is not possible to identify individual gas particles that were accreted by a BH in eagle, but the data do provide the cumulative accreted mass of each BH particle – meaning one can calculate the total gas mass that was lost via BH accretion. We find this to be |$\sim 10^6\, {\rm M}_\odot$| at max, implying that the H₂ component of the lost gas is negligible. Comparisons of particle IDs reveal that the H₂ mass lost through conversion of gas particles to BHs is also negligible. Additionally, it is possible that some new gas particles were captured by the subhalo and then converted to stars between successive snapshots. Although these would not be identified in our analysis, they do not cause any net change in the H₂ content and can therefore be ignored.

For each mode from (i) to (vii), Fig. 6 shows the median H₂ gain/loss rate divided by the H₂ mass (or the specific gain/loss rate) against the time since infall for the SF cluster galaxies with 10^10.5 < M_⋆/M_⊙ < 10^10.8, where each panel corresponds to a specific sample-redshift (as in previous figures). We have plotted the specific rates in particular because these are especially suited for assessing the efficacy of a process in governing the gas content (e.g. see Wright et al. 2022). The shaded regions encompass 20–80th percentiles. The vertical dashed line corresponds to the moment when the look-back time t = t_c.

The results indicate that, before a typical galaxy entered a cluster (i.e. leftwards of the vertical line), H₂ was mainly replenished via increase in |$m_\mathrm{H_2}$| of the ISM particles (solid purple), and star formation was the main mode of loss (solid yellow) for all redshifts. The two modes of gas accretion (solid cream and dashed-dot grey) caused replenishment rates that were lower than the contribution from ISM particles gaining H₂, by ≳ 0.75 dex, and therefore had a relatively weaker impact on the H₂ content. Similarly, the loss rate corresponding to the movement of gas from ISM to CGM (dashed-double-dot magenta) was about three times lower than that due to gas escape (dashed red).

As the galaxies approached infall, the replenishment rate due to particles gaining H₂ had decreased, whereas the star-formation-driven depletion rate remained constant and approached the former. The latter is particularly interesting, because it suggests that the reduction in SFR was caused by lowering of the H₂ content rather than more efficient star formation, echoing the result in Fig. 5. Meanwhile, the escape-driven rate increased by ≈0.25 dex relative to its value a Gyr prior to t_c, and thereby, became similar to the star-formation-driven rate.

Post infall (rightwards of the vertical line), the losses caused by star formation and escape dominated over the replenishment modes, and there was a sharp decline in the replenishment of H₂ via gas accretion, as is expected for galaxies in cluster environments (van de Voort et al. 2017; Wright et al. 2020, 2022). The curve corresponding to the introduction of fresh gas (dashed-dot grey) truncates ≈1 Gyr after infall, because gas supply is cut off for most galaxies. Overall, this shows that the steeper decline of |$M_\mathrm{H_2}$| of cluster galaxies post infall is a result of the dominance of star formation and gas escape over all the modes that replenish H₂.

There are some additional insights that can be derived from this analysis. For instance, it is worth noting that not all the H₂ gained as a result of increase in |$m_\mathrm{H_2}$| of ISM particles can be owed to cooling. In fact, as shown using the dotted blue curve, cooling accounts ≲50 per cent at any given moment. Heating (dashed orange) has similar contribution to the decrease in |$m_\mathrm{H_2}$| of ISM particles. On the other hand, we find that all the transport from CGM to ISM is accompanied by cooling, the movement of gas from ISM to CGM always involves heating, and almost all the escaped H₂ comes from ISM. Similar results are obtained for other mass bins.

4.2 Significance of stellar/AGN feedback (or a lack thereof)

While it is clear that gas escape was one of the major modes of H₂ loss post infall, whether this was primarily caused by stellar feedback, AGN feedback, or direct stripping of H₂, remains abstruse. In this regard, we assess the impact of the feedback mechanisms by comparing the injected energy against the binding energy of gas particles.

In order to calculate the cumulative energy released due to stellar feedback between successive snapshots, we first identify all the gas particles that were converted to stellar particles, and note the mass of each such stellar particle at the time of its formation (m_⋆). Then, the injected energy (in ergs) is determined as (Schaye et al. 2015),

where f_th is taken from the eagle data base and N is the number of particles. Likewise, the energy injected by the black hole is computed as

where c is the speed of light in vacuum and m_acc is the total mass accreted by the black hole between the successive snapshots.

We quantify the total binding energy of gas particles prior to their escape. Since the aim is to examine the impact of feedback on the H₂ loss caused by gas escape, we focus on the particles that were lost through this mode. Amongst these, we exclude the particles that did escape but were not representative of H₂. To this end, we rank order the escaped particles in decreasing order of their H₂ mass in the first snapshot, and determine the cumulative mass profile by summing up from lowest to highest rank. Then, we select all the particles with ranks below the rank corresponding to 90 per cent of their total H₂ mass, and add up their binding energies to obtain E_bind.

On examining the evolution of E_⋆/E_bind (not shown) for the galaxies pertaining to Fig. 6, we find that the E_⋆/E_bind remained fairly constant at ≈10⁻¹³ throughout the history of our galaxies. The E_agn/E_bind, on the other hand, had an erratic history where it was either 0 or ≲10⁻³. Therefore, neither of the feedback mechanisms were efficient enough to eject the molecular gas by themselves, before, or after infall into clusters. This also holds true for other stellar masses, including M_⋆ ≲ 10^9.5 where quenching is generally fast. Thus, it appears that the energy required for escape of H₂ from SF cluster galaxies was essentially provided by stripping mechanisms that they were subjected to after infall into clusters.

5 SUMMARY

Since the detection of H₂ requires observed galaxies to possess sufficient quantities of cold gas, most of the detected systems tend to be SF galaxies. The observational data for such galaxies in clusters, particularly at high redshifts, remains sparse, precluding any overall consensus regarding the relationship between the H₂ content of SF galaxies, the cluster environment, and redshift. In this paper, we investigated this relationship by leveraging the 100 Mpc ‘REFERENCE’ run from the eagle project (Crain et al. 2015; Schaye et al. 2015). We selected SF galaxies at four different redshifts ranging from z = 0 to z = 1, and obtained their H₂ content in post-processing. We divided each of these samples into cluster and non-cluster galaxies based on the mass of their groups, and examined differences in their H₂ content. We then investigated potential factors that modulate these differences, and their implications for galaxy evolution in clusters. Finally, we used particle tracking to reveal how the evolution of the H₂ content of SF cluster galaxies relates to star formation, molecular gas stripping, gas accretion, stellar/BH feedback, and other miscellaneous processes within the ISM.

The main results from this study are summarized as follows:

• At a fixed stellar mass, SF galaxies in clusters generally possess less H₂ than their non-cluster counterparts, but more H₂ than a typical cluster galaxy (Fig. 2). This provides theoretical support for observational studies that claim the same. We found that the offset in H₂ masses between the two samples increases at lower redshifts at a given stellar mass, and also varies with stellar mass at a given redshift: it is highest for |$M_\star \approx 10^{10.5}\, {\rm M}_\odot$| (upto ≈0.5 dex), but diminishes below |$M_\star \approx 10^{9.3}\, {\rm M}_\odot$|⁠.

• The |$M_\mathrm{H_2}$| histories of galaxies based on their main progenitors show that SF galaxies outside clusters did not possess more H₂ to begin with, and the H₂-deficit in SF cluster galaxies is a consequence of the steep decline of their |$M_\mathrm{H_2}$| relative to non-cluster counterparts around their time of first infall into a cluster (Fig. 3). This decline, and the infall, commenced earlier for most cluster galaxies. The offset between H₂ content of the two samples (shown in Fig. 2) is positively correlated to the time elapsed since infall for the cluster galaxies (see Fig. 4), though there are other factors that modulate the scatter in this trend – e.g. group mass, density at the virial radius, crossing time-scale (Table 2, Section 3.1).

• The evolution of H₂ content and SFR of cluster galaxies with respect to their first infall into a cluster showed that the two followed each other (Fig. 5), implying a fairly invariant star-formation efficiency before and after infall. We found that an average cluster galaxy in eagle with |$M_\star \gtrsim 10^{9.5}\, {\rm M}_\odot$| undergoes a ‘slow-then-rapid’ quenching, whereas those at lower masses quench rapidly. Most cluster galaxies were pre-processed in groups for ≈2 Gyr before entering a cluster, but the SF ones, being centrals, were not. This enabled better preservation of H₂ in SF cluster galaxies.

• By tracking the gas and stellar particles in the main progenitors of SF cluster galaxies, we were able to calculate the changes in H₂ mass that were caused by star formation, gas escape, gas accretion, and variations in the H₂ mass of gas particles in the ISM. We found that, before entering a cluster environment, the H₂ mass of an average SF galaxy was lost mostly to star formation and gas escape, but was compensated to a large extent by other processes; mainly by increase in |$m_{\rm H_2}$| of ISM particles (see Fig. 6). The sharp decline in the H₂ content since the time of infall occurred due to the domination of star formation and escape over the modes that replenished H₂. We compared the energy injected by feedback events against the binding energy of the escaped gas to infer that they were primarily removed via environmental stripping (Section 4.2).

We foresee improvements and exploration beyond the analysis presented in this paper. The limited resolution of eagle precludes self-consistent treatment of all the different phases of the ISM, but modelling these phases is crucial for properly capturing the physics germane to the formation and destruction of molecular gas. Although the post-processing prescription used in this work does account for such processes (Stevens et al. 2019), a self-consistent approach with better mass resolution would be more accurate. The analysis can be repeated with other cosmological hydrodynamical simulations [e.g. Horizon-AGN (Dubois et al. 2014), IllustrisTNG (Pillepich et al. 2018), FIRE-2 (Hopkins et al. 2018), Simba (Davé et al. 2019)], to investigate whether our conclusions are generic to galaxies that form in a ΛCDM cosmology, or if they are specific to eagle. Furthermore, the number statistics of SF cluster galaxies can be improved upon by using a box size larger than 100 Mpc, which would expand the scope to redshifts greater than z = 1 and more massive clusters, and also help make predictions for the high-redshift science that can be carried out with the JWST.

ACKNOWLEDGEMENTS

We would like to thank Aaron Ludlow for useful discussions and insights. AM acknowledges support from the Australian Government through a Research Training Program (RTP) Scholarship. ARHS is the recipient of the Jim Buckee Fellowship at The University of Western Australia. We acknowledge the Virgo Consortium for making their simulation data available. The eagle simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyères-le-Châtel. The data were processed on the OzSTAR supercomputer, which is managed through the Centre for Astrophysics & Supercomputing at the Swinburne University of Technology, and funded by Astronomy Australia Limited and the Australian Commonwealth Government. The analysis was performed with the help of astropy (Astropy Collaboration 2018), h5py, matplotlib (Hunter 2007), numpy (van der Walt, Colbert & Varoquaux 2011), scikit-learn (Pedregosa et al. 2011), and scipy (Virtanen et al. 2020) packages for python. The paper has been typeset using Overleaf.¹⁴

DATA AVAILABILITY

The eagle data used in this work is publicly available at http://icc.dur.ac.uk/Eagle/database.php, and the descriptions for halo catalogues and particle data are provided in McAlpine et al. (2016) and The EAGLE team (2017), respectively. The H₂ masses used in this work can be obtained through a reasonable request to the corresponding author. The observational data are accessible through the links provided in the respective papers.

Footnotes

References

APPENDIX A: Obtaining the H₂ content of gas particles

A1 Computing the neutral hydrogen content

To determine the neutral hydrogen content, we use the Rahmati et al. (2013) prescription that was calibrated using the traphic (Pawlik & Schaye 2008) radiative transfer simulations and accounts for photoionisation due to the extragalactic ultraviolet background (UVB), self-shielding, radiative recombination, and collisional ionisation. We begin by computing the total photoionisation rate Γ_phot using equation (A1) and table A1 of Rahmati et al. (2013). For this, the hydrogen number density (n_H) is determined using the hydrogen fraction (i.e. the fraction of gas mass in hydrogen, X) and the particle mass provided in the eagle data base. We assume the UVB as Haardt & Madau (2012), which is weaker by a factor of 3 compared to the UVB adopted for eagle (i.e. Haardt & Madau 2001), but none the less results in similar masses (Bahé et al. 2016; Crain et al. 2017). In eagle, the temperature of a star-forming gas particle is unrealistic and set by the polytropic equation of state mentioned in Section 2.1. For such particles, we adopt a temperature |$T=10^4\, {\rm K}$| in order to mimic the warm, diffuse ISM around young stellar populations (cf. Crain et al. 2017). Finally, we determine the fraction of hydrogen that is neutral (f_n) by inserting the values for Γ_phot, T and n_H in equation (A8) of Rahmati et al. (2013).

A2 Deriving the contribution from H₂ to the neutral hydrogen content

Next, we partition the neutral hydrogen into atomic and molecular components. We inform our choice of partitioning scheme – used for deriving the atomic and molecular contributions to gas particles in eagle – based on some comparative tests that were done on the five prescriptions mentioned in Section 2.1. In this section, we provide short descriptions for the prescriptions and their implementation.

• BR06: This is a scaling relation between molecular-to-atomic gas ratio and the mid-plane gas pressure (P) that was derived empirically using observations of 14 galaxies. The relation is given as

  $$\begin{eqnarray} R \equiv \frac{\Sigma _{\mathrm{H_2}}}{\Sigma _{\mathrm{H}\, \small {{\rm I}}}} = \left(\frac{P}{P_0}\right)^\alpha \, , \end{eqnarray}$$

  (A1)

  where Σ_x is the surface density of component x, P₀/k_B = 4.3 × 10⁴ cm⁻³ K and α = 0.92. We estimate the gas pressure in each gas particle based on its density, n_H, and temperature, T; i.e. P/k_B = n_HT, where k_B is the Boltzmann constant.

• L08: Similar to BR06, this is also an empirically derived relation between R and P but based on a relatively larger sample of 23 galaxies. It follows the same form as equation (A1) but differs in the value of constants, which are slightly lower at P₀/k_B = 1.7 × 10⁴ cm⁻³ K and α = 0.8.

• GK11: This prescription is based on various ‘fixed ISM’ hydrodynamical simulations that were run using adaptive mesh refinement (AMR), each with a fixed dust-to-gas ratio (same value for all) and interstellar radiation field at 1000 Å. The simulations included gas dynamics, radiative transfer of UV photons from stellar particles, a chemical network of hydrogen and helium with non-equilibrium cooling and heating rates, and molecular hydrogen formation on dust grains.

  There are two sets of fitting formulae in their paper that can be used to derive |$f_{\mathrm{H_2}}$|⁠. The first is their equation (6):

  $$\begin{eqnarray} f_{\mathrm{H_2}} = \frac{1}{1+\exp (-4x-3x^3)}\, , \end{eqnarray}$$

  (A2)

  where

  $$\begin{eqnarray} x \equiv \Lambda ^{3/7}\ln \left(D_{\rm MW} \frac{n_{\rm H}}{\Lambda n_*}\right)\, . \end{eqnarray}$$

  (A3)

   n_* = 25 cm⁻³, D_MW is the dust-to-gas ratio relative to the Milky Way value, and Λ is given by

  $$\begin{eqnarray} \Lambda \equiv \ln \left[1+gD_{\rm MW}^{3/7}\left(\frac{U_{\rm MW}}{15}\right)^{4/7}\right]\, . \end{eqnarray}$$

  (A4)

  Here, U_MW is the normalisation of the interstellar far-UV (FUV) flux at 1000 Å, g is a fitting formula used to account for the transition between the self-shielding regime (g ≈ 1) to dust shielding regime (⁠|$g\propto D_{\rm MW}^{-1}$|⁠):

  $$\begin{eqnarray} g \equiv \frac{1+\alpha s+s^2}{1+s}\, , \end{eqnarray}$$

  (A5)

  where

  $$\begin{eqnarray} \alpha \equiv \frac{2.5 U_{\rm MW}}{1+(0.5U_{\rm MW})^2}\, , \end{eqnarray}$$

  (A6)

   

  $$\begin{eqnarray} s \equiv \frac{0.04}{D_\star +D_{\rm MW}}\, . \end{eqnarray}$$

  (A7)

   D_⋆ denotes the transition point where the formation of H₂ mainly occurs through gas-phase reactions:

  $$\begin{eqnarray} D_\star \equiv 0.0015\ln [1+(3U_{\rm MW})^{1.7}]. \end{eqnarray}$$

  (A8)

  The second approach is given by their equation (10):

  $$\begin{eqnarray} f_{\mathrm{H_2}} = \left(1+\frac{\Sigma _{\mathrm c}}{\Sigma _{\mathrm{H}\, \small {{\rm I}}+\mathrm{H_2}}}\right)^{-2}, \end{eqnarray}$$

  (A9)

  where Σ_c is the characteristic neutral-gas surface density where the relationship between SFR surface density and molecular hydrogen surface density steepens. It is given by equation (14) in their paper, i.e.

  $$\begin{eqnarray} \Sigma _{\rm c} = \frac{20\Lambda ^{4/7}}{D_{\rm MW}\sqrt{1+U_{\rm MW}D^2_{\rm MW}}}. \end{eqnarray}$$

  (A10)

• K13: This is an extension of the theoretical model developed by Krumholz, McKee & Tumlinson (2008, 2009a, b) for accurate treatment of |$\mathrm{H}\, \small {{\rm I}}$| dominated regions. Under this model, a galactic gas disc is partitioned between |$\mathrm{H}\, \small {{\rm I}}$| (non star-forming) and H₂ (star-forming) based on the gas column density, metallicity and the ratio of the interstellar radiation field (ISRF) to the cold atomic ISM density. If the pressure implies co-existence of warm and cold atomic phases, the ratio of ISRF to density is nearly constant, and the conversion is only a function of column density and metallicity. For other regions with low ISRF and long depletion time-scales, the ratio is not fixed and there is a density floor imposed in order to maintain hydrostatic equilibrium.

  Similar to Garratt-Smithson et al. (2021), we apply equation (10) of Krumholz (2013) to determine molecular gas fractions of our gas particles:

  $$\begin{eqnarray} f_{\mathrm{H_2}} = \left\lbrace \begin{array}{lr}1 - 3\, S\, (4 + S)^{-1} & \forall ~S\lt 2\\ 0 & \forall ~S \ge 2 \end{array} \right. \, , \end{eqnarray}$$

  (A11)

  where

  $$\begin{eqnarray} S \equiv \frac{\ln (1 + 0.6\, \chi + 0.01\, \chi ^2)}{0.6\, \tau _{\rm c}}\, , \end{eqnarray}$$

  (A12)

   

  $$\begin{eqnarray} \tau _{\rm c} \equiv 0.066\, f_{\rm c}\, D_{\rm MW}\, \frac{\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}}{{\rm M}_{\odot }\, {\rm pc}^{-2}}\, , \end{eqnarray}$$

  (A13)

   

  $$\begin{eqnarray} \chi \equiv 72\, U_{\rm MW} \left(\frac{n_{\rm CNM}}{{\rm cm}^{-3}}\right)^{-1}. \end{eqnarray}$$

  (A14)

  Here, f_c is a clumping factor to account for the difference between the physical scale of an atomic–molecular complex, and the scale at which it is observed, i.e. spatial resolution of the simulation (for eagle, f_c = 5). n_CNM is cold neutral medium (CNM) density given by

  $$\begin{eqnarray} n_{\rm CNM} \equiv {\rm max}[n_{\rm CNM,2p},\, n_{\rm CNM,hydro}]\, , \end{eqnarray}$$

  (A15)

  where n_CNM,2p is the CNM density for a two-phase equilibrium between cold and warm neutral medium, i.e.

  $$\begin{eqnarray} n_{\rm CNM,2p} \equiv 23\, U_{\rm MW}\, \frac{4.1}{1 + 3.1\, D_{\rm MW}^{0.365}}\, {\rm cm}^{-3}\, , \end{eqnarray}$$

  (A16)

  and n_CNM,hydro is the lowest possible density for a CNM in hydrostatic equilibrium corresponding to the maximum temperature of T_CNM,max = 243 K (Wolfire et al. 2003); i.e.

  $$\begin{eqnarray} n_{\rm CNM,hydro} \equiv \frac{P_{\rm th}}{1.1\, k_{\rm B}\, T_{\rm CNM,max}}\, . \end{eqnarray}$$

  (A17)

  Here, the factor of 1.1 accounts for He, and P_th is the thermal pressure. The mid-plane pressure (P_mp) can be expressed as a sum of three contributions from (1) the self-gravity of the |$\mathrm{H}\, \small {{\rm I}}$| that is not gravitationally bound, (2) the weight of the |$\mathrm{H}\, \small {{\rm I}}$| within the gravitational potential from the bound H₂ clouds, and (3) the weight of the |$\mathrm{H}\, \small {{\rm I}}$| within the gravitational potential provided by stars and dark matter. P_th can be written as a function of the sound speed in the WNM (c). P_th is smaller than P_mp by a factor of α ≈ 5 due to additional support from magnetic fields, cosmic rays, and turbulence. Based on the expressions in Krumholz (2013), this implies

  $$\begin{eqnarray} P_{\rm th} = \frac{\pi \, G\, \Sigma _{\rm H\, {{\small {\rm I}}}}^2}{4\, A} \left(1 + 2\, R + \sqrt{(1 + 2\, R)^2+\mathcal {F}} \right)\, , \end{eqnarray}$$

  (A18)

  where

  $$\begin{eqnarray} \mathcal {F} \equiv \frac{32\, \zeta _d\, A\, f_w\, c_w^2\, \rho _{\rm sd}}{\pi \, G\, \Sigma _{\rm H\, {{\small {\rm I}}}}^2}\, , \end{eqnarray}$$

  (A19)

  and

  $$\begin{eqnarray} R \equiv \frac{f_{\rm H_2}}{1-f_{\rm H_2}}\, . \end{eqnarray}$$

  (A20)

  ρ_sd is the density of DM and stars within the gas disc, and ζ_d ≈ 0.33 is a numerical factor that depends on the shape of the gas surface density profile.
• GD14: Gnedin & Draine (2014) provides two equations for deriving the H₂ content [equations (6) and (8)], but only equation (8) is appropriate for the spatial resolution of eagle (≈0.7 kpc). It gives

  $$\begin{eqnarray} R = \left(\frac{\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}}{\Sigma _{R=1}} \right)^{\alpha }, \end{eqnarray}$$

  (A21)

  where |${\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}}$| is the neutral-hydrogen surface density, and Σ_{R = 1} is the surface density at which the atomic and molecular hydrogen fractions are equal:

  $$\begin{eqnarray} \Sigma _{R=1} = \frac{50\, \sqrt{0.001 + 0.1\, U_{\rm MW}}}{g\, (1 + 1.69\, \sqrt{0.001 + 0.1\, U_{\rm MW}})}\, {\rm M}_{\odot }\, {\rm pc}^{-2}. \end{eqnarray}$$

  (A22)

  Here, g is determined as

  $$\begin{eqnarray} g \equiv \sqrt{D_\star ^2 + D_{\rm MW}^2}\, , \end{eqnarray}$$

  (A23)

  where D_MW is the dust-to-gas ratio relative to the Milky Way value and

  $$\begin{eqnarray} D_\star \equiv 0.17\, \frac{2 + S^5}{1 + S^5}\, . \end{eqnarray}$$

  (A24)

   S is the sampling length of the simulation in units of 100 pc. α is calculated as

  $$\begin{eqnarray} \alpha \equiv 0.5 + \left(1 + \sqrt{\frac{U_{\rm MW}\, D_{\rm MW}^2}{600}} \right)^{-1}. \end{eqnarray}$$

  (A25)

  Thus, the value of R is essentially set by the values of |$\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}$| and S.

For implementing GK11, K13, and GD14, we follow the computational approach of Stevens et al. (2019). The neutral-hydrogen surface density can be expressed as

where Σ is the total-gas surface density. Since Σ is a two-dimensional quantity, whereas a gas particle in eagle possesses a three-dimensional density (ρ), we need to derive Σ by multiplying ρ with a measure of the particle’s size. To this end, we assume the gas to be thermalised (Schaye & Dalla Vecchia 2008) and approximate its size as the Jeans length for an ideal gas, L (Schaye 2001). This can be written as

where γ is the ratio of specific heat capacity at constant pressure to that at constant volume, k_B is the Boltzmann constant, T is the particle temperature, u is the specific internal energy, and G is the gravitational constant. Assuming that all gas is thermalised, we approximate

and

where f_mol is the ratio of H₂ to the total non-metallic component (H + He), i.e.

|$f_\mathrm{H_2}$| is the fraction of neutral hydrogen in molecular phase.

As both γ and u are dependent on |$f_\mathrm{H_2}$|⁠, solving equation (A26) for |$\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}$| would need prior knowledge of |$f_\mathrm{H_2}$|⁠. We compute |$f_\mathrm{H_2}$| iteratively, where it is initialised as 0 for the first iteration and used to determine f_mol [equation (A30)]. This is used to get u and γ from equations (A29) and (A28) which, along with equations (A27) and (A26), give |$\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}$|⁠. Then, we obtain |$f_\mathrm{H_2}$| using the equations in the prescription, where we take the |$\Sigma _{\rm H\, {{\small {\rm I}}}+H_2}$| just derived, S = 0.01L, compute U_MW based on Diemer et al. (2018), and assume D_MW to scale with metallicity and – similar to Lagos et al. (2015) – equal to the ratio of smoothed metallicity of the gas particle (Z) to the solar metallicity (Z_⊙ = 0.0127; Allende Prieto, Lambert & Asplund 2001). This is repeated until |$f_\mathrm{H_2}$| converges within 0.5 per cent of its value in the previous iteration, or else, if there have been 300 iterations.

APPENDIX B: IDENTIFYING VALID POST-PROCESSING PRESCRIPTIONS

Now, we elaborate on the tests that were performed on these prescriptions. This includes a comparison against molecular gas estimates for z = 0 galaxies, and the cosmic H₂ density across redshifts.

B1 Comparison against observations at z = 0

We compare the H₂ masses of z = 0 galaxies predicted by the five models against the observed values from the xCOLD GASS survey (Saintonge et al. 2017), which is an extension of the COLD GASS survey (Saintonge et al. 2011) to lower stellar masses and redshifts. It includes CO (1–0) observations of 532 galaxies at 0.01 < z < 0.05. Amongst these, 366 galaxies (with |$M_\star \gt 10^{10}\, {\rm M}_\odot$| and 0.025 < z < 0.050) are from COLD GASS, and the remaining 166 galaxies (with 10⁹ < M_⋆/M_⊙ < 10¹⁰ and 0.01 < z < 0.02) were randomly selected from the SDSS group sample. This data set serves as an ideal z = 0 benchmark for galaxy evolution studies, as it is an unbiased and representative selection based purely on redshift and stellar mass, with an even sampling of the stellar mass space.

Fig. A1 shows a scatter plot of |$M_{\mathrm{H_2}}$| against M_⋆ for z = 0 eagle galaxies with stellar masses above |$M_\star =10^9\, {\rm M}_\odot$|⁠, alongside galaxies in xCOLD GASS. The results from eagle are shown as a two-dimensional density plot, such that the density of galaxies increases from white to black. The black contours enclose 68 and 95 per cent of the eagle data. The detections and upper limits from xCOLD GASS are shown as filled and open circles, respectively.

The prescriptions of BR06 and L08 predict considerably lower H₂ masses compared to observations. This discrepancy is potentially due to simplifying assumptions used by these models – e.g. they do not account for relevant physical processes such as ionisation due to the local UV field, or the formation of molecular gas on dust grains. Also, the gas disc in the simulation has a larger scale height than expected (Bahé et al. 2016; Benítez-Llambay et al. 2018), which causes the thermal pressure in high-density gas and H₂ content to be lower than a galaxy with a resolved cold phase (Diemer et al. 2018). This also explains, at least in part, the good agreement of |$\mathrm{H}\, \small {{\rm I}}$| masses from these schemes with observations reported by Manuwal et al. (2022). Considering these factors, we choose to avoid these prescriptions.

We find that the two equations in GK11 [equations (A2) and (A9)] produce vastly different results, which indicates that this prescription may not be reliable for eagle. The K13 and GD14 prescriptions produce H₂ masses that are in reasonable agreement with the data from xCOLD GASS, and are also consistent with each other. These prescriptions also produce sensible |$\mathrm{H}\, \small {{\rm I}}$| masses (Manuwal et al. 2022). This is reasonable, as these schemes are based on detailed analytical/numerical treatment of key physical processes, e.g. formation of H₂ on dust grains, and ionisation due to interstellar radiation field.

B2 Evolution of the cosmic H₂ density

Additionally, we estimate the co-moving cosmic H₂ density in eagle for redshifts spanning z = 0 to z = 4.5¹⁵ and compare them against the compilation of observed cosmic H₂ densities from ALMACAL (Klitsch et al. 2019), ASPECS LP (Decarli et al. 2020; Magnelli et al. 2020), Berta et al. (2013), COLDz (Riechers et al. 2019), Garratt et al. (2021), Maeda, Ohta & Seko (2017), PHIBBS2 (Lenkić et al. 2020), Scoville et al. (2017), VLASPECS (Riechers et al. 2020), and xCOLD GASS (Fletcher et al. 2021).

Fig. B1 shows the H₂ density in eagle as a function of redshift as predicted by the K13 and GD14 models (blue and violet curves, respectively), along with the observational values. Note that observed values carry large uncertainties and span a wide density range at any given redshift. If we take observations at face value, a comparison with the results from eagle suggests that both the K13 and GD14 prescriptions produce H₂ densities that are below the majority of observed values – especially at z = 1, where the offset is ≈0.4 dex. We are unsure if this is primarily due to the limitations of the models or the observations, as observational estimates are likely to carry additional (systematic) uncertainties beyond the reported errors. For example, most of the surveys are biased towards gas-rich galaxies, especially at high redshifts. There is also uncertainty related to the CO-to-H₂ conversion factor (α_CO) used to estimate H₂ masses from CO flux. Nevertheless, we note that predicted H₂ densities from eagle lie within the range spanned by the observed values at all redshifts. The shape of the curves is also consistent with the general trend in observations: an increase in the cosmic density until z ≈ 1.6, and a decline at higher redshifts. We interpret this as an additional support for the reliability of these partitioning schemes.