-
PDF
- Split View
-
Views
-
Cite
Cite
Florencia Collacchioni, Claudia D P Lagos, Peter D Mitchell, Joop Schaye, Emily Wisnioski, Sofía A Cora, Camila A Correa, The effect of gas accretion on the radial gas metallicity profile of simulated galaxies, Monthly Notices of the Royal Astronomical Society, Volume 495, Issue 3, July 2020, Pages 2827–2843, https://doi.org/10.1093/mnras/staa1334
- Share Icon Share
ABSTRACT
We study the effect of the gas accretion rate (|$\dot{M}_{\rm accr}$|) on the radial gas metallicity profile (RMP) of galaxies using the eagle cosmological hydrodynamic simulations, focusing on central galaxies of stellar mass M⋆ ≳ 109 M⊙ at z ≤ 1. We find clear relations between |$\dot{M}_{\rm accr}$| and the slope of the RMP (measured within an effective radius), where higher |$\dot{M}_{\rm accr}$| are associated with more negative slopes. The slope of the RMPs depends more strongly on |$\dot{M}_{\rm accr}$| than on stellar mass, star formation rate (SFR), or gas fraction, suggesting |$\dot{M}_{\rm accr}$| to be a more fundamental driver of the RMP slope of galaxies. We find that eliminating the dependence on stellar mass is essential for pinning down the properties that shape the slope of the RMP. Although |$\dot{M}_{\rm accr}$| is the main property modulating the slope of the RMP, we find that it causes other correlations that are more easily testable observationally: At fixed stellar mass, galaxies with more negative RMP slopes tend to have higher gas fractions and SFRs, while galaxies with lower gas fractions and SFRs tend to have flatter metallicity profiles within an effective radius.
1 INTRODUCTION
The mass–metallicity relation (MZR) is the correlation between the gas-phase metallicity of galaxies and their stellar mass and has been established both in observations (e.g. Tremonti et al. 2004; Maiolino et al. 2008; Troncoso et al. 2014) and galaxy formation simulations (e.g. De Rossi et al. 2015, 2017; Ma et al. 2016; Collacchioni et al. 2018; Tissera et al. 2019) throughout a wide range in look-back time, 0 to ≈12 Gyr. The chemical enrichment of galaxies is expected to be the result of the interplay between inflows of pristine gas, outflows of enriched material and the star formation history of galaxies (Finlator & Davé 2008; Zahid et al. 2014; Bothwell et al. 2016), and hence the study of the MZR holds important information about the formation of galaxies.
Over the last decade, integral field spectroscopy (IFS) data from surveys such as CALIFA (Sánchez et al. 2012), KMOS (Stott et al. 2014), SAMI (Bryant et al. 2015; Poetrodjojo et al. 2018), and MUSE (Carton et al. 2018; Erroz-Ferrer et al. 2019) have greatly expanded the study of the MZR by providing spatially resolved information of the abundance of metals in galaxies, opening a new window on to the importance of chemical enrichment in the evolution of galaxies.
The study of radial metallicity profiles (RMPs), i.e. the abundance of metals locked in the gas component of galaxies, typically measured through nebular emission lines (see Maiolino & Mannucci 2019 for a recent review), as a function of the distance to the centre of their galaxy, can help to understand how chemical enrichment takes place in galaxies. This is usually quantified as |${\rm log}_{10}(Z_{\rm gas}/\mathrm{\it Z}_{\odot })={\alpha }\, r\, +\, b$|, with Zgas/Z⊙ being the gas metallicity in units of solar metallicity, and α being the slope of the RMP. Pioneering studies found that low-redshift galaxies (z ≤ 0.5) display a negative RMP gradient (α < 0), i.e. galaxies are more chemically enriched in the central regions than in the outskirts (Sánchez et al. 2012, 2013; Ho et al. 2015). This is generally interpreted through the inside-out formation scenario, in which stars at the centre of galaxies form earlier than those in the outskirts and hence have had more time to enrich the interstellar medium (ISM), naturally producing higher gas metallicities in the centre compared to the outskirts (Boissier & Prantzos 1999). However, further studies have found a variety of RMPs that complicate the picture. For example, Troncoso et al. (2014), based on the AMAZE project, found that 3 ≲ z ≲ 4 galaxies display positive RMPs (α > 0), which the authors interpret as low-metallicity gas flowing directly into the centres of galaxies (see also Cresci et al. 2010).
Instead of measuring a single α for the whole RMP, Sánchez et al. (2014) and Sánchez-Menguiano et al. (2016), using the CALIFA survey, and Belfiore et al. (2017), using the SSDS-IV MaNGA survey (among others authors), fit the RMP with two power laws, i.e. two different α for the inner and outer regions, finding puzzling results that cannot be easily interpreted within the inside-out scenario, such as a flattening in the outer regions, α ≈ 0, or relatively flat inner profiles (at least, less negative than expected) in high-mass galaxies (|$M_\star \gtrsim 10^{10} \, {\rm M_\odot }$|; Sánchez-Menguiano et al. 2016). This could be linked to feedback removing metals preferentially from the inner regions of the galaxy (e.g. Lagos, Lacey & Baugh 2013; Muratov et al. 2017), as well as significant star formation in the outer regions of the galaxy (e.g. due to mergers or close interactions) that can quickly enrich the gas.
Several studies have explored the possible physical causes behind the shape of the RMP. According to Sánchez et al. (2012), the change in the outer regions (beyond approximately two times the effective radius) could be due to different gas densities, the star formation history, or even the presence of bars, which can alter the flow of gas internal to galaxies. Yet, Zinchenko et al. (2019), also using the CALIFA survey, found no correlation between the RMP gradient and bars or spiral patterns, and Carton et al. (2018) did not find a correlation between α and the star formation rate (SFR) or stellar masses of galaxies when analysing MUSE data. This is in tension with the reports of weak correlations between α and the specific SFR (|$\rm sSFR ={\rm SFR}/{\it M}_{\star }$|) by Stott et al. (2014) and Wuyts et al. (2016) using KMOS surveys, and with stellar mass by Ho et al. (2015), who use a combined CALIFA and SDSS Data Release 7 galaxy sample. Moreover, Ho et al. (2015) found that measuring the RMP normalized by the galaxy’s effective radius does lead to a correlation between α and galaxy properties, implying a co-evolution of gas and stars.
It is fair to say that there is tension between the different observational results, which may, in part, be due to the different selection of galaxies, but also perhaps suggesting that the RMP is not well described by a single or double power law, and/or that other galaxy properties may do correlate more strongly with the RMP than sSFR, stellar mass, and morphology. In fact, Carton et al. (2015), measuring atomic hydrogen gas content (H i) with the WSRT of 50 SDSS galaxies, found that the RMP slopes show a strong correlation with the H i mass fraction (the H i-to-stellar mass ratio), such that galaxies with higher H i mass fractions also have more negative α.
From a theoretical perspective, a range of results from hydrodynamical simulations of galaxy formation have been reported that do not necessarily agree with each other. Tissera et al. (2016), using a simulated cubic volume of 14 Mpc on a side, argued that a range of slopes, from very negative to positive, can result from close galaxy encounters, presence of bars and/or low-metallicity gas accretion in the central regions of galaxies. This variety of outcomes is caused by these processes triggering star formation activity and SN-driven outflows in different ways. A similar result is reported by Sillero et al. (2017), who concluded that galaxy interactions generally lead to α > 0 for the case of major mergers of disc galaxies, but significant variations in the metallicity distribution and slopes are seen throughout the merger process. In addition, Sillero et al. (2017) found that the initial gas fraction of the galaxy merger and the SN feedback strength had an effect on α. An important limitation of the studies above is the poor statistics, which makes it difficult to identify the primary processes that lead to higher/lower α values. The latest generation of galaxy formation simulations has allowed a more thorough exploration of these trends with better statistics, thanks to the much larger cosmological boxes available. One such example is the eagle simulation (Crain et al. 2015; Schaye et al. 2015); its largest box is 100 Mpc on a side, ≈360 times larger volume than the studies above. Using eagle, Tissera et al. (2019) showed that very active merger histories are associated with galaxies displaying flatter RMPs, which naturally leads to lower mass galaxies having more negative RMPs. Tissera et al. (2016) and Sillero et al. (2017) agree in that galaxies with stellar masses M⋆ ≤ 1010 M⊙ show more negative RMPs than more massive ones. Other simulations, however, find contradictory results. Ma et al. (2017), using cosmological zoom-in simulations from the FIRE project (Hopkins et al. 2014), reported low-mass galaxies to have flatter RMPs compared to massive galaxies. Tissera et al. (2019) also found no correlation between α and the environment in which galaxies live, while the physical size of the gaseous disc of galaxies appeared as one of the only galaxy properties clearly correlated with α, in a way that more extended gas discs are associated with more negative α. Interestingly, Tissera et al. (2016) find that as redshift increases, high-mass galaxies (M⋆ ≥ 1010 M⊙) show less negative and sometimes even positive α.
Gas accretion is expected to shape both the stra formation and metallicity histories of galaxies and, hence, one would like to explore directly the correlation between the RMP and the gas accretion rate (|$\dot{M}_{\rm accr}$|). However, |$\dot{M}_{\rm accr}$| is not easily accessible in observations, motivating the simulation-based work to explore indirect tracers, such as star formation and stellar mass. Star formation is expected to be proportional to the gas supply, such that more infalling gas can trigger more star formation (e.g. Davé, Finlator & Oppenheimer 2011), while more massive galaxies live in more massive haloes (Behroozi, Wechsler & Wu 2013) and hence should be able to attract more gas into their gravitational potential. By studying directly the effect of gas supply on to galaxies, Perez, Michel-Dansac & Tissera (2011) found that more negative slopes are associated with low-metallicity gas accretion, while Sillero et al. (2017) predicted a correlation between α and sSFR only in cases where the gas inflow triggers star formation within the central regions on short time-scales. The proposed scenario is that significant |$\dot{M}_{\rm accr}$| on to the central regions of galaxies can increase the gas density, triggering high levels of star formation which rapidly enrich the inner regions of galaxies, leading to more negative slopes of the RMP. Other studies suggest extreme accretion, i.e. mergers, yields inverted metallicity gradients (Troncoso et al. 2014). Hence, the study of RMPs might be essential to constrain the properties of the gas that is being accreted on to the galaxies (i.e. pristine or pre-enriched), and the effect |$\dot{M}_{\rm accr}$| has on their chemical evolution (Finlator 2017).
The main limitation of the above studies is that they tend to focus on a small number of simulated galaxies, and hence it is difficult to separate effects related to gas accretion from other physical effects, such as galaxy mergers/interactions, depth of the gravitational potential well, and gas outflows. It is therefore essential to explore theoretically how |$\dot{M}_{\rm accr}$| affects the RMP of galaxies across cosmic time in a large sample of simulated galaxies over cosmologically representative volumes. This is the context in which our work develops. We use the eagle cosmological hydrodynamical simulations suite to study how |$\dot{M}_{\rm accr}$| correlates with the RMPs, breaking the RMP into different regions of the galaxies. Our objective is to understand which regions of galaxies are most affected by the gas supply, what the roles of stellar mass, SFR and gas content are in determining the RMP, and what the ‘smoking guns’ of gas accretion are, which can be seen in the RMP. We do this by combining simulated boxes at different resolutions to allow us to cover the stellar mass range |$M_{\star } \gt 10^9 \, \rm M_{\odot }$| at 0 ≤ z ≤ 1.
This paper is organized as follows. In Section 2, we briefly summarize the cosmological hydrodynamic simulations used in this work, as well as the galaxy sample selected, and introduce our definition of |$\dot{M}_{\rm accr}$|. We analyse how the RMP changes with |$\dot{M}_{\rm accr}$|, stellar mass, SFR, and gas content in Section 3, and show that |$\dot{M}_{\rm accr}$| is the property that is most strongly correlated with the slope of the RMP. In Section 4, we discuss the scope of our results and limitations, as well as the possibility of measuring the correlations found here in observations. Finally, our conclusions are presented in Section 5.
2 eagle SIMULATIONS
The eagle project1 is a set of cosmological hydrodynamical simulations that differ in box size, number of particles, numerical resolution, and subgrid physics. A brief description of the simulations and the underlying subgrid models is presented in this section. For more information on the model, the reader can refer to Schaye et al. (2015, hereafter S15) and Crain et al. (2015), while the public data release is documented in McAlpine et al. (2016). To run the eagle simulations, a modification of the N-body tree-PM smoothed particle hydrodynamics (SPH) code gadget-3 (last described by Springel 2005) is used, with the numerical methods referred to as anarchy (Schaller et al. 2015).
eagle follows the evolution of dark matter (DM) and gas particles, consistent with a flat Lambda cold dark matter cosmology characterized by the Planck Collaboration I (2014) parameters: matter density Ωm = 0.307, dark energy density |$\Omega _\Lambda = 0.693$|, baryon matter density Ωb = 0.04825, square root of linear variance of matter distribution σ8 = 0.8288, index of power spectrum of primordial adiabatic perturbations ns = 0.9611, Hubble parameter |$H_{\rm 0} = 100 \, h^{-1} \, {\rm km} \, {\rm s}^{-1} \, {\rm Mpc}^{-1}$| with h = 0.6777, and primordial helium abundance Y = 0.248.
A Friend-of-Friends (fof; Davis et al. 1985) algorithm is used to identify haloes of DM particles; while the subfind algorithm (Springel et al. 2001; Dolag et al. 2009) is used to identify the gravitationally bound subhaloes within a fof structure, linking the gas particles with their nearest DM particles. Theses subhaloes are then defined as the galaxies in the simulation. For each fof halo, the subhalo that contains the particle with the lowest gravitational potential value is defined as a central, and satellite galaxies are all of the remaining subhaloes.
eagle contains a range of subgrid models that are described below, and those have free parameters that are calibrated to the observed z ≈ 0 galaxy stellar mass function (GSMF), the relation between the galaxy stellar and black hole (BH) mass, and galaxy sizes of star-forming (SF) galaxies. For this work, we use two of the eagle simulations: the recalibrated high-resolution simulation of a 253 cMpc3 volume, which we will refer to as Recal-L25N752, and the reference simulation of a 1003 cMpc3 volume, called Ref-L100N1504. DM particle masses are 1.21 × 106 M⊙ for the former and 9.7 × 106 M⊙ for the latter, while the gas particle masses are 2.26 × 105 and 1.81 × 106 M⊙, respectively.
The subgrid physical models implemented in the eagle simulations are the following. Radiative cooling and photoheating on 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) are applied (Wiersma, Schaye & Smith 2009a), accounting not only for variations in metallicity but also in relative abundances. Hydrogen reionization is implemented as described in S15 at z = 11.55. Star formation follows the description in Schaye & Dalla Vecchia ( 2008), with a metallicity-dependent density threshold, in order to have a more realistic description of the gas transitions between warm, neutral, and cold phases (Schaye 2004). A Chabrier stellar mass function (IMF, Chabrier 2003) is adapted, with minimum and maximum masses of 0.1 and 100 M⊙, respectively. Stellar mass loss (Wiersma et al. 2009b), and mass and energy losses from Type II and Ia supernovae (SNe) are also implemented. Energy from SN is stochastically injected in nearby gas particles in the form of thermal feedback as in Dalla Vecchia & Schaye (2012). In a similar way, feedback from active galactic nuclei (AGN) is computed with a fixed efficiency, and the energy is injected into the surrounding gas medium in a thermally and stochastic way (for details, see S15).
eagle’s subgrid model imposes a gas temperature threshold that scales with the density of the gas (S15) to prevent spurious fragmentation due to the finite resolution of the simulation. Therefore, eagle does not model the cold gas phase, which makes the ISM of galaxies more homogeneous (or less clumpy) than we observe in real galaxies, leading to stars and dust, in principle, being more evenly mixed than it would be if a cold gas phase was included. Despite this limitation, as we increase the resolution of the simulation, though maintaining the temperature threshold above, the ISM displays increased levels of clumpiness that is not numerically driven and instead a physical effect (Trayford et al. 2017, 2020). In eagle, stars form in these clumps that have extensions of ≲1 kpc, which is larger than the typical molecular clouds in the local Universe. As a result, we expect metals to more easily mixed in the ISM compared to what we would expect if more clumpiness was produced. A consequence of this could be less metallicity variations across the galaxy, but given we are interested in metallicity profiles measured in radial bins of ≈1 kpc, this is likely not a big limitation. Current IFS observations used to constrain metallicity gradients have a similar spatial resolution to what we can achieve with the combination of eagle simulations of Recal-L25N752 for galaxies with M⋆ < 1010 M⊙ and Ref-L100N1504 for more massive galaxies (500–1000 pc for their median redshift; e.g. van de Sande et al. 2019).
eagle has proven to be very successful in a variety of comparisons. Several papers have shown that the reference eagle simulation reproduces observed galaxy properties such as the stellar mass function evolution from z = 0 to 4, the main sequence of galaxies (Furlong et al. 2015), the colour distribution of galaxies (Trayford et al. 2015; Wright et al. 2019), and the gas content of galaxies of a given stellar mass (Lagos et al. 2015, 2016; Bahé et al. 2016; Crain et al. 2017). Very relevant for this work, the mass–metallicity relation of stars and gas, the resolved mass–metallicity relation (Trayford & Schaye 2019), as well as the gas metallicity gradients (S15, Tissera et al. 2019), are well reproduced at M⋆ ≳ 1010 M⊙ for Ref-L100N1504. At lower stellar masses, deviations from the observations become important. The higher resolution, recalibrated run of the eagle suite alleviates this problem, enabling agreement down to stellar masses of 109 M⊙ (see S15; De Rossi et al. 2017). Hence, the simultaneous use of both the Ref-L100N1504 and the Recal-L25N752 simulations allows us to have very good statistics above 1010 M⊙ in stellar mass (mostly from the Ref-L100N1504 simulation), and to push down to 109 M⊙, thanks to the Recal-L25N752 simulation, with confidence that the observed mass–metallicity relation is well reproduced over the entire stellar mass range. However, we keep the analysis of the two simulated samples separated, unless otherwise stated. We remind the reader that the Recal-L25N752 simulation has slightly different parameters for the stellar and AGN feedback, as the aim of that run was to reproduce the z = 0.1 stellar mass function with a higher resolution (see discussion on weak and strong convergence in S15).
The stronger feedback of the Recal-L25N752 simulation tends to remove more metals from the ISM than what we see in Ref-L100N1504, which is partially why the metallicity of galaxies in the former is lower than in the latter at stellar masses <1010 M⊙. S15 discussed this in some length and showed that a high-resolution version of the Ref-L100N1504 (with exactly the same subgrid physics model and parameters as the Ref-L100N1504 run but with the volume and number of particles of the Recal-L25N752) produces gas metallicities that are ∼0.2 dex higher than Recal-L25N752 at stellar masses of ≈108 M⊙. Although the zero point is slightly different between these two simulations, the slopes are overall in agreement once we studied them relative to their stellar mass (as we will see in the following sections).
2.1 Computation of the gas accretion rate
We compute |$\dot{M}_{\rm accr}$| on to galaxies using a simple particle tracking methodology (following Neistein et al. 2012). For a given subhalo at a given simulation snapshot, we first identify gas particles that are classified as SF, as described in S15. We then define accreting gas particles as the subset of these that are
bound to the main progenitor subhalo at the previous (earlier) snapshot in the merger tree, and
not SF at that time.
Star particles that fulfil conditions (a) and (b) are also taken into account in the calculation of the accreting gas. These stellar particles are those that come from gas particles that were accreted as of the previous snapshot. With this criteria, we ensure the accreting gas came from the same subhalo studied and, at the same time, that the accretion itself changes the state of the gas (from non-SF to SF). In other words, we are interested in inflows that can trigger star formation. The latter is done to more easily connect with the observations, in which nebular emission lines are used as gas tracers and to measure the gas metallicity.
The choice of only counting the gas coming from the main progenitor is to compute the contribution of the accreted gas that comes from the diffuse circumgalactic medium (CGM) rather than galaxy mergers. We remark that previous hydrodynamical simulations show this ‘smooth accretion’ to be the main source of gas accretion on to galaxies (van de Voort et al. 2011). Wright et al. (in preparation) show that in eagle 80 per cent or more of the gas accreted on to haloes comes from smooth accretion (which includes both pristine accretion and pre-processed gas). |$\dot{M}_{\rm accr}$| is then simply defined as the mass of accreted gas particles that satisfy conditions (a) and (b), divided by the time interval between the snapshots. Progenitor subhaloes are identified using merger trees generated with the DHalo algorithm (Jiang et al. 2014; Qu et al. 2017).
For the Ref-L100N1504 simulation, we use 50 snapshots to compute |$\dot{M}_{\rm accr}$|, while for Recal-L25N752, we use 29 snapshots. The reason why we use different numbers of snapshots is that for Ref-L100N1504, we have a larger number of outputs available with a better time cadence but with a reduced number of gas particle properties, also called ‘snipshots’, compared to the default 29 snapshots that are publicly available from the eagle data base (McAlpine et al. 2017). For Recal-L25N752, however, we use the public data, which include 29 snapshots, with the complete list of particle properties. The time interval between outputs for Ref-L100N1504 is 0.42 Gyr from z = 0 → 0.03, and 0.4 Gyr from z = 1 → 1.05, while for Recal-L25N752, the time interval between outputs is 1.34 Gyr from z = 0 → 0.1, and 0.93 Gyr from z = 1 → 1.26. Although we are using different time-steps to compute accretion rates in the two simulations, this does not have an important impact on our results as accretion happens in long time-scales (see Mitchell et al., in preparation; Wright et al., in preparation). If we were studying outflows, however, this difference in time cadence would be important, as outflows happen in short time-scales and are most stochastic. Hereafter, we use the word ‘snapshots’ to refer to both these sets of outputs, making it clear which simulations we are referring to.
There is, however, an inherent weakness of the particle tracking method. If there is inflow that joins the SF ISM but leaves it before the snapshot in which the accreted gas is identified (with ‘leave’ we mean that the particles are transferred to another phase, e.g. non-SF ISM), then that inflow will not be counted. For more details about the method implemented here to measure the accreted gas, the reader can refer to Mitchell et al. (in preparation).
2.2 Simulated galaxy sample
Because of how sensitive the infalling gas is to environmental effects such as tidal stripping, ram pressure stripping, interactions, etc., which predominantly affect satellite galaxies, we limit this study to central SF galaxies only (i.e. the galaxy sitting in the deepest part of the potential well of haloes). We also apply the restriction of using only galaxies with M⋆ ≥ 109 M⊙ for the Recal-L25N752 simulation and M⋆ ≥ 1010 M⊙ for the Ref-L100N1504 one, due to the resolution of the respective simulations (which gives us a minimum star particles number of ≈6700 and ≈7700, respectively). We include those central galaxies whose |$\dot{M}_{\rm accr}$| corresponds to at least 10 SF gas particles, which translates into |$\dot{M}_{\rm accr}$| ≈9 × 10−2 M⊙ yr−1 at z = 0 and |$\dot{M}_{\rm accr}$| ≈5 × 10−2 M⊙ yr−1 at z = 1. In Appendix A, we test our main results against different choices of the latter minimum number of SF particles, and found them to be insensitive to thresholds up to 100 particles. Our results are therefore well converged against this limit.
Our final sample comprises a total of 1280 galaxies at z = 0 and 1642 at z = 1. Because of the size evolution of galaxies (galaxies at higher redshift are smaller at fixed stellar mass; Furlong et al. 2017), we limit our study to z ≤ 1, for which galaxies are better resolved. The median physical sizes of our selected galaxies2 in the Recal-L25N752 simulation are 3.98 and 3.17 kpc at redshifts z = 0 and 1, respectively. Similarly, the median galaxy sizes of our selected sample in the Ref-L100N1504 simulation are 4.22 and 3.11 kpc at redshifts z = 0 and 1, respectively. Recall that the latter galaxies are more massive than the ones selected from the Recal-L25N752 simulation. In both cases, these sizes are comfortably above the gravitational softening of the corresponding simulation, which are 0.35 and 0.70 physical kpc for Recal-L25N752 and Ref-L100N1504, respectively. Ludlow et al. (2019) showed that the spatial resolution of a simulation is approximately 0.05 times the mean particle spacing, which for the Ref-L100N1504 and Recal-L25N752 are 66.5 and 33.35 physical kpc at z = 0, respectively, leading to a spatial resolution of 3.3 and 1.6 kpc, respectively. This shows that the stellar mass selection applied here selects well-resolved galaxies in these two simulations.
2.3 Estimation of the gas metallicity profile
3 THE EFFECT OF |$\dot{M}_{\rm accr}$| ON THE GAS METALLICITY GRADIENTS
In this section, we present correlations between several galaxy properties obtained by the eagle simulations Ref-L100N1504 and Recal-L25N752, i.e. stellar mass, |$\dot{M}_{\rm accr}$|, SFR, and neutral gas fraction, and how they correlate with the RMPs. This is done at z = 0 and 1. The analysis has also been performed for redshifts in between, which produce consistent results. Hence, for the sake of conciseness, we only show z = 0, 1.
3.1 The connection between metallicity, stellar mass, and |$\dot{M}_{\rm accr}$|
Fig. 1 shows the gas metallicity versus stellar mass relation at redshift z = 0 (top panel) and 1 (bottom panel) for a combined sample of galaxies from the Recal-L25N752 and Ref-L100N1504 simulations, respectively, described in Section 2. In this case, the gas metallicity is estimated from all the SF gas particles within 30 kpc from the centre of potential of the galaxy (this is the typical radius in which the integral properties of galaxies are calculated in the eagle simulations; see S15 for details). Median values of the gas metallicity are shown as solid lines, while the dashed lines depict the 16th and 84th percentiles. Coloured pixels show the median |$\dot{M}_{\rm accr}$| of galaxies in bins of 0.25 dex in stellar mass and metallicity. In observations, gas metallicities are usually estimated from regions that can be smaller than 30 kpc. For this reason, we investigated the effect of choosing a smaller aperture (10, 15, and 20 kpc) and find that changes in the MZR are negligible, except at the very high mass end, stellar masses |${\gt}10^{11} \, \rm M_{\odot }$|, where we see an increase in the scatter. These differences, however, are small, and hence we decide to continue the analysis with an aperture of 30 kpc to ease the comparison with other simulation papers.

Stellar mass–gas metallicity relation (MZR) at z = 0 (top panel) and 1 (bottom panel) for galaxies in a combined simulated sample (combining the galaxy samples of the Recal-L25N752 and Ref-L100N1504 simulations; see details in Section 2). The gas metallicity is estimated using all the SF gas particles within 30 kpc from the centre of potential of the galaxy. In each panel, the black solid and thin dashed lines represent the median and 16th–84th percentiles, respectively. Black symbols with error bars in the top panel show the observational measurements of the MZR from Tremonti et al. (2004). At fixed redshift, galaxies with higher stellar mass show higher gas metallicity. At fixed stellar mass, galaxies show a higher gas metallicity at lower redshift. Coloured pixels represent the median |$\dot{M}_{\rm accr}$| of the galaxies in bins of stellar mass and metallicity (containing at least 10 galaxies). Galaxies with higher stellar mass show higher |$\dot{M}_{\rm accr}$|, with |$\dot{M}_{\rm accr}$| values decreasing as redshift decreases. Thick grey dotted and dashed lines show the median of galaxies in two bins of |$\dot{M}_{\rm accr}$| as labelled in the panels (presented in units of |$\rm M_{\odot }\, yr^{-1}$|). At fixed stellar mass and redshift, there is a tendency for the gas metallicity to be anti-correlated with |$\dot{M}_{\rm accr}$|.
As expected, the gas metallicity increases with stellar mass until a threshold mass above which the metallicity saturates or even decreases. This threshold mass increases with redshift from |${\approx}10^{10}\rm \, M_{\odot }$| at z = 0 to |$10^{11}\, \rm M_{\odot }$| at z = 1. This is related to the transition from galaxies that are preferentially SF to preferentially passive. At high redshift, this happens at a higher stellar mass than at z = 0, as seen from the development of the red sequence in eagle (Trayford et al. 2016; Wright et al. 2019). At z = 1, Trayford et al. (2016) showed that the red sequence in eagle is populated at the very high mass end (M⋆|$\gt 10^{11} \, \rm M_{\odot }$|) by central galaxies and at low masses (M⋆|$\lt 10^{9.5} \, \rm M_{\odot }$|) by satellite galaxies, with a dearth of galaxies at around |$10^{10} \, \rm M_{\odot }$|. This dearth region gets increasingly populated towards z = 0. Trayford et al. (2016) and Bower et al. (2017) argued that this is due to AGN feedback being effective only at the highest masses at z = 1, which is where BHs have grown enough to a regime where they can be efficient at driving outflows. The combination of the AGN effect and the lower number of SF galaxies leads to a flattening of the MZR moving from higher to lower stellar masses as time goes on, and can even appear as a decrease of the metallicity at the high mass end. The latter has been reported in observations (Sánchez et al. 2017). This may in part may be due to systematic uncertainties. In these very metal-rich galaxies, it can happen that the electron temperature decreases, complicating the metallicity estimations (Marino et al. 2013). If we add to this scenario the fact that massive galaxies also have the higher |$\dot{M}_{\rm accr}$| at each redshift, then this inflow can be diluting the metallicity of the gas, contributing to the decrease seen in Fig. 1. At the low-mass end, satellites only make a small contribution to the total number of galaxies and hence they do not impact the shape of the MZR.
The normalization of the MZR also increases with decreasing redshift, and its shape (in terms of a slope of the relation) changes with redshift as well. Thus, the MZR evolves in normalization and shape, implying that a galaxy’s stellar mass and chemical abundance evolve at different rates. This has been reported by Guo et al. (2016) for the Ref-L100N1504 simulation and by De Rossi et al. (2017) for Recal-L25N752. This evolution in normalization is a natural success of the eagle simulations and generally not easily obtained in other simulations (see Collacchioni et al. 2018 for a discussion). It is interesting to note that, at fixed redshift, galaxies with higher stellar mass also have a higher |$\dot{M}_{\rm accr}$|, on average. This positive correlation is tight (about ∼0.5 dex of dispersion) at all redshifts analysed in this work, and it is in agreement with previous theoretical works (e.g. van de Voort et al. 2011; Mollá et al. 2016; Correa et al. 2018, see also Appendix B).
From Fig. 1, it can also be seen that, at fixed stellar mass, the gas metallicity increases as |$\dot{M}_{\rm accr}$| decreases, regardless of the redshift considered. Thick grey lines denote the MZR for galaxies in different ranges of |$\dot{M}_{\rm accr}$| (as it is described in the panels). Galaxies with lower values of |$\dot{M}_{\rm accr}$| have a MZR with a higher normalization, i.e. with higher values of metallicity. This trend resembles the observed fundamental mass–metallicity relation plane (Lara-López et al. 2010; Mannucci et al. 2010), which correlates the metallicity, the stellar mass, and the SFR in the sense that, at fixed stellar mass, the metallicity increases as the SFR decreases. The existence of this relation in eagle was demonstrated and discussed by De Rossi et al. (2017) and Matthee & Schaye (2018). In this scenario, the accretion of pristine gas has the dual effect of diluting the gas metallicity and triggering star formation (Troncoso et al. 2014). However, the fact that De Rossi et al. (2017) found a stronger anti-correlation between sSFR and metallicity in simulations with stronger stellar feedback indicates that outflows, which preferentially eject metals, also play an important role.
Galaxies have lower |$\dot{M}_{\rm accr}$| as redshift decreases (at fixed stellar mass), showing an evolution with time. This is expected as the overall gas accretion rate on to haloes, and therefore on to galaxies, decreases with time (Correa et al. 2018). However, it is important to emphasize that large differences are expected between |$\dot{M}_{\rm accr}$| on to haloes and on to galaxies (e.g. van de Voort et al. 2011; van de Voort 2017; Nelson et al. 2013, 2019). The interplay between |$\dot{M}_{\rm accr}$| and SFR will be analysed in Section 3.5.
3.2 The relation between the RMP and |$\dot{M}_{\rm accr}$|
The aim of this work is to get a better understanding of how |$\dot{M}_{\rm accr}$| and other galaxy properties affect the RMPs. As was discussed in Section 1, the slope of the RMP measured over the whole galaxy does not depend on stellar mass or the environment in which the galaxy is immersed, particularly in the eagle simulations (Tissera et al. 2019). However, since we are interested in finding out whether there is a preferential radius below/above which the profiles are significantly affected by gas accretion, we divide them into several regions and analyse the corresponding slopes instead of analysing a single slope of the RMP. We measure the power-law slope of the RMP, α, in regions defined by the radii that contain half of the stellar mass of the galaxy, referred to as r50.
Tissera et al. (2019), also using eagle, found that the slope of the gas metallicity profile correlates strongly with the half-neutral gas mass size radius (i.e. the radius that contains half of the neutral, i.e. atomic plus molecular, gas mass), but that the half-stellar mass radius plays a small role. This suggests that the effective radius of the neutral gas content would be a more appropriate property to study the RMP. However, this property is only rarely available in observations and more difficult to define robustly. For example, observations of H i have shown that integrating longer, which allows to push towards lower H i column densities, continues to reveal H i, making the half-gas mass radius sensitivity-dependent (e.g. Oosterloo et al. 2007; Heald et al. 2011; Kamphuis et al. 2013). The stellar r50 is a well-defined quantity that can be measured robustly for hundreds of thousands of galaxies (see Lange et al. 2016 for an example in the nearby Universe).
Fig. 2 shows the relation between the slope of the RMP in the inner (left-hand panels) and outer (right-hand panels) regions of eagle galaxies as a function of |$\dot{M}_{\rm accr}$| at z ≤ 1. Solid and dashed lines show the median values and 16th–84th percentiles, respectively. Results for both simulations are shown separately using different colours. There is a small offset between the two simulations in this figure, which is not surprising as for Recal-L25N752 we are including all galaxies with stellar masses |$M_{\star } \ge 10^9 \, \rm M_{\odot }$|, while for the Ref-L100N1504 simulation, we limit the sample to galaxies with stellar masses |$M_{\star } \ge 10^{10} \, \rm M_{\odot }$|. We show later that eliminating the dependence on stellar mass brings the two simulations into agreement.

Median values of the RMP slopes, α, for each simulation used in this work as a function of |$\dot{M}_{\rm accr}$|. Top and bottom panels show galaxies at z = 0 and 1, respectively. The left- and right-hand panels depict different regions of the galaxies, as labelled. Galaxies are colour-coded by the simulation they are generated with, being red for Recal-L25N752 and blue for Ref-L100N1504. Dashed lines depict the 16th–84th percentiles. We show individual galaxies where bins have fewer than 10 objects (circles). The RMP slope in the inner region (r/r50 ≤ 1) is negative and it is tightly correlated with |$\dot{M}_{\rm accr}$| for both simulations at all redshifts analysed, becoming more negative for higher |$\dot{M}_{\rm accr}$|. The Recal-L25N752 simulation shows a more negative slope, especially for higher values of |$\dot{M}_{\rm accr}$|. At larger radii (1 < r/r50 ≤ 5), we find a weaker correlation between α and |$\dot{M}_{\rm accr}$| with much larger scatter, with values of α ≈ −0.08. As explained in Section 2.2, the cuts in stellar mass for our galaxy samples are M⋆ ≥ 109 M⊙ for Recal-L25N752 and M⋆ ≥ 1010 M⊙ for Ref-L100N1504. The difference in mass accounts for most of the difference between the simulations (see Section 3.3). Note that the y-axis range changes from the left- to the right-hand panels, which is done to better highlight the values spanned by the data.
For the inner region, it can be seen that galaxies with a higher |$\dot{M}_{\rm accr}$| display more negative RMP slopes α, with 1σ scatter of ≈0.1 |${\rm dex} \, r_{\rm 50}^{-1}$|. The anti-correlation between α and |$\dot{M}_{\rm accr}$| is present in both simulations, although the trend is stronger for Recal-L25N752. This trend is seen at all redshifts studied. A simple physical picture to interpret this trend would be the accreted gas having lower metallicity than the ISM in the galaxy, and hence diluting the gas as it falls on to the disc. This would steepen the RMP because conservation of specific angular momentum will cause the accreted gas to settle in the outskirts of galaxies (see Tissera et al. 2019 for an analysis of this scenario in eagle). Another explanation is that the correlation in Fig. 2 may be caused by these two galaxy properties correlating with a third, more fundamental one. We explore this in upcoming sections.
In the outer regions, 1 < r/r50 < 5, we find that α stays close to ≈−0.08 |${\rm dex} \, r_{\rm 50}^{-1}$|, with a scatter of about ≈0.25 dex. It is interesting to note that the α values at 1 < r/r50 < 5 are significantly less negative than the RMP slope in the inner parts for galaxies that have large |$\dot{M}_{\rm accr}$|, while galaxies with low |$\dot{M}_{\rm accr}$| have a similar slope at r/r50 < 1 and 1 < r/r50 < 5. There is a weak trend of |$\dot{M}_{\rm accr}$| being non-monotonically related to α in these outer regions, with |$0.5 \le \dot{M}_{\rm accr} {\rm [M_\odot \, yr^{-1}]) \le 1.5}$| being associated with more negative values of the outer slope. This effect is however, weak, with differences of ≲0.05 dex in α and only appreciable in the simulation Ref-L100N1504. The large scatter obtained in eagle at 1 < r/r50 < 5 is likely driven by the low numbers of SF gas particles galaxies have at their outskirts, making the measurement of α very noisy. We cannot therefore conclusively say that the relation between α and |$\dot{M}_{\rm accr}$| at 1 < r/r50 < 5 is really non-monotonic, but we can assert that the outer parts of galaxies have only mildly negative RMPs. A possibility leading to the flat RMP at r/r50 > 1 is that the gas metallicity of the ISM reaches that of the halo gas, which may act as a metallicity floor. To confirm this a more detailed analysis of the gas distribution should be made, which is beyond the scope of this work and will be left for future research. In addition, higher resolution simulations would be required to increase the number of SF particles in the outer regions to better constrain the values of α there.
Previous simulation work have found galaxy mergers and close encounters to have a significant effect on their chemical enrichment patterns in a way that the inflowing gas from the interaction dilutes the metallicity of the central regions of the galaxy (Montuori et al. 2010). Observations of local ultraluminous infrared galaxies suggest similar trends (Rupke, Veilleux & Baker 2008). We explore this possibility by analysing the slopes of the galaxies that suffered a major or minor merger since the previous snapshot (i.e. galaxies that had more than one progenitor of M⋆ ≳ 108 M⊙). We find that only ≈10 per cent of our sample has experienced a recent major or minor merger event, and that the slopes α in their inner and outer regions are very similar to the overall galaxy sample. Hence, we do not find any evidence that mergers affect the RMPs of galaxies in significant ways (at least at z ≤ 1). This, in part, may be due to the majority of galaxy mergers at z ≲ 0.5 in eagle having low gas fractions (neutral gas-to-stellar mass ratios ≲0.1; see fig. 2 in Lagos et al. 2018) and hence having a reduced effect on the gas metallicity profiles. These are certainly different in nature from the objects studied in Rupke et al. (2008), which were selected to be highly starbursting galaxies, and different from the simulations of Montuori et al. (2010), which are gas-rich major mergers. Therefore, it is not surprising that we find a weaker effect of galaxy mergers.
Given the strong (anti-)correlation we obtain between α at r < r50 with |$\dot{M}_{\rm accr}$|, we analyse the latter in detail in the upcoming sections. From now on, unless otherwise specified, we only focus on the RMP at the inner region.
3.3 The relation between the RMP and stellar mass
Fig. 3 shows the RMP inner slope, α, as a function of stellar mass. Both simulations display an anti-correlation so that more negative α are associated with more massive galaxies. This may seem at first sight in contradiction to the results of Tissera et al. (2019), who also used eagle, but the main difference is that here we show α measured within r50, while Tissera et al. (2019) measured α between 0.5 × r50 and 2 × r50. As a result, the RMPs in Tissera et al. (2019) are slightly flatter than the ones we calculate for r < r50. To our knowledge, there is no other result in the literature with the RMP slope calculated using the radial cuts that we adopt. As a consequence, we do not compare our results with observations and other theoretical findings due to the lack of consistency. We, however, encourage other works, specially observational ones, to use the radial cuts we propose here.

The median relation of the slope of the inner RMP, α, as a function of the stellar mass, M⋆. Top and bottom panels show galaxies at z = 0 and 1, respectively. The Recal-L25N752 simulation is shown with thin red lines, while the Ref-L100N1504 simulation is shown with thick blue lines. In all cases, dashed lines depict the 16th–84th percentiles, and individual galaxies are shown with circles where bins have fewer than 10 objects.
The two simulations analysed here, Recal-L25N752 and Ref-L100N1504, in the stellar mass region when they overlap, are offset from each other, suggesting that stellar mass is not the primary property determining α. This is not necessarily surprising, as the environments of galaxies with a stellar mass of |${\approx}10^{10} \, \rm M_{\odot }$| in the Recal-L25N752 are not necessarily the same as in the Ref-L100N1504 simulation. Furlong et al. (2015) showed that the SFRs of galaxies of stellar masses |$10^9{-}10^{10} \, \rm M_{\odot }$| are higher in Recal-L25N752 than in Ref-L100N1504, while becoming more similar at |${\gt}10^{10} \, \rm M_{\odot }$|. This tells us that |$\dot{M}_{\rm accr}$| is higher in galaxies of stellar masses |${\approx}10^{10} \, \rm M_{\odot }$| in Recal-L25N752 than Ref-L100N1504 (see Appendix B for details regarding the |$\dot{M}_{\rm accr}$|−M⋆ relation). We demonstrate below that |$\dot{M}_{\rm accr}$| is indeed the primary driver of α and that the offset seen in Fig. 3 is driven by the different |$\dot{M}_{\rm accr}$| Ref-L100N1504 and Recal-L25N752 galaxies of stellar masses |${\approx}10^{10} \, \rm M_{\odot }$| have.
We quantify which property, |$\dot{M}_{\rm accr}$| or stellar mass, shows a stronger correlation with the inner slope α, via the Spearman’s rank-order correlation coefficient (Rs), finding that |$\dot{M}_{\rm accr}$| gives absolute values that are similar (Recal-L25N752) or larger (Ref-L100N1504, above ≈0.1) than those obtained with stellar mass for all redshifts. Hence, in terms of scatter, α is better correlated with |$\dot{M}_{\rm accr}$| than with stellar mass.
3.4 The relation between the RMP and |$\dot{M}_{\rm accr}$| at fixed mass

The residuals of the α–M⋆ relation (Δα) as a function of the residuals of the |$\dot{M}_{\rm accr}$|–M⋆ relation (Δlog (|$\dot{M}_{\rm accr}$|)), with α being the inner slope (r/r50 ≤ 1) of the RMP. Top and bottom panels show galaxies at z = 0 and 1, respectively. Red and blue represent the Recal-L25N752 and Ref-L100N1504 simulations, respectively. The 16th–84th percentiles are depicted as dashed lines, and solid lines show medians. We show individual galaxies where bins have fewer than 10 objects (red symbols). Having eliminated the dependence on stellar mass, we still see a very strong anti-correlation between α and |$\dot{M}_{\rm accr}$|, indicating this to be independent of stellar mass. The two simulations are in good agreement. Ref-L100N1504 covers a wider range of |$\dot{M}_{\rm accr}$| due to the larger cosmological volume compared to Recal-L25N752.
By doing this calculation, we eliminate the stellar mass dependence from both axes. We find that the anti-correlation between Δα and |$\Delta \dot{M}_{\rm accr}$| remains. We conclude that, at least for the inner regions, galaxies of fixed mass with a higher |$\dot{M}_{\rm accr}$| display steeper (negative) RMPs. Although there is a general expectation that higher |$\dot{M}_{\rm accr}$| can lead to changes in the RMP, to our knowledge, this is the first time such results are quantified in cosmological simulations. Also note that for the residuals in Fig. 4, we do not see any differences between the standard and higher resolution eagle simulations, suggesting that the differences seen in Fig. 2 were due to differences in their stellar mass distributions.
3.5 The relation between the SFR, |$\dot{M}_{\rm accr}$|, and the RMP
Observationally, measuring the gas accretion rate is very difficult as the accreted gas is expected to be low density and hence be very faint (e.g. Fox & Davé 2017). In addition, this inflowing gas is expected to have a relatively weak kinematic signature (compared to outflowing gas), which makes it challenging to separate spectroscopically from the intrinsic velocity dispersion of the resident ISM (Bouché 2017). On the other hand, the integrated SFR of a galaxy is much easier to infer and is expected to be closely correlated to the gas accretion rate. The latter is the basis for the ‘equilibrium models’ of galaxy evolution. In these models, the gas accretion rate is perfectly balanced by the combination of SFR and outflow rates, and |$\dot{M}_{\rm accr}$| can then be inferred from the SFR and outflow rate. It follows then that |$\dot{M}_{\rm accr}$| is linearly correlated with the SFR (e.g. Davé, Finlator & Oppenheimer 2012; Lilly et al. 2013), modulated by the mass-loading parameter (the ratio between the outflow rate and SFR). As a consequence, one would expect the correlations discussed in Section 3.2 to also extend to the SFRs of galaxies, i.e. the RMP slope, α, decreases with increasing SFR.
A correlation between α and the SFR of galaxies can have other interpretations beside the one provided by the ‘equilibrium models’. For example, it is possible that the primary process controlling α is the metal enrichment of the ISM due to recently formed stars. In the latter case, one could imagine that the SFR may be the more fundamental parameter causing the change of the slope α rather than the gas accretion rate. This could be an interesting outcome since there still is debate as to whether there is a positive, negative, or even null correlation between α and the SFR in galaxies (see Section 1).
Both of these two scenarios would result in an anti-correlation between α and the SFR, but in the former case, one expects the scatter to be larger than for the relation between α and |$\dot{M}_{\rm accr}$|, while in the latter, the opposite would be a more natural outcome. To disentangle between these two cases, we study the same correlations of Figs 2 and 4, but for the SFRs of galaxies.
Before doing so, we verify the existence of a tight correlation between |$\dot{M}_{\rm accr}$| and SFR in Fig. 5. In agreement with the expectation of the equilibrium models described above, we find a strong correlation between |$\dot{M}_{\rm accr}$| and the SFR at all redshifts studied here, with a 1σ scatter ≲0.5 dex. The scatter significantly decreases with SFR, from ≈0.5 dex at |$\rm SFR\lesssim 0.45\, \rm M_{\odot }\, yr^{-1}$| to ≈0.1 dex at |$\rm SFR \gtrsim 6 \, \rm M_{\odot } \, yr^{-1}$|. Although we do expect an overall correlation between the SFR and |$\dot{M}_{\rm accr}$|, the scatter here may be artificially small considering our definition of gas accretion, which is based on gas accretion that leads to star formation (see Section 2.1 and Mitchell et al., in preparation, for details). Interestingly, the predicted median relation evolves only weakly with redshift, at least at z ≤ 1, with differences of ≲0.3 dex between z = 1 and 0 at fixed SFR.

|$\dot{M}_{\rm accr}$| as a function of the SFR of galaxies in the Recal-L25N752 (red) and Ref-L100N1504 (blue) simulations. Top and bottom panels show eagle galaxies at z = 0 and 1, respectively. The 16th–84th percentiles are shown as dashed lines. We show individual galaxies in bins with fewer than 10 objects as symbols. The dot–dashed lines show the 1:1 relation. We remind the reader that Recal-L25N752 galaxies were selected to have M⋆ ≥ 109 M⊙, while Ref-L100N1504 galaxies were selected to have M⋆ ≥ 1010 M⊙. Both simulations display a positive and tight correlation between |$\dot{M}_{\rm accr}$| and SFR. Also, at fixed SFR, galaxies display a weak decrease in their |$\dot{M}_{\rm accr}$| of ≈0.2–0.3 dex from z = 1 to 0 (only appreciable for galaxies with |${\rm SFR} \lt 10^{0.5} \, {\rm M_\odot yr}^{-1}$|).
The left-hand panel of Fig. 6 shows the slope of the RMP in the inner regions (r/r50 < 1) as a function of the SFR. This relation shows a similar trend and scatter as the one with |$\dot{M}_{\rm accr}$| (see Fig. 2). As we did in Section 3.3, we quantify the correlations, finding the α–|$\dot{M}_{\rm accr}$| correlation to give a similar (Recal-L25N752) or higher (Ref-L100N1504) Rs by ∼0.1 than the α–SFR relation (again, |$\dot{M}_{\rm accr}$| reaching Rs values more negative than −0.3), suggesting that the more fundamental correlation is that with |$\dot{M}_{\rm accr}$|. Tissera et al. (2019) showed that eagle galaxies with more negative slopes tend to have a larger fraction of their stars formed recently, consistent with the fact that there has been more gas accretion leading to such star formation activity, and with the clear anti-correlation we find here between α and the SFR. We caution, however, that Tissera et al. (2019) measured a single slope of the RMP rather than separating it into different regions of the disc.

Left-hand panels: solid lines show median values of inner (r/r50 ≤ 1) RMP slope as a function of the SFR at redshift z = 0 (top panel) and 1 (bottom panel), for the Ref-L100N1504 (blue) and Recal-L25N752 (red) simulations. A similar anti-correlation is obtained as for |$\alpha\!-\!\dot{M}_{\rm accr}$| (left-hand panel of Fig. 2). Right-hand panel: solid lines show median values of the residuals of the relation between the inner (r/r50 ≤ 1) RMP slope and stellar mass, as a function of the residuals of the SFR and the stellar mass. Redshifts are as labelled. Once we control for stellar mass (right-hand panels), the anti-correlation is weaker than that between Δ α and |$\Delta \log (\dot{M}_{\rm accr})$| shown in Fig. 4. This suggests that the anti-correlation between the slope of the RMP and the SFR is only a byproduct of the anti-correlation between the RMP slope and |$\dot{M}_{\rm accr}$|. All panels: the 16th–84th percentiles are depicted with dashed lines. We show individual galaxies where bins have fewer than 10 objects (circles).
As was done for |$\dot{M}_{\rm accr}$|, we remove the dependence on stellar mass by studying the residuals of the α–M⋆ relation as a function of the residuals of the SFR–M⋆ relation (see equation 3) in the right-hand panel of Fig. 6 for z = 0 and 1. As explained in Section 3.2, the residuals are constructed as the difference between the property (i.e. slopes or SFR) and the median at the stellar mass of the galaxy. As was the case for Fig. 4, controlling for stellar mass brings the two simulations into agreement. The simulation Recal-L25N752 displays a weak anti-correlation between |$\Delta \, \alpha$| and |$\Delta \, \rm SFR$| (Rs, the Spearman’s rank-order correlation coefficient, of ≈−0.2), similar to the anti-correlation between |$\Delta \, \alpha$| and |$\Delta \, \dot{M}_{\rm accr}$| (Fig. 4). However, the simulation Ref-L100N1504 shows a stronger anti-correlation with Rs ≈ −0.3. These results suggest that the α–SFR relation is a byproduct of the other properties coming into play, such as that of |$\dot{M}_{\rm accr}$| and stellar mass. Therefore, we conclude that the inner RMP slope is primarily set by the gas accretion rate rather than by the SFR.
3.6 The relation between the RMP and the gas fraction
In the so-called ‘equilibrium model’, |$\dot{M}_{\rm accr}$| regulates the gas content, SFR and metallicity of galaxies. The gas fraction of a galaxy is therefore expected to be a tracer of gas accretion (modulo the time-scale to convert gas into stars and outflows). The gas fraction is attractive as is more readily available in observations than |$\dot{M}_{\rm accr}$|, and can be linked back to the latter (although under some assumptions). Hence, one would expect that a galaxy formation simulation, which overall captures the nature of SF galaxies, produces a relation between a galaxy’s gas metallicity and gas fraction. Here, we explore the relation between the RMP and the gas fraction of galaxies.
We use the neutral gas fraction measurement of Lagos et al. (2016), defined as the neutral gas mass within 30 kpc, |$M_{\rm H\,{\small I}}+M_{\rm H_2}$|, over the sum of the former and baryon mass in the same aperture, |$M_{\star }+M_{\rm H\,{\small I}}+M_{\rm H_2}$|. The neutral gas mass is found by applying the ionized-to-neutral gas separation in post-processing following Rahmati et al. (2013) (see Lagos et al. 2015 for details of the subgrid phase partition of gas particles).
The left-hand panels of Fig. 7 show the RMP slope at r < r50 as a function of the neutral gas fraction at z = 0 and 1 for the Ref-L100N1504 (blue) and Recal-L25N752 (red) simulations. If we first look at the Ref-L100N1504 simulation, we see that there is a tendency for the RMP to decrease with radius more steeply as the gas fraction increases for all redshifts studied. This is because higher gas fractions are associated with higher |$\dot{M}_{\rm accr}$|. However, in the Recal-L25N752 simulation, we see a reversal of that relation. This is because the gas fraction is strongly anti-correlated with stellar mass (see fig. 1 in Lagos et al. 2016). Hence, it is necessary to remove the stellar mass dependence to unveil a possible correlation between the RMP and gas fraction.

Left-hand panels: solid curves show median values of the inner (r/r50 ≤ 1) RMP slope as a function of the neutral gas fraction fgas at redshift z = 0 (top panel) and 1 (bottom panel), for the Ref-L100N1504 (blue) and Recal-L25N752 (red) simulations. A similar anti-correlation is observed as in Fig. 2 for the inner part. Right-hand panels: solid lines show median values of the residuals of the relation between the inner (r/r50 ≤ 1) RMP slope and stellar mass, as a function of the residuals of the fgas and the stellar mass. Redshifts are as labelled. As in Fig. 4, the inner region RMP slope displays an anti-correlation with fgas even after eliminating the dependence on stellar mass. All panels: the 16th–84th percentiles are depicted with dashed lines. We show the individual galaxies where bins have fewer than 10 objects (red symbols).
The right-hand panel of Fig. 7 shows the relation between the residuals of the RMP slope and the residuals of the gas fraction (Δfgas) as previously defined in equation (3) (Section 3.2). The correlation is shown at z = 0 (top panel) and 1 (bottom panel). It is interesting to note that once the stellar mass dependence is removed, a strong anti-correlation is found, which is reminiscent of Fig. 4. This could be related to the fact that all the galaxy properties mentioned in this work (i.e. |$\dot{M}_{\rm accr}$|, SFR, and fgas) modulate, but to different degrees, the RMP slope at fixed stellar mass.
We compute the Spearman’s rank-order correlation coefficient Rs for the Δα–Δfgas relation (Rs,gas) and compare it with that obtained for the Δα–Δ |$\dot{M}_{\rm accr}$| relation (Rs,accr). For both relations, Rs is similar for the simulations Ref-L100N1504 and Recal-L25N752, with values of ≈−0.2 and ≈−0.3, respectively. It is worth mentioning though that the absolute values of Rs indicate the α–|$\dot{M}_{\rm accr}$| relation to be stronger than the α–fgas relation at all redshifts considered.
It is curious that the relation of the RMP slope with fgas changes so much going from the left- to the right-hand panels of Fig. 7, compared to what happens with |$\dot{M}_{\rm accr}$| (Figs 2 and 4) and the SFR (Fig. 6). These large differences arise because the correlation between |$\dot{M}_{\rm accr}$| and gas mass (H i+H2) changes with stellar mass. The relation is positive at low stellar mass, but flattens or even inverts as the stellar mass increases. Hence, at low stellar masses (M⋆ ≲ 1010 M⊙), gas accretion mostly increases the gas fraction of the galaxy (which is the case of dwarf galaxies), while at high stellar masses (M⋆ ≳ 1010.5 M⊙), gas accretion may trigger AGN feedback, which then reduces the gas fraction. Thus, even though the relation between the RMP slope and fgas gives us important constraints on the effect of gas accretion, it is neither as simple nor as direct a relation as the one we obtain between the RMP slope and |$\dot{M}_{\rm accr}$|.See Appendix C for the results when only considering the molecular gas component (H2) for the gas fraction, which can be of use since atomic gas (i.e, |${\rm H\,{\small I}}$|) is currently not accessible in observations at z ≳ 0.4.
4 DISCUSSION
We showed in the previous Sections 3.2 and 3.5 that the gas accretion plays an important role in shaping the inner slope of the RMP. We have also concluded that other galaxy properties, such as SFR, stellar mass, and gas fraction, have a secondary role in the process of altering the slope, and are most likely driven by how these correlate with |$\dot{M}_{\rm accr}$|. To better visualize the effect of gas accretion on the RMP, we show in Fig. 8 the RMP of galaxies in bins of stellar mass and |$\dot{M}_{\rm accr}$|, for the Ref-L100N1504 simulation, at z = 0 (left-hand panels) and 1 (right-hand panels). The top panels show galaxies with stellar masses in the range of M⋆/M⊙ ∈ [1010, 1010.3], while bottom panels show a stellar mass range of M⋆/M⊙ ∈ [1010.3, 1011]. We rank the |$\dot{M}_{\rm accr}$| values and present the median RMPs of galaxies in three |$\dot{M}_{\rm accr}$| percentiles: the bottom 25th, between 40th and 60th, and the top 25th. This allows us to compare galaxies in a way that is independent of the overall |$\dot{M}_{\rm accr}$| evolution.
![RMP at z = 0 (left-hand panels) and 1 (right-hand panels) for galaxies in the Ref-L100N1504 simulation. The upper panels show galaxies with a stellar mass of M⋆/M⊙ ∈ [1010, 1010.3], while the bottom panels show the range M⋆/M⊙ ∈ [1010.3, 1011]. Different colours represent the percentiles of $\dot{M}_{\rm accr}$, as labelled. The lines and shaded regions show the median RMP and its 16th–84th percentiles, respectively. Galaxies with higher $\dot{M}_{\rm accr}$ (dashed lines) show a steeper slope at r < r50. Since the softening for the Ref-L100N1504 simulation is 0.7 kpc and we find average values of r50 ≈ 4.22 kpc at z = 0 and r50 ≈ 3.11 kpc at z = 1, the profiles will be resolved for r/r50 ≳ 0.17 and r/r50 ≳ 0.23, respectively. It is interesting to note that, at fixed stellar mass, there is a vertical offset of metallicity according to different values of $\dot{M}_{\rm accr}$, as was discussed in Section 3.1.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/mnras/495/3/10.1093_mnras_staa1334/1/m_staa1334fig8.jpeg?Expires=1749333158&Signature=TQxWTVjdRGkwGcA8PgLT3~1Vi2yE-JGDjVq~PI46npzPD3H33v3T9f1SKZ6OZ06HPitayHYPaVnc3Js8F9Y28TsGhKO2voFAmYO4u2D6FNXVbJJ46yFw-W4p1ZuQfrROJ7ruiBY2shWEtfTZeNn6IUlGZvTm~6~HUDo1X45~WmG9N1zcRiSSaM4sljoA7v-ej8zgof~9Gyux75i4uPPhc9Y9~9hLT24CqzxJtKMWFlMji86UcV5FI6cyVIRyRtW6X9T5ydDDb-uplNL6o3Ex~wRedo-l9XViqLJz5VFNH1F2z5TgZ0oN7sOp9vDmdJBYTfJoWOekNi2oCGxre0y67g__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
RMP at z = 0 (left-hand panels) and 1 (right-hand panels) for galaxies in the Ref-L100N1504 simulation. The upper panels show galaxies with a stellar mass of M⋆/M⊙ ∈ [1010, 1010.3], while the bottom panels show the range M⋆/M⊙ ∈ [1010.3, 1011]. Different colours represent the percentiles of |$\dot{M}_{\rm accr}$|, as labelled. The lines and shaded regions show the median RMP and its 16th–84th percentiles, respectively. Galaxies with higher |$\dot{M}_{\rm accr}$| (dashed lines) show a steeper slope at r < r50. Since the softening for the Ref-L100N1504 simulation is 0.7 kpc and we find average values of r50 ≈ 4.22 kpc at z = 0 and r50 ≈ 3.11 kpc at z = 1, the profiles will be resolved for r/r50 ≳ 0.17 and r/r50 ≳ 0.23, respectively. It is interesting to note that, at fixed stellar mass, there is a vertical offset of metallicity according to different values of |$\dot{M}_{\rm accr}$|, as was discussed in Section 3.1.
Galaxies with higher |$\dot{M}_{\rm accr}$| show a steeper inner slope, independently of their stellar mass and redshift. In the stellar mass range |$10^{10}\!-\!10^{10.3}\rm \, M_{\odot }$| at z = 0 (top left-hand panel in Fig. 8), galaxies have similar gas metallicities at the centre; however, galaxies with low |$\dot{M}_{\rm accr}$| have RMPs with a flattened core that can extend to quite large radii (r ≈ 0.5 × r50), while galaxies with higher |$\dot{M}_{\rm accr}$| do not show a core in their RMP, and instead fall off sharply. This behaviour, however, is not universal (the top right-hand panel does not show such a sharp fall). Note that at large radii, r ≈ 10 × r50, galaxies display a flattening of their RMP, probably due to the ISM reaching the CGM metallicity. The latter happens at systematically smaller radii for galaxies with high |$\dot{M}_{\rm accr}$|. In more massive galaxies (bottom left-hand panel in Fig. 8), the flattened RMP core in galaxies of low |$\dot{M}_{\rm accr}$| is even more prominent, reaching out to r ≈ r50, while the higher |$\dot{M}_{\rm accr}$| galaxies do not display a core. At z = 1, the main difference is that the gas metallicities are overall lower than at z = 0, and the differences in the normalization of the RMPs at fixed stellar mass for different |$\dot{M}_{\rm accr}$| is larger.
Curiously, profiles with the lower values of |$\dot{M}_{\rm accr}$| in massive galaxies (i.e. with stellar masses in the range M⋆/M⊙ ∈ [1010.3, 1011]) show a dip in the very internal parts of the RMP (r < 0.5 × r50). To understand this behaviour in more detail, we analyse the ratio between central BH-to-stellar mass ratio (MBH/M⋆) of the galaxies in Fig. 8. We find that, independent of redshift, MBH/M⋆ decreases with increasing |$\dot{M}_{\rm accr}$| at fixed stellar mass, indicating that MBH is the most massive for the lowest |$\dot{M}_{\rm accr}$|. In addition, we find that the most massive galaxies at z = 1 of the bottom 25th percentile |$\dot{M}_{\rm accr}$| have a median MBH/M⋆ ≈ 0.0007, which is about twice the ratio of the other |$\dot{M}_{\rm accr}$| bins and stellar masses. At z = 0, the trend between MBH/M⋆–|$\dot{M}_{\rm accr}$| is still present but with smaller differences between the different |$\dot{M}_{\rm accr}$| percentiles. This means that BHs in the subsample of the z = 1 most massive, low-|$\dot{M}_{\rm accr}$| galaxies have the most potential to affect their host galaxies. This strong AGN activity is expected to remove the gas from the centre of the galaxy, specially the heavy elements, as the internal parts tend to be more metal-rich. Following the arguments in Bower et al. (2017), the gas in the centre of these galaxies forms stars and ignites a rapid BH growth phase, which produces an outflow, which, in turn, expels significant significant amounts of gas and metals from the galaxy before these metals have time to mix with the surrounding gas. The dip seen in galaxies with M⋆/M⊙∈ [1010.3, 1011] for the lowest |$\dot{M}_{\rm accr}$| at z = 1 is therefore likely a product of the strong AGN feedback phase these galaxies are going through, which removes metals from their centres before these metals have time to mix in the ISM.
Our analysis can also give insight into where the gas recently accreted on to the galaxy is. Galaxies with higher |$\dot{M}_{\rm accr}$| display steeper gradients in the internal parts. Thus, gas accretion appears to be diluting the metals in an outside-in fashion, which translates into more negative RMPs.
An important consideration to interpret the observed trends above is the way in which we measure |$\dot{M}_{\rm accr}$|. Because we are specifically tracking particles that at the time of interest are SF (those that have SFR > 0), but were not in the galaxy in the previous time-step, it is natural to expect our measure to be biased towards gas accretion, which leads to star formation and hence leads to more important changes in the RMP in the centre (where most star formation takes place3). With this technique, we are therefore likely seeing the most effect the gas accretion can have on the internal RMP. However, our measurement of gas accretion is a lower limit as we are, by construction, ignoring the gas particles that join the ISM but are not related to star formation, as well as particles that were ejected in between snapshots.
Another important caveat to keep in mind is that we are radially averaging the effect of gas accretion, which simulations show is not necessarily axisymmetric. This may be washing out some of the more localized effects the gas accretion can have, particularly if it is strongly filamentary (e.g. Kacprzak 2017; Putman 2017; Ho; Martin & Turner 2019). In the future, we will study the gas accretion effects in the two-dimensional gas metallicity distribution of galaxies with the aim of quantifying how much variation is expected in individual galaxies due to the generally asymmetric nature of gas accretion.
Despite these limitations, our work clearly establishes the important role of gas accretion in shaping the RMP of galaxies, which we are able to control for other effects, such as stellar mass and SFR variations, thanks to the statistics of the eagle simulations, finding that gas accretion is the primary responsible for the slope of the RMP. Note that metallicities, SFRs, and gas fractions were not used in the process of parameter tuning in eagle and, hence, the result we present in this paper represents a true prediction of the simulation.
Our results establish some clear correlations that we would expect at fixed stellar mass: galaxies with steeper (more negative α) RMPs at r < r50 tend to have higher neutral gas fractions and |$\dot{M}_{\rm accr}$| (the latter being more difficult to test).
Deep measurements of the H i content of galaxies that push down to low column densities, typical of the circumgalactic medium (Popping et al. 2009), together with metallicity gradient estimates from IFS surveys, will be an ideal combination to test our predictions. Instruments such as the Australian Square Kilometre Array (SKA) Pathfinder (ASKAP; Johnston et al. 2008) and the Karoo Array Telescope (MeerKAT; Booth et al. 2009), and, in the future, the SKA will allow the former measurements, while instruments such as SAMI (Bryant et al. 2015), MUSE (Carton et al. 2018), and other IFS surveys, allow for the latter. Note that our results cannot be easily extrapolated to z > 1 as galaxy mergers are expected to become more common (as shown by Qu et al. 2017; Lagos et al. 2018 for eagle) and because our |$\dot{M}_{\rm accr}$| ignores the gas that comes from galaxy mergers (as we are aiming to quantify ‘smooth’ accretion), this may become an important shortfall. In addition, the |$\dot{M}_{\rm accr}$| at z ≳ 1 is expected to be much more collimated and to penetrate down to the galaxy more easily than at z ≲ 1 (e.g. Correa et al. 2018). This may have the effect of directly feeding the central regions of the galaxy rather than from outside-in, possibly driving the inverted RMP seen in observations at z ≈ 3–4 (Cresci et al. 2010; Troncoso et al. 2014)
Recently, Patrício et al. (2019) presented observational measurements of the RMP of three strongly lensed galaxies at |$z=0.6, \, 0.8$|, and 1, and reported these to have more negative α than z = 0 galaxies. This is consistent with our findings as galaxies at z ≈ 1 have higher |$\dot{M}_{\rm accr}$|, therefore leading to more negative α. However, to confirm this observationally, a larger sample of galaxies is required to study the RMP at fixed stellar mass across cosmic time.
5 CONCLUSIONS
In this work, we use two simulations from the eagle project, the reference large-volume, standard-resolution simulation of |$100 \, {\rm Mpc}$| on a side (Ref-L100N1504), and the recalibrated high-resolution simulation of |$25 \, {\rm Mpc}$| on a side (Recal-L25N752), to study how the slope of the RMP of galaxies, α, changes with the gas accretion rate (|$\dot{M}_{\rm accr}$|) and other galaxy properties. Because of the resolution limits of these simulations, we focus on galaxies with |$M_\star \ge 10^{10} \, {\rm M_\odot }$| for the Ref-L100N1504 simulation, and |$M_\star \ge 10^9 \, \rm M_{\odot }$| for the Recal-L25N752 simulation. We also limit our study to central SF galaxies and to redshifts z ≤ 1.
We use a particle tracking method to find the gas particles that are being accreted on to galaxies, including those that are converted into stars. We are especially interested in smooth accretion rather than accretion in the form of galaxy mergers, and hence we select only those gas particles that were not part of another galaxy in the previous simulation snapshots. Our aim is to find whether the |$\dot{M}_{\rm accr}$| can be robustly connected to changes in the RMP as a primary driver, and hence controlling for differences in stellar mass and SFR is essential. The eagle simulations are an ideal tool for this purpose, as its combination of volume and resolution allows us to look into the internal structure of galaxies as well as providing us with enough statistics to explore galaxy properties at fixed stellar mass. Here, we focused only on central SF galaxies with at least 10 accreted gas particles in the last ≈100 Myr, which come from sources other than galaxy mergers, at redshifts z = 0 and 1.
We summarize our conclusions as follows:
The gas accretion rate is positively correlated with stellar mass (Fig. 1) and SFR (Fig. 5), at all redshifts studied. We also find that, at fixed stellar mass, higher gas metallicity galaxies are associated with lower |$\dot{M}_{\rm accr}$| (Fig. 1). Together, these results are consistent with the MZR–SFR relation (also called fundamental metallicity relation in the literature; e.g. Mannucci et al. 2010).
A tight negative correlation is found between the inner RMP slopes (measured within r/r50 ≤ 1) and |$\dot{M}_{\rm accr}$| at all redshifts studied (Fig. 2). Galaxies with higher gas accretion rates tend to have steeper RMPs. Even though galaxies change their |$\dot{M}_{\rm accr}$| with time at fixed stellar mass, this anti-correlation does not seem to evolve. At large radii, r/r50 > 1, we find a weak trend for the RMPs of galaxies to become flatter as |$\dot{M}_{\rm accr}$| increases. However, this trend is characterized by a very large scatter caused by noise, which prevents us from drawing any strong conclusion. A higher resolution simulation would be required to confirm this trend.
A clear anti-correlation between the inner RMP slope (r/r50 ≤ 1) and |$\dot{M}_{\rm accr}$| remains even when eliminating the dependence on stellar mass (Fig. 4).
The SFR is not as strongly anti-correlated with the slope α as |$\dot{M}_{\rm accr}$| is (Fig. 6), indicating that the latter is more fundamental in shaping the RMP.
We also obtain a relation between the neutral gas fraction and the slope of the inner RMP (at r/r50 ≤ 1; Fig. 7), but it is less clear than with |$\dot{M}_{\rm accr}$|. However, the gas fraction is a useful proxy as it is more readily accessible observationally than |$\dot{M}_{\rm accr}$|.
When analysing the RMP binned by redshift and stellar mass, we see that galaxies with the lowest |$\dot{M}_{\rm accr}$| show flatter inner slopes (even cored RMPs), while galaxies with the highest |$\dot{M}_{\rm accr}$| display steeper negative slopes (Fig. 8).
In the future, we will aim to reveal what properties are causing the high scatter at r/r50 > 1, as well as studying galaxies in two dimensions rather than radially averaged (Marino et al. 2016; Trayford & Schaye 2019). Even though the assumptions made in the calculation of |$\dot{M}_{\rm accr}$| are simple, it is clear that the latter plays a primary role in altering the RMPs of galaxies to a degree that we hope will be testable with a combination of sensitive IFS instruments, absorption-line studies, and deep H i observations, e.g. from local surveys such as SAMI (Bryant et al. 2015) and CALIFA (Sánchez et al. 2012), and, in the near future, the MUSE MAGPI4 survey (MAGPI collaboration, in preparation) at z ≈ 0.3.
ACKNOWLEDGEMENTS
We thank the referee for their insightful and constructive report, which helped the clarity of this paper. We thank Yannick Bahé for fruitful discussions and suggestions. 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èresle-Châtel. FC acknowledges Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina), and the Endeavour Scholarships and Fellowships for their supporting fellowships. CL is funded by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. The Cosmic Dawn Center is funded by the Danish National Research Foundation. SAC acknowledges funding from CONICET (PIP-0387), Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, PICT-2013-0317), and Universidad Nacional de La Plata (11-G150), Argentina. EW acknowledges support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.
Footnotes
We compute the median value and dispersion of the half-stellar mass radius, r50, for the whole sample at each redshift and simulation.
We remind the reader that this is done for ease of comparison with observational measurements, which typically use nebular emission lines that trace H ii regions.
REFERENCES
APPENDIX A: CONVERGENCE OF THE |$\dot{M}_{\rm accr}$|–RMP RELATION
Throughout this work, we consider central, star-forming (SF) galaxies with at least 10 SF gas particles that have been accreted (which were outside the galaxy ≈100 Myr ago). This minimum of particles translates to values of |$\dot{M}_{\rm accr}$||${\approx}9 \times 10^{-2} \, {\rm M_\odot\,yr^{-1}}$|. Here, we assess the robustness of our findings to the minimum number of gas accreted particles allowed.
We show in Fig. A1 the relation between the RMP slope residuals, Δα, and the |$\dot{M}_{\rm accr}$| residuals, |$\Delta \dot{M}_{\rm accr}$|, for four SF gas particles minimum number cuts, from at least 1 particle to 100 particles, as labelled. We remind the reader that we define these residuals in Section 3.2. For the sake of clarity, we only show the 1σ scatter for our fiducial cut of 10 particles for both simulations (the magenta shaded region depicts the Recal-L25N752 simulation, while cyan represents the Ref-L100N1504 simulation).

Same as Fig. 4, but now comparing different minimum number of gas particles of the accreted gas above which we include galaxies in our analysis. Solid, dashed, and dotted lines show |$\dot{M}_{\rm accr}$| in galaxies with ≥10, ≥50, and ≥100 gas particles of accreted gas, respectively. The error bars show the 1σ percentiles and are shown only for the case of ≥10 gas particles, for clarity.
We can see that for all redshifts analysed, there is almost no difference between the four SF gas particles number cuts. It is interesting to further analyse the cases where galaxies have fewer than 10 accreted SF gas particles. For simulation Ref-L100N1504, 17 galaxies at z = 0 (≈10 per cent of the sample at this redshift) have fewer than 10 particles and we find them to follow the same trends reported in this paper. At z = 1, only one galaxy from this simulation is found to have fewer than 10 gas particles (it has |$\dot{M}_{\rm accr}$||${\approx}10^{-2.5} \, {\rm M_\odot\,yr^{-1}}$| and α ≈ 0.015). In contrast, simulation Recal-L25N752 does not have any central SF galaxy with fewer than 10 accreted SF gas particles (at least until z ≤ 1). We hence conclude that our results are robust to the chosen minimum number of gas particles for the accreted gas, above which we include galaxies in our analysis.
APPENDIX B: GAS ACCRETION RATE AS A FUNCTION OF STELLAR MASS
In order to understand in more detail how |$\dot{M}_{\rm accr}$| affects the chemical evolution of galaxies, it is relevant to analyse how |$\dot{M}_{\rm accr}$| depends on stellar mass, M⋆. This is the purpose of this appendix.
Fig. B1 shows the median |$\dot{M}_{\rm accr}$|−M⋆ relation at redshifts z = 0 (top panel) and z = 1 (bottom panel). Solid lines show the medians, while dashed lines represent the 16th–84th percentile ranges. The black dot–dashed line is a linear function with slope ≡ 1 to use as a reference.

|$\dot{M}_{\rm accr}$| as a function of M⋆ at redshift z = 0 (top panel) and 1 (bottom panel). Simulations Recal-L25N752 and Ref-L100N1504 are shown in red and blue, respectively. Solid lines shown the median values of the relation, while dashed lines represent the 16th–84th percentile ranges. Individual galaxies are shown where bins have fewer than 10 objects (symbols). The dot–dashed line shows the function |$\rm log_{10}(\dot{\it M}_{\rm accr}) = log_{10}({{\it M}_\star }) \, + \, 9.5$|, as a reference.
As expected, the |$\dot{M}_{\rm accr}$|−M⋆ relation has a positive slope, i.e. galaxies with increasing M⋆ also show an increasing |$\dot{M}_{\rm accr}$|, and at fixed stellar mass, |$\dot{M}_{\rm accr}$| decreases with time. The relation is quite tight (about ≈0.5 dex of scatter, as quantified by the 16th–84th percentile range). A small offset between the two simulations is seen for galaxies with M⋆ ≈ 1010 M⊙, though not as pronounced as the one from Fig. 3. This offset is not unexpected due to the different environments in both simulations (Furlong et al. 2015). Despite this, we see that once stellar mass is accounted for, a clear relation emerges between |$\dot{M}_{\rm accr}$| and the slope of the RMP.
APPENDIX C: THE RELATION BETWEEN THE MOLECULAR HYDROGEN FRACTION AND THE GAS ACCRETION RATE
Since, observationally, it is easier to measure molecular gas H2 than atomic hydrogen |${\rm H\,{\small I}}$| at z ≳ 0.5, we test if our conclusions from Section 3.6 might change if only considering the former.
Fig. C1 is a replica of Fig. 7, but with a twist in the definition of the gas fraction. Now, we consider only the contribution of the H2 and so we have that gas fraction is defined as |$f_{\rm gas} = M_{\rm H_2} / (M_\star + M_{\rm H_2})$|. At first sight, one of the differences is that the values of the fgas drop considerably compare to our previous definition. This is due to the fact that |$M_{\rm H\,{\small I}}$| is more abundant than the molecular component. Another interesting thing to stand out is that the maximum values of fgas increases with redshift, showing why through observations is easier to measure the |$M_{\rm H_2}$| at redshift z ≳ 0.5.

Same as Fig 7, but now taking only the contribution of the molecular gas, H2, in the definition of the gas fraction, i.e. fgas = H2/(M⋆ + H2).
With this, we conclude that whether the definition of gas fraction involves H2 or |${\rm H\,{\small I}}$| (not shown), the results we found in Section 3.6 remain the same.