Exploring the angular momentum – atomic gas content connection with eagle and IllustrisTNG

ABSTRACT

We use the Evolution and Assembly of GaLaxies and their Environments (eagle) and IllustrisTNG (The Next Generation) cosmological simulations to investigate the properties of the baryonic specific angular momentum (j), baryonic mass (M), and atomic gas fraction (f_atm) plane for nearby galaxies. We find EAGLE and TNG to be in excellent agreement with each other. These simulations are also consistent with the results obtained with eXtended GALEX Arecibo SDSS Survey (xGASS) for gas fractions greater than 0.01. This implies that the disagreements previously identified between xGASS and predictions from simple analytical disc stability arguments also holds true for eagle and tng. For lower gas fraction (the regime currently unconstrained by observations), both simulations deviate from the plane but still maintain good agreement with each other. Despite the challenges posed by resolution limits at low gas fractions, our findings suggest a potential disconnect between angular momentum and gas fraction in the gas-poor regime, implying that not all gas-poor galaxies have low specific angular momentum.

1 INTRODUCTION

Angular momentum is a key property of galaxies, as it is linked to their formation and evolutionary history. It is now known that the stellar angular momentum scales with mass (the so-called Fall relation; Fall 1983; Fall & Romanowsky 2013; Lapi, Salucci & Danese 2018; Posti et al. 2018a, b; Sweet et al. 2018; Stone, Courteau & Arora 2021; Du et al. 2022; Di Teodoro et al. 2023), and that the scatter in the relation is correlated with morphology and stellar structure (e.g. Fall 1983; Romanowsky & Fall 2012; Cortese et al. 2016; Pulsoni et al. 2023). However, it is still unclear whether stellar structure is the primary physical driver of the scatter in the Fall relation or simply a proxy of the overall accretion history of galaxies. Indeed, from a theoretical point of view, the growth of angular momentum of discs should be tightly connected to their ability to accrete gas (e.g. Mo, Mao & White 1998; Boissier & Prantzos 2000), potentially implying a more fundamental role of cold gas in driving the scatter of the relation. This is also consistent with theoretical work focused on the link between gas content and disc stability (e.g. Obreschkow et al. 2016; Stevens et al. 2018; Romeo 2020). The well-motivated theory behind the connection of a galaxy’s gas content and specific angular momentum has been empirically tested with modest observational samples, but both the sample size and manner in which those data are analysed can be improved.

The advent of new data sets allows more detailed investigations of the role of gas content in the build-up of angular momentum in galaxies. In particular, Hardwick et al. (2022a, hereafter Paper 1) used the eXtended GALEX Arecibo SDSS Survey (xGASS; Catinella et al. 2010, 2018) sample, a deep H i survey which is representative of the local Universe, to study the stellar Fall relation. We found that the most strongly correlated parameter with the scatter of the Fall relation is H i gas fraction, not bulge-to-total ratio, particularly for low stellar masses and when isolating the disc component of galaxies. In Hardwick et al. (2022b, hereafter Paper 2), we expanded this work and investigated the connection between a galaxy’s baryonic angular momentum, baryonic mass and atomic gas fraction, and found that a tight plane exists between these three parameters. However, this plane deviates from predictions from simple analytical models of disc stability (e.g. Obreschkow et al. 2016; Romeo 2020), which predict a steeper slope and, most importantly, a stronger dependence on gas fraction. A similar result was also found by Mancera Piña et al. (2021a, b) studying a sample of ∼100 disc galaxies with resolved H i rotation curves, suggesting that more detailed modelling is needed to fully unveil the physical connection between mass, angular momentum, and cold gas content.

The natural next step is to extend such a comparison between xGASS and theoretical models to cosmological hydrodynamical simulations, which are not limited by simplifying assumptions of the models above. While recent years have seen a dramatic increase in the number of studies focused on the origin and drivers of the scatter of the Fall relation (Genel et al. 2015; Teklu et al. 2015; Lagos et al. 2017; Zoldan et al. 2018; Marshall et al. 2019; Stevens et al. 2019; Wang et al. 2019; Elson, Głowacki & Davé 2023). In particular, qualitative results from Rodriguez-Gomez et al. (2022) suggest that observational trends from Paper 1 are reproduced in IllustrisTNG. However, no work to date has explicitly investigated the baryonic specific angular momentum (j_bar) – baryonic mass (M_bar) – atomic gas fraction (f_atm) plane (hereafter, the JMG plane), quantified its slope and scatter, and carried out a detailed comparison with observations. This is critical not only to obtain some insights into the physics driving this relation, but also to investigate how universal this JMG plane really is.

Despite the xGASS sample showing a greater diversity of galaxies than surveys of just disc galaxies, the subsample for which angular momentum can be estimated from H i line widths is still biased towards galaxies that have H i contents above the detection limit of the survey. Additionally, although xGASS is a large sample in the context of similar observational samples, the number of galaxies is quite small when compared to large simulation volumes. Therefore, in this work, we wish to further test the robustness of this JMG plane with even greater statistics and galaxy diversity by using large hydrodynamical cosmological simulations. First, we aim to test how well the observations and simulations agree in this parameter space. Then we can use the simulation data to better explore the physical connection between these three parameters. As these simulations are not limited to observational constraints, they allow us to probe down to low gas fractions and weaker rotational support than is possible with observational surveys.

This paper is set out as follows. We start by describing the archival observational and simulation data used in this work from xGASS, eagle (Evolution and Assembly of GaLaxies and their Environments; Crain et al. 2015; Schaye et al. 2015), and IllustrisTNG (Illustris The Next Generation, hereafter, tng; Nelson et al. 2018; Pillepich et al. 2018) in Section 2. We then explain how we create mock detection samples for the simulations in Section 3. The results of comparing eagle and tng to xGASS in the M_bar–j_bar–f_atm parameter space is presented in Section 4. We then discuss the implications of these results in Section 5 and conclude in Section 6.

2 DATA DESCRIPTION

2.1 xGASS

Our observational data set comes from xGASS (Catinella et al. 2010, 2018), which includes galaxies in the stellar mass range of 10⁹ to |$10^{11.5}\,\rm {M}_{\odot }$| across the redshift range 0.01 < z < 0.05 and was selected from SDSS DR6 (Adelman-McCarthy et al. 2008). Galaxies were observed with the Arecibo telescope until detected in H i or until a gas fraction limit of 2–10 per cent was reached. The sample, in particular at high stellar masses, was selected to have a nearly flat stellar mass distribution. Overall, xGASS represents arguably the best representative sample of integrated gas properties in the local Universe. In this work, we use the same subsample of 564 xGASS galaxies that were used in Papers 1 and 2. Briefly, this includes only galaxies detected in H i and for which we could accurately determine kinematics; i.e. an inclination greater than 30 deg and not affected by confusion within the radio beam. We use the stellar and baryonic specific AM that was calculated in Papers 1 and 2, respectively, which are publicly available.¹ These were estimated by combining stellar mass surface density profiles (within a 10R_e aperture) with H i widths. AM was calculated within a large aperture to ensure convergence for galaxies with high Sérsic indices (Sérsic 1963).

2.2 eagle and tng

In this work, we compare our observational results to simulation data from both the eagle (Crain et al. 2015; Schaye et al. 2015) cosmological hydrodynamical simulation and tng (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018) magnetohydrodynamical cosmological simulation. For comparisons to eagle, we use the ‘reference’ eagle model and AM values determined in Lagos et al. (2017), and for tng we use the quantities obtained by Stevens et al. (2019).

For consistency, we use the 100 comoving Mpc simulation box for both eagle and tng, as they have a comparable number of dark matter particles (1504³ and 1820³, respectively). eagle is simulated with smoothed particle hydrodynamics using the gadget-3 code (Springel 2005; Springel et al. 2008), while tng has discretized gas elements within a moving Voronoi mesh implemented using the arepo code (Springel 2010).

There are various advantages of simultaneously comparing both works to xGASS. First, while both simulations include subgrid models that calculate feedback (from stars and accreting black holes), gas cooling, star formation and black hole growth, the details of the modelling are dramatically different with, for example, the AGN feedback models have significant differences in their physical bases as well as numerical implementations (see Crain et al. 2015; Schaye et al. 2015; Nelson et al. 2018; Pillepich et al. 2018 for further details). Secondly, the way the subgrid models are calibrated is different. eagle calibrates its subgrid models on two scaling relations; the z = 0.1 galaxy stellar mass function and the size–mass relation of disc galaxies. In contrast, the primary scaling relations used to calibrate tng subgrid models are the cosmic SFR density history, z = 0 galaxy stellar mass function and the z = 0 stellar–halo mass relation, with additional scaling relations used as secondary constraints (the black hole–bulge mass relation, the gas fraction of haloes within R_500c and the stellar size–mass relation, all at z = 0). Thirdly, the way Lagos et al. (2017) and Stevens et al. (2019) calculated AM is different, as we describe below.

Both Lagos et al. (2017) and Stevens et al. (2019) calculated AM and associated quantities, for the same selection of galaxies; z = 0 and |$M_{\star } \gt 10^9 \rm {M}_{\odot }$| (10 803 galaxies for eagle and 20 876 for tng). eagle and tng do not model atomic gas directly; instead, they determine it in post-processing. We use the H i gas mass from Lagos et al. (2015) and Stevens et al. (2019) that were calculated using the Krumholz (2013) theoretical model. This model determines molecular hydrogen fractions as a function of the total column density of neutral hydrogen, metallicity, and the density of the stellar disc, and then uses this to infer atomic hydrogen content (for more details of this model see Krumholz 2013). f_atm is then defined to be 1.35 M_H i/M_bar, the same as our definition for xGASS (the factor of 1.35 adds the approximate contribution from helium). We calculate the baryonic AM as the sum of the AM from stars, H i and H₂. In both Lagos et al. (2017) and Stevens et al. (2019), AM is calculated within apertures to be comparable to observations. From Lagos et al. (2017), we use the AM calculated within 5R_e (half-mass radius of stars). In Stevens et al. (2019), they calculate AM within what they define as the ‘BaryMP’ radius, which is the radius where the gradient of the cumulative baryonic mass profile converges (Stevens et al. 2014). This BaryMP radius is 6.9R_e on average.

In Paper 1, for xGASS, we calculated AM within an aperture of 10R_e to ensure all of our galaxies had their AM converged. As R_BaryMP varies for each galaxy to a radius where the baryonic mass converges, the AM values determined for tng will also likely be converged and comparable with xGASS. Additionally, if the equations of AM are solved analytically with a single Sérsic profile (Sérsic 1963), a galaxy with a Sérsic index less than 2 will have its AM converge by 5R_e (Paper 1). Therefore, for the majority of galaxies, the AM determined within 5R_e or 10R_e will be comparable and we assume that the values determined in Lagos et al. (2017) are appropriate to compare with our xGASS values.

As we will show in Section 4, despite these differences, the agreement between the two simulations and xGASS is striking, highlighting how none of the differences in the way the key parameters investigated here significantly affect our analysis. To reiterate, these results would remain unchanged even if eagle AM were calculated in the same way as tng. This is because the mock detection samples are dominated by galaxies with a significant disc component, for which AM has already converged to the total value at 5R_e (see also Paper 1).

3 MOCK DETECTION SAMPLE

The left column of Fig. 1 shows the H i gas fraction of all galaxies in eagle (top row) and tng (bottom row).

The medians for each simulation are represented by blue and magenta dashed lines (for eagle and tng, respectively). These can be compared to the weighted median of xGASS in black (these are taken from Catinella et al. 2018 and in this figure are weighted to recover a volume-limited sample and include H i non-detections). Both eagle and tng agree with xGASS for intermediate stellar masses (∼10¹⁰ M_⊙ to 10¹¹ M_⊙), but eagle galaxies at low stellar masses are more gas–poor than xGASS, while high–stellar–mass tng galaxies are slightly more gas–rich than xGASS (consistent with what was found by Davé et al. 2020). However, it should be noted that the xGASS sample selection is very different from that of the simulations. The xGASS sample was deliberately selected from Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009) to oversample high stellar masses, which is different to the volume-limited samples of both eagle and tng. This can be seen in the stellar mass distributions shown in the left column of Fig. 2, where the full xGASS sample is shown as the hollow black histograms, while eagle and tng are blue and magenta, respectively.

To carry out a more fair comparison between xGASS, eagle, and tng, we extract from the simulations a sample that has the same gas fraction limit and stellar mass distribution as xGASS. This can be seen in the right column of Figs 1 and 2. This reduces our eagle sample from 10 803 to 7037 galaxies and our tng sample from 20 876 to 14 919 galaxies (although it should be noted that 26 eagle galaxies and 3201 tng galaxies had no mass in H i, so are not shown in the left panel of Fig. 1). Once the H i detection cut is applied, unsurprisingly the agreement between xGASS and the simulations improves. In the right-hand column of Fig. 1 both samples follow approximately the same distribution with differences now reduced to 0.4 dex or less, primarily at low stellar masses. The right-hand column of Fig. 2, eagle and tng now overlap with the xGASS subsample distribution by construction.

Throughout this paper when comparing xGASS to eagle/tng, we will show both the full sample and the ‘mock detection sample’ (H i detection cut and stellar mass weighting). This allows us to distinguish between physical differences and selection effects.

4 RESULTS

4.1 Fall relation

We first compare tng and eagle to xGASS observations in the mass–specific AM relation (Fall relation; Fall 1983). In Fig. 3, we show the stellar Fall relations (panel a) and baryonic Fall relations (panel b) for both simulations.

The left columns show all of the simulated galaxies without a selection cut, while the right columns are only the galaxies in the simulation above the xGASS detection limit (mock detection sample). Each axis shows the simulation running median as the coloured dashed line, and a comparison to the xGASS observational sample in black.

We present both columns, to show the effect the sample selection has in this parameter space. In both the stellar and baryonic cases, the mock detection sample has fewer low AM galaxies than the full samples for both simulations. This results in a tighter distribution of galaxies around the median for the mock detection sample. The median AM for low stellar mass galaxies is also higher for the mock detection sample than for the full sample.

In the right column, both simulations show good agreement with the xGASS median for M_⋆ > 10^10.3 M_⊙ (i.e. less than 0.1 dex difference), which then increases to a maximum discrepancy at M_⋆ = 10⁹ M_⊙ (0.35 dex for tng and 0.42 dex for eagle). We note that at stellar masses below ∼10¹⁰ M_⊙ the resolution of both simulations has an impact on the AM measurements determined (Wilkinson et al. 2023). Therefore, the disagreement at low stellar masses should not be overinterpreted, and we conclude that the agreement between xGASS and the simulation data are reasonable for the mock detection sample.

We also show the baryonic Fall relation in Fig. 3(b). We note a better agreement between xGASS and the simulation data for the baryonic Fall relation than the stellar Fall relation. Without any selection cuts (left column) for eagle, there is a maximum discrepancy of 0.42 dex at 10¹¹ M_⊙, and less than 0.2 dex difference for M_bar < 10^10.6 M_⊙. tng has less than 0.2 dex discrepancy at all baryonic masses. When the mock detection sample is considered (right column), any disparities between either of the simulations and the xGASS observations are effectively eliminated. Specifically, tng and eagle have discrepancies of less than 0.1 dex and less than 0.2 dex at all baryonic masses, respectively. This better agreement for the baryonic Fall relation is especially important for this work, as the remainder of the analysis will focus on the baryonic j.

4.2 j_bar–M_bar–f_atm plane

In this subsection, we now compare both eagle and tng to the j_bar–M_bar–f_atm plane, which was found for xGASS in Paper 2. The top row of Fig. 4 shows j_bar against M_bar in four evenly log-spaced f_atm bins. Each column shows medians (in 0.3 dex M_bar bins) of the eagle and tng galaxies in that atomic gas fraction bin. The median in each bin is shown by a blue (eagle) or magenta (tng) line, and the shaded regions show the range of the 16th to 84th percentile for each bin. For comparison, the black lines show the xGASS JMG plane at fixed gas fractions. We show these lines for the gas fractions at the bin edges of each of the columns (f_atm = 0.01, 0.03, 0.1, 0.3 is solid, dashed, dotted and dot–dashed, respectively). In the bottom row of Fig. 4, we show the residual in dex of these medians with respect to the mid-point of the xGASS JMG plane lines.

In Fig. 4, we only show the mock detection sample, as, in this projection and gas fraction range, there is very little difference between the full simulation sample and the mock detection sample. The mock detection sample is simply a cut in H i gas fraction with a stellar mass weighting applied. Therefore, the only difference that can be seen between the two samples, is more galaxies in the first column for the full simulation samples that extend to lower baryonic masses. However, for completeness in appendix  A we also show the spread of the galaxies in this parameter space as two-dimensional (2D) histograms for both the full sample and the mock detection sample (eagle and tng are shown as Figs A1 and A2, respectively).

Regardless of whether the full sample or mock detection samples are used, for gas fractions greater than 0.01, the JMG plane is in excellent agreement for both the simulations and xGASS. The two middle panels of Fig. 4 show simulation medians with the same slope as the xGASS plane. The left-hand and right panels have slopes that are slightly shallower than the xGASS plane but are consistent within their scatter. In the left panel, this slight difference in slope is simply due to low statistics in this f_atm bin, which can be seen by the large shaded region. In the right panel, different values of f_atm dominate at various baryonic masses, leading to a change in slope. Specifically, galaxies with f_atm ≈ 1 are predominantly within the range of 9.5 < log₁₀(M_bar/M_⊙) < 10. The peak at log₁₀(M_bar/M_⊙) ≈ 9.8 is where these galaxies dominate the median and drive up the average j_bar at this location. In all panels, the simulation medians have an offset of less than 0.2 dex from the allowed region of the xGASS JMG plane.

There are small offsets in normalization, but these are consistent with the offsets observed in gas fraction and j in Figs 1 and 3. In many ways, the agreement between eagle, TNG, and xGASS in this parameter space is remarkable. First, it is interesting that both simulations are in good agreement with each other (maximum difference of 0.23 dex between the two simulations’ medians) given that they rely on different codes and subgrid physics prescriptions. Often these differences will result in the simulations having slightly different predictions, such as the H i fraction – stellar mass relation (Fig. 1), where low stellar mass eagle galaxies are more gas poor than tng galaxies. However, when gas fraction, baryonic mass, and baryonic specific AM are all considered together, as we have in Fig. 4, there is very little difference between the two simulations’ predictions. Secondly, it is also intriguing how well these simulations agree with the xGASS JMG plane. In particular, this agreement is strongest when considering j_bar, M_bar, and f_atm together, rather than when isolating either the H i– stellar mass relation (Fig. 1) or Fall relation (Fig. 3) separately.

It should be emphasized that, despite eagle and tng having their subgrid models calibrated against many observational scaling relations (Crain et al. 2015; Schaye et al. 2015; Nelson et al. 2018; Pillepich et al. 2018), this is likely not the cause of the tight agreement seen between the simulations and observations in Fig. 4. First, although both simulations are calibrated to reproduce the observed stellar mass–size relation, which is closely linked to the Fall relation, Figs 3(a) and (b) show that the simulated Fall relations have larger discrepancies between the simulations and observations than is seen for the JMG plane in Fig. 4. Therefore, the simulations calibration to reproduce the stellar mass–size relation is not the sole cause of the tight agreement between the JMG plane and eagle/tng. Secondly, neither of these simulations are calibrated to reproduce cold gas content, as the majority of the subgrid models calibrate for stellar content. Therefore, the tight agreement seen between simulations and observations in the j_bar, M_bar, and f_atm plane is not due only to the simulations calibration, and is instead a prediction of the models.

We also chose to fit a JMG plane directly to the simulation data using the hyper-fit (Robotham & Obreschkow 2015) Bayesian hyperplane fitting tool, as we did for xGASS data in Paper 2. As the full simulation data are heavily skewed in their baryonic mass and AM distribution, we choose to only fit a JMG plane to the mock detection sample (we explore this skewness more in Section 4.3). The best-fitting parameters of a JMG plane with the form

are shown in Table 1.

These values can be compared to those found for Paper 2 in the bottom row of the table. We see that the tng JMG plane parameters are within errors of the xGASS JMG plane (except for the β parameter), while the eagle simulation parameters deviate by more than 3σ. Although, it should be noted that the errors provided on these parameters are uncertainties from the MCMC (Markov Chain Monte Carlo) chain and do not incorporate any uncertainties on individual galaxy values, so will likely be an underestimate of the true error. The projection of both of these fits are shown in Figs A1 and A2, where it is clear that the best fit is not always an accurate representation of the data distribution, in particular at low gas fraction, so that median values provide a more fair comparison. As we will see in the next section, this is also because at low gas fractions the JMG plane may no longer be able to properly describe the distribution of galaxies in both EAGLE and TNG. We can also compare the standard deviation of galaxies from the JMG plane in the j direction, which is given in the right column of Table 1. The spread of galaxies is the smallest in xGASS, with tng being marginally larger and eagle being ∼50 per cent larger. In an attempt to better understand these differences, in the next subsection, we look into the scatter around the JMG plane in more depth.

4.3 Scatter from the JMG plane

In addition to studying the shape and slope of the JMG plane, it is also important to investigate the scatter around it. The offsets of galaxies from the xGASS JMG plane in the j direction are shown in the left and middle panels of Fig. 5. The only difference between the two panels is that the left panel shows all galaxies within either the eagle or tng simulations (blue and magenta histograms, respectively), whereas the middle panel shows only galaxies that are in the mock detection samples.

When all galaxies in the simulations are considered, the offsets are strongly positively skewed. This is due to the simulations containing galaxies that have extremely low gas fractions, while still maintaining a moderate baryonic AM. In fact, the majority of the galaxies in this long extended tail have gas fractions that are at or below the limit of where they would be considered accurate. One way to determine if a gas fraction is ‘accurate’ is to count the number of gas particles each contributing at least 5 per cent of their mass to neutral gas. The blue hollow histogram shows the result of excluding all galaxies with less than 10 gas particles that reach this criterion for EAGLE. This illustrates that only a small fraction of galaxies in this extended tail have a sufficient number of gas particles to be reliable. Therefore, when considering the tail of the entire sample, the precise position of a galaxy with respect to the plane should be regarded as an approximation. We will further explore this result in Section 5.

Although this extended tail is dominated by galaxies with uncertain gas masses, there is still a small number of galaxies with reliable gas masses and high offsets. This could indicate that despite gas fraction being a strong predictor of a galaxy’s baryonic AM for galaxies with a gas fraction greater than ∼0.01, this relationship breaks down for galaxies with lower gas fractions. In practice, this is unsurprising as, by construction, the dominant mass component will set the baryonic AM of a galaxy. In the gas-poor regime, the gas is no longer the dominant mass component; therefore, the JMG plane will be inaccurate. However, the interesting point to note is that not all galaxies in eagle and tng will be slow rotators once they have depleted their gaseous reservoir. We will explore this more in Section 5, but to summarize, this has two implications; first, the JMG plane is only valid for galaxies that have a significant fraction of their mass in H i and second, the physical connection between baryonic AM and gas fraction is not universal for all galaxies. In summary, the xGASS JMG plane is only applicable to galaxies with a gas fraction greater than ∼0.01.

This limit is close to the limit applied to create our mock detection sample. In the middle panel of Fig. 5, we show that once galaxies with low-gas fractions are removed from the simulation samples, then the long positive skewed tail is removed. The offsets now have a distribution that is much closer to a normal distribution but it is not centred on zero. This can be seen by the values printed in the top left corner, which show the mean (μ) and standard deviation (σ) of a Gaussian fit to these offsets (in blue and magenta for eagle and tng, respectively). This shows that both simulations on average lie ∼0.1 dex below the JMG plane. The standard deviation of these offsets also implies that eagle has a larger spread around the JMG plane (σ = 0.23) than tng (σ = 0.18).

In the right column of Fig. 5, instead of showing both simulations’ offsets from the best-fitting JMG plane to xGASS, we now compare the offset of the simulations from the best-fitting JMG plane to their own data. Although this does not significantly affect the standard deviation of these offsets, they are now centred on zero.

We also compare the offsets of the xGASS galaxies from the xGASS JMG plane in the middle and right panels of Fig. 5 with the black hollow histograms. The xGASS observational vertical scatter is σ = 0.15. This is similar to the spread found for tng (σ = 0.17) but is smaller than the eagle spread (σ = 0.23). When we calculated the xGASS JMG planes’ σ, we did not attempt to calculate an intrinsic scatter (which takes into account the observational errors of the values), meaning this should instead reflect a combination of measurement errors as well as intrinsic variations of galaxies from the JMG plane. This was due to us not being confident in the exact values of our errors, as it was difficult to combine observational errors, methodology uncertainty and errors introduced from assumptions for individual galaxies (see Papers 1 and 2 for more details). For the simulations, it is also hard to determine the exact errors on each galaxy’s AM, as there is error associated with splitting the gas particles into phases (as we will elaborate on in Section 5), the uncertainty in the assumptions used to determine AM and particle shot noise. Therefore, it is unclear if the difference in scatter seen between the observations and simulations is statistically significant.

4.4 JMG plane offsets and their relationship to star formation rates

The next science question that simulations and their increased statistics allow us to explore is; ‘Is there any residual dependence on SFR that the JMG plane does not encapsulate?’ To address that question we look at the distribution of galaxies with respect to the JMG plane and their ΔMS (offset from the star-forming main sequence). We attempted to do this analysis with xGASS in Paper 2 but were limited in our statistics, so could not draw any strong conclusions. Given that these simulations have a factor of 20 higher statistics, this should not be an issue with these simulation data.

In Fig. 6, we investigate whether variations above and below the star-forming main-sequence influence a galaxy’s position relative to the JMG plane. In this figure, we show eagle (top row) and tng (bottom row) offsets for the mock detection sample (same as those in the right panel of Fig. 5). This shows offsets from either the eagle or tng JMG planes (for the top and bottom rows, respectively). Overlaid on these background histograms are subsamples of galaxies that are either 0.5 dex above (blue) or below (red) the star-forming main sequence. The main sequence is re-defined for each simulation’s mock detection sample using the curved main sequence from Leslie et al. (2020).² In tng, the mock detection sample, above the MS subsample and below the MS subsample have a similar distribution around the JMG plane. For eagle, galaxies below the MS have a slightly larger spread around the plane than galaxies above but have a similar mean. These result qualitatively does not change if we consider offsets from the JMG plane in the f_atm direction. This implies that variations in SFR around the main sequence do not seem to be mirrored by variations around the JMG plane, and vice versa, and that structure and SFR are not directly influencing each other.

5 DISCUSSION

Our analysis has shown that in both eagle and tng, galaxies with gas fractions greater than 0.01 lie on a M_bar – j_bar–f_atm plane that is remarkably similar to the empirical one found in Paper 2. Therefore, the tension we found in Paper 2 between the Obreschkow et al. (2016) gravitational stability model (hereafter O16 model) is also true for eagle and tng. The O16 model has the same qualitative trends as our simulation data but the exact exponents of the O16 model differ from what we find for eagle and tng. As we already discussed in Paper 2, the O16 model makes several simplifying assumptions for the sake of an analytical argument. These assumptions result in the model predicting a single H i profile shape given a fixed q ≔ j_barσ/GM_bar. By using the increased statistics of the simulation data, and its ability to accurately resolve cold gas profiles, we find that the H i profiles have additional baryonic mass dependence which is not encapsulated in this model. A similar result was found in Stevens et al. (2018) with halo mass. An example of tng H i profiles within a small q-interval is shown in Fig. 7. This shows that different baryonic masses exhibit differently shaped H i profiles.The discrepancy between O16 and the cosmological simulations might be due to this non-trivial scale dependence.

The agreement between xGASS, tng, and eagle suggested that we can use the simulations to gain a deeper understanding of the properties and shape of the JMG plane and its implications for galaxy evolution.

As we have shown in Fig. 4, galaxies in the M_bar–j_bar–f_atm parameter space are well described by a plane (in log space) in both simulations for gas fractions above 0.01. However, below this threshold there is an indication that galaxies may deviate from the JMG plane, as illustrated by the long tail in the distribution shown in the left panel of Fig. 4. Galaxies with very low gas fractions have higher baryonic specific AM than predicted by the JMG plane. Although, it should be noted that these galaxies have gas fractions that are at the limit of what would be considered reliable (see results section). Since the simulations do not directly model gas phases, the H i mass of a galaxy is determined in post-processing, with each gas particle assigned a percentage for neutral and then atomic gas. For galaxies to have such low gas fractions, they have either a small number of gas particles contributing to the atomic mass and/or each gas particle contributes a very small percentage of its mass-to-atomic gas. In the very gas–poor regime, potential errors in the gas phase splitting and particle shot-noise add significant statistical and systematic uncertainties to the gas fraction. To address this, in Fig. 4 we also show the blue hollow histogram for a subsample of galaxies in eagle, where uncertain H i masses are removed. In this figure, we define this subsample as galaxies that have at least 10 gas particles, each with at least 5 per cent of their mass in neutral gas. This approach provides a more conservative representation of the distribution of eagle galaxies around the plane, now showing only a faint indication of a tail with large offsets.³ Despite the majority of galaxies within this long extended tail having uncertain H i masses, there are still 38 galaxies with offsets greater than 1 dex and reliable gas fractions. This provides speculative evidence that gas fraction may not accurately predict a galaxy’s baryonic AM in the gas–poor regime. These reliable large offset galaxies all have low gas fractions (f_atm < 10⁻³), whilst possessing moderate baryonic AM. Therefore, whatever mechanism caused the galaxies to deplete their gaseous reservoirs did not significantly reduce their baryonic AM. We inspected these outlier galaxies, and they all appear to be quiescent disc galaxies with the majority of their velocity fields showing regular rotation. Only two have undergone a merger (i.e. stellar mass ratio greater than 0.1) in the last gigayear.

It is worth noting that although the gas fractions of most gas-poor galaxies are uncertain (and consequently, their exact offsets from the JMG plane are uncertain), the galaxies as a whole are still resolved, making it unlikely for them to become gas–rich if a higher resolution simulation was conducted. In addition, as their j_bar and M_bar are dominated by stars, these values can be considered accurate. We find that gas–poor galaxies maintain a similar j_bar–M_bar distribution, once they fall below a gas fraction of 0.01. Therefore, excluding all of these gas–poor galaxies from our analysis limits the conclusions that we can draw and is somewhat unnecessary. To address this, in Fig. 8, we adopt a slightly different reliability measure; we assume that all galaxies with a total H i mass greater than the mass of one gas particle to be accurate (for eagle|$M_{\rm {gas\ particle}} = 1.81 \times 10^{6} \rm {M}_{\odot }$| and for tng|$M_{\rm {gas\ particle}} = 1.4 \times 10^{6} \rm {M}_{\odot }$|⁠). Galaxies with H i masses below this threshold are assigned an upper limit of M_H i = M_{gas particle}. This alternative criterion obtains a similar result to setting a limit based on the number of particles, despite not explicitly checking for an adequate number of gas particles, and allows us to easily apply the same condition to both eagle and tng. Fig. 8 shows galaxies with f_atm > 0.01, with the left panel showing all galaxies with M_HI > M_{gas particle} and the right panel the galaxies with H i masses set to the upper limit.This figure shows that the distribution of f_atm < 0.01 galaxies is preferentially above the JMG plane. M_HI > M_{gas particle} galaxies are, on average, 0.12 and 0.24 dex above the JMG plane, while galaxies set to the upper limits are 0.53 and 0.62 dex above the plane (for EAGLE and TNG, respectively). In other words, galaxies with low gas fractions possess higher AM at a fixed mass and gas fraction than expected for that gas fraction. We expect that the true distribution of gas–poor galaxies around the plane would be somewhere in between the distributions seen in Figs 5 and 8. Once upcoming cosmological simulations that model gas phases directly become available, we will be able to determine the exact relationship between AM and gas content in the gas–poor regime.

This result is noteworthy, as we see that a considerable number of galaxies have moderate angular momentum values despite having little-to-no gas. If all gas-poor galaxies were slow rotators, the JMG plane would still hold in the gas-poor regime, which is not the case. Therefore, despite gas being a strong indicator of a galaxy’s angular momentum in the gas-normal regime, it breaks down in the gas-poor regime. A similar result was found in Paper 1 where we showed that the scatter of the stellar Fall relation was strongly correlated with H i gas fraction at low masses (M_⋆ < 10^10.25 M_⊙), even more so than bulge-to-total ratio. However, when we looked at the high-mass regime (M_⋆ > 10^10.25 M_⊙) then H i gas fraction became less dominant, with bulge-to-total ratio having a slightly stronger correlation with scatter than gas fraction.

It is not surprising that both eagle and tng indicate the presence of a substantial population of passive galaxies with minimal gas but significant angular momentum. Seminal studies of the Virgo cluster already showed that a large fraction of passive galaxies are structurally more similar to discs than ellipticals (Binggeli, Sandage & Tammann 1988; Lisker, Grebel & Binggeli 2006). More recently, the advent of integral field spectroscopic surveys has firmly established that the vast majority of passive galaxies show stellar angular momenta not too dissimilar from those observed in star-forming galaxies (e.g. Cortese et al. 2019; Wang et al. 2020; Cortese et al. 2022; Fraser-McKelvie & Cortese 2022). This clearly highlights how star formation quenching and major structural transformation are two separate (and not always associated) processes in the evolution of galaxies.

The fact that the JMG plane is not universal and valid only for a sub-sample of the local galaxy population does not reduce its importance for galaxy evolution studies. As both eagle and tng implement their star formation and feedback processes differently, the fact that both simulations agree so well with each other could imply that the JMG plane is primarily set by gravitational processes and how cold gas settles in galaxies and reaches equilibrium, rather than any major connection with the way H i is used for star formation.

The idea that the shape of the JMG plane is disconnected from the gas–star formation cycle in galaxies (i.e. both quenching and star-forming stages) is further supported by the lack of any correlation between a galaxy’s offset from the j_bar–M_bar–f_atm plane and its position with respect to the star-forming main sequence. In other words, galaxies that have higher SFR with respect to the main sequence, are not preferentially above or below the JMG plane and vice versa. The trends shown in Fig. 6 are qualitatively the same if offsets are calculated in the f_atm direction. This implies that the physical processes causing an increase (or decrease) in SFR are not driven by processes that cause an increase (or decrease) in atomic gas fraction with respect to the j_bar–M_bar–f_atm plane. This is interesting given that ΔMS is strongly correlated with offsets from the M_⋆–f_atm relation (e.g. Saintonge & Catinella 2022) and we see a similar correlation for the offsets from the M_bar–f_atm relation for both eagle and tng. We can speculate on potential scenarios that could cause the excess cold gas (with respect to the JMG plane) to not be correlated with SFR, such as a large ring of stable H i in the outskirts of a galaxy, which could result in an increase in H i gas fraction, without triggering a starburst event. However, more work is needed to determine the process (or processes) driving the scatter of this JMG plane, and why it is disconnected from galaxies star formation rates.

6 CONCLUSIONS

This study presents a comprehensive comparison between the j_bar–M_bar–f_atm plane for xGASS data presented in Paper 2, and cosmological simulation data from eagle and tng. We compared all the galaxies in each simulation volume, and mock detection samples, to determine how sample selection could be affecting our results. We summarize our main conclusions as follows:

• The j_bar–M_bar–f_atm plane found for the xGASS sample is consistent in both orientation and scatter with the eagle and tng mock detection samples and full simulation samples, for f_atm > 0.01.

• There is moderate evidence that for gas fractions below f_atm ∼ 0.01, the simulations deviate from the empirical JMG plane, asymptoting towards a constant j_bar–M_bar distribution that no longer depends on f_atm.

• The scatter in this JMG plane is independent of ΔMS (for f_atm > 0.01), suggesting that the processes causing deviations from the star-forming main sequence do not affect the processes causing deviations from the JMG plane.

It would be interesting for future works to investigate tracking these simulated galaxies through different redshift snapshots to see if this gives insights into the drivers of scatter in this JMG plane and the factors contributing to its deviation at low gas fractions.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their comments that improved the quality of this manuscript. JAH and LC acknowledge support from the Australian Research Council (FT180100066). Parts of this research were conducted by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. DO is a recipient of an Australian Research Council Future Fellowship (FT190100083), funded by the Australian Government. ARHS is funded through the Jim Buckee Fellowship at UWA. 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, Bruyeres-le-Chatel. The IllustrisTNG simulations were undertaken with compute time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

DATA AVAILABILITY

All of the xGASS data used in this work are publicly available at www.xgass.icrar.org. The eagle simulations are publicly available [see McAlpine et al. (2015) and The EAGLE team (2017) for how to access eagle data]. For access to the tng data used here, please contact ARHS. Otherwise, the public-facing tng data base has similar – although, not identically calculated – galaxy properties available at www.tng-project.org/data/

Footnotes

References

APPENDIX: DETAILED INVESTIGATION OF THE SIMULATION PLANES

For conciseness, in Fig. 4 we present only the comparison between the xGASS JMG plane and simulations for the binned medians of the mock detection sample. For completeness, in this section, we also show the full eagle and tng samples, as well as showing the 2D histogram distribution of the galaxies. This is shown in Figs A1 and A2 for eagle and tng, respectively.

We chose to only show the mock detection sample in Fig. 4 because, in the gas fraction range that is greater than 0.03, the binned medians of the full sample and the mock detection sample are almost identical, which can be seen when comparing the top and bottom rows of Figs A1 and A2. The only differences are seen for 0.01 ≤ f_atm ≤ 0.03, where there is a smaller range in M_bar for the mock detection sample. However, the qualitative agreement with the JMG plane is similar for both the mock and full samples.

Figs A1 and A2 also highlight that the eagle simulation has a much larger spread in j_bar values (at fixed baryonic mass and atomic gas fraction) than tng. This is most evident when comparing the full samples of both simulations, but is also seen when comparing the mock detection samples. This large spread in the data for the full simulation samples, and in particular, the non-Gaussianity of the scatter, which, when combined with the asymptotic behaviour of galaxies approaching a fixed j_bar at extremely low gas fractions, meant that we could not fit a JMG plane directly to the full simulation data. hyper-fit (Robotham & Obreschkow 2015) is a Bayesian fitting tool that assumes data to be normally distributed around the JMG plane. Although the code will give a mathematically correct solution even when the data is not normally distributed, this solution is not physically meaningful. Therefore, we do not show the fits to the full simulation data in this work. As the mock detection sample is closer to a Gaussian distribution, we fit a JMG plane to these samples. The parameters of this fit are given in Table 1 and shown in the bottom rows of Figs A1 and A2 as blue and magenta hashed regions (for eagle and tng, respectively).