Effective yields as tracers of feedback effects on metallicity scaling relations in the EAGLE cosmological simulations

Effective yields, $y_{\rm eff}$, are defined by fundamental galaxy properties (i.e., stellar mass -$M_{\star}$-, gas mass -$M_{\rm gas}$- and gas-phase metallicity). For a closed-box model, $y_{\rm eff}$ is constant and equivalent to the mass in metals returned to the gas per unit mass locked in long-lived stars. Deviations from such behaviour have been often considered observational signatures of past feedback events. By analysing EAGLE simulations with different feedback models, we evaluate the impact of supernovae (SN) and active galactic nuclei (AGN) feedback on $y_{\rm eff}$ at redshift $z=0$. When removing supermassive black holes (BH) and, hence, AGN effects, in simulations, galaxies are located around a plane in the $M_{\star} - M_{\rm gas} - {\rm O/H}$ parameter space (being O/H a proxy for gas metallicity, as usual), with such a plane roughly describing a surface of constant $y_{\rm eff}$. As the ratio between BH mass and $M_{\star}$ increases, galaxies deviate from that plane towards lower $y_{\rm eff}$ as a consequence of AGN feedback. For galaxies not strongly affected by AGN feedback, a stronger SN feedback efficiency generates deviations towards lower $y_{\rm eff}$, while galaxies move towards the opposite side of the plane (i.e., towards higher values of $y_{\rm eff}$) as SN feedback becomes weaker. Star-forming galaxies observed in the Local Universe are located around a similar 3D plane. Our results suggest that the features of the scatter around the observed plane are related to the different feedback histories of galaxies, which might be traced by $y_{\rm eff}$.


INTRODUCTION
The simplest model proposed for the metal enrichment of a galaxy is the so-called 'closed-box', which assumes that the system is isolated from its environment.Clearly, there is no mass flow towards or outwards the system and, as a consequence, the baryonic mass remains constant with time.Following the prescriptions of a closed-box model, the evolution of a system comes down to the synergy between the stellar and gas components, whose variations are regulated only by star formation and stellar yields (Pagel & Patchett 1975).On the one hand, stars are created after the collapse of molecular clouds at a certain rate that depends on the availability of cold gas reservoirs.On the other hand, as a fair number of stars in each generation reaches the end of their lifetime, successive events of supernovae (SN) and stellar winds enrich the interstellar medium (ISM) with heavier elements.If we further assume that the ejected material from stars is homogeneously and instantaneously mixed, the evolution of the metallicity of the gas obeys a simple analytical expression: ★ E-mail: candelazerbo@gmail.com† E-mail: mariaemilia.dr@gmail.comwhere  Z denotes the true stellar yield, defined as the mass of newly produced metals (via nucleosynthesis) expelled by a generation of stars with regard to the total mass that remains inside of long-lived stars and compact remnants.Zgas () and μ() = Mgas ()/( M★ () + Mgas ()) are the gas metallicity and gas mass fraction at time , respectively, where M★ is the stellar mass.
There is abundant observational evidence that a closed-box behaviour does not constitute a suitable formation scenario for most galaxies.Instead, a more appropriate description can be achieved considering the interaction between them and the surrounding environment, through inflows and outflows of mass (e.g.Edmunds & Pagel 1984;Edmunds 1990;Tremonti et al. 2004;Dalcanton 2007;Erb 2008;Tortora et al. 2022).In this context, by comparing the predictions of a closed-box model with the behaviour of real galaxies, different works have tried to address the relative impact of feedback effects (e.g.outflows) on galaxies of different masses (e.g.Lara-López et al. 2019).Nevertheless, such approach relies on the validity of other closed-box model approximations.Regarding the instantaneous mixing assumption, it is expected to be accurate when the metallicity of galaxies is estimated by means of the oxygen abundance in HII regions (e.g.Dalcanton 2007).Since the production of oxygen is predominantly related to winds of massive stars and SN events, a short time scale for recycling mixing is a reasonable ap-proximation.As O is the most common gas-phase metallicity tracer, such assumption is generally valid.With respect to perfect mixing, the existence of metallicity gradients within galaxies indicate that mixing scales seem to be long.Nevertheless, it is expected that the latter issues lead to smaller deviations from the closed-box model than those generated by gas flows (e.g.Tremonti et al. 2004).We notice, however, that some works suggest that gas-rich dwarfs might require a more careful inspection (e.g.Werk et al. 2011).
If we invert equation ( 1) and evaluate the expression using the gas metallicity ( gas ) and gas mass fraction ( =  gas /( ★ +  gas ), being  ★ the stellar mass) of observed galaxies, we can define the effective yield as usual (e.g., Dalcanton 2007): gas ln(1/) . (2) The effective yield is constant and equal to the true stellar  Z for galaxies that evolve as closed-box systems.On the contrary, Edmunds (1990) demonstrated that inflows and outflows of material lead to  eff ≤  Z as a result of metal-enriched outflows and/or the accretion of primordial gas.The only exception that can generate the opposite trend is the accretion of gas with metallicity similar to or higher than the system.Nevertheless, this situation is highly unlikely and it is usually not considered.Furthermore, Dalcanton (2007) showed that in gas-rich galaxies the most effective mechanism for lowering the value of effective yields are metal-rich outflows.Conversely, metalpoor accreted gas does not have a significant impact on the effective yield of galaxies with high gas mass fraction (Dalcanton 2007).In this case, even though the accretion will lower the gas metallicity,  gas , the gas mass fraction will increase as well, leaving  eff almost unaffected.This particular scenario is known as a pseudo-closed-box equilibrium.
Hydrodynamical simulations have proved to be powerful tools for explaining observations and gaining more insight into the intertwined processes that occur during the evolution of real galaxies (e.g.Ma et al. 2016;Torrey et al. 2019).In particular, it has been shown that feedback processes powered up by SN events and AGN are essential in shaping galaxy metallicity scaling relations (e.g.De Rossi et al. 2017).In general, different works based on simulations show that SN feedback plays an important role on the regulation of the star formation activity of low-mass galaxies, generating a decrease in their star formation rate (SFR) and, hence, in their metallicity (e.g.Brooks et al. 2007).
State-of-the art cosmological simulations predict also a critical role of AGN feedback on the determination of metallicity scaling relations of massive galaxies.For example, De Rossi et al. (2017) carried out a detailed analysis of the 'Evolution and Assembly of GaLaxies and their Environments' (eagle, Schaye et al. 2015;Crain et al. 2015) suite of cosmological hydrodynamical simulations, showing that they are able to broadly describe the observed flattening of the mass-metallicity relation (M ★ Z g R) for massive galaxies in the Local Universe. 1 These authors reported that AGN feedback plays a central role in regulating the chemical evolution of massive galaxies, driving such behaviour.The primary cause of these effects appears to be the energy and momentum released by AGN, leading to a depletion of cold gas reservoirs in their host galaxies by heating and/or the ejection are still debated.Different observational methods for inferring the key galaxy properties yield significantly different answers (Telford et al. 2016) and, hence, the comparison between simulations and observations is not straightforward.
of metal-enriched material.Consequently, the process inhibits both star formation and chemical evolution.
In the last decades, a key relationship between  eff and the baryonic mass,  bar =  ★ +  gas , has been widely studied.At the low-mass end, several works have reported that effective yields increase with baryonic mass (Garnett 2002;Tremonti et al. 2004;Lee et al. 2006;Ekta & Chengalur 2010).Three different plausible channels have been often proposed to explain the decrease of  eff towards lower  bar : the more efficient removal of metals, via galactic winds, from shallower potential wells (Garnett 2002;Tremonti et al. 2004;Silich & Tenorio-Tagle 2001); the infall of pristine gas (Sánchez Almeida et al. 2014, 2015); and, a higher ISM mixing efficiency, driven by the migration of metal-poor gas from the galaxy outskirts towards the more metal-enriched central region (Ekta & Chengalur 2010).As higher masses are considered, Tremonti et al. (2004) first reported a change in the behaviour of observed galaxies at  bar ≳ 10 10 M ⊙ , which show a flatter  bar −  eff relation, on average.According to De Rossi et al. (2017), eagle simulations predict that  eff tends to decrease above a similar characteristic mass.In the simulations, such a trend appears to be the result of the cumulative effects of previous AGN feedback, which operate through the following mechanisms: heating the star-forming gas component, thereby suppressing star formation and leading to a passive galaxy, as well as expelling metalenriched material through galactic outflows.
LL19 performed a detailed comparison between the observed  bar −  eff relation and that obtained from eagle at  ★ > 10 9 M ⊙ , which is consistent with the mass range studied by De Rossi et al. (2017).LL19 reported, for the first time, an anti-correlation at the high-mass end for both observed and simulated galaxies.In addition, their results showed a clear bimodal behaviour when galaxies are separated by stellar age.On the one hand, younger galaxies present higher values of  eff as we consider higher masses.They also exhibit higher gas mass fractions, specific star formation rates sSFR and lower star formation efficiencies (SFE = SFR/ gas ).On the other hand, old galaxies tend to have high  bar and, as the considered stellar age increases, the  bar −  eff relation becomes flatter until an anti-correlation appears.This population is characterised by low gas fractions, sSFR and SFE, that is to say, is composed of passive galaxies whose star formation has been quenched.By using eagle simulations, in this article, we provide new insights about the origin of the  bar −  eff relation, considering different feedback models.We also try to assess the capability of  eff to diagnose the accumulated effects of SN and/or AGN feedback processes on 2D and 3D metallicity scaling relations.Given that the eagle set of simulations are publicly available and offer results from different models of SN and AGN feedback efficiencies, they are suitable for our study.Furthermore, eagle simulations were used in the work upon which we based our study, Lara-López et al. (2019) (hereafter, LL19), and show very good agreement with the observed values of  eff .Our paper is divided into the following sections.In Section 2, we present a brief description of eagle simulations and our galaxy sample.In Section 3, we explore the role of SN and AGN feedback processes on  eff and other related quantities, such as oxygen abundance, stellar mass and gas fraction.In Section 4, we analyse the 3D scaling relation defined by  ★ , oxygen abundance and gas mass, exploring the impact of feedback on its features.We discuss results from simulations and perform a comparison with observations in Section 5. Finally, our conclusions are summarised in Section 6.

THE EAGLE SIMULATIONS
The eagle cosmological hydrodynamical simulations were run with a modified version of the treepm-sph gadget3 code (Springel 2005).The joint evolution of dark matter and baryons are tracked within cosmological representative volumes, considering different periodic co-moving boxes, mass resolutions, and sub-grid physics models.The sub-grid prescriptions take into account unresolved processes, such as star formation, stellar evolution, radiative cooling, photoionization heating, metal enrichment and feedback associated with massive stars.Most eagle models include seeding, growth and merging of supermassive black holes (BH), and AGN feedback; see Schaye et al. (2015) and Crain et al. (2015), for full details.
A ΛCDM cosmology is adopted, with parameters consistent with Planck Collaboration (2015), namely Ω Λ = 0.693, Ω m = 0.307, Ω b = 0.0483,  s = 0.961,  = 0.248, and ℎ = 0.677, where symbols have their usual meaning.A glass-like particle initial configuration was implemented as initial condition, with a second-order Lagrangian perturbation following Jenkins (2010), using the public Panphasia Gaussian white noise field (Jenkins & Booth 2013).Full details about the generation of the initial conditions can be found in appendix B of Schaye et al. (2015).
Dark matter haloes are identified using the Friends-of-Friends algorithm (FoF, Davis et al. 1985), and baryonic particles are assigned to the same FoF halo as their nearest dark matter neighbour.Galaxies containing baryons and dark matter are then identified with the subfind algorithm (Springel 2005;Dolag et al. 2009).The galaxy that hosts the most bound dark matter particle in a halo is defined as the central galaxy of that halo, and the remaining subhalos are classified as satellite galaxies. 2ithin the eagle suite, different simulations are identified according to the linear co-moving extent  of the simulated cubic volume and the particle count .For example, a simulation with label L0025N0376 corresponds to a cubic box with 25 co-moving megaparsecs (cMpc) of side length, performed with 376 3 dark matter particles and an equal initial number of gas particles.Also, the complete label of a given eagle simulation includes a prefix that indicates the sub-grid model adopted.The reference model ('Ref' prefix) was calibrated to reproduce  ≈ 0 observations (see Section 2.1).A recalibrated model ('Recal' prefix) was also considered, which improves the agreement with observational data for the highest resolution simulations within eagle suite.In this article, we particularly analyse the simulations so-called RefL0025N0376 and RefL0050N0752, in order to compare them with runs corresponding to different feedback parameters but similar  −  combinations.For the sake of simplicity, the former simulations will be regarded as 'RefL25' and 'RefL50', respectively.Variations with respect to the reference model were tested by modifying one or more parameters of the sub-grid modules (Crain et al. 2015).In this work, in addition to the reference case, we analyse models that assume: a weaker and a stronger efficiency of SN feedback ('WeakFB' and 'StrongFB' models, respectively), no AGN feedback ('NoAGN' model), and an enhancement of the gas temperature increase due to AGN feedback (Δ AGN = 10 9 K; i.e. 'AGNdT9' model); see Section 2.1, for details.Simulations that apply variations of SN and AGN feedback parameters were run with a similar numerical resolution, but using simulated boxes of  = 25 cMpc and  = 50 cMpc, respectively.Hence, given the smaller volume, the former simulations include a lower number of galaxies.In the following sections, we summarise the main sub-grid models and physical parameters included in eagle that are relevant for our study.For a more complete explanation, the reader is referred to Schaye et al. (2015) and Crain et al. (2015).

Summary of EAGLE subgrid implementation
In this section, we briefly describe the most relevant aspects of subgrid physics in eagle simulations.In particular, we provide a detailed description of the key subgrid parameters involved in the different feedback models examined in this study, which are summarised in Table 1.Parameters corresponding to the reference model were calibrated considering the following  ≈ 0 observables: the galaxy stellar mass function (GSMF), galaxy sizes, and the relation between BH mass and stellar mass.Other models used in this work test the predictions of single-parameter variations with respect to the reference implementation.

Star formation and chemical enrichment
Radiative cooling, photoheating and chemical enrichment are implemented element by element following Wiersma et al. (2009a,b), tracking individually 11 chemical elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe).The model also assumes an optically thin gas component in ionisation equilibrium, which is exposed to both, an ionizing UV/X-ray background (Haardt & Madau 2001) and the cosmic microwave background.The simulations track stellar masslosses of the aforementioned elements considering three channels: stellar winds and type II supernovae (SNII) from  ★ > 6 M ⊙ stars, type Ia supernovae (SNIa) originated in catastrophic mass transfer between close binary stars, and winds from stars belonging to the asymptotic giant branch (AGB).Stellar yields that depend on the initial metal abundance are adopted: yields from Portinari et al. (1998) that consider mass-loss from massive stars were used for SNII, while yields of Marigo (2001) and Wiersma et al. (2009b) were implemented for AGB stars and SNIa, respectively.
Following Schaye & Dalla Vecchia (2008), star formation in eagle is implemented stochastically, assuming a volume density threshold of total hydrogen  * H that depends on metallicity  (Schaye et al. 2015;Crain et al. 2015).The hydrogen number density,  H , is related to the overall gas density,  g , considering  H ≡   g / H , where  is the hydrogen mass fraction and  H is the mass of a hydrogen atom (Crain et al. 2015).A temperature floor  eos , associated with the equation of state  eos ∝  4/3 g , is also applied, which is normalised to  eos = 8 × 10 3 K at  H = 10 −1 cm −3 .When gas particles fulfil the conditions  H >  * H and log 10 (/K) < log 10 ( eos /K) + 0.5, they are considered star-forming (SF) gas particles, and are assigned a star formation rate (SFR) that follows the Kennicutt-Schmidt relation (Kennicutt 1998).

Energy feedback from star formation
Thermal feedback from star formation is applied stochastically, adopting the feedback model of Dalla Vecchia & Schaye (2012).This model carries out a stochastic selection of neighbouring gas particles that are heated by a temperature increment of 10 7.5 K.In particular, a fraction  th of energy from core-collapse supernovae (SNII) is injected, taking into account the local metallicity and gas density.Such energy is released into the ISM 30 Myr after the birth of a stellar population (Schaye et al. 2015;Crain et al. 2015), and is Table 1.Main physical parameters of the eagle simulations used in this work.From left to right, the columns show the simulation identifier, side length of the co-moving volume (), number of dark matter particles ( , which is similar to the initial number of gas particles), initial baryonic particle mass ( g ), dark matter particle mass ( DM ), asymptotic maximum (  th,max ) and minimum (  th,min ) values of  th (see equation 3), the parameters that control the characteristic density and the power-law slope of the density dependence of the energy feedback from star formation ( H,0 and  n , respectively), the (sub-grid) accretion disc viscosity parameter ( visc ), the temperature increment of stochastic AGN heating (Δ AGN ), and the number of galaxies extracted from the simulation ( gal , see Section 2.2).The upper section corresponds to simulations run with the reference model, which has been calibrated to reproduce  ≈ 0 observations (RefL25 and RefL50 simulations).The lower section comprises simulations run with models featuring single-parameter variations with respect to the reference one.given by: where  ⊙ = 0.0127 is the solar metallicity,  H,birth is the density inherited by the star particle from its parent gas particle,  is the metallicity,  th,min and  th,max are the asymptotic values of  th , and   ,  n and  H,0 are free parameters that were chosen to reproduce the  ≈ 0 GSMF and the galaxy mass-size relation (see Schaye et al. 2015, for details).In the model,   =  n is assumed.
In eagle, galaxy formation and evolution is governed primarily by the supply of gas into the ISM, being this process regulated by feedback.Changing the efficiency of star formation feedback has a significant impact on many galaxy properties (e.g.Crain et al. 2015 and references therein;De Rossi et al. 2017).The reference eagle model adopts  th,max = 3 and  th,min = 0.3, while, in the WeakFB and StrongFB models, these values are scaled by a factor of 0.5 and 2, respectively (see Table 1).We note that the StrongFB model predicts significant changes in the stellar mass fraction ( ★ / h ) as a function of halo mass ( h ) when compared with the reference model (Crain et al. 2015), which is the one calibrated against observations.In particular, the StrongFB model predicts the formation of galaxies with significantly lower  ★ and metallicities, at a given  h , than the reference prescription.In this article, the StrongFB and WeakFB models are considered to test the effects of varying SN feedback efficiencies on metallicity scaling relations but, given that they are not adjusted to reproduce observations, caution should be taken if using them to interpret the behaviour of real galaxy populations.

AGN feedback
In eagle simulations (except for the NoAGN model), feedback from AGN quenches star formation in massive galaxies, shapes the gas profiles in the inner parts of their host haloes, and regulates the growth of BH.When the host halo mass of a given galaxy increases above 10 10 ℎ −1 M ⊙ , a seed black hole of mass 10 5 ℎ −1 M ⊙ is placed inside it following Springel et al. (2005).The black hole then grows as a result of mergers and gas accretion.The model assumes a modified Bondi-Hoyle accretion rate (Rosas-Guevara et al. 2015;Schaye et al. 2015) that is regulated with a viscosity parameter  visc .AGN feedback is implemented thermally and stochastically, choosing random neighbouring particles and heating them with a temperature increment Δ AGN .A higher value of Δ AGN implies more energetic and intermittent feedback events, leading to reduced radiative losses within the ISM.As shown by De Rossi et al. (2017), a higher Δ AGN drives a stronger AGN feedback impact on metallicity scaling relations of massive galaxies.
With the exception of the NoAGN model, the eagle simulations analysed here adopt a value of  visc = 2.In the reference model, Δ AGN = 10 8.5 K is assumed, while, in the AGNdT9 model (higher AGN feedback temperature increment), Δ AGN = 10 9 K.In the NoAGN model, the BH growth and AGN feedback implementations are entirely turned off.
Finally, it is worth mentioning two caveats.Throughout this study, we assess the characteristics of black holes associated with each galaxy by utilising the galaxy catalogues publicly available through the eagle database.As mentioned in McAlpine et al. (2016), the variable representing the black hole mass (denoted  BH , in this article) does not correspond to the mass of the central BH in a galaxy, but rather represents the cumulative value of all BH assigned to that subhalo.Nonetheless, for  BH > 10 6 M ⊙ , this closely approximates the mass of the most massive BH.Additionally, care should be taken when interpreting the variable used to quantify the accretion rate of BH (  BH,acc , in this paper).The time sampling of the simulation outputs may not capture the high temporal variability of BH accretion rates accurately.
As previously indicated, Table 1 summarises the values of the main subgrid parameters implemented in the eagle simulations used in this paper.

Galaxy sample and definitions
Following De Rossi et al. (2017), we selected simulated galaxies with  ★ ⩾ 10 9 M ⊙ to avoid resolution issues.The number of systems resulting from this selection criteria in our different simulations is indicated in the last column of Table 1.Note that both, central and satellite galaxies are included in our galaxy samples.
As different feedback models predict galaxy populations with different galaxy size distributions (see, e.g.Crain et al. 2015), our main analysis is based on integrated galaxy properties calculated without imposing aperture limits (i.e.all particles identified as belonging to the galaxy by subfind are considered).In this way, we can focus on the general effects of feedback on the global properties of galaxies and not on their spatial distribution.Thus, unless stated otherwise, we do not apply aperture corrections to galaxy properties.However, for the purpose of comparing our results with observations presented in LL19 (Section 5.3), we recompute simulated quantities within similar apertures as those employed in their fiducial sample, to which we will refer to as 'LL19 sample'.As shown in De Rossi et al. (2017) and LL19, aperture effects do not significantly affect the main features of  eff and metallicity scaling relations, but can generate moderate changes in their slopes and absolute normalisations.
In this work, the effective yield  eff of a simulated galaxy is defined by equation ( 2).Since gas-phase metal abundances are usually inferred from SF regions, unless otherwise specified and for the sake of consistency, we derive  eff ,  and  bar considering only the SF gas component of our simulated galaxies.Besides, throughout this article,  gas and  gas stand for the mass and metallicity of the SF gas component of a given galaxy, respectively.In addition, for the sake of comparison with previous works, we characterise the metallicity in terms of oxygen abundance (O/H), since oxygen is generally the most abundant heavy element in mass (Maiolino & Mannucci 2019).We estimate O/H as usual: 12 + log(O/H) = 12 + log( O / H ), where   represents the number density of the chemical element  in the SF gas component of a galaxy.A comparison with similar quantities obtained from the total gas component is carried out in Section 5.2.

FEEDBACK EFFECTS ON EFFECTIVE YIELDS AND ASSOCIATED 2D SCALING RELATIONS
The effective yields of galaxies combine information regarding relevant properties of these systems, such as their gas and stellar content, and the corresponding metallicities.In this section, we study the effects of SN and AGN feedback on 2D fundamental scaling relations that involve such quantities.We try to determine how different feedback models affect  eff .

Effective yield-baryonic mass relation
As mentioned in the Introduction, the well-known  eff −  bar relation has been frequently used in the literature to evaluate the relative accumulated effects of feedback in observed galaxies of different masses.
In Fig. 1, we analyse the  eff −  bar relation for our set of simulations.Symbols are colour-coded with  BH / ★ , which quantifies the dominance of the BH in each galaxy. BH / ★ can be regarded as a rough measure of the accumulated effects of AGN feedback during the evolution of a galaxy, in the sense that, at a fixed total stellar mass, galaxies with larger (smaller)  BH are expected to have been more (less) significantly affected by AGN feedback along their formation histories.
According to the top panels in Fig. 1 and consistently with results from LL19, variations of the AGN feedback model affect mostly  eff of massive galaxies ( bar ≳ 10 10 M ⊙ ).At low masses, where  BH / ★ ≲ 10 −3.5 for the bulk population, all simulations predict similar weak  eff −  bar correlations with low scatter and high  eff .At higher  bar ≳ 10 10 M ⊙ , the NoAGN model predicts the highest  eff at a given mass.On the other hand, the RefL50 and AGNdT9 models predict a strong median decrease of  eff with  bar , caused by the presence of massive galaxies with dominant BH ( BH / ★ ≳ 10 −3.5 ) and very low  eff .The scatter at high masses also increases with Δ AGN .For the RefL50 and AGNdT9 simulations, the low  eff of massive galaxies can be associated with the accumulated impact of AGN feedback driven by their dominant BH, which can heat SF gas, suppress star formation and chemical evolution, and drive outflows of metal-enriched material out of galaxies (see, De Rossi et al. 2017 and LL19, for a discussion).Thus, our findings suggest a close connection between  eff and  BH / ★ .
The effects of varying the SN feedback efficiency are presented in the bottom panels of Fig. 1.In the case of the strong SN feedback model, all galaxies show low  BH / ★ and  eff , which are below those corresponding to systems of similar masses in the RefL25 model.On the other hand, the WeakFB model predicts a higher percentage of galaxies with high  BH / ★ at all masses.At low masses, the WeakFB and RefL25 models predict similar median  eff −  bar relations, whereas, at higher masses, the WeakFB simulation shows a stronger  eff −  bar anti-correlation, as a consequence of the higher number of galaxies with dominant BH.The WeakFB model also predicts a larger scatter.Hence, it is important to highlight that the impact of SN feedback extends beyond the direct removal of metal-enriched material, as observed in the StrongFB simulation.SN feedback also appears to affect  eff through the regulation of the BH-to-stellar mass growth: interestingly, a weaker SN feedback tends to drive a more pronounced impact from AGN feedback, whereas a stronger SN feedback diminishes the relevance of AGN feedback.As already discussed by Henriques et al. (2019), SN feedback has a critical role in regulating the efficiency of AGN feedback; when the stellar mass reached by the galaxy is large enough to avoid mass ejection by SN feedback, more cold gas becomes available for star formation and BH growth, thus giving place to AGN feedback.
According to the aforementioned findings, both SN and AGN feedback seem to affect the shape, normalisation and scatter of the  eff −  bar relation in different ways, which suggests that its features can provide clues about the relative accumulated effects of feedback in real galaxies of different masses.Considering the predictions of our simulations: 1) at a given  bar > ∼ 10 10 M ⊙ , lower  eff can be associated with a stronger impact of AGN feedback; 2) at a given  bar < ∼ 10 10 M ⊙ , galaxies with higher than the median  eff can be associated with a weaker SN feedback efficiency and a weaker AGN feedback impact; 3) at a given  bar < ∼ 10 10 M ⊙ , galaxies with  eff well below the median, correspond generally to systems with a strong AGN impact and a weak SN feedback efficiency.
Finally, we note that some previous works that analyse EAGLE galaxies in the context of analytical models adopt a 'true stellar yield' parameter   = 0.04 (log 10 (  ) ≈ −1.4; e.g.Sharma & Theuns 2020;Zenocratti et al. 2022), which seems to be suitable for describing the behaviour of most simulated systems.Such value is consistent with the net yields expected for simple stellar populations (SSPs) with a Chabrier IMF (e.g.Madau & Dickinson 2014).The effective yields shown in Fig. 1 are below the aforementioned true yield parameter, a fact that is expected as simulated galaxies are not closed boxes.In addition to SN and AGN feedback, they are affected by the accumulated effects of past mergers and gas flows.Besides, simulated galaxies are composed by a mix of SSPs with different metallicities and, thus, different net yields.

Metallicities, stellar masses and gas fractions
As effective yields are calculated from metallicities, stellar masses and gas fractions (equation 2), the study of feedback effects on such quantities is relevant to explain their impact on  eff .In this section, we start exploring well-known 2D galaxy scaling relations involving O/H,  ★ and .Consequences of varying AGN and SN feedback models on such relations are analysed in Fig. 2 and Fig. 3, respec- tively, where scatter plots are colour-coded by  eff , and symbols are scaled to the value of  BH / ★ .

AGN feedback effects
Left panels of Fig. 2 show the M ★ Z g R for different AGN feedback models.At the low-mass end, different feedback prescriptions predict high  eff and a slightly similar positive M ★ Z g R slope.We notice a break-point at  ★ ∼ 10 10.3 M ⊙ , where the RefL50 and AGNdT9 models predict a transition from a correlation to an anti-correlation for higher masses.On the other hand, a steeper correlation is obtained for the NoAGN model at such intermediate masses, with a gradual flattening towards  ★ ≳ 10 11.5 M ⊙ .The steeper negative slope of the high-mass end of the M ★ Z g R with increasing Δ AGN is accompanied by an increasing number of metal-poor galaxies with low  eff .In fact, comparing the three simulations, galaxies with the lowest  eff are obtained in the RefL50 and AGNdT9 implementations, and correspond to systems with high  ★ and low O/H.In addition, the scatter plots symbols in Fig. 2 are scaled with  BH / ★ .It is clear that massive galaxies with lower metallicities have more dominant BH, indicating that they could have been more affected by AGN feedback.On the contrary, we note that the highest values of  eff are obtained in the simulation NoAGN and correspond to systems with low  ★ , high O/H and null  BH / ★ .Our findings are in agreement with De Rossi et al. (2017), who claimed that the negative slope of the M ★ Z g R at high  ★ can be driven by AGN feedback through two different main channels: i) the heating of SF gas and its subsequent change to a non-star-forming (NSF) phase, which drives the quenching of star formation and, consequently, the suppression of chemical evolution; ii) the ejection of metal-enriched gas out of galaxies. eff seem to be good tracers of this situation, reaching lower and lower values as the global metallicity decreases due to the increasing AGN feedback influence.
As seen in the middle panels of Fig. 2, all AGN feedback models generate similar median  ★ −  anti-correlations, with  ranging between 0 and 0.5 for smaller galaxies, and showing almost negligible values, for massive ones.Although weaker, there is a trend for massive galaxies of decreasing their SF gas fraction with increasing Δ AGN .Such a behaviour is expected given that AGN feedback heats the gas, fostering its transition to a NSF phase.It is also worth noting that, as Δ AGN increases, there is a higher number of galaxies with low  eff and very low , specially for high-mass galaxies.Therefore, considering equation (2) and our previous analysis of the M ★ Z g R, the lower  eff at the high-mass end are a consequence of both, the lower  and lower metallicity of the corresponding galaxies.
Finally, right panels of Fig. 2 show the  − O/H relation.Consistently with the findings of De Rossi et al. (2017) for the eagle 'Recal' model, our simulations predict a tight decrease of metallicity with increasing SF gas fraction, showing a larger scatter towards lower values of .Typical galaxies with  > ∼ 0.1 depict a median relation that follows roughly a curve corresponding to a constant log 10 ( eff ) ≈ −1.76.Interestingly, such tight relation represents well the behaviour of galaxies derived from our three different AGN feedback models, with the only exception of gas-poor systems.In particular, when AGN feedback is turned off, a higher number of gas-poor systems show higher metallicities and  eff .On the other hand, the RefL50 model predicts a significant population of gas-poor systems with lower than average metallicity and lower than average  eff .And, the latter population is even larger for the AGNdT9 model.As shown before, very low  eff are related to the presence of galaxies with dominant BH (i.e.high  BH / ★ values) and, hence, probably more affected by AGN feedback.In particular, three ranges of  eff are clearly distinguished in the  − O/H plane: 1) galaxies with high  eff follow the log 10 ( eff ) ≈ −1.76-curve, showing slightly higher than average metallicity at a given ; 2) galaxies with intermediate  eff also follows the same curve, showing slightly lower than average metallicity at a given ; 3) low  eff correspond to galaxies with low metallicities and low gas-fractions, significantly departing from the log 10 ( eff ) ≈ −1.76-curve.The latter region is only significantly populated when AGN feedback is turned on and correspond to massive galaxies ( ★ > ∼ 10 10 M ⊙ ) with the most dominant BH.

SN feedback effects
In this section, we explore the effects of varying SN feedback efficiency on O/H,  ★ and .Given the smaller box of the simulations analysed here (RefL25, WeakFB, StrongFB; see Section 2), the number of galaxies in our selected samples is smaller than those studied in the previous section (RefL50, NoAGN, AGNdT9; see Table 1).We also note that, in the case of the StrongFB simulation, the number of galaxies is a factor of ≈ 2 lower than for the WeakFB run (Table 1).This is related to the stronger efficiency of feedback in the former case, which leads to a decrease in the star formation of galaxies due to gas heating and SN-driven winds.These effects prevent that many galaxies surpass our imposed lower mass limit of  ★ = 10 9 M ⊙ .
In Fig. 3, left panels, the effects of different SN feedback prescriptions on the M ★ Z g R are compared.We clearly see distinct patterns for the different simulations.As reported in previous works, the RefL25 model predicts a flat median relation with moderate scatter.Regarding the WeakFB model, it leads to a higher normalisation for the median M ★ Z g R, which shows slightly higher positive (negative) slope at low (high)  ★ .The WeakFB model also produces a larger scatter in metallicity at all masses.Such scatter is mostly generated by the appearance of a population of galaxies with low metallicities, low  eff and dominant BH (larger symbols, see also Fig. 1).But, in contrast with RefL50 and AGNdT9 simulations (see Fig. 3), a significant number of dominant BH are located in low-mass galaxies in the WeakFB simulation.This could be caused by the weaker SN efficiencies in the latter case, which allows more gas to cool down and reach the galaxy centre, enhancing the growth of the central BH (see, also, the discussion in Section 3.1).According to these results, a weak feedback efficiency can lead to significant AGN feedback impact in galaxies of all masses, contributing to lower the metallicity even for low-mass systems.On the other hand, the StrongFB model do not lead to the formation of dominant BH, driving a M ★ Z g R with a positive slope along the whole mass range.These trends support again a scenario where dominant BH are required for reproducing the flattening of the M ★ Z g R at high masses (De Rossi et al. 2017).In addition, the normalisation of the M ★ Z g R is lower for the StrongFB model, which is expected as a stronger SN feedback fosters the ejection of metal enriched material out of galaxies.In the middle panels of Fig. 3, we compare the  −  ★ relation for our three SN feedback scenarios.In the case of the reference model,  decreases with increasing  ★ , with a large scatter around the median relation.Low values of  eff can be associated with gaspoor galaxies.A similar behaviour is obtained for the WeakFB model but, with a lower normalisation for , and a higher number of gaspoor galaxies with low  eff .This is consistent with a more efficient SF process and gas consumption in the WeakFB simulation: given the less significant influence of SN feedback, galaxies tend to form their stellar component earlier.In addition, as previously discussed, low- eff galaxies in these simulations also have dominant BH, which prevent further gas cooling, quenching star formation and chemical evolution.It is interesting to see that only systems with dominant BH reach negligible .In the case of the StrongFB simulations and contrary to expectations, low-mass galaxies show higher median  than in the other simulations.This can be due to a stronger SN feedback heating during the peak of the SF histories of these systems, which could caused the ejection of gas from the shallower potential well of these galaxies, with the ejected material remaining in the outer part of them.As the potential wells become deeper at later epochs, the ejected material could have been re-accreted, driving an increase in  and a decrease in metallicity; this explains the presence of a more significant population of galaxies with low O/H and high  in the StrongFB simulation.
With respect to the  − O/H relation (Fig. 3, right panels), the RefL25 model predicts similar trends to those discussed before for the RefL50 simulation: O/H decreases with increasing , a larger scatter is obtained at lower  and, once more, the median  − O/H relation follows a curve of constant log 10 ( eff ) ≈ −1.76.But, in contrast to the results obtained from variations of Δ AGN , variations in SN feedback efficiencies generate more significant departures from the log 10 ( eff ) ≈ −1.76 curve.For the WeakFB model, the  − O/H relation is flatter, shows larger scatter and exhibits a small offset towards higher metallicities (and, hence, higher  eff ) with respect to the RefL25 model.On the other hand, for the StrongFB model, the  − O/H relation presents a similar slope to that associated to the RefL25 model, but departs around 0.2 dex towards lower metallicities (following roughly the curve of constant log 10 ( eff ) ∼ −2).

Effective yield-BH connection
The results obtained so far suggest the existence of an anti-correlation between  eff and  BH / ★ .This is verified in Fig. 4, which presents a scatter plot of the dependence of  eff on  BH / ★ for all models considered.Note that symbols are colour-coded according to  ★ .
It is clear that, for galaxies with high  BH / ★ ≳ 10 −3.5 , all simulations predict a strong anti-correlation between  eff and  BH / ★ .Although the median  eff − BH / ★ relation seems not to depend on the feedback model, the number of galaxies with higher  BH / ★ increases for higher AGN heating temperatures.Besides, galaxies with  BH / ★ ≳ 10 −3.5 tend to have  ★ ≳ 10 10 M ⊙ in all simulations, with the only exception of the WeakFB run.For the latter model, lower mass galaxies can have dominant BH, as we showed before.Our results indicate that AGN feedback seems to be the main responsible for the decrease of  eff as  BH / ★ increases at  BH / ★ > ∼ 10 −3.5 .On the other hand, SN feedback has influence on the BH growth through the regulation of the gas reservoir.Clear signatures of SN feedback impact on  eff are evident at low  BH / ★ , where lower  ★ and an almost constant  eff value are obtained for all simulations.At  BH / ★ ≲ 10 −3.5 , the StrongFB and WeakFB models predict a median  eff of ∼ 10 −1.9 and ∼ 10 −1.7 , respectively.On the other hand, the reference and AGNdT9 models show, for similar  BH / ★ ,  eff ∼ 10 −1.76 , which corresponds to the median  eff for the NoAGN case.Additionally, it is interesting to note that the StrongFB model predicts  eff within a narrow range of intermediate values (≈ 10 −2.2 − 10 −1.8 ) for almost all galaxies, compared with the wider range of  eff values covered by other feedback models (≈ 10 −2.4 − 10 −1.6 ).Our trends could be explained considering that a strong SN feedback efficiency prevents from reaching very high  eff , whereas the lack of very massive BH (required for a strong AGN feedback impact), in the case of a strong SN feedback, prevents from reaching very low  eff .To sum up, Fig. 4 suggests that, in the case of galaxies with no dominant BH,  eff seem to be mostly determined by the efficiency of SN feedback, while, for galaxies with dominant BH,  eff show a strong anti-correlation with  BH / ★ due, mainly, to AGN feedback.
Finally, we note that the  BH / ★ range obtained from our analysis is consistent with recent observations of massive galaxies (e.g.Graham & Sahu 2023).This is not a surprise as the reference model in eagle simulations was calibrated using the relation between BH mass and stellar mass in a similar mass range (Section 2).
Summarising the results from Section 3, AGN feedback tends to favour the simultaneous decrease of the SF gas fraction and metallicity of high-mass galaxies, leading to lower  eff .On the other hand, a stronger SN feedback efficiency leads to lower  eff at all stellar masses as a consequence of the decrease of metallicity at a given SF gas fraction.And, a weak SN feedback efficiency drives a more complex behaviour, fostering the formation of dominant BH for a significant number of galaxies even at low  ★ : galaxies with dominant BH show lower O/H, lower SF gas fractions and, hence, lower  eff , while galaxies with low  BH / ★ exhibit higher metallicities at a given SF gas fraction, and, hence, higher  eff .

THE STELLAR MASS-METALLICITY-GAS MASS PARAMETER SPACE
Lara As  eff are defined from these key galaxy properties and taking into account the results discussed in the previous section, it is reasonable to expect that  eff could trace the impact of feedback processes on the features of the aforementioned 3D galaxy metallicity scaling relations.In this section, we examine the relation between  ★ , 12 + log 10 (O/H) and  gas (hereafter, M ★,g Z g R), for our complete set of eagle simulations.We aim at evaluating how feedback processes affect the relations between these three fundamental properties and try to assess the role of  eff as a feedback indicator.Fig. 5 and Fig. 6 show the M ★,g Z g R for simulations with different AGN and SN feedback prescriptions, respectively.Different symbols are colour-coded according to  BH / ★ .For the NoAGN simulation,  BH is always zero so galaxies are shown with an uniform grey colour.Comparing results from different panels, we see that AGNdT9 and WeakFB simulations present the largest dispersion.On the other hand, the smallest dispersion is obtained for the NoAGN simulation, which, indeed, predicts an M ★,g Z g R that can be well represented by a plane.Such NoAGN plane also seems to be a good representation of the bulk of galaxies with no dominant BH (i.e.low  BH / ★ ) in other simulations.
We determined the characteristic plane associated with the NoAGN model by performing a least-square fit, of a two order polynomial, with the R hyperplane fitting package hyper-fit (Robotham & Obreschkow 2015).The plane is well described by the following expression: where the parameters  1 ,  2 and  3 are shown in Table 2 for different data samples.The first row in the table corresponds to the plane plotted in Fig. 5 and Fig. 6, and discussed here.We will refer to this plane as the NoAGN fitting plane (hereafter, 'NoAGN-fp').
The departures (residuals) from the NoAGN-fp can be quantified by calculating the orthogonal deviation of each galaxy  from such a plane as: where the subscript  indicates quantities corresponding to the given galaxy .
For the sake of clarity, the following convention will be adopted in this work: at a fixed  ★ and O/H, galaxies with higher (lower)  gas than the value corresponding to the NoAGN-fp are assumed to be located 'over' ('under') the plane, presenting positive (negative) deviations with respect to it.
A comparison between the NoAGN-fp and the features of the M ★,g Z g R for different feedback models is carried out in Section 4.1.In Section 4.2, we analyse the connection between the latter findings and the distribution of  eff values for galaxy populations in different simulations.In Section 4.3, we evaluate how different feedback scenarios can affect the deviations of the M ★,g Z g R from the NoAGN-fp.

Feedback impact on the M
In Fig. 5 and Fig. 6, we compare the best fitting plane obtained for the NoAGN simulation (grey surface) with the M ★,g Z g R associated to different feedback models.We see that the NoAGN-fp can roughly describe the behaviour of the vast majority of galaxies with no dominant BH.In addition, when varying SN or AGN feedback prescriptions, clear distinct trends are obtained for the location of galaxies with respect to the NoAGN-fp.Fig. 5 evaluates the impact of AGN feedback on the M ★,g Z g R. As previously discussed, the increase of Δ AGN leads to a more significant number of galaxies with higher  BH / ★ .At the same time, the AGNdT9 model predicts the highest number of galaxies spread below the NoAGN-fp, displaying also the most substantial negative deviations.Another key feature is the dependence of the distance to the plane on the dominance of the BH, which is quantified by the parameter  BH / ★ .Galaxies with lower  BH / ★ tend to be located closer to the plane, with almost null deviations with respect to it or even slightly positive ones.Conversely, systems with higher  BH / ★ show larger negative deviations.Therefore, our findings suggest that galaxies tend to deviate down from the NoAGN-fp as they are more affected by the heating of gas by AGN feedback.
The consequences of varying SN feedback are analysed in Fig. 6.In principle, we can clearly see that a change in the SN feedback efficiency affects not only the scatter of the M ★,g Z g R but also the displacement of the bulk of the galaxy population with respect to the NoAGN-fp.For the StrongFB model, all galaxies tend to be located below the NoAGN-fp.As discussed before, a stronger SN feedback leads to a moderate average decrease in O/H, at a given , for most of the galaxies (Section 3), which explains their moderate displacement down the NoAGN-fp.For such galaxies,  BH / ★ ∼ 10 −4 − 10 −3 , which are values well below the maximum ones reached by galaxies in other simulations ( BH / ★ ∼ 10 −2.5 − 10 −2 ).Thus, in the case of the StrongFB model, larger deviations from the NoAGN-fp are not possible given that AGN feedback effects seem to be limited.As SN feedback efficiency decreases (see plots corresponding to the simulation WeakFB), most of the galaxies move, on average, closer to the NoAGN-fp, some of them reaching positive deviations with respect to it.It is interesting to highlight that, in the WeakFB simulation, galaxies above the NoAGN-fp (purple galaxies) are those characterised by low values of  BH / ★ ( < ∼ 10 −3.5 ).Thus, these galaxies more closely resemble a closed box model for two main reasons.Firstly, the low SN feedback efficiency implemented in this model ensures a weak effect of stellar winds and energy injection due to SN feedback.Secondly, the low values of  BH / ★ indicate that they should not have been strongly affected by gas outflows caused by AGN feedback.On the other hand, the WeakFB model also predicts a significant number of galaxies with high  BH / ★ , which tend to be located at larger distances below the NoAGN-fp.As mentioned before, SN and AGN feedback processes seem to be indirectly intertwined for such galaxies: due to their lower SN feedback efficiency, they could have evolved faster due to an early efficient consumption of their cold gas reservoir, reaching  = 0 with higher  BH / ★ and, thus, having been more affected by AGN feedback.
A comprehensive analysis of the differences between galaxies in proximity to the NoAGN-fp and those situated at greater distances is undertaken in greater depth within Section 5.1.

Effective yields as feedback tracers
To evaluate the connection between  eff and the feedback-driven features of the M ★,g Z g R in Fig. 5 and Fig. 6, we can compare the M ★,g Z g R predicted by different feedback models with the  eff distribution of simulated galaxies.In this sense, we note that, since  eff is given by gas metallicity and gas mass fraction (equation 2), a unique  eff value is associated to each point within the 3D parameter space defined by  ★ , O/H and  gas .Interestingly, simulated galaxies located around the NoAGN-fp show a roughly constant  eff ≈ 10 −1.76 , being these trends consistent with the flatter  bar −  eff relation described by galaxies in the NoAGN simulation (see Fig. 1).Such characteristic value  eff ≈ 10 −1.76 approximates the average  eff of galaxies in the NoAGN model.In addition, the bulk of the galaxy population in all simulations is located within a region where the NoAGN-fp seems to be a good local representation of a surface defined by  eff ≈ 10 −1.76 .In this context, orthogonal deviations from the NoAGN-fp in our 3D parameter space can be associated, locally, with variations in the values of  eff .As we will show in next section, the highest  eff are reached by those galaxies with the largest positive deviations with respect to the plane, which correspond to systems with lower  BH / ★ .On the other hand, galaxies that are located under the NoAGN-fp plane with larger deviations from it show higher  BH / ★ and lower  eff (see, also, Fig. 4).
It is important to remember that our galaxy samples include both central and satellite galaxies of haloes (Section 2.2).We verify that similar general trends would be obtained if only central galaxies were considered for our analysis; specially, a similar NoAGN-fp would emerge since few satellites lie within the considered mass range in the NoAGN simulation.Satellite galaxies only tend to increase the dispersion towards lower masses as a consequence of the presence of systems with high  eff .This is consistent with the enhanced metallicity of eagle satellite galaxies reported by Bahé et al. (2017).These authors studied the eagle reference model and concluded that satellite galaxies tend to be more metal enriched than equally massive centrals.They demonstrated that this phenomenon is primarily driven by the removal of metal-poor SF gas from the outer regions of galaxies, accompanied by the suppression of metal-poor inflows resulting from the removal of gas from the galaxy halo.Consequently, satellite galaxies in eagle simulations present higher metallicities than central galaxies of the same mass, showing positive deviations from the NoAGN-fp.

Residual Analysis
In this section, we quantify the deviations from the NoAGN-fp by using the residuals  defined in equation ( 5).We remind that, at a given  ★ and metallicity,  is taken to be positive (negative) for higher (lower) values of  gas with respect to the NoAGN-fp.
Fig. 7 shows  as a function of  BH / ★ for different feedback models.Symbols are colour-coded according to  eff .For models including the reference SN feedback efficiency (RefL25, RefL50, AG-NdT9), galaxies with less dominant BH ( BH / ★ < ∼ 10 −3.5 ) show negligible median residuals.This is expected since the NoAGN feedback model, from which the NoAGN-fp is derived, adopts a reference SN feedback efficiency and considers no AGN feedback effects; hence, such a behaviour should be recovered for galaxies with low  BH / ★ in simulations with the default SN feedback model.On the other hand, for galaxies with no dominant BH, the residuals become negative (positive) when implementing a stronger (weaker) SN feedback efficiency relative to the reference model.In the case of galaxies with more dominant BH ( BH / ★ ≳ 10 −3.5 ),  decreases with  BH / ★ in all simulations, as it is also evident from Fig. 5-6; this is a consequence of the increasing relevance of AGN feedback.
The tight relation between  eff and  is clear in Fig. 7 (compare it, also, with Fig. 4), which can be well described by a linear function, as we checked.By performing a linear least squares fit to the data obtained in the NoAGN model, this relation takes the form: being valid for galaxies in all simulations.As discussed before, the NoAGN-fp locally represents a surface of constant  eff in our 3D parameter space. measures the length of an orthogonal vector from the plane to a given galaxy, so that this vector should be parallel to the gradient of the scalar field defined by  eff .Hence, increasing , implies continuously crossing isosurfaces associated with different  eff values.

DISCUSSION
Our results suggest a close connection between the accumulated effects of AGN and SN feedback, and the features of scaling relations involving stellar mass, metallicity and gas mass.In addition, given that  eff is defined from those quantities, it seems to be a good tracer to constrain different feedback scenarios.In particular, we detect a characteristic plane in the 3D parameter space defined by the aforementioned properties, which is surrounded by galaxies with a reference SN feedback and a negligible impact of AGN feedback.Interestingly, the plane locally describes a surface of constant  eff , so that orthogonal departures from the plane can be directly associated to variations in  eff (i.e.different  eff isosurfaces are crossed when moving away from the plane).An increasing influence of SN or AGN feedback, generates deviations towards lower  eff .On the other hand, a weaker SN feedback can only drive positive variations of  eff for galaxies not significantly affected by AGN (i.e. with low  BH / ★ ); in other case, AGN feedback causes the decrease of  eff .
In this section, we discuss about the nature of our findings and their implications.We also explore our results in the context of observational data.

The nature of deviations from the NoAGN-fp
In this section, we try to get more insight into the nature of the residuals (, see equation 5) from the NoAGN-fp (Section 4), which, in principle, seem to be strongly dependent on the dominance of the BH inside galaxies (Fig. 7).Our aim is to explore plausible important mechanisms that can enhance the growth of the BH mass in our selected galaxies and lead to such deviations.
We address the impact of merger events by analysing the formation  1, for details).Horizontal dotted lines highlight the null distance to the NoAGN-fp.Note that for the NoAGN model,  BH is null, so no galaxies are plotted in this case.
histories of galaxies studied in the preceding sections.For this analysis, we consider mergers experienced by the main progenitor branch3 of each galaxy selected at  = 0.The results are summarised in Fig. 8, where  is plotted against the total baryonic mass (including the stellar, SF gas and NSF gas components) accumulated through merger events (i.e. tot,accr =   bar,tot, , where  considers all satellites that merge onto the main progenitor along the galaxy evolution).The upper panel corresponds to simulations which adopt different SN feedback efficiencies, while the lower panel examines the impact of implementing different temperature increases due to AGN feedback.With the exception of the StrongFB model, a consistent trend appears across the majority of models.Notably, we observe that galaxies with larger negative  have acquired a more significant amount of baryonic mass through merger events.In addition, a more close inspection suggests also that these galaxy mergers have driven the accretion of gas and BH growth in simulated systems, which reach  = 0 with more than 80 per cent of the total baryonic mass and more than 40 per cent of the total gas in the form of NSF gas.These findings are in agreement with the higher  BH / ★ shown by galaxies with  < 0: merger events can supply significant amounts of gas and potentially generate instabilities which can lead to the migration of material towards the inner galaxy regions, hence boosting the gas accretion rate onto BH and the injection of energy via AGN feedback.Moreover, mergers between BH can take place during galaxy mergers, leading to an increase in the BH mass of the remnant systems.This behaviour is not seen in the case of the StrongFB model, likely due to the lower baryonic masses reached by galaxies in this sample (see Fig. 1); merger events are expected to have had a more significant role on the formation of more massive galaxies.
We also analysed the BH accretion rates in our  = 0 galaxies,  BH,accr (z=0) , and the median of this quantity along the galaxy lifetime calculated along the main branch of progenitors of each galaxy, ⟨  BH,accr (z) ⟩. 4 The former provides information regarding the current state of the BH, while the latter gives us details about the historical value of the BH accretion rate of each galaxy (see Section 2.1.3,for further details).This information is plotted in Fig. 9.It is clear that, for all models, galaxies with higher  BH / ★ tend to show a stronger present and past average BH activity, confirming that such systems have been more significantly affected by AGN feedback.

Gas heating impact on effective yields
AGN and SN feedback can affect the gas component of galaxies through two main channels: 1) gas heating, which could quench the star formation and, hence, the chemical evolution of galaxies; 2) metal-enriched gas outflows, which deprive the galaxies from fuel for forming new stars and lower the metallicity of the systems (see De Rossi et al. 2017, for a discussion about these effects in eagle).So far, we have defined all our gas quantities considering only the SF gas component of galaxies (Section 2.2), which can vary through the two aforementioned mechanisms.But, a close inspection  1, for details).The curves represent the median relations, and the error bars denote the 25th and 75th percentiles.Symbols corresponding to the NoAGN simulation are coloured grey, since the variable  BH / ★ is null for all galaxies in this model.
of our galaxy sample in different simulations reveals that galaxies below the NoAGN-fp tend to have significantly high percentages of NSF gas (≳ 85 per cent, for the majority of galaxies).On the other hand, galaxies above the plane, which show the highest  eff , tend to have lower amounts of NSF gas.These findings suggest that feedback-driven gas heating is probably the main process responsible for  eff negative variations and deviations down the NoAGN-fp.In this context, it is interesting to explore the behaviour of galaxies if re-defining gas properties in order to include the NSF phase.Thus, we re-estimated all our galaxy properties considering the whole gas component (including SF and NSF gas) within each system.Gas heating can foster the transition from the SF to the NSF gas phase, affecting quantities derived from SF gas, while, by definition, the total gas component is not affected by such transition.Thus, the analysis of the variations of the latter component with feedback allows to probe the relevance of gas heating against that of outflows.
In Fig. 10, we analyse the 3D distribution of galaxies in the parameter space given by the total gas mass ( gas,tot =  gas,SF +  gas,NSF , where  gas,NSF stands for the NSF gas component),  ★ , and oxygen abundance derived from the total gas component (O/H) tot .Symbols are colour-coded with  eff,tot =  gas,tot /ln(1/ tot ), where  gas,tot is the metallicity corresponding to the whole gas phase and  tot =  gas,tot /( ★ +  gas,tot ).The left and right panels show results for the NoAGN and RefL50 models, respectively.We see that, when taking into account the total gas component, galaxies in the NoAGN simulation are again located around a plane (NoAGN-fptot-gas, hereafter), showing even lower dispersion around the plane than that obtained in Fig. 5 and 6.The new output parameters from the fitting of equation ( 4) to the NoAGN simulated sample and the corresponding scatter can be found in Table 2 (second row).
The right panel of Fig. 10 shows the M ★,g Z g R for the RefL50 simulation.We see that, when using the whole gas component for defining our properties, negative deviations from the plane decrease significantly.In fact, the NoAGN-fp-tot-gas seems to constitute a good representation for almost all the galaxy population even for the RefL50 model, with very few systems exhibiting very negative  1, for details), colour-coded with  eff,tot .All quantities are calculated using the total gas component.The grey shadow represents the fitting plane obtained from the NoAGN simulation.More viewing angles are available in the Supplementary Material.
departures.A more in-depth examination of the latter simulation reveals that, for galaxies exhibiting the highest  BH / ★ ratios, the median deviation from the NoAGN-fp is  ≈ −0.5 when calculating properties solely from the star-forming gas component.Conversely, when calculating properties from the whole gas component, the median departure from the NoAGN-fp-tot-gas is  ≈ −0.25, consistent with the smaller scatter around the plane (see Table 2).These results are in line with feedback-driven gas heating (which fosters the transition from the SF to the NSF gas phase) being the main responsible for the negative deviations from the NoAGN-fp (and the associated decrease of effective yields) in Fig. 5 and 6.
Finally, we emphasise that findings from this section do not imply that galaxies roughly behave as closed-boxes when taking into account their whole gas component.All simulated galaxies are affected by gas inflows, outflows and mergers during their evolution.Our results point out the importance of the criteria used for characterising the gas phase for the correct use of  eff as a tracer of feedback impact on metallicity scaling relations.

Comparison with observations
LL19 performed a detailed comparison between effective yields in eagle simulations and observations.Regarding the observational sample, they combined optical and radio wavelengths to calculate combine properties, such as gas fractions, baryonic masses and effective yields.LL19 employed spectroscopic data from the Sloan Digital Sky Survey (SDSS), HI information from the Arecibo Legacy Fast Arecibo L-band Feed Array (ALFALFA) survey, data from the GASS and COLD GASS surveys, and a sample of star forming galaxies from the Virgo cluster.For deriving metallicities, they use the approaches in Pilyugin & Grebel (2016); Pilyugin et al. (2018), applying a correction to account for oxygen locked up in dust.For more details about the estimates of observed properties, the reader is referred to LL19.In the case of simulations, LL19 studied different AGN feedback models.Additionally, these authors analysed the recalibrated high-resolution run (RecalL25) within the eagle suite, which predicts a M ★ Z g R slope that, at low  ★ , agrees better with some observational works.5However, as shown by LL19, the main general trends obtained for  eff (and its associated scaling relations) are preserved when using intermediate resolution simulations run with the reference model, as those used in this article (Section 2).
In order to compare LL19 observational data with the simulation predictions presented in this article, we followed the methodology described in LL19 to re-calculate our simulated properties in such a way that they are consistent with the instrument apertures and definitions associated with LL19 observations.Very briefly, global quantities typically measured from optical diagnostics (e.g. ★ , O/H,  gas ) were re-estimated taking into account the mass enclosed by a spherical aperture of radius 30 kpc.As discussed in Schaye et al. (2015), stellar masses obtained in this way are comparable to those derived from a projected circular aperture of the Petrosian radius.Besides, given that gas metallicities of the observed galaxy sample are inferred from HII SF regions, we evaluate chemical abundances by studying the SF gas component.Regarding gas mass calculations, we apply a larger aperture of 70 kpc (Crain et al. 2017), which roughly corresponds to the Arecibo L-Band Feed Array (ALFA) FWHM beam size of ∼ 3.5 arcmin (Giovanelli et al. 2005) at the median redshift of the GASS sample,  = 0.037 (Catinella et al. 2010).Besides, as in LL19, we quantify the simulated gas mass by using the total hydrogen mass enclosed by 70 kpc ( H,tot , hereafter).
Fig. 11 shows the M ★,g Z g R corresponding to simulations NoAGN (left-hand panel) and RefL50 (right-hand panel), obtained when using quantities calculated consistently with the LL19 observations, as explained above.Again, galaxies corresponding to the NoAGN model tend to be located around a plane, which seems to roughly describe a surface with constant  eff .In Table 2 (third row), we report the output parameters from the fitting of equation ( 4) to the latter plane (NoAGN-fp-LL19, hereafter).Once more, we obtain that larger negative (positive) deviations from the plane can be associated with lower (higher)  eff .The bulk of galaxies in the RefL50 simulation also aligns closely with the latter plane, with larger deviations observed in galaxies with the lowest  eff .Considering the scatter around the NoAGN-fp-LL19 ( ≈ 0.3, Table 2), we can predict a characteristic value  0,obs = −0.3( eff ≈ 10 −2.4 ), such that observed ing the exact value of the slope and zero point of the M ★ Z g R. The use of different selection criteria and methods for the estimates of stellar masses and metallicities affects the comparison between different observational data.  1, for details), colour-coded with  eff .All quantities are calculated using similar apertures and definitions as in LL19.The grey surface represents the fitting plane (edge-on view) obtained for the NoAGN simulation.More viewing angles are available in the Supplementary Material.galaxies with residuals  < ∼  0,obs ( eff < ∼ 10 −2.4 ) have been probably affected by a strong AGN feedback.
For the sake of comparison with observational samples constituted by SF galaxies, we re-analysed the M ★,g Z g R corresponding to the RefL50 simulation (Fig. 11, right panel) for systems with sSFR > 10 −11 [yr −1 ]. Results are shown in Fig. 12.It is clear that very large deviations from the NoAGN-fp-LL19 (edge-on grey surface) are not expected for star-forming galaxies.These trends align with the notion that such deviations are primarily generated by galaxies that are more affected by AGN feedback, as extensively discussed in previous sections.Nevertheless, we would like to notice that, although host galaxies currently showing AGN activity are usually removed when studying metallicity scaling relations, AGN feedback might have occurred in cyclic episodes that affected the progenitors of observed galaxy populations.In eagle, AGN activity also evolves with time and all BH are actively accreting gas some part of the time, giving place to the AGN phenomenon.In other words: although the accretion rate onto the BH could be low today, the BH activity could have been more significant in the past.This can be seen in Fig. 9, where it is clear that the BH accretion rates at  = 0 (  BH,accr (z=0) ) can be much lower than the median of this quantity along the galaxy lifetime (⟨  BH,accr (z) ⟩).In the particular case of galaxies shown in Fig. 12, we find that, for ≈ 50% of the sample,  BH,accr (z=0) /⟨  BH,accr (z) ⟩ < 1.
Finally, in Fig. 13, we analyse the LL19 observational sample in the 3D parameter space given by gas mass, stellar mass and O/H. 6Encouragingly, the observed galaxy sample is located around our previously defined NoAGN-fp-LL19, which is shown as a yellow surface in the figure.The good agreement between observations and eagle simulations regarding the behaviour of effective yields has already been noted by LL19.The highest (lowest) observed  eff can be associated to galaxies above (below) the plane, with the residuals  from such a plane, tracing variations in  eff .And, according to findings in this paper, the behaviour of  eff (and, hence, of ) is significantly modulated by the joint accumulated effects of SN and AGN feedback.Although disentangling the impact of the latter processes could not be plausible in many cases, the location of galaxies 6 The reader is refereed to LL19 for a description of the procedures used for deriving galaxy quantities in observations.with respect to the plane can provide some hints about their feedback histories, as we discuss in the next section.

Connecting the features of the M ★,g Z g R to the plausible feedback histories of galaxies
Our findings suggest a close connection between feedback,  eff and the location of galaxies in the 3D  ★ -O/H- gas parameter space, which could be useful for proposing feedback scenarios for real galaxy populations.As extensively discussed, simulated galaxies that were not significantly affected by AGN feedback ( BH / ★ ≈ 0), are well represented by a plane in the aforementioned 3D space, with this plane located along a region of an almost constant high  eff .Interestingly, the existence of such a plane seems to be robust against the different variations that we implemented for defining our main properties (i.e.gas mass, stellar mass and O/H), as explained in previous sections.But, our different definitions lead to moderate changes in the normalisation and orientation of the plane (Table 2).In particular, when defining the properties of simulated galaxies according to those of SF galaxies in the Local Universe, galaxies in the latter sample are located around the plane derived from simulations.
Our findings so far suggest that the location of galaxies around such characteristic plane could provide clues about the feedback histories of real galaxies, as indicated below: • Galaxies that are well represented by the plane show an almost constant high  eff and nearly null residuals .Such galaxies are not expected to have been significantly affected by AGN feedback ( BH / ★ ≈ 0).
• Considering the scatter obtained for the plane (e.g. in Table 2), galaxies with residuals  below such characteristic value can be considered more affected by the accumulated effects of AGN feedback.Such galaxies are expected to display very high  BH / ★ and very low  eff , compared with the median population located closer to the plane.
• For a reference AGN feedback efficiency, galaxies more strongly affected by SN feedback tend to lie slightly below the NoAGN-fp ( ≲ 0) and does not reach very high  eff values.Thus, galaxies with the highest  eff (those above the plane) cannot be regarded as systems strongly affected by SN feedback.
• For a reference AGN feedback efficiency but a weak SN feedback, we detect two different trends: 1) galaxies with higher  BH / ★ show lower  eff and negative distances to the plane ( < 0); 2) galaxies with lower  BH / ★ have higher  eff and tend to be located above the NoAGN-fp ( ≳ 0).Hence, galaxies with the highest  eff (those above the plane) can be associated to systems with both, a weak SN feedback efficiency and a weak AGN feedback impact.
We note that we do not have information regarding black hole masses for our observational sample.Future determination of these masses would be very useful to test the predictions of our simulations.
Finally, it is worth acknowledging that the comparison between simulations and observations could not be so straightforward.As discussed in De Rossi et al. (2017), the normalisation of the observed M ★ Z g R is still a matter of discussion, with different observational works reporting a variety of results due to the use of distinct metallicity calibrators.In addition, uncertainties in the nucleosynthetic yields implemented in simulations can also affect the determination of absolute metallicity values.This issue regarding the calibration of absolute metallicities affects also the estimate of  eff (equation 2).In this context, the comparison of global metallicities and  eff between simulations and observations should be addressed with care.

SUMMARY AND CONCLUSIONS
We have analysed the effective yields ( eff ) of galaxies in eagle cosmological hydrodynamical simulations that implement different SN and AGN feedback prescriptions: reference models calibrated against some observations (simulations RefL50 and RefL25); a model with weaker SN feedback (WeakFB simulation), and another with stronger SN feedback (StrongFB simulation); a simulation without BH and, hence, without AGN feedback (NoAGN simulation); and a simulation which applies a higher gas temperature increment (Δ AGN ) associated to AGN feedback (AGNdT9 simulation), which drives a stronger AGN feedback impact according to previous works.In particular, we try to evaluate the role of SN and AGN accumulated feedback effects on the determination of  eff and the features of its associated scaling relations.For our main analysis, galaxy properties whose definitions involve gas (e.g. chemical abundances, gas fractions and  eff ) were derived from the SF gas-phase.But, we also explored the effects of using other plausible definitions.Our most relevant findings and conclusions can be summarised as follows: • In agreement with previous works, the reference model predicts an average positive (negative) slope for the  eff −  bar relation (Fig. 1) at low  bar < ∼ 10 10 M ⊙ (high  bar > ∼ 10 10 M ⊙ ).An increase of Δ AGN generates a decrease of such a slope at high masses but has a no significant impact for low-mass systems.Interestingly, similar effects are obtained when applying a weaker SN feedback efficiency, given that it fosters the formation of a significant number of galaxies with dominant BH (i.e., with high BH mass-stellar mass ratio,  BH / ★ ).On the other hand, no dominant BH are formed in the case of a strong SN feedback scenario.In comparison with the reference model, a stronger SN feedback leads to an overall decrease of  eff at all masses and, also, predicts a positive slope for the  eff −  bar for the whole simulation mass range.
• To understand the origin of  eff behaviour, we analysed the feedback impact on the key galaxy properties involved in the  eff definition: gas-phase metallicity, stellar mass and gas mass (Fig. 2 and 3).AGN feedback seems to drive the simultaneous decrease of the SF gas fraction and metallicity of massive galaxies (which have more dominant BH), which explains their lower  eff .On the other hand, a stronger SN feedback efficiency favours the decrease of metallicity at fixed SF gas fraction, also leading to lower  eff .In the case of a weak SN feedback efficiency, a more complex behaviour arises because of the formation of BH in a significant number of galaxies along our whole mass range: as a consequence of AGN feedback, galaxies with dominant BH show lower O/H, lower SF gas fractions and, hence, lower  eff , while galaxies with low  BH / ★ tend to have higher metallicities at a given SF gas fraction and, thus, higher  eff .
• We found a clear anti-correlation between  eff and the  BH / ★ for galaxies in all studied feedback models (Fig. 4): galaxies with  BH / ★ < ∼ 10 −3.5 show a higher and roughly constant  eff while, for galaxies with  BH / ★ > ∼ 10 −3.5 ,  eff strongly decreases with  BH / ★ .
• We also studied the distribution of galaxies in the 3D parameter space defined by gas-phase metallicity,  ★ and gas mass (M ★,g Z g R) for different SN and AGN feedback models (Fig. 5 -6).Interestingly, we found that, when AGN feedback is turned off, the M ★,g Z g R is well described by a plane (NoAGN-fp), which is locally situated along a surface of constant  eff .Such a plane can roughly describe the behaviour of galaxies with no dominant BH (low  BH / ★ ) in all simulations, regardless of the implemented feedback model.Orthogonal departures from the plane () can be directly associated to variations in  eff (i.e.different  eff isosurfaces are crossed when moving away from the plane).Deviations from the plane that generates an increase of  eff are regarded as 'positive' (i.e.galaxies are considered to be above the plane) in this work, and 'negative' (i.e.galaxies are considered to be below the plane), otherwise.
• Galaxies highly affected by the accumulated effects of AGN feedback (i.e.systems with high  BH / ★ ) show very low  eff and present large negative deviations from the NoAGN-fp (Fig. 4, 7).If galaxies are not significantly affected by AGN feedback ( BH / ★ ≈ 0), they are located close to the NoAGN-fp, having an almost constant high value of  eff .
• For a reference AGN feedback efficiency, a strong SN feedback generates a moderate displacement of the bulk galaxy population downward the NoAGN-fp, with such galaxies showing intermediate  eff values.For a reference AGN feedback efficiency but a weak SN feedback, we detect two different regimes: 1) galaxies with higher  BH / ★ show lower  eff and negative deviations from the plane; 2) galaxies with lower  BH / ★ have higher  eff and tend to be located above the NoAGN-fp (Fig. 4, 7).
• A deeper analysis suggests that systems with larger negative  present higher  BH / ★ as a consequence of their more significant amount of mass accreted via galaxy mergers (Fig. 8).This is because galaxy mergers tend to boost the gas accretion rate of the BH and, also, can drive BH mergers.The present and past average BH accretion rates of galaxies with higher  BH / ★ are also higher, indicating that such systems have been more significantly affected by AGN feedback during their evolution (Fig. 9).
• Feedback-driven gas heating seems to be the main process responsible for larger  eff negative variations and higher deviations down the NoAGN-fp.Gas heating foster the transition from the SF to the non-SF (NSF) gas phase, affecting all our quantities derived from SF gas.If re-estimating these galaxy properties by using the total gas component (including SF and NSF gas), almost the whole galaxy population in all studied simulations is located around a new plane in the 3D space given by gas-phase metallicity,  ★ and gas mass (Fig. 10).The scatter around this plane is lower and no significant negative deviations from it are obtained.
• The M ★,g Z g R of observed galaxies in the Local Universe also seems to be located around a plane, with the highest (lowest) effective yields associated to galaxies above (below) the plane.If re-calculating our simulated quantities mimicking the observational definitions, the observed and simulated trends show very good agreement.In particular, the new plane derived from the NoAGN model (NoAGNfp-LL19, Fig. 11) by using such quantities represents again a surface of constant  eff and, encouragingly, this plane can characterise well the behaviour of observed galaxies with similar  eff (Fig. 13).In addition, if restricting our simulated galaxy sample to systems with sSFR > 10 −11 [yr −1 ] as in observations of star-forming galaxies, large negative  are not obtained (Fig. 12).
To sum up, according to the findings in this paper, the evolution of  eff (and, hence, of ) is significantly affected by the joint accumulated effects of SN and AGN feedback.Thus, the determination of  eff (or, equivalently, ) and its associated scaling relations could provide relevant hints regarding the different feedback histories of real galaxy populations.Nevertheless, as discussed, caution should be taken when estimating the quantities involved in the definition of  eff , for a correct interpretation of observational data.In a forthcoming article, we will extend the present study, addressing the effects of feedback on  eff at  > 0 (for preliminary results, see Zerbo et al. 2022).

Figure 1 .
Figure 1. bar −  eff relation colour-coded according to the ratio between BH and stellar mass ( BH / ★ ) for  = 0 eagle galaxies.Top (bottom) panels show results for models with different AGN (SN) feedback parameters (see Table 1, for details), as indicated in the legend.Solid lines represent the  = 0 median relation, dotted lines in middle and left panels depict the  = 0 median relation for the reference model.Error bars denote the 25th and 75th percentiles.Mass bins populated with less than 10 elements ( bin < 10) are marked with white circles.

Figure 2 .
Figure 2. Simulations that evaluate the effects of varying the AGN feedback model (see Table 1, for details).Left panels:  ★ − O/H relation.Middle panels:  ★ −  relation.Right panels:  − O/H relation.The connection between each of these relations and the effective yields is highlighted by colour-coding symbols according to the values of  eff .Symbols are also scaled to the value of  BH / ★ .Solid lines represent the  = 0 median relation, dotted lines in lower panels depict the  = 0 median relation for the reference model.Error bars denote the 25th and 75th percentiles.Mass bins populated with less than 10 elements ( bin < 10) are marked with white circles.

Figure 3 .
Figure 3. Simulations that evaluate the effects of varying supernova (SN) efficiency (see Table 1, for details).Left panels:  ★ − O/H relation.Middle panels:  ★ −  relation.Right panels:  − O/H relation.The connection between each of these relations and the effective yields is highlighted by colour-coding symbols according to the values of  eff .Symbols are also scaled to the value of  BH / ★ .Solid lines represent the z = 0 median relation, dotted lines in lower panels depict the  = 0 median relation for the reference model.Error bars denote the 25th and 75th percentiles.Mass bins populated with less than 10 elements (N bin < 10) are marked with white circles.

Figure 4 .
Figure 4.  eff vs  BH / ★ for different AGN (top panels) and SN (bottom panels) feedback models (see Table1, for details).Symbols are colour-coded according to  ★ .Solid curves correspond to the median relations and error bars depict the 25th and 75th percentiles.For the sake of comparison, the median relation corresponding to the reference model for each simulation set is shown in each panel.

Figure 5 .
Figure 5. Scatter plot of galaxies in the  ★ − O/H −  gas space for simulations that evaluate the effects of varying AGN feedback (see Table 1, for details).Left panel: 3D visualisation.Right panel: 2D projection with  = −1 and  =  2 = 2.49.Symbols are colour-coded with  BH / ★ as an indicator of the dominance of the BH mass in each galaxy.Note that, for the NoAGN simulation,  BH is null for all galaxies, so they are shown with an uniform grey colour.For more viewing angles, see the Supplementary Material.

Figure 6 .
Figure 6.Similar to Fig. 5, but for the simulations that evaluate the effects of varying SN efficiency (seeTable 1, for details).For more viewing angles, see the Supplementary Material.

Figure 7 .
Figure 7. Residuals from the NoAGN-fp, , vs log 10 (  BH / ★ ), colour-coded with  eff .Top (bottom) panels compare different AGN (SN) feedback models (see Table1, for details).Horizontal dotted lines highlight the null distance to the NoAGN-fp.Note that for the NoAGN model,  BH is null, so no galaxies are plotted in this case.

Figure 8 .
Figure 8. Residuals ( ) from the NoAGN-fp as a function of the total mass (stars + SF gas + NSF gas) accreted via merger events, colour-coded by the dominance of BH, for models that varies the efficiency of SN feedback (top panel) and AGN feedback (bottom panel) (see Table1, for details).The curves represent the median relations, and the error bars denote the 25th and 75th percentiles.Symbols corresponding to the NoAGN simulation are coloured grey, since the variable  BH / ★ is null for all galaxies in this model.

Figure 9 .
Figure9.BH accretion rate at  = 0 vs.  BH / ★ , colour-coded by the median BH accretion rate along the galaxy lifetime calculated for the main branch of each galaxy, for models that varies the efficiency of SN feedback (top panel) and AGN feedback (bottom panel) (see Table1, for details).The curves represent the median relations, and the error bars denote the 25th and 75th percentiles.No data is plotted for the NoAGN simulation due to the null values of  BH for all galaxies in this model.

Figure 10 .
Figure 10.3D scatter plot of galaxies in the  ★ − (O/H) tot −  gas,tot space for the NoAGN (left-hand panel) and RefL50 (right-hand panel) simulations (see Table1, for details), colour-coded with  eff,tot .All quantities are calculated using the total gas component.The grey shadow represents the fitting plane obtained from the NoAGN simulation.More viewing angles are available in the Supplementary Material.

Figure 11 .
Figure 11.3D scatter plot in the  ★ − O/H −  H,tot space for the NoAGN (left-hand panel) and RefL50 (right-hand panel) simulations (see Table1, for details), colour-coded with  eff .All quantities are calculated using similar apertures and definitions as in LL19.The grey surface represents the fitting plane (edge-on view) obtained for the NoAGN simulation.More viewing angles are available in the Supplementary Material.

Figure 12 .
Figure 12. 3D scatter plot in the  ★ − O/H −  gas space for the RefL50 simulation.Only galaxies with sSFR > 10 −11 [yr −1 ] are plotted.The colour bar indicates the  BH / ★ of each galaxy.All quantities are calculated using similar apertures and definitions as in LL19.The grey surface represents the best fitting plane (edge-on view) obtained for the NoAGN simulation.To explore other viewing angles, check the Supplementary Material.

Figure 13 .
Figure13.3D scatter plot in the  ★ − O/H −  gas space for the LL19 observational data, colour-coded with  eff .The yellow surface represents the best fitting plane obtained for the NoAGN simulation, when quantities are calculated using similar apertures and definitions as in LL19 (see text for details).We encourage the reader to see the Supplementary Material for more viewing angles.
Numbers in bold indicate such variations.

Table 1
, for details).For more viewing angles, see the Supplementary Material.

Table 2 .
(Robotham & Obreschkow 2015)nd  3 ; fifth to seventh columns, respectively) from the fitting of equation (4) to the NoAGN simulated sample, as indicated in the first column.The 2D plane fit of the 3D data was performed using the R package hyper-fit(Robotham & Obreschkow 2015).The eighth column indicates the intrinsic scatter orthogonal to the hyperplane ( ).The second column indicates the stellar component used in the calculations: total stellar mass or stellar mass used in LL19 (i.e. the component enclosed within a radius of 30 kpc, for simulations).The third column shows the considered gas component: total gas component, total SF gas phase or the gas component used in LL19 (i.e. total hydrogen mass enclosed within a radius of 70 kpc, for simulations).The fourth column indicates the gas component used for estimates of O/H: total gas component, total SF gas phase or that used in LL19 (i.e.SF gas within a radius of 30 kpc, for simulations).