# Synthetic nebular emission from massive galaxies – I: origin of the cosmic evolution of optical emission-line ratios

Michaela Hirschmann, Stephane Charlot, Anna Feltre, Thorsten Naab, Ena Choi, Jeremiah P. Ostriker, Rachel S. Somerville

<figure><img src="image_1.png" /></figure>

## To cite this version:

Michaela Hirschmann, Stephane Charlot, Anna Feltre, Thorsten Naab, Ena Choi, et al.. Synthetic nebular emission from massive galaxies – I: origin of the cosmic evolution of optical emission-line ratios. Monthly Notices of the Royal Astronomical Society, 2017, 472 (2), pp.2468-2495. 10.1093/mn-ras/stx2180. hal-01645688

HAL Id: hal-01645688

https://hal.science/hal-01645688v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Synthetic nebular emission from massive galaxies – I: origin of the cosmic evolution of optical emission-line ratios

Michaela Hirschmann,<fnref n="1" />★ Stephane Charlot,<fnref n="1" />★ Anna Feltre,<fnref n="1" />,<fnref n="2" />★ Thorsten Naab,<fnref n="3" /> Ena Choi,<fnref n="4" /> Jeremiah P. Ostriker<fnref n="5" />,<fnref n="6" /> and Rachel S. Somerville<fnref n="4" />,<fnref n="7" />

Accepted 2017 August 21. Received 2017 August 21; in original form 2017 May 31

## ABSTRACT

Galaxies occupy different regions of the [O III]λ5007/H β-versus-[N II]λ6584/H α emission-line ratio diagram in the distant and local Universe. We investigate the origin of this intriguing result by modelling self-consistently, for the first time, nebular emission from young stars, accreting black holes (BHs) and older, post-asymptotic giant branch (post-AGB) stellar populations in galaxy formation simulations in a full cosmological context. In post-processing, we couple new-generation nebular-emission models with high-resolution, cosmological zoom-in simulations of massive galaxies to explore which galaxy physical properties drive the redshift evolution of the optical-line ratios [O III]λ5007/H β, [N II]λ6584/H α, [S II]λλ6717, 6731/H α and [O I]λ6300/H α. The line ratios of simulated galaxies agree well with observations of both star-forming and active local Sloan Digital Sky Survey galaxies. Towards higher redshifts, at fixed galaxy stellar mass, the average [O III]/H β is predicted to increase and [N II]/H α, [S II]/H α and [O I]/H α to decrease – widely consistent with observations. At fixed stellar mass, we identify star formation history, which controls nebular emission from young stars via the ionization parameter, as the primary driver of the cosmic evolution of [O III]/H β and [N II]/H α. For [S II]/H α and [O I]/H α, this applies only to redshifts greater than z = 1.5, the evolution at lower redshift being driven in roughly equal parts by nebular emission from active galactic nuclei and post-AGB stellar populations. Instead, changes in the hardness of ionizing radiation, ionized-gas density, the prevalence of BH accretion relative to star formation and the dust-to-metal mass ratio (whose impact on the gas-phase N/O ratio we model at fixed O/H) play at most a minor role in the cosmic evolution of simulated galaxy line ratios.

**Key words:** methods: numerical – galaxies: abundances – galaxies: evolution – galaxies: formation – galaxies: general.

## 1 INTRODUCTION

The emission from ionized interstellar gas contains valuable information about the nature of the ionizing radiation and the physical conditions in the interstellar medium (ISM) in a galaxy. In fact, prominent optical emission lines are routinely used to estimate the

density, chemical abundances and dust content of the ISM and whether ionization is dominated by young massive stars (tracing the star formation rate – hereafter SFR), an active galactic nucleus (hereafter AGN) or evolved, post-asymptotic giant branch (hereafter post-AGB) stars (e.g. Izotov & Thuan 1999; Kobulnicky, Kennicutt & Pizagno 1999; Kauffmann et al. 2003; Nagao, Maiolino & Marconi 2006; Kewley & Ellison 2008; Morisset et al. 2016). These different types of ionizing sources produce distinct, well-defined correlations between the intensity ratios of strong lines, such as H α, H β, [O I]λ6300, [O II]λλ3726, 3729, [O III]λ5007, [N II]λ6584

<footnote n="1">¹Sorbonne Universités, UPMC-CNRS, UMR7095, Institut d' Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="2">²Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon, UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote n="3">³Max-Planck-Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 85741 Garching, Germany</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, Rutgers, The State University of New Jersey, New Brunswick, NJ 08854, USA</footnote>

<footnote n="5">⁵Department of Astronomy, Columbia University, New York, NY 10027, USA</footnote>

<footnote n="6">⁶Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="7">⁷Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010</footnote>

<footnote>*E-mail: hirschma@iap.fr (MH); charlot@iap.fr (SC); anna.feltre@univ-lyon1.fr (AF)</footnote>
and [S II]λλ6717, 6731 (herafter simply [O I], [O II], [O III], [N II] and [S II]). Three of the most widely used line-ratio diagnostic diagrams, originally defined by Baldwin, Phillips & Terlevich (1981, hereafter BPT) and Veilleux & Osterbrock (1987), relate the [O III]/H β ratio to the [N II]/H α, [S II]/H α and [O I]/H α ratios. These diagrams have proven useful to identify the nature of the ionizing radiation in large samples of galaxies in the local Universe (e.g. Kewley et al. 2001; Kauffmann et al. 2003).

Over the past decade, rest-frame optical spectra have become available for increasingly large samples of more distant galaxies, at redshifts $z \sim 0.5-3$, through near-infrared (NIR) spectroscopy (e.g. Pettini & Pagel 2004; Hainline et al. 2009; Steidel et al. 2014; Shapley et al. 2015), in particular with the NIR multi-object spectrographs MOSFIRE (McLean et al. 2010) and FMOS (Kimura 2010). Interestingly, all these observations indicate that star-forming (SF) galaxies at $z \gtrsim 1$ have systematically larger [O III]/H β ratio, at fixed [N II]/H α ratio, than their present-day counterparts from the Sloan Digital Sky Survey (SDSS; see e.g. Shapley et al. 2005; Lehnert et al. 2009; Yabe et al. 2012; Steidel et al. 2014; Shapley et al. 2015; Strom et al. 2017).

The physical origin of this intriguing observational feature is being heavily debated and several explanations have been proposed: high-redshift galaxies could have typically higher ionization parameters than local ones because of higher typical electron densities, higher SFRs or a higher volume-filling factor of the ionized gas, originating from a higher average gas density/pressure (e.g. Brinchmann, Pettini & Charlot 2008; Hainline et al. 2009; Lehnert et al. 2009; Steidel et al. 2014; Hayashi et al. 2015; Kashino et al. 2017), although a large electron density by itself does not seem to explain the offset in all cases (e.g. Rigby et al. 2011; Hayashi et al. 2015). Other studies appeal to the evolution of the gas-phase metallicity and an enhanced N/O abundance ratio in high-redshift galaxies (e.g. Masters 2014; Shapley et al. 2015; Masters, Faisst & Capak 2016). Instead, Steidel et al. (2014, 2016, see also Strom et al. 2017) argue that the observed offset originates primarily from a harder stellar ionizing-radiation field in distant galaxies. Additional explanations include the contribution by weak, unresolved AGN emission concurrent with stellar emission (e.g. Wright et al. 2010) and observational selection effects (Juneau et al. 2014).

For a large part, this diversity of explanations arises from the intrinsic degeneracies affecting photoionization models when adopting different prescriptions of ionizing radiation, combined with the difficulty of distinguishing evolution from selection effects when observing different samples of galaxies at different redshifts. In this context, theoretical models of the nebular emission from galaxies in a full cosmological framework could provide valuable insight into the connection between observed emission lines and the underlying ISM and ionizing-source properties as a function of cosmic time. Yet, fully self-consistent models of this kind are currently limited by the performance of cosmological radiation-hydrodynamic simulations and insufficient spatial resolution on the scales of individual ionized regions around stars and active nuclei. To circumvent these limitations, some pioneer studies proposed the post-processing of cosmological hydrodynamic simulations and semi-analytic models with photoionization models to compute the cosmic evolution of nebular emission (Kewley et al. 2013; Orsi et al. 2014; Shimizu et al. 2016). Only Kewley et al. (2013) investigate the evolution of emission-line ratios in cosmic time, combining chemical enrichment histories from cosmological hydrodynamic simulations with photoionization models of SF galaxies. These authors explore the influence of ISM conditions on the SF sequence in the [O III]/H β

versus [N II]/H α BPT diagram, as well as the potential influence of an AGN.<fnref n="1" /> Kewley et al. (2013) find that the SF sequence can be shifted to higher [O III]/H β by 'extreme' ISM conditions in high-redshift galaxies, such as large ionization parameters, high gas densities and/or hard ionizing radiation, but in unknown relative proportions. To reach more specific conclusions requires the self-consistent modelling of nebular emission from different gas components ionized by different sources in simulated galaxies.

We achieve this in this study by modelling, for the first time in a largely self-consistent way, the nebular emission from galaxies in a full cosmological context. We account for the integrated nebular emission from not only young stars (as in Orsi et al. 2014; Shimizu et al. 2016), but also AGN and post-AGB stars, based on the star formation and chemical enrichment histories of the simulated galaxies. Specifically, we post-process high-resolution, cosmological zoom-in simulations of massive galaxies with recent nebular-emission models of galaxies and AGN. The simulations include modern prescriptions for star formation, chemical enrichment, stellar feedback, black hole (hereafter BH) growth and AGN feedback (Choi et al. 2017; Núñez et al. 2017). The nebular-emission models of star-forming galaxies include improved prescriptions for the stellar ionizing radiation and a self-consistent treatment of metal depletion on to dust grains (Gutkin, Charlot & Bruzual 2016; but note that dust evolution is not followed explicitly in the simulations). We extend here these models to include the nebular emission from post-AGB stars. For the emission from AGN narrow-line regions, we appeal to the models of Feltre, Charlot & Gutkin (2016). The integrated nebular emission of a model galaxy is then the sum of the star-forming, post-AGB and AGN components. This set-up provides an ideal basis to answer the questions we wish to address in this study: can we account for the observed evolution of optical emission-line ratios, in particular the systematically larger [O III]/H β ratio of high-redshift galaxies at fixed [N II]/H α ratio? If yes, what role do the different sources of ionizing radiation and ISM properties play in the origin of this trend?

The paper is structured as follows. In Section 2, we present the general theoretical framework of this study, including the zoom-in simulations of massive galaxies, the nebular-emission models and the way in which we combine the former with the latter. Sections 3 and 4 describe our main results about the cosmic evolution of galaxies in standard optical line-ratio diagnostic diagrams and the potential physical origin of this evolution. We discuss our findings in the context of previous theoretical and observational studies and address possible caveats of our method in Section 5. Finally, Section 6 summarizes our results.

## 2 THEORETICAL FRAMEWORK

### 2.1 High-resolution simulations of massive haloes

To achieve the analysis presented in this paper, we performed a set of 20 high-resolution, cosmological zoom-in simulations of massive galaxies based on initial conditions from Oser et al. (2010, 2012), who computed the evolution in a full cosmological context of 39 galaxies with present-day halo masses between $7 \times 10^{11} M_{\odot} h^{-1}$ and $2.7 \times 10^{13} M_{\odot} h^{-1}$ ($H_0 = 100 \text{ h km s}^{-1}$). We performed these simulations with a modified version of the highly parallel,

<footnote n="1">¹ Kewley et al. (2013) consider different metal enrichments and ISM conditions in the narrow-line regions around AGN, but they do not rely on any simulation predictions for these quantities.</footnote>
smoothed particle hydrodynamics (SPH) code GADGET3 (Springel, Di Matteo & Hernquist 2005) and SPHGal (Hu et al. 2014, see also Choi et al. 2017; Núñez et al. 2017), as described in the next paragraphs. We note that our simulations differ slightly from those presented in Choi et al. (2017), particularly in the prescriptions for AGN and stellar feedback. These changes hardly affect the properties of simulated galaxies, and hence, they have a negligible impact on the results presented in this paper.

### 2.1.1 The hydrodynamic simulation code SPHGal

To overcome traditional fluid-mixing problems encountered in classical SPH codes (Agertz et al. 2007), our ‘modern’ simulation code SPHGal (Hu et al. 2014) includes a density-independent pressure-entropy SPH formulation (Ritchie & Thomas 2001; Saitoh & Makino 2013; Hopkins 2013), a Wendland C⁴ kernel with 200 neighbouring particles (Dehnen & Aly 2012), an improved artificial viscosity (Cullen & Dehnen 2010) and an artificial thermal conductivity (Read & Hayfield 2012). Moreover, to guarantee a proper treatment of shock propagation and energy feedback, a limiter of the adaptive time-step scheme of SPH ensures that neighbouring particles have similar time steps (Saitoh & Makino 2009; Durier & Dalla Vecchia 2012). For further details on these numerical schemes and their performance in test runs, we refer the reader to Hu et al. (2014).

SPHGal also follows baryonic processes, such as star formation, chemical enrichment, metal-line cooling, stellar and AGN feedback and ultraviolet photoionization background. Specifically, star formation and chemical evolution are modelled as described in Aumer et al. (2013), assuming chemical enrichment via type Ia and type II supernovae and AGB stars, with chemical yields from Woosley & Weaver (1995), Iwamoto et al. (1999) and Karakas (2010), respectively. We trace 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) in both gas and star particles. For gas particles, we include metal diffusion to allow a more realistic mixing of metals released into the ambient (possibly more metal-poor) gas. The net cooling rates are calculated from the individual element abundances, gas temperatures and densities, accounting for a redshift-dependent ultraviolet background (Haardt & Madau 2001).

Stars are assumed to form stochastically out of gas particles if the gas density exceeds a threshold value $n_{\text{th}} = n_0(T_{\text{gas}}/T_0)^3(M_0/M_{\text{gas}})^2$, where $T_{\text{gas}}$ and $M_{\text{gas}}$ are the temperature and mass of the gas particle, and $n_0 = 2 \text{ cm}^{-3}$ and $T_0 = 30\,000$ K (see Section 2.1.2). Gas particles with densities above $n_{\text{th}}$ are Jeans unstable. Their SFR is calculated as $d\rho_*/dt = \eta\rho_{\text{gas}}/t_{\text{dyn}}$, where $\rho_*$, $\rho_{\text{gas}}$ and $t_{\text{dyn}}$ are the stellar and gas densities and gas dynamical time-scale. The star formation efficiency, $\eta$, is set to a value of 0.025 reproducing the observed Schmidt-Kennicutt relation.

Star formation is regulated by both stellar and AGN feedback. We adopt the approach outlined in Núñez et al. (2017) for early stellar and SN feedback. Early feedback from young, massive stars includes ultraviolet radiative heating (within a Stroemgren sphere) and mass, energy, momentum and metal injection by stellar winds. SN feedback includes mass and metal release into the ambient gaseous medium, together with energy and momentum input during the momentum-conserving free-expansion phase of type I and type II SN blast waves ($v_{\text{out, SN}} = 4500 \text{ km s}^{-1}$).<fnref n="2" /> Mass, metals, momentum and energy from low- and intermediate-mass stars are also

transferred to surrounding gas particles in the form of slow winds ($v_{\text{out, AGB}} = 10 \text{ km s}^{-1}$), mimicking an AGB phase with mass loss. Finally, AGN feedback is tied to the prescription for BH growth. BHs are represented by collisionless sink particles, a BH seed of $10^5 M_{\odot}$ being placed at the density minimum of any dark-matter halo whose mass exceeds $10^{11} M_{\odot}$<fnref n="3" /> BHs can further grow via two channels: gas accretion and merger events with other BHs. Gas accretion is assumed to follow a statistical Bondi-Hoyle approach (Bondi 1952), whereby a gas particle is accreted on to a BH with a probability given by the volume fraction of the gas particle lying within the (unresolved) Bondi radius of the BH (e.g. Choi et al. 2012).

To compute AGN feedback from this prescription, we do not make the widely used assumption of considering *only* (spherical) thermal energy release into the ambient medium (as is the case in, e.g. the Illustris, Magneticum and EAGLE simulations; see Genel et al. 2014; Hirschmann et al. 2014; Schaye 2015, and the recent reviews by Somerville & Davé 2015; Naab & Ostriker 2017). Instead, we rely on a more physically motivated approach including both mechanical and radiative feedback (Ostriker et al. 2010; Choi et al. 2017). Specifically, we incorporate the effect of AGN-driven winds (motivated by observed broad-absorption-line winds) by randomly selecting gas particles in the vicinity of the BH (with a probability given by the feedback efficiency), which are given a velocity kick of $10\,000 \text{ km s}^{-1}$ perpendicular to the gaseous disc. Kicked particles share momentum with their two nearest neighbours, the residual energy being deposited into the gas particles as thermal energy. This allows us to roughly capture the Sedov-Taylor expansion phase of a blast wave (roughly 70 per cent in thermal energy, 30 per cent in kinetic energy). Radiative feedback from Compton and photoionization heating due to X-ray radiation from the accreting BH, radiation pressure associated with X-ray heating and the Eddington force are also included. Coupling between X-ray radiation and the surroundings follows detailed small-scale simulations by Sazonov et al. (2005). Accretion is not limited by the Eddington rate, but the Eddington force acting on electrons is self-consistently included. We refer the reader to Choi et al. (2015, 2017) for more details about AGN feedback modelling. We note that our set of 20 zoom-in simulations does not include any metallicity-dependent heating prescription. This is justified by the fact that, as shown by Choi et al. (2017), such refinements are not found to have any significant impact on basic properties of massive galaxies.

Choi et al. (2017) show how the hydrodynamic simulation described above, and in particular the sophisticated prescription for AGN feedback, can generate many realistic properties of massive galaxies, such as star formation histories, baryon conversion efficiencies, sizes, gas fractions and hot-gas X-ray luminosities. It is worth mentioning that these last two quantities are often over- or underestimated when adopting ‘traditional’ prescriptions for AGN feedback.

### 2.1.2 The simulation set-up

The dark matter haloes chosen for zoom-in re-simulations were selected from a dark matter-only N-body simulation with a comoving periodic box length $L = 100$ Mpc and $512^3$ particles (Moster et al. 2010). The cosmological parameters, based on WMAP3 measurements, are taken to be $\sigma_8 = 0.77$, $\Omega_m = 0.26$, $\Omega_{\Lambda} = 0.74$

<footnote n="2">² This is a simplified version of the full, 3-phase blast-wave model adopted in Núñez et al. (2017) and Choi et al. (2017).</footnote>

<footnote n="3">³ These halo-threshold and BH-seed masses were chosen to roughly reproduce the Magorrian et al. (1998) relation and follow theoretical calculations of BH formation by Stone et al. (2017).</footnote>
and $h = 0.72$ (see e.g. Spergel et al. 2003). The simulation was started at $z = 43$ and run to $z = 0$, with a dark-matter particle mass $M_{DM} = 2 \times 10^8 M_{\odot} h^{-1}$ and a fixed comoving gravitational softening length of $2.52 h^{-1}$ kpc. We refer the reader to the original papers of Oser et al. (2010, 2012) for more details about the simulation set-up.

From this simulation, Oser et al. (2010) selected 39 haloes with masses in the range $7 \times 10^{11} - 2.7 \times 10^{13} M_{\odot} h^{-1}$ at $z = 0$ for re-simulation. To construct the initial conditions for the high-resolution re-simulations, individual haloes are traced back in time, and all particles closer to the halo centre than twice the radius where the mean density drops below 200 times the critical density of the universe at any given snapshot are identified. These dark matter particles are replaced with dark matter and gas particles at higher resolution ($\Omega_b = 0.044$, $\Omega_{dm} = 0.216$). The new dark matter particles have a mass $m_{dm} = 2.5 \times 10^7 M_{\odot} h^{-1}$, i.e. 8 times smaller than the original ones, while the gas particles have a mass $m_{gas} = 4.2 \times 10^6 M_{\odot} h^{-1}$, equal to that of star particles. The comoving gravitational softening length of the dark matter particles is $890 h^{-1}$ pc, and that of the gas and star particles $400 h^{-1}$ pc. Here, we select for re-simulation 20 of the most massive haloes identified by Oser et al. (2010), with $z = 0$ virial masses between $3 \times 10^{12} M_{\odot}$ and $3 \times 10^{13} M_{\odot}$, and associated central galaxy masses (computed as the stellar mass within a tenth of the virial radius) between $3 \times 10^{10} M_{\odot}$ and $3 \times 10^{11} M_{\odot}$ ($h = 0.72$).

### 2.1.3 Mass assembly histories

To investigate the redshift evolution of different galaxy properties (including emission-line ratios), we construct stellar merger trees for the sample of 20 model galaxies described in the previous section. As in Oser et al. (2012), we start by using a friends-of-friends algorithm to identify, at any simulation snapshot, a central galaxy – the host (i.e. the most massive galaxy sitting at the minimum of the halo potential well) – and its surrounding satellite (less massive) galaxies. We require a minimum of 20 stellar particles (i.e. a minimum mass of about $1.2 \times 10^8 M_{\odot}$) to identify a galaxy. At $z = 2$, all galaxies in our sample are more massive than about $10^{10} M_{\odot}$, implying that, at $z < 2$, we resolve mergers down to a mass ratio of at least 1 : 100. In the analysis presented in the remainder of this paper, we trace back at every time step only the most massive progenitor of a present-day galaxy, i.e. we focus on central galaxies.

## 2.2 Modeling of nebular emission

We post-process the re-simulations of 20 galaxies presented in Section 2.1 to include nebular emission. To achieve this, we adopt the recent prescriptions of Gutkin et al. (2016) and Feltre et al. (2016) to compute the nebular emission arising from young massive stars (Section 2.2.1) and narrow-line regions of AGN (Section 2.2.2). We also design a prescription to account for the contribution by post-AGB stars to nebular emission (Section 2.2.3). All emission-line models presented in this paper were computed using version c13.03 of the photoionization code CLOUDY (Ferland et al. 2013), always accounting for the interaction of photons, electrons and atomic ions with dust grains in H II regions.

### 2.2.1 Nebular emission from star-forming galaxies

We adopt the grid of nebular-emission models of star-forming galaxies computed by Gutkin et al. (2016). These calculations combine the latest version of the Bruzual & Charlot (2003) stellar population

synthesis model (Charlot & Bruzual, in preparation) with CLOUDY, following the method outlined by Charlot & Longhetti (2001). In brief, the approach consists in convolving the spectral evolution of a typical, ionization-bounded H II region powered by a new-born star cluster with a star formation history, to compute the nebular emission of a whole galaxy. In this context, the parameters of the photoionization model of this typical H II region should be interpreted as effective (i.e. galaxy-wide) ones, describing the ensemble of H II regions and the diffuse gas ionized by stars throughout the galaxy (see section 2.3 of Charlot & Longhetti 2001 and section 2.2 of Gutkin et al. 2016 for details). A ‘closed geometry’ is used in CLOUDY to perform these calculations, as appropriate for spherical H II regions.

The grid computed by Gutkin et al. (2016, see their table 3) encompasses models in wide ranges of interstellar (i.e. gas+dust-phase) metallicities, $Z_*$, ionization parameters, $U_*$, dust-to-metal mass ratios, $\xi_d$, H II-region densities, $n_{H,*}$ and carbon-to-oxygen abundance ratios, (C/O)$_*$.

A main feature of these models (which incorporate secondary nitrogen production) is the self-consistent treatment of metal depletion on to dust grains, which allows one to relate the gas-phase metallicity measured from nebular emission to the total interstellar metallicity $Z_*$. The ionization parameter is defined as the dimensionless ratio of the number density of H ionizing photons to that of hydrogen. Gutkin et al. (2016) adopt the same metallicity for the ionizing stars as for the ISM. This is consistent with our simulations, in which the metallicity of newly born stars is (by construction) very similar to that of the warm/cold gas. The parameters of this model grid are summarized in Table 1.<fnref n="4" /> It is worth noting that, as pointed out by Charlot & Longhetti (2001, see also Gutkin et al. 2016), the actual mass of the effective star cluster (namely the ionization photon rate, Q) has no influence on the predicted nebular emission at fixed $U_*$ and $n_{H,*}$, due to a degeneracy between this mass and the gas filling factor entering the definition of the ionization parameter (see Section 2.3.1 below). We adopt here the emission-line predictions of Gutkin et al. (2016) for 10-Myr-old stellar populations with constant SFR (sufficient to reach a steady population of H II regions) and a standard Chabrier (2003) IMF (consistent with the IMF adopted in the simulations), truncated at 0.1 and 100 $M_{\odot}$. As we shall see in Section 5, increasing the upper mass cut-off of the IMF to 300 $M_{\odot}$ hardly affects our results.

### 2.2.2 Nebular emission from AGN

For the narrow-line regions of AGN, we adopt the grid of nebular-emission models of Feltre et al. (2016). In this prescription, the spectrum of an AGN is approximated by a broken power law of adjustable index $\alpha$ in the frequency range of ionizing photons (equation 5 of Feltre et al. 2016). In the narrow-line region, gas is assumed to be distributed in clouds of a single type. These models include dust and radiation pressure. An ‘open geometry’ is used in CLOUDY, as appropriate for gas with a small covering factor.

The grid of AGN nebular-emission models is parametrized in terms of the interstellar metallicity in the narrow-line region, $Z_*$, the ionization parameter of this gas, $U_*$, the dust-to-metal mass ratio, $\xi_d$, the density of gas clouds, $n_{H,*}$, and the carbon-to-oxygen abundance ratio (C/O)$_*$ (see Table 1 for details). We note that the

<footnote n="4">⁴ The values of the ionization parameter reported in Table 1 are a factor of 9/4 larger than those in table 3 of Gutkin et al. (2016). This is because we label models here in terms of the volume-averaged ionization parameter (see equation 1 of Section 2.3.1).</footnote>
<table><caption>**Table 1.** Overview of the parameter space of the nebular-emission models for young stars, AGN and post-AGB stars. To select the SF, AGN and PAGB models appropriate for each galaxy at each simulation time step, we adopt a fixed spectral slope of AGN ionizing radiation ($\alpha = -1.7$), fixed dust-to-metal mass ratio ($\xi_d = 0.3$) and fixed ionized-gas density ($n_{H,*} = 10^2 \text{ cm}^{-3}$, $n_{H,\bullet} = 10^3 \text{ cm}^{-3}$ and $n_{H,\odot} = 10 \text{ cm}^{-3}$). We further match the AGN, SF and PAGB ionization parameters, interstellar (i.e. gas+dust-phase) metallicity, C/O abundance ratio and age and metallicity of post-AGB stars to those of the simulated galaxy. The emission-line luminosities are scaled to the SFR, AGN luminosity and mass of post-AGB stellar population of the galaxy.</caption><thead><tr><th>Parameter space</th><th>SF models (Gutkin et al. 2016)</th><th>AGN models (Feltre et al. 2016)</th><th>PAGB models (this work)</th></tr></thead><tbody><tr><td><strong>Ionizing spectrum</strong><br>(matched/fixed)</td><td><strong>10-Myr-old stellar</strong><br>population with <strong>const</strong><br><strong>SFR(fixed)</strong>, stellar metallicity<br>same as that of gas (matched)</td><td>UV slope α = -1.2, -1.4, -1.7,<br>-2.0 (fixed)</td><td>3-, 5-, 7-, 9-Gyr-old stellar<br>populations (matched) Z<sub>☉</sub>, stars =<br>0.008, 0.014, 0.017, 0.02 (matched)</td></tr><tr><td><strong>Interstellar metallicity Z</strong><br>(matched)</td><td>Z<sub>☉</sub> = 0.0001, 0.0002, 0.0005,<br>0.001, 0.002, 0.004, 0.006,<br>0.008, 0.014, 0.017, 0.02,<br>0.03</td><td>Z<sub>☉</sub> = 0.0001, 0.0002, 0.0005,<br>0.001, 0.002, 0.004, 0.006, 0.008,<br>0.014, 0.017, 0.02, 0.03, 0.04,<br>0.05, 0.06, 0.07</td><td>Z<sub>☉</sub> = 0.0001, 0.0002, 0.0005, 0.001,<br>0.002, 0.004, 0.006, 0.008, 0.014,<br>0.017, 0.02, 0.03, 0.04, 0.05, 0.06,<br>0.07</td></tr><tr><td><strong>Ionization parameter</strong><br>log U, function of the<br>average gas density<br>(matched)</td><td>log U<sub>*</sub> = -0.65, -1.15,<br>-1.65, -2.15, -2.65, -3.15,<br>-3.65</td><td>log U<sub>•</sub> = -0.65, -1.15, -1.65,<br>-2.15, -2.65, -3.15, -3.65,<br>-4.65</td><td>log U<sub>☉</sub> = -2.15, -2.65, -3.15,<br>-3.65, -4.15, -4.65</td></tr><tr><td><strong>Dust/metal mass ratio</strong><br>ξ<sub>d</sub> (fixed)</td><td>0.1, 0.3, 0.5</td><td>0.1, 0.3, 0.5</td><td>0.1, 0.3, 0.5</td></tr><tr><td><strong>Ionized-gas density</strong><br>log (n<sub>H</sub>/cm<sup>3</sup>) (fixed)</td><td>log (n<sub>H,*</sub>) = 2.0, 3.0, 4.0</td><td>log (n<sub>H, •</sub>) = 2.0, 3.0, 4.0</td><td>log (n<sub>H, ☉</sub>) = 1.0, 2.0, 3.0</td></tr><tr><td><strong>C/O abundance ratio in</strong><br>solar units (matched)</td><td>(C/O)<sub>*</sub>/(C/O)<sub>☉</sub> = 0.1, 0.2,<br>0.27, 0.38, 0.52, 0.72, 1.0</td><td>(C/O)<sub>•</sub>/(C/O)<sub>☉</sub> = 0.1, 0.2, 0.27,<br>0.38, 0.52, 0.72, 1.0</td><td>(C/O)<sub>☉</sub>/(C/O)<sub>☉</sub> = 1.0</td></tr><tr><td><strong>Model normalization</strong><br>(matched)</td><td>Star formation rate SFR</td><td>AGN luminosity L<sub>AGN</sub></td><td>Mass of evolved stars M<sub>☉, stellar</sub></td></tr></tbody></table>

introduction of non-solar (C/O) ratios is an improvement over Feltre et al. (2016), who considered only models with the solar value. This refinement hardly affects the predictions for optical emission-line ratios presented in this paper, but it has a major influence on the predicted ultraviolet emission-line ratios (Hirschmann et al. in preparation).

### 2.2.3 Nebular emission from post-AGB stars

To describe the nebular emission from quiescent, passively evolving galaxies, we must also account for the emission from diffuse gas heated by evolved, post-AGB stars. We build a grid of such ‘PAGB’ models using spectra of single-age, evolved stellar populations (computed with the same version of the Bruzual & Charlot 2003 stellar population synthesis code as used by Gutkin et al. 2016) as input to the photoionization code CLOUDY. We compute models for a range of stellar populations ages (between 3 and 9 Gyr) and metallicities, Z<sub>☉, stars</sub> (see Table 1). In contrast to the SF models (Section 2.2.1), we do not impose the interstellar metallicity to be the same as that of the ionizing stars in the PAGB models. Instead, we treat it as an independent parameter, Z<sub>☉</sub>. This is motivated by the fact that old stars do not necessarily have the same metallicity as the diffuse ISM in which they evolve (we are interested here in the emission from diffuse gas ionized by all hot post-AGB stars, including beyond the short-lived planetary-nebula phase). In fact, our simulations show that, particularly at redshifts z < 2, old stellar populations can be significantly more metal-rich than cold/warm gas, whose metallicity is often diluted by late infall of metal-poor gas. The other model parameters controlling the nebular emission from post-AGB stars, summarized in Table 1, are consistent with those of the SF models: the gas ionization parameter, U<sub>☉</sub>, dust-to-metal mass ratio, ξ<sub>d</sub>, and hydrogen density, n<sub>H, ☉</sub>.

## 2.3 Coupling nebular-emission models with zoom-in simulations

We couple the extensive grid of nebular-emission models described in Section 2.2, with the simulations of massive galaxies described in Section 2.1, by selecting an SF, an AGN and a PAGB emission-line model for each simulated galaxy at each redshift step.<fnref n="5" /> The sum of these three components makes up the integrated nebular emission of a model galaxy. In practice, we select the SF/AGN/PAGB models appropriate for each galaxy by self-consistently matching all model parameters possibly available from the simulations (e.g. metallicity of the star-forming gas). The quantities derived from the simulated galaxies are calculated within 1/10 of the halo virial radius for the SF and PAGB models, and within 1 kpc of the BH for the AGN models. Those parameters that cannot be retrieved from the simulation are set to standard values (indicated in bold characters in Table 1). This is the case for, e.g. the dust-to-metal mass ratio, ξ<sub>d</sub>, and the hydrogen gas density in individual ionized regions, n<sub>H</sub>, since we do not model dust physics in the simulation and also cannot resolve individual H II regions, AGN narrow-line regions and PAGB-star environments. We adopt ξ<sub>d</sub> = 0.3 in all SF, AGN and PAGB models of nebular emission in this paper. This value, close to that of ξ<sub>d</sub>, ☉ = 0.36 in the Solar neighbourhood (Gutkin et al. 2016), lies in the middle of the parameter range. Moreover, we adopt n<sub>H,*</sub> = 10<sup>2</sup> cm<sup>-3</sup> and n<sub>H, •</sub> = 10<sup>3</sup> cm<sup>-3</sup>, as typical gas densities estimated from optical line-doublet analyses of H II regions and AGN (see e.g. sections 5.6

<footnote n="5">⁵ Note that in this study, we select only one emission-line model per galaxy since we are focusing on integrated spectral properties. Future studies will, instead, focus on spatially resolved emission properties requiring a more refined coupling procedure.</footnote>
and 13.4 of Osterbrock & Ferland 2006).<fnref n="6" /> For the PAGB models, we adopt $n_{H, \odot} = 10^{-3}$ cm$^{-3}$, which should be more appropriate for the diffuse environment of old, post-AGB stars in the ISM of mature galaxies (adopting $n_{H, \odot} = 10^{-3}$ cm$^{-3}$ would not change significantly our results). The potential impact of adopting different values for these fixed parameters is discussed in Section 5.

In the next paragraphs, we describe in more detail the way in which we couple the SF, AGN and PAGB nebular models with our zoom-in galaxy simulations.

### 2.3.1 Matching SF models to simulated galaxies

With each galaxy at each simulation time step, we associate the SF emission-line model from the Gutkin et al. (2016) grid described in Section 2.2.1 with closest star and gas parameters. As mentioned earlier, the parameters of the Gutkin et al. (2016) models are effective (i.e. galaxy-wide) ones, describing the ensemble of H II regions and the diffuse gas ionized by young stars throughout the galaxy. We select the grid metallicity $Z_*$ and carbon-to-oxygen ratio (C/O), closest to the simulated global (i.e. galaxy-wide) metallicity $Z_{gas, glob}$ and abundance ratio (C/O)$_{gas, glob}$ of the warm-gas phase. Here, we consider gas particles with temperatures between 10 000 K, typical of giant H II regions and narrow-line regions of AGN, and 100 K, the lowest temperature achievable via atomic gas cooling in the simulations (the remainder of the gas in the simulations is in the diffuse, hot ionized phase). We compute the volume-averaged ionization parameter of the simulated galaxy as (see e.g. equation B.6 of Panuzzo et al. 2003)<fnref n="7" />

$$U_{\text{sim},*} = \frac{3\alpha_{\text{B}}^{2/3}}{4c} \left( \frac{3Q_{\text{sim},*} \epsilon^2 n_{\text{H},*}}{4\pi} \right)^{1/3} . \qquad (1)$$

Here $Q_{\text{sim},*}$ is the rate of ionizing photons (obtained by multiplying the SFR of the simulated galaxy by the rate output by a 10-Myr-old stellar population with unit SFR and metallicity $Z_*$), $n_{H,*} = 10^2 \text{ cm}^{-3}$ (see above), $\alpha_{\text{B}}$ is the case-B hydrogen recombination coefficient and $\epsilon$ is the volume-filling factor of the gas, defined by

$$\epsilon = \frac{n_{\text{gas, glob}}}{n_{\text{H,*}}}, \qquad (2)$$

where $n_{\text{gas, glob}}$ is the volume-averaged, global (hydrogen) gas density, again considering gas particles with temperatures between 10 000 K and 100 K. In the rare cases where the volume-averaged gas density exceeds the adopted hydrogen density in the ionized regions, i.e. $\epsilon > 1$, we set the filling factor to unity, i.e. $\epsilon = 1$. Note that a larger $n_{\text{gas, glob}}$ at fixed $n_{H,*}$ implies a larger filling factor $\epsilon$, i.e. a more compact arrangement of individual gas clumps around the ionizing source, and hence, a larger incident flux of ionizing photons per unit gas area. In our approach, therefore, $U_{\text{sim},*}$ depends on the simulated SFR via $Q_{\text{sim},*}$. This is justified by the observed relation between ionization parameter and metallicity (fig. 2 of Carton et al. 2017), together with that between metallicity and specific SFR (fig. 7 of Mannucci et al. 2010), for SDSS star-forming galaxies, which imply a positive correlation between $\log U_*$ and SFR at

fixed galaxy stellar mass. Combined with the filling factor derived from the simulation (equation 2), this uniquely defines the quantity $Q_{\text{sim},*} \epsilon^2$ entering the definition of $U_{\text{sim},*}$ at fixed $n_{H,*}$ (equation 1). We select the SF model with ionization parameter $\log U_*$ closest to $\log U_{\text{sim},*}$ (computed by Gutkin et al. 2016 for an arbitrary combination of effective star-cluster mass and gas filling factor yielding the same $Q_* \epsilon^2$; see Section 2.2.1). This uniquely defines the Gutkin et al. (2016) model associated with each simulated galaxy at each time step.

### 2.3.2 Matching AGN models to simulated galaxies

We adopt a procedure similar to that outlined in the previous section to associate nuclear activity of a galaxy at any simulation time step with an AGN emission-line model from the Feltre et al. (2016) grid described in Section 2.2.2. The ISM conditions for the AGN model are taken to be the central (and not global) ones of the simulated galaxy, i.e. in the vicinity of the BH. Specifically, we select gas particles in a comoving sphere of 1-kpc radius around the BH to compute the central warm-gas metallicity, $Z_{\text{gas}, 1\text{kpc}}$, central carbon-to-oxygen ratio, (C/O)$_{\text{gas}, 1\text{kpc}}$, and central volume-averaged gas density, $n_{\text{gas}, 1\text{kpc}}$ (we have checked that adopting radii in the range 0.4–3 kpc instead of 1 kpc hardly affects our results). This should be roughly appropriate to probe the narrow-line regions around AGN with luminosities in the range found in our simulations (see fig. 3 of Hainline et al. 2014 and the model AGN luminosities in Fig. 6 below). As before (equation 2), a larger $n_{\text{gas}, 1\text{kpc}}$ at fixed $n_{H,*}$ implies a larger volume-filling factor of the narrow-line region. We compute the AGN luminosity from the simulated BH accretion rate (BHAR) as in Hirschmann et al. (2014). For the AGN ionizing spectrum, we adopt a fixed ultraviolet slope $\alpha = -1.7$ of the flux per unit frequency ($\propto \nu^{\alpha}$; indicated in bold face in Table 1) and an amplitude scaled to the simulated AGN luminosity. This allows us to compute the rate of ionizing photons produced by the AGN, noted $Q_{\text{sim},*}$, and the corresponding ionization parameter, noted $U_{\text{sim},*}$, via equation (1). Then, we select the Feltre et al. (2016) model with closest $Z_*$, $\log U_*$ and (C/O)$_*$.

### 2.3.3 Matching PAGB models to simulated galaxies

To select a PAGB emission-line model from the grid presented in Section 2.2.3 for each galaxy at each simulation time step, we start by computing the average age and metallicity of all star particles older than 3 Gyr (when a significant population of post-AGB stars starts to build up in a simple stellar population). We identify the available grid age and metallicity ($Z_*$, stars) closest to these values in Table 1 and compute the rate of ionizing photons from PAGB stars, noted $Q_{\text{sim},*}$, based on the mass in stars older than 3 Gyr in the simulation. Since our simulation does not allow us to distinguish gas in star-forming clouds from that in the diffuse ISM (which would require to resolve scales of a few tens of parsecs), for the PAGB model we adopt the same global interstellar metallicity $Z_{\text{gas, glob}}$, abundance ratio (C/O)$_{\text{gas, glob}}$ and volume-averaged gas density $n_{\text{gas, glob}}$ as for the SF model above. Assuming $n_{H,*} = 10^{-3}$ cm$^{-3}$ accounts for the fact that the gas seen by old stars is more diffuse (equation 2) than that seen by young stars ($n_{H,*} = 10^2 \text{ cm}^{-3}$). Then, using equation (1), we compute the ionization parameter of the gas ionized by post-AGB stars, noted $U_{\text{sim},*}$, and select the model with closest $Z_*$, $\log U_*$ and (C/O)$_*$ in Table 1.

<footnote n="6">⁶ Densities measured from emission-line doublets are those of the emitting gas clumps.</footnote>

<footnote n="7">⁷ The definition of the volume-averaged ionization parameter in equation (1) differs from that of the quantity $\langle U \rangle$ in equation (7) of Charlot & Longhetti (2001) by a factor of 3/4, and from that of $U_S$ in equation (7) of Gutkin et al. (2016) by a factor of 9/4. These different model-labelling choices have no influence on the actual CLOUDY calculations, for which the input parameter is the rate of ionizing photons (equation 8 of Charlot & Longhetti 2001).</footnote>
### 2.3.4 Total emission-line luminosities and line ratios

The procedure described in the previous paragraphs allows us to compute the contributions of young stars, an AGN and post-AGB stars to the luminosities of various emission lines (such as $L_{H\alpha}$, $L_{H\beta}$, $L_{OIII}$, etc.) in a simulated galaxy. The total emission-line luminosities of the galaxy can then be calculated by summing over these three contributions. For line luminosity ratios, we adopt for simplicity the notation $L_{OIII}/L_{H\beta} = [O\,III]/H\beta$. In this study, we focus on exploring four line ratios in the optical regime: $[O\,III]\lambda5007/H\beta$, $[N\,II]\lambda6584/H\alpha$, $[S\,II]\lambda6717$, 6731/H$\alpha$ and $[O\,I]\lambda6300/H\alpha$. The strengths of metal lines reflect a combination of photoionization processes of the corresponding elements and excitation processes of the resulting ions via collisions with photo-electrons. We note that, in this paper, we do not consider attenuation by dust outside H\,II regions and compare our predictions with observed emission-line ratios corrected for this effect (if not provided by the original studies, we applied a correction based on the $H\alpha/H\beta$ ratio and the Calzetti et al. 2000 attenuation curve). By design, the above optical-line ratios are anyway little sensitive to dust, as they are defined by lines close in wavelength (for reference, the corrections are less than ~0.015 dex for a V-band attenuation of $A_V \sim 1$ mag).

## 3 COSMIC EVOLUTION OF OPTICAL EMISSION-LINE RATIOS

Knowing the most massive (i.e. main) progenitor of every (central) galaxy at any simulation time step allows us to explore the cosmic evolution of emission-line ratios. We start by investigating the dependence of optical emission-line ratios on the evolution of stars, BHs and the ISM in three galaxies illustrating different physical properties (Section 3.1). Then, we explore the predicted evolution of optical emission-line ratios for the full set of 20 simulated massive galaxies described in Section 2.1 and compare this with observations at various redshifts (Section 3.2). Our sample of zoom-in simulations of massive, mostly quiescent, present-day galaxies is well-suited for this analysis, since it allows us to probe a mass range of star-forming galaxies at high redshift (corresponding to the main progenitors) similar to that sampled by available observations (Section 3.2.6). A natural consequence of following the evolution of the main progenitors of present-day massive galaxies is that the average galaxy mass increases from high to low redshift in our simulated sample. To understand how this feature can affect our results, we also investigate emission-line evolutionary trends for galaxies in a fixed stellar mass range at all redshift.

### 3.1 Three case studies

Fig. 1 illustrates the evolution of various quantities pertaining to stars, BHs and the ISM, for three galaxies of our simulation set, M0175, M0190 and M0215 (from left to right). These quantities are: mass assembly history (first row); SFR and BH accretion rate (second row); global and central warm-gas metallicities (third row); global and central C/O ratios (fourth row); and global and central volume-averaged gas densities (fifth row). The sixth row shows, as an example, the redshift evolution of two emission-line ratios, $[O\,III]/H\beta$ (thick line) and $[N\,II]/H\alpha$ (thin line). In the last row, we show the evolution of the same models in the $[O\,III]/H\beta$ versus $[N\,II]/H\alpha$ diagram (colour-coded by redshift as in the sixth row). Different symbols refer to three ranges of BAHR-to-SFR ratio, which we associate to: SF galaxies, with $BHAR/SFR < 10^{-4}$ (stars); composite galaxies, with $10^{-4} \le BHAR/SFR \le 10^{-2}$

(circles); and active galaxies, with BHAR/SFR > $10^{-2}$ (triangles). In addition, galaxies in which $H\beta$ emission from post-AGB stars exceeds that from both AGN and young stars are indicated by squares. Such galaxies are thought to constitute a sub-group of the classical population of low-ionization (nuclear) emission-line regions (LINER/LIER; e.g. Singh 2013; Belfiore et al. 2016). For reference, we also indicate the location of local ($z \sim 0.1$) SDSS galaxies in these line-ratio diagrams (grey shaded areas and contours), together with standard observational criteria to distinguish SF galaxies (below the dashed line) from composites (between the dashed and dotted lines), AGN (above the dotted line) and LI(N)ER (in the bottom-right quadrant defined by dot-dashed lines), according to Kewley et al. (2001, dotted line) and Kauffmann et al. (2003, dashed and dot-dashed lines). We now describe in more detail our findings for this illustrative sample of three galaxies.

#### 3.1.1 Galaxy M0175

Galaxy M0175 (the left-most column in Fig. 1) has a fairly quiet mass assembly history dominated by 'in situ' star formation without any major merger (top panel). Yet, around $z=2$, an infall of gas clouds drives a peak in BH accretion (second panel). AGN feedback from this event stops further BH accretion and also makes star formation drop on a longer time-scale. Only at low redshift does new gas infall (not merger-driven) trigger a second peak in BH accretion. These peaks in BH accretion induce peaks in the otherwise gradually rising $[N\,II]/H\alpha$ ratio. This is because the harder ionizing radiation of an AGN relative to stars makes the electronic temperature higher, and hence, the collisionally excited $[N\,II]$ line stronger (see fig. 1 of Feltre et al. 2016). During these events, the galaxy makes excursions in the regions of the $[O\,III]/H\beta$ versus $[N\,II]/H\alpha$ diagram populated by AGN and composites (blue and orange triangles in the bottom panel). We note that, unlike the $[N\,II]/H\alpha$ ratio, the $[O\,III]/H\beta$ ratio starts to decline below $z=3$, and drops further by almost an order of magnitude around $z \sim 1$. This behaviour results from the combination of two factors: the drop in SFR at additionally slightly decreasing global average gas density, which makes the ionization parameter also drop (equations 1–2); and the rise in metallicity at nearly constant C/O ratio. This rise makes cooling more efficient in H\,II regions, reducing the electronic temperature and in turn the luminosity of $[O\,III]\lambda5007$, whose excitation potential is higher than that of $[N\,II]\lambda6584$ (Stasińska 1980; Gutkin et al. 2016). We note that the rise in $[N\,II]/H\alpha$ towards low redshift is boosted by the inclusion of secondary nitrogen production in the nebular models (Section 2.2.1). The drop in average central gas density towards low redshift further contributes to the above trends in $[O\,III]/H\beta$ and $[N\,II]/H\alpha$ by making the AGN ionizing parameter smaller.

As a result of these trends, galaxy M1075 'moves' from the top-left to bottom part of the SF branch in the $[O\,III]/H\beta$ versus $[N\,II]/H\alpha$ diagram in Fig. 1 (area below the dashed line, along the observed sequence of SDSS star-forming galaxies). At later times, when the radiation from young stars and the AGN is very weak, post-AGB stars become the dominant source of ionizing photons, causing the galaxy to appear as a LI(N)ER (yellow and green squares). At $z=0$, M1075 disappears from line-ratio diagram because no warm gas is left.

#### 3.1.2 Galaxy M0190

The mass assembly history of galaxy M0190 (top panel, middle column in Fig. 1) is very different from that of M0175. Galaxy M0190
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Redshift evolution of different properties of three re-simulated galaxies with (middle column) and without (left and right columns) major-merger events: stellar mass (first row); BH accretion rate (green dashed lines, second row); SFR (blue solid lines, second row); global (blue solid) and central (green dashed) interstellar metallicity (third row); global (blue solid) and central (green dashed) C/O ratio (fourth row); and global and central gas density (fifth row). The sixth and seventh rows show the implied redshift evolution of the integrated [O iii]/H β and [N ii]/H α ratios and the corresponding [O iii]/H β versus [N ii]/H α diagram (where symbols are colour-coded according to redshift as along the x-axis of the sixth row). Different symbols refer to different galaxy types (stars: SF; circles: composite; triangles: AGN-dominated; and squares: post-AGB-dominated galaxies; see Section 3.1 for details). Also shown for reference are observations of local SDSS galaxies (grey shaded areas and contours), together with standard observational criteria to distinguish SF galaxies (below the dashed line) from composites (between the dashed and dotted lines), AGN (above the dotted line) and LI(N)ER (in the bottom-right quadrant defined by dot-dashed lines), according to Kewley et al. (2001, dotted line) and Kauffmann et al. (2003, dashed and dot-dashed lines).</figcaption></figure>
experiences three major mergers with mass ratio above 1:4 – one at $z = 1.8$, one at $z = 0.7$ and one at $z = 0.2$ (indicated by orange dashed lines), through which roughly 60 per cent of the final galaxy stellar mass is accreted. Interestingly, the BH accretion history (second panel) reveals peaks of AGN activity associated with the two mergers at $z < 1$, but not the high-redshift one. This is presumably because a turbulent environment with torques and radial gas flows at high redshift maintains the BHAR at a high fraction of the Eddington rate. At $z < 1$, instead, spikes of AGN activity (resulting in AGN luminosities of $10^{45}-10^{46}$ erg s⁻¹) can be triggered even by small events during a ‘smooth’ accretion phase. The impact of AGN feedback on star formation for M0190 is also moderated by the repeated inflow of gas, which maintains the SFR above $1\ M_{\odot}\ yr^{-1}$ over most of the galaxy’s history. Only after the last major merger at $z = 0.2$ do both star formation and BH accretion shut down, as a consequence of feedback-driven gas heating and outflows. High levels of (pristine) gas accretion also affect the central and global interstellar metallicities, as can be seen from the dip in metallicity evolution between $z \gtrsim 1$ and $z = 0.2$ (third panel).

The above histories of star formation, BH accretion and chemical enrichment have consequences for emission-line properties: despite several pronounced nuclear-activity peaks at $0.2 < z < 1$, M0190 remains in the star-forming and composite regions of the [O III]/H β versus [N II]/H α diagram without strong excursions in the AGN region during this evolutionary phase (green and turquoise triangles in the bottom panel). This is primarily because sustained star formation (and hence, high ionization parameter) and low metallicity contribute to maintaining high [O III]/H β and low [N II]/H α at redshifts down to $z \sim 0.2$ (Section 3.1.1; see also fig. 2 of Feltre et al. 2016). Only below $z = 0.2$, when star formation drops and metallicity rises, can the galaxy move to the AGN region of the diagram (yellow symbols).

### 3.1.3 Galaxy M0215

As a third example, galaxy M0215 (top panel, right-most column in Fig. 1) exhibits a smooth mass-assembly history with only three minor mergers at redshifts $z < 2$. Star formation remains elevated, at a rate of 10–100 $M_{\odot}\ yr^{-1}$, down $z = 0.7$. BH accretion is also smooth down to $z = 0.3$, when spikes of AGN activity set in (second panel). Below $z = 0.7$, the SFR starts to decline as new gas supply does not compensate consumption through star formation and BH accretion, until AGN-driven winds shut down entirely star formation at $z < 0.2$. The lack of new supply of metal-poor gas, together with the cooling of a hot enriched halo, also makes the average interstellar metallicity rise sharply from sub-solar to super-solar values at $z < 0.7$ (third panel). This sudden change in metallicity triggers a sharp drop by nearly an order of magnitude in [O III]/H β ratio and a rise in [N II]/H α ratio, which translate into a ‘bi-modal’ feature in the corresponding line-ratio diagram (bottom panel). During the SF phase at $z \gtrsim 1$, the galaxy remains on the top left of the SF branch, with $log([O III]/H\beta) > 0.8$ and $log([N II]/H\alpha) < -1$, while around $z = 0.7$, it jumps from the SF-galaxy to composite regions in the bottom right area, with $log([O III]/H\beta) < 0.5$ and $log([N II]/H\alpha) > -1$. Later on, enhancement in [N II]/H α from nuclear activity (Section 3.1.1) triggers excursions of the galaxy in the AGN region of the diagram. This illustrates how BH accretion causing low-redshift galaxies to populate the AGN region of the [O III]/H β versus [N II]/H α diagram does not necessarily require merger events (Li et al. 2008). At $z < 0.2$, after star formation and BH accretion have been suppressed, post-AGB stars take over

the production of ionizing photons, making the galaxy appear as a LI(N)ER (red squares).

In summary, therefore, a generic feature of all three examples of massive-galaxy evolution shown above is the drop in [O III]/H β and the rise in [N II]/H α from high to low redshift. These general trends arise from a combination of overall increasing metallicity and decreasing SFR. Peaks in the BHAR typically trigger peaks in [N II]/H α, provided that the (central) metallicity is large enough and not diluted by infall of (metal-poor) gas. Around these trends, the exact mass assembly and merger history of a galaxy can strongly affect nebular emission on a case-by-case basis, which leads to substantial scatter in the predicted line ratios at given redshift and mass, complicating the interpretation. The C/O abundance ratio and average gas density, instead, appear to have a minor influence on the evolution of [O III]/H β and [N II]/H α. In Section 4 below, we explore in more detail the potential role of various physical parameters in driving the observed evolution of optical-line ratios in galaxies at different redshifts.

## 3.2 Population study

We now investigate the evolution of the full set of 20 zoom-in simulations of massive galaxies and their main progenitors in various optical emission-line ratio diagrams at different redshifts. We consider redshift bins including several simulation snapshots, and hence, potentially several tens of emission-line galaxies. As noted earlier in this section, such a sample is naturally characterized by an increase in average galaxy stellar mass from high to low redshift, accompanied by an increase in metallicity. In what follows, we start by investigating the properties of the full set of simulated galaxies in optical line-ratio diagrams at different redshifts, including the evolution of the SF sequence and of average emission-line ratios (Sections 3.2.1–3.2.4). Then, we assess the potential influence of observational selection effects on these properties (Sections 3.2.5) and, by considering only galaxies and progenitors in a fixed stellar-mass bin of (0.3–1.0) × $10^{11} M_{\odot}$ at all redshifts, that of the mass distribution of our sample (Section 3.2.6).

### 3.2.1 Line-ratio diagrams at redshifts $z < 0.5$

In the top row of Fig. 2, we show the locations of galaxies and their main progenitors extracted from all simulation snapshots at redshifts $z < 0.5$, in three line-ratio diagrams defined by [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α. The symbol types have the same meaning as in the bottom panels of Fig. 1 (stars: SF; circles: composite; triangles: AGN; squares: PAGB), but are now colour-coded according to global interstellar metallicity (as indicated in the middle panel). The grey shaded areas and contours indicate the location of SDSS galaxies in each diagram. To perform a meaningful comparison between models and observations, we show the effect of requiring a typical flux limit of $5 \times 10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2}$ for all simulated emission lines (e.g. table 1 of Juneau et al. 2014). The filled symbols correspond to galaxies satisfying this criterion, and the open symbols to those too faint to be detected. Also shown in these diagrams are standard observational criteria to distinguish SF galaxies from composites, AGN and LI(N)ER. The dashed, dotted and dot-dashed lines in the [O III]/H β versus [N II]/H α diagram have the same meaning as in Fig. 1, while the dotted lines in the other two diagrams distinguish SF galaxies (bottom left) from AGN (top) and LI(N)ER (bottom right; Kewley et al. 2001).

Fig. 2 shows that, at $z < 0.5$, simulated galaxies satisfying our conservative flux detection limit occupy the same areas as SDSS
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Optical emission-line ratio diagrams, [O III]/H β versus [N II]/H α (left column), [S II]/H α (middle column) and [O I]/H α (right column), for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors (symbols, colour-coded according to global interstellar metallicity, as indicated on the top scale), extracted from all simulation snapshots in different redshift intervals (different rows), as described in Section 3.2. The symbols and different lines have the same meaning as in Fig. 1. The simulations are compared to observations of local SDSS galaxies (grey shaded area, first row; also reported in light grey in higher-redshift bins for reference) as well as of distant galaxies (black diamonds) by Yabe et al. (2012, z ~ 1.4) and Steidel et al. (2014, z ~ 2.3). The small open symbols show the synthetic line ratios of all galaxies, regardless of luminosity, while the large filled symbols show galaxies above a flux limit of 5 × 10⁻¹⁷ erg s⁻¹ cm⁻² in all lines.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Average [O III]/H β emission-line ratio in bins of [N II]/H α (left-hand panel), [S II]/H α (middle panel) and [O I]/H α (right-hand panel) for the star-forming subset [i.e. with log (BHAR/SFR) < −4] of the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors, in different redshift ranges (connected stars of different colours). The SDSS data (in grey) and dashed and dotted lines are the same as in Fig. 2.</figcaption></figure>

galaxies in all three line-ratio diagrams (galaxies with line luminosities below the flux limit are typically very quiescent). Moreover, in general, simulated galaxies of SF, composite, AGN and PAGB types appear to fall in regions of the diagrams corresponding to the observationally defined SF, composite, AGN and LI(N)ER categories. This is remarkable in that, in our approach, the different types are connected to physical parameters, such as the BHAR/SFR ratio (for SF galaxies, composites and AGN) and the contribution to total H β luminosity (for PAGB-dominated galaxies). An in-depth investigation of the usefulness of this connection for the interpretation of nebular emission from distant galaxies will be the subject of a future study. For the remainder of the present paper, we conclude on the basis of this comparison with SDSS galaxies in the nearby Universe that our sample of zoom-in simulations of massive galaxies is well anchored at $z < 0.5$ to investigate the evolution of nebular emission with cosmic time.

### 3.2.2 Line-ratio diagrams at redshifts $z > 0.5$

Rows 2–6 of Fig. 2 show the analogue of the first row for the redshift bins $z = 0.5–1.0$, $1.0–1.5$, $1.5–2.0$, $2.0–3.0$ and $3.0–5.0$ (from top to bottom). As expected from the evolution of the three example galaxies in Fig. 1, Fig. 2 confirms that, for the whole sample of simulated galaxies, [O III]/H β tends to increase and [N II]/H α to decrease from low to high redshift, and that part of this trends at least is attributable to the lower global interstellar metallicity of high-redshift galaxies. The [S II]/H α and [O I]/H α ratios also tend to decrease towards high redshift. A most notable result from Fig. 2 is the consistency between the SDSS line properties of simulated galaxies brighter than the flux-detection limit (filled symbols) and the observed properties of star-forming galaxies from the sample of Yabe et al. (2012) at $z \sim 1.4$ (black diamonds with error bars in the third row) and that of Steidel et al. (2014) at $z \sim 2.3$ (black diamonds with error bars in the fifth row). In Section 4, we will exploit our fully self-consistent simulations to gain insight into the physical parameters that, along with global interstellar metallicity, are likely to contribute to the observed evolution of galaxies in these optical line-ratio diagrams.

### 3.2.3 Evolution of the SF-galaxy sequence in optical line-ratio diagrams

Motivated by a number of observational studies (e.g. Steidel et al. 2014; Shapley et al. 2015; Kashino et al. 2017; Strom et al. 2017), we quantify in Fig. 3 the redshift evolution of the average ‘star-forming-galaxy sequence’ in optical line-ratio diagrams. We define this sequence as the average [O III]/H β of SF galaxies (i.e. those with BHAR/SFR < $10^{-4}$; see Section 3.1) in bins of [N II]/H α (left panel), [S II]/H α (middle panel) and [O I]/H α (right panel), as shown by stars and solid lines in the three diagrams of Fig. 3. The different colours correspond to different redshift bins, as indicated. As in Fig. 2, the shaded areas show the distributions of local ($z \sim 0.1$) SDSS galaxies in these diagrams, while the dashed and dotted curves mark the upper boundary of the SF-galaxy region according to Kauffmann et al. (2003) and Kewley et al. (2001).

Fig. 3 shows that the rise in [O III]/H β as a function of redshift predicted by our zoom-in simulations amounts to roughly 0.4–0.5 dex at fixed [N II]/H α = 0.1, [S II]/H α = 0.1 and [O I]/H α = 0.01. These values are widely consistent with observed trends for the relation between [O III]/H β and [N II]/H α (e.g. Steidel et al. 2014; Kashino et al. 2017; Strom et al. 2017), while the existence of similar observational evidence for the relation between [O III]/H β and [S II]/H α is still debated (e.g. Shapley et al. 2015; Kashino et al. 2017, see also Section 5). We note that the predicted evolutionary trend of [O III]/H β is a specific outcome from our zoom-in simulations including AGN feedback. In the runs without AGN feedback, instead, [O III]/H β almost does not evolve with redshift. We discuss this in more detail in Appendix A.

### 3.2.4 Evolution of the global galaxy population in optical line-ratio diagrams

It is also of interest to examine the evolution of the emission-line properties of the global galaxy population, beyond that of purely SF galaxies. In Fig. 4, we show the redshift evolution of the average [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α for our full sample
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Redshift evolution of the average [O iii]/H β (red, top left), [N ii]/H α (beige, top right), [S ii]/H α (green, bottom left) and [O i]/H α (blue, bottom right) ratios of the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors (large filled circles and solid lines), together with the ±1σ scatter around the mean relations (shaded areas). The small filled circles and dashed lines show the results obtained when considering only galaxies above a flux limit of 5 × 10⁻¹⁷ erg s⁻¹ cm⁻² in all lines. Dot-dashed and dotted lines show the mean relations for SF- and composite/AGN-dominated galaxies, respectively. Also shown for references are observed mean ratios of SDSS emission-line galaxies in the local Universe (black open square: all galaxies; grey filled square: SF galaxies only) and of different samples of galaxies at various redshifts, as indicated (Yabe et al. 2014; Steidel et al. 2014, 2016; Vanzella 2016).</figcaption></figure>

of simulated galaxies (large filled circles and solid lines), together with the ±1σ scatter around this mean evolution (shaded areas). The small filled circles and dashed lines show the (very similar) evolution obtained when considering only galaxies above the flux-detection limit of 5 × 10⁻¹⁷ erg s⁻¹ cm⁻² (as in Fig. 2). As expected from our previous findings, Fig. 4 shows that the trends of increasing [O iii]/H β and decreasing [N ii]/H α, [S ii]/H α and [O i]/H α from low to high redshift predicted by our simulations are in good agreement with various observational results at z ≲ 3 (from SDSS; Steidel et al. 2014, 2016; Yabe et al. 2014; Vanzella 2016, black symbols, as indicated in the top-right panel), despite the sparse statistics of the simulations.

In the lowest redshift bin, the predicted average [O iii]/H β, [S ii]/H α and [O i]/H α agree reasonably well with the observed ones of SDSS galaxies at z ~ 0.1, but the predicted [N ii]/H α ratio appears significantly lower than observed out to z ~ 2.5. A mismatch in gas-phase metallicity is not likely to account for this discrepancy, given the good general agreement between models and observations for the other three ratios. At z ~ 0.1, part of the discrepancy could arise from the larger fraction of AGN-dominated galaxies in the SDSS sample relative to the simulations (which, because of low statistics, do not represent the same galaxy populations as in the SDSS and other observational samples), as the

average [N ii]/H α of AGN-dominated galaxies (dotted line in the top-right panel of Fig. 4) is significantly larger than that of SF galaxies (dot-dashed line). In fact, selecting only SF galaxies – i.e. those below the Kauffmann et al. (2003) criterion in the [O iii]/H β versus [N ii]/H α diagram – in the SDSS sample (grey filled square in each panel of Fig. 4) brings down [N ii]/H α far more significantly than [O iii]/H β, [S ii]/H α and [O i]/H α. The same argument does not hold for the mismatch with the Yabe et al. (2012) measurement of [N ii]/H α at z ~ 1.4 and the Steidel et al. (2014) one at z ~ 2.5, since these samples include, respectively, 0 and 4 per cent of AGN. The higher-than-predicted [N ii]/H α (and marginally lower-than-predicted [O iii]/H β) of these samples could arise from a difference in global interstellar metallicity. For reference, we checked that the evolution of the galaxy mass-metallicity relation out to z ~ 3 in our simulations is broadly consistent with that in fig. 8 of Maiolino et al. (2008).

### 3.2.5 Influence of flux-selection effects

The small differences between the dashed and solid lines in Fig. 4 indicate that, for our sample of simulated massive galaxies and their main high-redshift progenitors, evolution effects dominate
<figure><img src="image_8.png" /><figcaption>**Figure 5.** [O III]/H β (red, top row), [N II]/H α (beige, second row), [S II]/H α (green, third row) and [O I]/H α (blue, bottom row) as a function of stellar mass, in different redshift bins (different columns), for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors. In each panel, individual (pastel circles) and average (bright circles with error bars) line ratios are shown with and without including flux limit of $5 \times 10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2}$ in all lines (filled/big and open/small circles, respectively). The grey shaded area highlights the stellar-mass range $10.5 < \log (M_{\text{stellar}}/M_{\odot}) < 11.0$. Also shown are observations of local SDSS (black open square: all galaxies; grey filled square: SF galaxies only) and of galaxies at $z \sim 2.3$ from Steidel et al. (2014, black open triangles).</figcaption></figure>

over flux-selection effects in determining the redshift trends of [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α. It is important to stress the significance of this finding in the context of the recent study by Juneau et al. (2014), who showed the difficulty in disentangling evolution from selection effects in analyses of emission-line properties of distant galaxy samples. While a more robust conclusion would require a larger, statistically complete sample of simulated galaxies (not achievable with our current limited set of zoom-in simulations; see Section 5), the results of Fig. 4 already demonstrate the usefulness of our fully consistent modelling to interpret observations of the nebular emission from distant galaxy populations.

### 3.2.6 Influence of the stellar-mass evolution of simulated galaxies

Despite the generally encouraging consistency between model and observed trends in Fig. 4, we must remember that our predictions rely on a statistically small set of simulated massive galaxies, whose progenitor masses at high redshift are systematically smaller than

today. In this regard, it is important to check the influence of the global evolution of stellar mass (and metallicity) in the simulations on the predicted redshift trends of emission-line ratios. To this goal, we show in Fig. 5 the individual (pastel circles) and average (bright circles) line ratios as a function of stellar mass for galaxies and their main progenitors in different redshift bins (different columns), with and without including flux selection (filled/big and open/small circles, respectively). The grey shaded area highlights the stellar-mass range $10.5 < \log (M_{\text{stellar}}/M_{\odot}) < 11.0$ in all panels. Fig. 5 shows that, even in a fixed stellar-mass range, [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α for galaxies in our simulations exhibit strong redshift evolution, in good agreement with that inferred from the comparison between SDSS and the Steidel et al. (2014) data for [O III]/H β and [N II]/H α (both samples being available over a range of stellar masses; we note in passing the good overlap between the observed and model stellar-mass ranges at $z = 2–3$). Moreover, at a given redshift, the average line ratios show only a modest dependence on stellar mass, suggesting that any stellar-mass bias would have a limited influence on population-wide line-ratio evolution. Such evolution may instead depend more sensitively on the
evolution of other ISM and radiation properties, even at fixed stellar mass. This will be the topic of the next section.

# 4 PHYSICAL ORIGIN OF THE EVOLUTION OF OPTICAL EMISSION-LINE RATIOS

The physical origin of the observed evolution of optical emission-line ratios in galaxies is a heavily debated issue (Section 1). In this section, we use our models to investigate the origin of the predicted redshift evolution of [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α in terms of ISM and radiation properties of simulated galaxies. We focus on those physical quantities that control the evolution of the SF, AGN and PAGB models (Section 2.2). Our approach offers the unique advantage of exploring separately the relative influence of different physical parameters on the evolution of emission-line properties. In the next subsections, we carry out such an analysis for the whole set of simulated galaxies and their main progenitors, as well as for galaxies and progenitors in a fixed stellar-mass range.

## 4.1 Redshift evolution of the parameters controlling nebular emission

To identify the physical processes responsible for the drop in [O III]/H β and the rise in [N II]/H α, [S II]/H α and [O I]/H α from high to low redshift, we start by exploring how the parameters controlling the SF, AGN and PAGB nebular-emission models change with redshift. These parameters are potential drivers of the cosmic evolution of optical emission-line ratios. We show the redshift evolution of these parameters in the different panels of Fig. 6, for the SF (first and second rows), AGN (third and fourth rows) and PAGB (bottom row) models, and for the 20 simulated galaxies and their main progenitors in our sample (thin lines, colour-coded according to final stellar mass, as indicated). In each panel, we also show the corresponding average evolution (thick solid line) and associated ±1σ scatter (grey shaded area).

Some parameters exhibit little variation with redshift, such as the average global carbon-to-oxygen abundance ratio, (C/O)gas, glob, ISM density, $n_{gas, glob}$, AGN Eddington ratio, $f_{fedd}$, central interstellar metallicity $Z_{gas, 1kpc}$, and PAGB ionization parameter, $U_{sim, \odot}$. These parameters are not likely to drive the evolutionary trends in emission-line ratios seen in the simulations. Other parameters rise significantly from high to low redshift, such as the global interstellar metallicity, $Z_{gas, glob}$, AGN luminosity, $L_{AGN}$, and mass and age of PAGB stellar populations. Instead, the SF ionization parameter, $U_{sim, \star}$ (driven by the SFR), and the AGN ionization parameter, $U_{sim, \bullet}$ (driven by the average central ISM density, $n_{gas, 1kpc}$) show a marked decline. We note that, since $U_{sim, \bullet}$ scales as $L_{AGN}^{1/3}$ and $n_{gas, 1kpc}^{2/3}$ (equations 1 and 2), the drop in $U_{sim, \bullet}$ by itself indicates that $L_{AGN}$ has a negligible influence on this parameter (and hence on optical-line ratios; see Feltre et al. 2016). The rise in BHAR/SFR ratio further suggests a growing relevance of AGN nebular emission towards low redshift, and the drop in specific SFR an increasing importance of post-AGB stellar populations. These different evolving quantities, aside from $L_{AGN}$, can potentially drive the trends seen in emission-line ratios.

## 4.2 Relation between line-ratio and parameter evolution

To reveal which of the potential driving parameters identified in Section 4.1 effectively contribute to the redshift evolution of emission-line ratios predicted by the models (and seen in the observations), we plot in Fig. 7 the average [O III]/H β (red lines),

[N II]/H α (beige lines), [S II]/H α (green lines) and [O I]/H α (blue lines) ratios against the physical parameters controlling the SF (top row), AGN (middle and bottom rows) and PAGB (top and bottom rows) models, together with the ±1σ scatter about these mean relations (pastel shaded areas). We note that, while the average trends in Fig. 7 were derived by considering all simulated galaxies and their main progenitors at redshift $z < 5$, we have checked that the trends hardly depend on the exact redshift range chosen within this interval. We now examine the dependence of each optical-line ratio in Fig. 7 on those model parameters identified above as potential drivers of the cosmic evolution of galaxy nebular emission, separating between the influence of the SF, AGN and PAGB components.

### 4.2.1 [O III]/H β ratio

(i) *Influence of SF-related parameters:* The top-left panel of Fig. 7 shows that [O III]/H β (in red) decreases with increasing global interstellar metallicity. This is primarily because a rise in metallicity makes cooling more efficient (and oxygen is a major coolant), which reduces the electronic temperature in the ionized gas, and thus, the probability for collisional excitation of [O III] optical transitions (e.g. Stasińska 1980; Gutkin et al. 2016). A secondary effect is the softer ionizing radiation of metal-rich stars relative to metal-poor ones, which lowers the number density of O$^{++}$ ions, new-born stars having the same metallicity as star-forming gas in our simulations (Section 2.2.1; see e.g. Gutkin et al. 2016). The other panels in the top row of Fig. 7 further show that [O III]/H β rises with increasing (specific) SFR and volume-averaged gas density, and hence, by virtue of equations (1) and (2), with increasing SF ionization parameter. The rise in [O III]/H β with $n_{gas, glob}$ follows from the implied larger volume-filling factor at fixed $n_{H,\star}$ (equation 2), which increases the probability for oxygen to be photoionized twice. Double ionization is also favoured when a rise in the rate of ionizing photons is induced by a larger SFR at fixed $n_{gas, glob}$.

Combining these findings with the results of Fig. 6 suggests that both a rising global interstellar metallicity and declining strength of the stellar radiation field (i.e. SFR) can contribute to the drop in [O III]/H β from high to low redshift in Fig. 4. Instead, the roughly constant $n_{gas, glob}$ with redshift argues against a primary influence of this parameter on cosmic [O III]/H β evolution.

(ii) *Influence of AGN-related parameters:* The middle row of Fig. 7 shows that [O III]/H β strongly increases with increasing central gas density, $n_{gas, 1kpc}$ (third panel), and hence, increasing AGN ionization parameter, $U_{sim, \bullet}$ (second panel), for reasons analogous to those described above for the dependence of [O III]/H β on $n_{gas, glob}$ and $U_{sim, \star}$. Based on the results from Fig. 6, we conclude that the drop in $n_{gas, 1kpc}$ from high to low redshift can also contribute to the cosmic evolution of [O III]/H β.

(iii) *Influence of PAGB-related parameters:* We already noted above that global interstellar metallicity, which enters PAGB as well as SF nebular-emission models, could contribute to the drop in [O III]/H β from high to low redshift. The bottom row of Fig. 7 further shows that, among the other potential driving parameters identified in Section 4.1, [O III]/H β appears to vary significantly with only the mass – and not the age – of post-AGB stellar populations. We note that the drop in [O III]/H β at increasing $M_{PAGB}$ follows from the greater average metallicity of the most massive PAGB stellar populations, hosted by the most massive galaxies. These stellar populations produce softer ionizing spectra, and hence fewer O$^{++}$ ions, than their more metal-poor counterparts in less massive galaxies, an effect which dominates over the rise in the total number
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Redshift evolution of the different physical quantities used to select SF (first and second rows), AGN (third and fourth rows) and PAGB (bottom row) nebular-emission models, for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors (thin lines, colour-coded according to final stellar mass, as indicated): global interstellar metallicity ($Z_{\text{gas, glob}}$); SF ionization parameter ($U_{\text{sim,\*}}$); global gas density($n_{\text{gas, glob}}$); SFR; specific SFR; global C/O ratio; central interstellar metallicity ($Z_{\text{gas, 1kpc}}$); AGN ionization parameter ($U_{\text{sim,\bullet}}$); central gas density ($n_{\text{gas, 1kpc}}$); AGN luminosity ($L_{\text{AGN}}$); AGN Eddington ratio log ($f_{\text{edd}}$); BHAR/SFR ratio; PAGB ionization parameter ($U_{\text{sim, \circ}}$); mass ($M_{\text{PAGB}}$) and average age (AgePAGB) of stars older than 3 Gyr. In each panel, the thick black line and shaded area show the average evolution and $\pm 1\sigma$ scatter around it, including all galaxies and progenitors, while the thick dashed line shows the average evolution of galaxies and progenitors in the fixed stellar-mass range $10.5 < \log(M_{\text{stellar}}/M_{\odot}) < 11.0$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Average [O III]/H β (red), [N II]/H α (beige), [S II]/H α (green) and [O I]/H α (blue) ratios of the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors (thick lines), plotted against the same galaxy physical parameters as in Fig. 6, plus the mean metallicity of stars older than 3 Gyr ($Z_{\text{PAGB}}$). Shaded areas show the ±1σ scatter around the mean relations.</figcaption></figure>

of PAGB stars from low- to high-mass galaxies. The dependence of $M_{\text{PAGB}}$ on redshift in Fig. 6 confirms that this parameter can also contribute to the cosmic evolution of [O III]/H β.

### 4.2.2 [N II]/H α ratio

(i) *Influence of SF-related parameters:* The top row of Fig. 7 shows that [N II]/H α (in beige) rises significantly with increasing global interstellar metallicity and decreasing (specific) SFR, and in turn, decreasing SF ionization parameter. The rise of [N II]/H α with $Z_{\text{gas, glob}}$ follows primarily from the implied increased abundance of nitrogen, for which secondary production is included in the models (Section 2.2.1). Additionally, the softer ionizing radiation of metal-rich stars relative to metal-poor ones lowers the production probability of multiply ionized nitrogen relative to N⁺. A smaller probability of multiply ionized nitrogen is also the reason for the rise of [N II]/H α with decreasing (specific) SFR and $U_{\text{sim, *}}$ (by analogy with our discussion of [O III]/H β above). In light of the dependence of $Z_{\text{gas, glob}}$ and (specific) SFR on redshift in Fig. 6, we conclude that both parameters can contribute to the rise in [N II]/H α from high to low redshift in Fig. 4.

(ii) *Influence of AGN-related parameters:* Fig. 7 shows that [N II]/H α rises significantly with decreasing central gas density, and hence, decreasing AGN ionization parameter (middle row). This is because a lower $n_{\text{gas, 1kpc}}$ at fixed $n_{\text{H, *}}$ (lower $U_{\text{sim, *}}$) leads to a lower probability of multiply ionizing nitrogen at the expense of N⁺. Also, [N II]/H α rises with increasing BHAR/SFR ratio (bottom row), because the harder ionizing spectrum of an AGN relative

to young stars makes the electronic temperature – and collisional excitation of [N II] – larger. The redshift dependence of $n_{\text{gas, 1kpc}}$ and the BHAR/SFR ratio in Fig. 6 further suggest that these parameters are likely to contribute to the cosmic evolution of [N II]/H α.

(iii) *Influence of PAGB-related parameters:* As noted above, global interstellar metallicity, which enters PAGB as well as SF nebular-emission models, can contribute to the rise in [N II]/H α from high to low redshift. The bottom row of Fig. 7 further shows that [N II]/H α rises with increasing mass and increasing metallicity of post-AGB stellar population. This is because metal-rich stellar populations in massive galaxies produce softer ionizing spectra (and hence, less multiply-ionized nitrogen) than their more metal-poor counterparts in lower-mass galaxies. The dependence of $M_{\text{PAGB}}$ on redshift in Fig. 6 indicates that this parameter (linked to $Z_{\text{PAGB}}$) can also contribute to the cosmic evolution of [N II]/H α.

### 4.2.3 [S II]/H α and [O I]/H α ratios

(i) *Influence of SF-related parameters:* The top row of Fig. 7 shows that the dependence of [S II]/H α (green lines) and [O I]/H α (blue lines) on SF-related parameters is similar to that of [N II]/H α, for the same reasons as outlined for that ratio in Section 4.2.2. We note that [S II]/H α and [O I]/H α depend less strongly than [N II]/H α on global interstellar metallicity, because of the influence of secondary production in the case of nitrogen. We conclude that, as for [N II]/H α, global interstellar metallicity and (specific) SFR are likely to contribute to the rise of [S II]/H α and [O I]/H α from high to low redshift in Fig. 4.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Relative contributions by different components to the redshift evolution of the average [O iii]/H β (red, top left), [N ii]/H α (beige, top right), [S ii]/H α (green, bottom left) and [O I]/H α (blue, bottom right) ratios of the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors. In each panel, black stars (joined by a solid line) show the relation obtained using a simplified base model of pure-SF nebular emission, with fixed SF ionization parameter, log (U<sub>sim,★</sub>) = −3.5, and global interstellar metallicity, log (Z<sub>gas, glob</sub>) = log (Z<sub>☉</sub>). The other symbols/lines show the effect of adding to this model, step by step: AGN nebular emission (triangles/dotted line); PAGB nebular emission (squares/dashed line); and cosmic evolution of Z<sub>gas, glob</sub> (upside-down triangles/triple-dot-dashed line) and U<sub>sim,★</sub> (circles/solid line and shaded ±1σ scatter around the mean relation; identical to Fig. 4).</figcaption></figure>

(ii) *Influence of AGN- and PAGB-related parameters:* [S II]/H α and [O I]/H α depend in a similar way to [N II]/H α on AGN- and PAGB-related parameters in Fig. 7. Thus, central gas density (AGN ionization parameter), BHAR/SFR ratio and, to a lesser extent than for [N II]/H α (because of the shallower dependence on M<sub>PAGB</sub> in Fig. 7) mass of post-AGB stellar populations can also contribute to the cosmic evolution of [S II]/H α and [O I]/H α.

**4.3 Relative influence of the different drivers of cosmic line-ratio evolution**

In Section 4.2, we have seen that global and central interstellar metallicities, (specific) SFR, BHAR/SFR ratio and mass of post-AGB stellar population can all contribute to the cosmic evolution of the optical line ratios in Fig. 4. To identify the relative influence of these different parameters on the evolution of nebular emission, we now proceed with a 'cumulative' approach: we start by combining our zoom-in simulations of galaxy formation with a simplified model of pure-SF nebular emission, for which we adopt fixed U\* and Z\*, (Sections 2.2.1 and 2.3.1). Then, step by step, we examine how the predicted line-ratio evolution changes after introducing AGN and PAGB nebular emission and accounting for cosmic

evolution of $U_{\star}$ and $Z_{\star}$. For simplicity, we explore in this way only the overall effect of adding the AGN and PAGB components, without distinguishing between, e.g. the influence of central gas density and BHAR/SFR ratio. In practice, we start from an SF nebular-emission model with fixed log $U_{\text{sim},\star} = -3.5$ and $Z_{\star} = Z_{\text{gas, glob}} = Z_{\odot}$. We show in Fig. 8 the implied (flat) redshift evolution of the mean [O iii]/H β, [N ii]/H α, [S ii]/H α and [O I]/H α ratios obtained with this base model for our sample of simulated galaxies (black stars and solid lines).

4.3.1 [O III]/H β ratio

The top-left panel of Fig. 8 shows that adding the contribution from AGN nebular emission to the base model (red triangles and dotted line) slightly raises [O III]/H β relative to the black line, but without introducing any redshift evolution. Thus, central gas density and, in turn, AGN ionization parameter do not appear to affect significantly the cosmic evolution of [O III]/H β. The inclusion of PAGB nebular emission (red squares and dashed line) reduces slightly [O III]/H β at z < 1, but the effect is still negligible at higher redshift. Instead, the effect of incorporating in $Z_{\star}$ the evolution of the global interstellar metallicity ($Z_{\text{gas, glob}}$) predicted by the simulations (red upside-down
triangles and triple-dot-dashed line) is drastic: this causes $[O\text{iii}]/H\beta$ to drop sharply from high to low redshift. At $z < 1$, this trend is accentuated by the inclusion, through $U_*$, of the drop in average SFR predicted by the simulations (red circles and solid line, which, together with the shaded area, are the same as in Fig. 4).

Hence, the cosmic evolution of $[O\text{iii}]/H\beta$ predicted by our simulations results from combined effects of different parameters governing SF and PAGB (but not AGN) nebular emission. Global interstellar metallicity emerges as the main driver of this evolution, with minor contributions from SFR evolution and post-AGB stellar populations at redshift $z < 1$.

### 4.3.2 $[N\text{i}]/H\alpha$ ratio

The top-right panel of Fig. 8 shows that adding the contribution from AGN nebular emission to the base model (beige triangles and dotted line) makes $[N\text{i}]/H\alpha$ increase significantly at redshift $z < 1$. Including PAGB nebular emission reinforces slightly this trend (beige squares and dashed line). As in the case of $[O\text{iii}]/H\beta$, accounting for the evolution of $Z_{\text{gas, glob}}$ predicted by the simulations is crucial and makes $[N\text{i}]/H\alpha$ increase significantly from high to low redshift (beige upside-down triangles and triple-dot-dashed line). This rise is intensified at $z < 1$ by the inclusion of SFR evolution (beige circles and solid line).

Our simulations therefore predict that the cosmic evolution of $[N\text{i}]/H\alpha$ is regulated by different parameters governing SF, AGN and PAGB nebular emission. At $z > 1$, the evolution is dominated by that of global interstellar metallicity, and at $z < 1$, by the drops in central gas density (driving the AGN models) and SFR, together with the build-up of metal-rich populations of post-AGB stars.

### 4.3.3 $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ ratios

The bottom panels of Fig. 8 show that adding the contribution from AGN nebular emission to the base model makes both $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ (green/blue triangles and dotted lines) increase sharply at redshifts below unity. Including PAGB nebular emission reinforces significantly these trends (green/blue squares and dashed lines). In strong contrast with $[O\text{iii}]/H\beta$ and $[N\text{ii}]/H\alpha$, accounting for the evolution of global interstellar metallicity predicted by the simulations has only a modest effect on the evolution of $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$, and only at redshift $z > 2$ (green/blue upside-down triangles and triple-dot-dashed lines). Instead, the inclusion of SFR evolution strengthens the rise in $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ over the whole range from high to low redshift (green/blue circles and solid lines).

Hence, as for $[O\text{iii}]/H\beta$ and $[N\text{ii}]/H\alpha$, our simulations predict that the cosmic evolution of $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ arises from combined effects of different parameters governing SF, AGN and PAGB nebular emission. In this case, global interstellar metallicity plays only a minor role, and the evolution of $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ appears to be driven primarily by that of central gas density and SFR and the build-up of metal-rich populations of post-AGB stars.

## 4.4 Influence of the intrinsic stellar-mass evolution of simulated galaxies

It is important to check how the global evolution of stellar mass in our simulations (Section 3.2.6) affects our conclusions regarding the physical origin of cosmic evolution of optical line ratios in Section 4.3. We already demonstrated that the predicted cosmic evolution of $[O\text{iii}]/H\beta$, $[N\text{ii}]/H\alpha$, $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ persists

even when considering galaxies (progenitors) in a fixed stellar-mass range at all redshifts (Fig. 5). However, SFR and global interstellar metallicity, identified above as primary drivers of the evolution of nebular emission, correlate with stellar mass in the simulations, reflecting the galaxy mass-metallicity relation and SF main sequence (see the model predictions by Hirschmann et al. 2013; Hirschmann, De Lucia & Fontanot 2016). In fact, the thick dashed lines in Fig. 6 show the average redshift evolution of the parameters controlling the SF, AGN and PAGB nebular-emission models, when considering only galaxies in the fixed stellar-mass range $10.5 < \log (M_{\text{stellar}}/M_{\odot}) < 11.0$, independent of redshift. The evolutionary trends are similar to those obtained when including the full sample (thick solid lines), except for the loss of redshift dependence of the average interstellar metallicity (top-left panel) and the average mass of post-AGB stellar populations (second panel of bottom row). We also checked that the results of Fig. 7 do not change when considering galaxies in a fixed mass range.

In Fig. 9, we show the analogue of Fig. 8 obtained when considering only galaxies with masses in the range $10.5 < \log (M_{\text{stellar}}/M_{\odot}) < 11.0$ at all redshifts. As anticipated in the discussion of Fig. 5 (Section 3.2.6), the predicted global evolution of optical-line ratios remains very strong in this case (circles and solid lines in all panels). Moreover, the influence of different physical quantities on this evolution is qualitatively similar to that in Fig. 8, although the relative contributions by global interstellar metallicity and SFR change significantly. At fixed stellar mass, global interstellar metallicity now has only a weak impact on the cosmic evolution of $[O\text{iii}]/H\beta$ and $[N\text{ii}]/H\alpha$ down to $z = 0$ (red/beige upside-down triangles and triple-dot-dashed lines), and SFR is the primary driver of the evolution of any line ratio at $z > 1.5$ (circles and solid lines). This is because, at fixed stellar mass, SFR evolves more strongly with redshift than global interstellar metallicity (thick dashed lines in Fig. 6). The bottom two panels of Fig. 9 further show than nebular emission from AGN and post-AGB stellar populations still contributes significantly to the increase in $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ at $z < 1.5$. We note that, in this case, the rising influence of the PAGB contribution at fixed stellar mass does not come from an increase in $M_{\text{PAGB}}$ (Fig. 6), but simply from the drop in specific SFR.

Remarkably, therefore, we find that identifying the physical origin of cosmic evolution of optical-line ratios depends sensitively on the stellar-mass properties of the considered sample. When considering our full set of simulated galaxies, for which mass evolves with redshift, the associated evolution of global interstellar metallicity accounts for most of the predicted evolution of line ratios. Instead, when considering galaxies in a fixed stellar-mass range of $10.5 < \log (M_{\text{stellar}}/M_{\odot}) < 11.0$, global interstellar metallicity has a weaker influence on the predicted evolution, which is primarily driven by SFR (through the ionization parameter), at all redshifts for $[O\text{iii}]/H\beta$ and $[N\text{ii}]/H\alpha$ and at $z > 1.5$ for $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$.

To summarize the above trends, we show graphically with arrows in Fig. 10 how redshift, $[O\text{iii}]/H\beta$, $[N\text{ii}]/H\alpha$, $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ (different rows) change qualitatively as different radiation and ISM properties increase (different columns), for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors. Grey cells refer to properties varying with redshift, grey/white-gradient cells to properties whose dependence on redshift vanishes when selecting a fixed stellar-mass range, and white cells to properties independent of redshift. Coloured arrows highlight those physical quantities expected to drive the cosmic evolution of $[O\text{iii}]/H\beta$, $[N\text{ii}]/H\alpha$, $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$, according to our analysis of Figs 8 and 9.
<figure><img src="image_12.png" /><figcaption>Figure 9. Same as Fig. 8, but for the subsample of galaxies and progenitors with stellar masses in the range 10.5 < log (M<sub>stellar</sub>/M<sub>☉</sub>) < 11.0.</figcaption></figure>

## 4.5 Physical origin of the evolution of SF galaxies in the [O III]/H β versus [N II]/H α diagram

It is of interest to examine the evolution of the SF-galaxy population predicted by our simulations in the [O III]/H β versus [N II]/H α diagram. This is shown in Fig. 11, where stars indicate the average line ratios of all SF-dominated galaxies at different redshifts (as indicated). The predicted evolution is primarily driven by global interstellar metallicity and (specific) SFR (ionization parameter). SF galaxies in the top-left quadrant of Fig. 11 tend to have low global interstellar metallicities, log (Z<sub>gas, glob</sub>/Z<sub>☉</sub>) < -0.2, and high specific SFR, log (sSFR/yr<sup>-1</sup>) > 10.5, while those in the bottom half have on average high global interstellar metallicities, log (Z<sub>gas, glob</sub>/Z<sub>☉</sub>) > -0.2, and low specific SFR, log (sSFR/yr<sup>-1</sup>) < 10.5. The rise in [N II]/H α with redshift is also influenced by an increase in BHAR/SFR ratio. Considering only galaxies (progenitors) in a fixed stellar-mass range of 10.5 < log (M<sub>stellar</sub>/M<sub>☉</sub>) < 11.0 hardly affects the predicted evolution in Fig. 11, as shown by the yellow dashed line. In this case, global interstellar metallicity has a negligible influence on the evolution, which is mostly driven by the drop in (specific) SFR from high to low redshift. It is worth pointing out that, if AGN feedback is not included in the simulations, the predicted star formation histories are much flatter and do not show such significant drop in (specific) SFR (see Choi et al. 2017). Hence, we find that AGN feedback plays a key role in the predicted [O III]/H β evolution of SF-dominated galaxies. We discuss this in more detail in Appendix A.

## 5 DISCUSSION

In the last two sections, we demonstrated that the redshift evolution of optical emission-line ratios predicted by our simulations of massive galaxies is widely consistent with available observations. We investigated the physical origin of this evolution in terms of ISM and ionizing-source parameters (as summarized by Fig. 10). In this section, we discuss the potential influence of other parameters, so far fixed at standard values in the simulations (Section 2.3), on these results. We also discuss some caveats of our nebular-emission models and zoom-in simulations. Finally, we put our results on the cosmic evolution of emission-line ratios in the context of previous observational and theoretical studies.

### 5.1 Influence of fixed model parameters

When coupling nebular-emission models with galaxy simulations in Section 2.3, we adopted standard values of parameters, such as the IMF upper mass cut-off, $m_{up} = 100~M_{☉}$; the hydrogen density in H II regions, $n_{H, *} = 100~cm^{-3}$; the slope of AGN ionizing radiation, $\alpha = -1.7$; and the dust-to-metal mass ratio in ionized gas, $\xi_{d} = 0.3$ (Table 1). We now investigate the influence of these parameters of the predicted [O III]/H β and [N II]/H α ratios of simulated galaxies (the influence on [S II]/H α and [O I]/H α being similar to that on [N II]/H α). By analogy with Fig. 2, we show in Fig. 12 the locations of galaxies and their main progenitors in the [O III]/H β versus [N II]/H α diagram in three redshift bins, z = 0–0.5 (left column),
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Schematic summary of Figs 6, 7, 8 and 9: arrows visualize qualitatively how redshift, $[O\text{ III}]/H\beta$, $[N\text{ II}]/H\alpha$, $[S\text{ II}]/H\alpha$ and $[O\text{ I}]/H\alpha$ (different rows) change as different radiation and ISM properties increase (different columns, referring to the same physical quantities as shown in Fig. 6), for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors. Grey cells refer to properties found to vary with redshift, grey/white-gradient cells to properties whose dependence on redshift vanishes when selecting a fixed stellar-mass range, and white cells to properties independent of redshift. Coloured arrows highlight those physical quantities expected to drive the cosmic evolution of $[O\text{ III}]/H\beta$, $[N\text{ II}]/H\alpha$, $[S\text{ II}]/H\alpha$ and $[O\text{ I}]/H\alpha$, according to our analysis of Figs 8 and 9.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Redshift evolution in the $[O\text{ III}]/H\beta$ versus $[N\text{ II}]/H\alpha$ diagram of the average location of the star-forming subset [i.e. with $\log(\text{BHAR/SFR}) < -4$] of the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors (stars colour-coded according to redshift, joined by the red line). The yellow dashed line shows the result obtained when including only galaxies in the fixed stellar-mass range $10.5 < \log(M_{\text{stellar}}/M_{\odot}) < 11.0$. The reference SDSS data (in grey) and black dashed line are the same as in Fig. 3. SF galaxies in the top-left quadrant tend to have low global interstellar metallicities and high specific SFR (and ionization parameter), while those in the bottom half have on average high global interstellar metallicities and low specific SFR (ionization parameter).</figcaption></figure>

1–1.5 (middle column) and 2–3 (right column). The top row shows the results obtained when using standard values of all fixed parameters (identical to the first, third and fifth panels in the left-hand column of Fig. 2). The other rows show the results obtained when adopting $m_{\text{up}} = 300 M_{\odot}$ (second row), $n_{\text{H,\*}} = 10^3 \text{ cm}^{-3}$ (third row), $\alpha = -1.2$ and $-2.0$ (fourth and fifth rows) and $\xi_d = 0.1$ and $0.5$ (sixth and seventh rows). The average SF-galaxy sequence in the top panels (yellow line) is reported in other panels (yellow dashed line) for easy comparison with the SF sequence obtained using other parameters (red lines).

### 5.1.1 Upper mass cut-off of the IMF

Increasing the upper mass cut-off of the IMF from $m_{\text{up}} = 100–300 M_{\odot}$ leads to a harder ionizing spectrum of the stellar populations, since stars with initial masses greater than 100 $M_{\odot}$ evolve at higher effective temperatures than lower-mass stars. As shown by Gutkin et al. (2016), this causes only a slight increase in $[O\text{ III}]/H\beta$ and $[N\text{ II}]/H\alpha$, illustrated by the marginal difference between the red and yellow lines in the second row of Fig. 12. Thus, we do not expect that changes in $m_{\text{up}}$ over time would strongly affect the evolutionary trends in these optical line ratios.

### 5.1.2 Hydrogen density in H II regions

As outlined in Gutkin et al. (2016, see their section 3.3), a rise in the density of gas clumps in H II regions from $n_{\text{H,\*}} = 10^2–10^3 \text{ cm}^{-3}$ increases the probability of collisional de-excitation over radiative cooling, especially for infrared fine-structure transitions, resulting in a slight increase of cooling through optical transitions, and hence $[O\text{ III}]/H\beta$ and $[N\text{ II}]/H\alpha$. The effect can be significant at large interstellar metallicities, but vanishes at low metallicities (Fig. 5 of
<figure><img src="image_15.png" /><figcaption>**Figure 12.** [O iii]/H β versus [N ii]/H α diagram for the 20 simulated massive galaxies of Section 2 and their main high-redshift progenitors, in different redshift intervals (different columns), and for different values of four parameters otherwise fixed by default (at $m_{up} = 100 M_{\odot}$, $n_{H,*} = 10^2 \text{ cm}^{-3}$, $\alpha = -1.7$ and $\xi_d = 0.3$, top row) in the simulations: adopting an upper IMF cut-off of $m_{up} = 300 M_{\odot}$ (second row); an ionized-gas density of $n_{H,*} = 10^3 \text{ cm}^{-3}$ (third row), a slope of AGN ionizing radiation $\alpha = -1.2$ and $-2.0$ (fourth and fifth rows) and a dust-to-metal mass ratio $\xi_d = 0.1$ and $0.5$ (sixth and seventh rows). All symbols, lines and observational data points have the same meaning as in Fig. 2. The average SF-galaxy sequence in the top panels (yellow line) is reported in other panels (yellow dashed line) for easy comparison with the SF sequence obtained using other parameters (red lines).</figcaption></figure>
Gutkin et al. 2016). This is illustrated by the small difference between the red and yellow lines in the third row of Fig. 12, which decreases from low to high redshift together with mean global interstellar metallicity (Fig. 6). We find a similarly weak effect when increasing $n_{H,\bullet}$ and $n_{H,\odot}$ in the AGN and PAGB models. Hence, even if the ionized-gas density in high-redshift galaxies were typically larger than that in local ones, as is plausible according to some observational studies (e.g. Steidel et al. 2014; Sanders et al. 2016; Kashino et al. 2017), our conclusions about the main drivers of the cosmic evolution of optical-line ratios would not be altered.

### 5.1.3 Hardness of AGN ionizing radiation

Stellar ionizing radiation in our models is controlled by the physical properties of simulated galaxies. For SF models, the ionizing radiation depends on the properties of young stars, whose metallicity is assumed to be the same as that of star-forming gas (as justified in Section 2.2.1). The impact of stellar metallicity on optical-line ratios is in any case negligible compared to that of interstellar metallicity (e.g. Gutkin et al. 2016). For PAGB models, the age and metallicity of old, post-AGB stellar populations predicted by our simulations define unambiguously the corresponding ionizing radiation. In contrast, we set the hardness of AGN ionizing radiation via a fixed power-law index $\alpha = -1.7$ (Section 2.3.2). A larger (lower) $\alpha$ would imply a harder (weaker) spectrum, resulting in enhanced (reduced) $[O\,III]/H\beta$ and $[N\,II]/H\alpha$ (Feltre et al. 2016). This is confirmed by Fig. 12, which shows at all redshifts a marked decline in these ratios for AGN and composite galaxies as $\alpha$ drops from -1.2 (fourth row) to -2.0 (fifth row). Fixing $\alpha$ at a value different from -1.7 in the AGN models would therefore affect the emission-line properties of AGN and composite galaxies, but not their cosmic evolution. Current observations provide hardly any constraints on the dependence of $\alpha$ on redshift and/or BH accretion rate and suggest simply a range of plausible values between -1.2 and -2.0 (e.g. Feltre et al. 2016). Even if $\alpha$ dropped systematically from high to low redshift, this would induce a corresponding drop not only in $[O\,III]/H\beta$, but also in $[N\,II]/H\alpha$, which would be unlikely to contribute to the observed evolutionary trends in Fig. 4.

### 5.1.4 Dust-to-metal mass ratio (and relation to N/O)

Increasing the dust-to-metal mass ratio depletes metal coolants from the gas phase. The electronic temperature rises, as does cooling through collisionally excited optical transitions. This implies a strong rise in $[N\,II]/H\alpha$ and a less strong rise in $[O\,III]/H\beta$, since oxygen is a refractory element strongly depleted from the gas phase. This is illustrated by the fifth and sixth rows of Fig. 12, showing predictions for $\xi_d = 0.1$ and 0.5, and where the average SF galaxy sequence (red line) lies, respectively, below and above that for $\xi_d = 0.3$ (yellow line). A systematic drop in $\xi_d$ from high to low redshift would therefore induce a corresponding drop in $[O\,III]/H\beta$. However, this would be accompanied by an even stronger drop in $[N\,II]/H\alpha$, which argues against a significant influence of the dust-to-metal mass ratio on the observed evolutionary trends in Fig. 4. Moreover, this would go against recent support for a lower dust-to-metal mass ratio in metal-poor galaxies compared to metal-rich galaxies (Rémy-Ruyer 2014), also predicted by some theoretical studies (Inoue 2003; Asano et al. 2013).

It is worth mentioning that, O being a refractory element and N not, our models account naturally for the rise in gas-phase N/O abundance ratio, $(N/O)_{gas}$, implied by a rise in $\xi_d$ (see Fig. 1 of

Gutkin et al. 2016).<fnref n="8" /> An elevated $(N/O)_{gas}$ in star-forming galaxies at high redshift has been suggested by some observational studies (e.g. Shapley et al. 2015, see Section 5.3.1 below). Our results indicate that, while this could contribute for a modest part to the observed drop in $[O\,III]/H\beta$ from high to low redshift (by less than 0.2 dex, according to the difference between the red and yellow lines in Fig. 12), most of the trend in $[O\,III]/H\beta$, and all of that in $[N\,II]/H\alpha$, are likely to be dominated by the drop in (specific) SFR and rise in global interstellar metallicity (which make $[O\,III]/H\beta$ decline by ~0.8 dex; see Fig. 9).

Further support for a negligible role of the dust-to-metal mass ratio (and hence N/O in our models) in the cosmic evolution of $[O\,III]/H\beta$ comes from recent observations (e.g. Rémy-Ruyer 2014; De Cia et al. 2016; Wiseman et al. 2017) and predictions from a semi-analytic model accounting for dust formation and destruction (Popping, Somerville & Galametz 2017). Both types of studies suggest that, at given gas metallicity, the dust-to-metal mass ratio in massive galaxies hardly changes with redshift out to $z=6$, the models further predicting that this is also the case at fixed stellar mass (see figs 5 and 6 of Popping et al. 2017). In these models, galaxies more massive than $\log(M_{\text{stellar}}/M_{\odot}) \sim 9.5$ are predicted to have dust-to-metal mass ratios around 0.2, reasonably close to the typical observed values of 0.3–0.4 in present-day galaxies (with a scatter sampling the full range from 0.1 to 1.0). This reinforces our choice of adopting a fixed $\xi_d = 0.3$ at all redshifts in our simulations (Table 1).

## 5.2 Caveats

### 5.2.1 Nebular-emission models: shocks, leakage of ionizing photons and dusty tori

(i) *Shock models:* The SF, AGN and PAGB nebular-emission models used in this paper include several important features allowing accurate comparisons with observations, such as improved prescriptions for stellar interiors and atmospheres, the self-consistent treatment of depletion of metals on to dust grains and the inclusion of secondary nitrogen production. This results in robust predictions of emission-line fluxes for given combinations of adjustable model parameters. In contrast to other theoretical studies (Orsi et al. 2014; Shimizu et al. 2016), we also account self-consistently for nebular emission from AGN and post-AGB stars, in addition to that from young stars. Yet, we do not include nebular emission from radiative shocks potentially produced by, e.g. starburst- and AGN-driven outflows and cloud–cloud collisions in galaxy interactions (e.g. Rich et al. 2010; Sharp & Bland-Hawthorn 2010; Rich, Kewley & Dopita 2011; Soto et al. 2012; Weistrop et al. 2012).

Since many observations indicate that galactic winds exist in a significant fraction of galaxies out to high redshifts (e.g. Steidel et al. 2010; Genzel et al. 2011; Kornei et al. 2012; Newman 2012), shock excitation may be an important contributor to optical-line ratios in high-redshift galaxies. Theoretically, shock models were developed by Dopita (2003), Allen et al. (2008, fast shocks) and Rich et al. (2010, 2011, slow shocks). Based on these models,

<footnote n="8">⁸ A rise in $(N/O)_{gas}$ in our models could also be induced by a change in C/O ratio (Fig. 1 of Gutkin et al. 2016). However, this is unlikely for the simulations of massive galaxies presented here, which show hardly any evolution of C/O over cosmic time (Fig. 6). We note that other scenarios not considered in this work, such as chemically differential galactic winds (Vincenzo et al. 2016) or changes in the nitrogen yield, may also influence $(N/O)_{gas}$.</footnote>
Kewley et al. (2013, their Fig. 6) show that metallicity and shock velocity are the primary parameters defining the location of pure-shock models in the $[O\text{III}]/H\beta$ versus $[N\text{II}]/H\alpha$ diagram. In fact, shock models behave very similarly to AGN models, with $[N\text{II}]/H\alpha$ strongly increasing with increasing metallicity at roughly constant $[O\text{III}]/H\beta$. For shock models of any velocity and metallicity, the results of Kewley et al. (2013) indicate that $[O\text{III}]/H\beta$ is always predicted to be larger than unity ($\log [O\text{III}]/H\beta > 0$), which could, when assuming that shocks primarily arise in high-redshift galaxies, contribute to the cosmic drop in $[O\text{III}]/H\beta$ from high to low redshift. However, reality is more complex: shocks are generally not expected to dominate line emission from a galaxy and are usually associated with star formation and AGN activity (see e.g. Kewley et al. 2013). Exploring self-consistently the contribution from shocks to the cosmic evolution of emission-line ratios would therefore require that we identify shocks in our simulations and couple their properties (velocity, metallicity) with, for example, the fast- and slow-shock models of Allen et al. (2008) and Rich et al. (2010, 2011). We postpone such an analysis to a future study.

(ii) *Leakage of ionizing photons:* All SF, AGN and PAGB models adopted in this paper are ionization bounded, i.e. they do not include any leakage of ionizing photons into the circumgalactic medium. The escape of ionizing radiation from high-redshift galaxies is a heavily debated issue, with so far, at redshifts around 3, only few galaxies identified as genuine ionizing-photon leakers (de Barros 2016; Shapley et al. 2016; Bian et al. 2017). The escape of substantial ionizing radiation from density-bounded HII regions is expected to decrease the emission contribution of low-ionization species, affecting emission-line ratios of high-to-low ionization species, such as $[O\text{III}]/[O\text{II}]$ (e.g. Jaskot & Oey 2013, and references therein). Recent calculations by Izotov, Thuan & Gusova (2017, their fig. 6) further suggest that $[O\text{III}]/H\beta$ is hardly affected by the escape of ionizing photons from density-bounded HII regions. We might expect $[N\text{II}]/H\alpha$, which involves a lower-ionization line, to be more strongly affected. A detailed analysis of this issue is postponed to future work.

(iii) *Dust obscuration on torus level:* The AGN nebular-emission models used in this paper neglect the potential influence of a dust torus on emission-line luminosities. In principle, accounting for the presence of a dust torus would reduce the amount of ionizing photons emitted by the central AGN capable of reaching the narrow-line region, hence lowering all emission-line luminosities by a same factor (corresponding to the fraction of solid angle occupied by the torus). For an AGN-dominated galaxy, the predicted integrated line ratios would not change, but in galaxies with substantial SF and PAGB components, the relative contribution to these ratios by the AGN would be altered. Thus, some galaxies identified as composite systems in the simulations could instead move to the SF and LINER regions in the $[O\text{III}]/H\beta$ versus $[N\text{II}]/H\alpha$ diagram (Figs 2 and 12). We can speculate that such changes might become more significant towards high redshift, where torus obscuration could become more relevant (e.g. Hasinger 2008). While worthwhile, a more detailed investigation of this effect is beyond the scope of this study.

### 5.2.2 Cosmological zoom-in simulations: statistical completeness

A potential source of inaccuracy of all current cosmological simulations is related to the necessary assumption of often rather simplified sub-resolution models for baryonic processes, which cannot accurately capture the properties of a multiphase ISM (see e.g. Naab & Ostriker 2017). Specifically, different models for stellar (e.g. Guedes

et al. 2011; Hirschmann et al. 2013, 2016; Hopkins et al. 2013; Stinson et al. 2013) and AGN feedback (e.g. Dubois et al. 2013; Steinborn et al. 2015, 2016; Choi et al. 2017; Weinberger et al. 2017) have been shown to be capable of significantly affecting various galaxy properties at all times. In this context, it is important to stress that the simulations presented in this paper rely on physically motivated models for stellar and AGN feedback, constrained in part by small-scale simulations of the ISM (see Section 2.1 above and Choi et al. 2017, Hirschmann et al., in preparation). This provides fairly realistic massive galaxies in terms of chemical enrichment, SFR and BHAR histories, stellar populations, etc. Specifically, we checked that the predictions of our simulations are widely consistent with observed relations between SFR, gas-phase metallicity and galaxy stellar mass at redshifts out to $z=2$ (using observational constraints from Daddi et al. 2007; Elbaz et al. 2007; Maiolino et al. 2008; Andrews & Martini 2013; Maier et al. 2015). This success puts our predictions for emission-line ratios – controlled by these physical quantities – on firm grounds.

A drawback of our suite of cosmological zoom-in simulations is their sparse statistics, since we analyse only 20 massive galaxies, whose most massive progenitors we follow back in time. As a result, our sample may be missing a contribution by less massive galaxies with emission-line fluxes potentially above the adopted flux-detection limit. On the one hand, this implies an intrinsic stellar-mass evolution of simulated galaxies, whose effect we addressed in Sections 3.2.6 and 4.4. On the other hand, the fraction of AGN/composite galaxies in our limited set of zoom-in simulations may not be realistic and comparable to that of observationally sampled galaxy populations. The fraction of AGN, and hence, their relevance to the cosmic evolution of emission-line ratios, may change when considering a statistically complete population of galaxies and AGN. This would affect more $[N\text{II}]/H\alpha$, $[S\text{II}]/H\alpha$ and $[O\text{I}]/H\alpha$ than $[O\text{III}]/H\beta$, which is insensitive to nebular emission from AGN (Figs 8 and 9). Increasing our set of zoom-in simulations or appealing to statistically complete samples of galaxy simulated using semi-analytic models and large cosmological boxes could help us improve this limitation, but lies beyond the scope of this study.

## 5.3 Comparison with previous explanations for the cosmic evolution of $[O\text{III}]/H\beta$

In Section 4, our combination of zoom-in galaxy simulations with versatile nebular-emission models allowed us to show that a decline in (specific) SFR (ionization parameter) and a rise in global interstellar metallicity are the main potential drivers of the observed drop in $[O\text{III}]/H\beta$ of star-forming galaxies from high to low redshifts. At fixed stellar mass, the higher $[O\text{III}]/H\beta$ of distant galaxies relative to local ones is primarily attributable to an elevated (specific) SFR. In contrast, a larger IMF upper-mass cut-off, higher ionized-gas density, harder ionizing radiation and larger depletion of metals on to dust grains (and associated elevated gas-phase N/O ratio) are not expected to play any significant role in the cosmic evolution of optical-line ratios. A number of different and sometimes contradictory explanations have been proposed in previous, mostly observational studies. We now replace our model predictions in the context of these studies.

### 5.3.1 Observational studies

Different observational investigations have favoured different physical origins of the cosmic evolution of $[O\text{III}]/H\beta$, with no unique
conclusion having been drawn so far. Possible explanations include the prevalence at high redshift of increased SFR (e.g. Kashino et al. 2017), higher hydrogen/electron densities (e.g. Brinchmann et al. 2008; Lehnert et al. 2009), increased contribution by an AGN to nebular emission (e.g. Wright et al. 2010), elevated (N/O)gas ratio (e.g. Shapley et al. 2015) and harder stellar ionizing radiation (e.g. Steidel et al. 2014; Strom et al. 2017) compared to low redshift. We now examine the conclusions from these studies in light of our theoretical predictions.

(i) Gas-phase metallicity: Based on a sample of 701 star-forming galaxies at $z = 1.4–1.7$ from the Fibre Multi-Object Spectrograph (FMOS)-COSMOS survey, Kashino et al. (2017) find that the empirically determined gas-phase metallicities of galaxies with masses $log(M_{\text{stellar}}/M_{\odot}) \gtrsim 11$ are similar to those of local SDSS counterparts, while less massive galaxies exhibit a rise in gas-phase metallicity from high to low redshift. Kashino et al. (2017) show that gas-phase metallicity can affect [O III]/H β and [N II]/H α, but they further note that a rise in metallicity at fixed ionization parameter cannot account for the observed cosmic evolution of [O III]/H β and [N II]/H α. They conclude that a change in ionization parameter is likely to be the primary cause of this evolution (their Fig. 12 and section 4.1). These results are consistent with our conclusion in Section 4.4 that, at fixed stellar mass in the range $10.5 < \log(M_{\text{stellar}}/M_{\odot}) < 11.0$, the cosmic evolution of optical-line ratios is primarily driven by (specific) SFR – i.e. ionization parameter – rather than global interstellar metallicity.

(ii) Ionization parameter and SFR: Brinchmann et al. (2008), Steidel et al. (2014, from the analysis of 251 star-forming galaxies at $z \sim 2.3$ with the Keck/MOSFIRE spectrometer), Hayashi et al. (2015, from the analysis of 118 star-forming galaxies at $z \sim 1.5$ with the Subaru Fibre Multi Object Spectrograph) and Kashino et al. (2017, see above) attribute at least part of the enhanced [O III]/H β ratio of high-redshift galaxies relative to low-redshift ones to a higher ionization parameter (see also Cullen et al. 2016). So far, no consensus has been reached on whether this quantity plays a primary (Kashino et al. 2017) or only minor (Steidel et al. 2014) role. Hayashi et al. (2015) and Kashino et al. (2017) discuss different possible origins of an elevated ionization parameter (higher star-formation efficiency, top-heavy IMF, harder ionizing spectra, etc.) without drawing any final conclusion. Our results in Section 4, based on the self-consistent modelling of galaxy formation and nebular emission, confirm the importance of an evolving ionization parameter in the cosmic evolution of [O III]/H β, especially at fixed galaxy stellar mass. In these models, (specific) SFR is the main driver of the change in ionizing parameter over cosmic time.

(iii) Hydrogen/electron density in H II regions: Densities larger than in typical local H II regions have been repeatedly measured in high-redshift star-forming galaxies and put forward as a possible cause for enhanced [O III]/H β (e.g. Brinchmann et al. 2008; Lehnert et al. 2009; Steidel et al. 2014; Sanders et al. 2016). Yet, most observational studies tend to agree that changes in ionized-gas density alone cannot account for the cosmic evolution of [O III]/H β (e.g. Rigby et al. 2011; Hayashi et al. 2015; Kashino et al. 2017; Strom et al. 2017). Specifically, based on the analysis of 380 star-forming galaxies at $z = 2–3$ with the Keck/MOSFIRE spectrometer, Strom et al. (2017) argue that measurements of the density-sensitive [O II]λλ3727, 3729 doublet for galaxies at small and large offsets from the local SF sequence in the [O III]/H β versus [N II]/H α diagram are nearly identical. This observational conclusion is consistent with our theoretical finding that changes in $n_{H,*}$ hardly affect the cosmic evolution of these optical-line ratios (third row in Fig. 12).

(iv) Additional contribution from (weak) AGN: Most recent observational studies rule out enhanced contribution by a central accreting BH as the main driver for the cosmic evolution of [O III]/H β, since the high-redshift star-forming galaxies showing this evolution do not exhibit any spectral signature of an AGN (such as strong high-ionization lines; Steidel et al. 2014; Strom et al. 2017; Kashino et al. 2017). This is consistent with our theoretical prediction that the cosmic evolution of [O III]/H β in star-forming galaxies (and composite and AGN galaxies) is not driven by nebular emission from an AGN, which can cause at most a slight offset at any redshifts (Figs 8 and 9).

(v) Enhanced (N/O)gas ratio at fixed (O/H)gas ratio: Some studies favour an enhanced (N/O)gas ratio at fixed (O/H)gas as the primary driver of the cosmic evolution of [O III]/H β and [N II]/H α (Masters 2014; Shapley et al. 2015; Yabe et al. 2015; Cowie, Barger & Songaila 2016; Masters et al. 2016; Sanders et al. 2016).<fnref n="9" /> In fact, these studies attribute the offset of the ‘SF branch’ in the line-ratio diagram at high redshift to a rise in [N II]/H α at fixed [O III]/H β, rather than to a rise in [O III]/H β at fixed [N II]/H α. This interpretation is consistent with the absence of any offset between high- and low-redshift galaxies in the [O III]/H β versus [S II]/H α diagram noted by Shapley et al. (2015, from the analysis of 133 star-forming galaxies at $z \sim 2.3$ with the Keck/MOSFIRE spectrometer), who also find that only low-mass galaxies, with $log(M_{\text{stellar}}/M_{\odot}) < 10$, exhibit higher [N II]/H α at high redshift than at low redshift. The fact that (N/O)gas could be enhanced at fixed (O/H)gas at high redshift has been challenged by Strom et al. (2017), who find that the relation between these two abundance ratios is redshift invariant. Moreover, Kashino et al. (2017) do find an offset in [O III]/H β at fixed [S II]/H α in their high-redshift sample relative to local galaxies. Our theoretical predictions are more in line with these recent studies, since we also find that, at fixed [S II]/H α, massive galaxies (with $M_{\text{stellar}} > 3 \times 10^{10} M_{\odot}$) have on average larger [O III]/H β at high redshift than at low redshift (Fig. 2). Even if enhanced (N/O)gas at fixed (O/H)gas at high redshift (corresponding to higher dust-to-metal mass ratio in our models) can contribute to the cosmic evolution of [O III]/H β and [N II]/H α, our simulations strongly favour this as the dominant process for the offset of the SF branch in this line-ratio diagram at high redshift (Section 5.1.4 and bottom row of Fig. 9), at least for massive galaxies. For less massive galaxies (with $M_{\text{stellar}} < 3 \times 10^{10} M_{\odot}$), we cannot draw any robust conclusion on this point.

(vi) Harder stellar ionizing radiation: Steidel et al. (2014, 2016) and Strom et al. (2017) favour harder stellar ionizing radiation (presumably from metal-poor, massive binaries) as the main driver of the higher typical [O III]/H β ratio of high-redshift star-forming galaxies relative to local ones. Since the Gutkin et al. (2016) SF models used in our analysis do not include any prescription for metal-poor massive binary stars, we cannot draw any robust conclusion regarding the importance of this particular stellar component. Interestingly, Gutkin et al. (2016) reproduce remarkably well the observed ultraviolet and optical emission-line properties of the composite Steidel et al. (2016) spectrum while ignoring any extra component of metal-poor, massive binary stars.<fnref n="10" /> They also show that stellar metallicity, which controls the hardness of stellar ionizing radiation, has only a

<footnote n="9">⁹ These studies find hardly any evidence at high redshift for harder ionizing radiation and enhanced ionization parameter at fixed metallicity.</footnote>

<footnote n="10">¹⁰ Gutkin et al. (2016) did not consider the observed far-ultraviolet stellar-emission continuum spectrum, which Steidel et al. (2016) used to constrain their models, in addition to emission line properties.</footnote>
minor influence on [O III]/H β relative to gas-phase metallicity. The fact that at fixed stellar mass, global interstellar metallicity (which is coupled to stellar metallicity in our approach) has only a negligible impact on [O III]/H β (Section 4.4) suggests that, in our analysis, harder stellar ionizing radiation is unlikely to play a major role in the cosmic evolution of [O III]/H β.

In summary, our theoretical explanation for the cosmic evolution of [O III]/H β is closest to that recently proposed by Kashino et al. (2017), who favour the ionization parameter as the primary driver of this evolution (a minor part of which could come from harder ionizing radiation at high redshift), with potential additional contributions by interstellar metallicity and ionized-gas density.

### 5.3.2 Theoretical studies

So far, only Kewley et al. (2013) investigated theoretically possible reasons for an evolving [O III]/H β ratio using galaxy properties extracted from cosmological hydrodynamic simulations. However, these authors focused exclusively on chemical enrichment histories extracted from simulations, which they injected into nebular-emission models, exploring independently the space of other ISM parameters (hydrogen density, ionization parameter) and an AGN component (central gas metallicity, ionization parameter). The results of Kewley et al. (2013) indicate that the location of the star-forming sequence in the [O III]/H β versus [N II]/H α diagram at any redshift depends mainly on ISM conditions. Specifically, ‘extreme’ ISM conditions at high redshift, driven by (a combination of) higher ionization parameter, harder ionizing radiation and higher electron density, can shift the sequence towards larger [O III]/H β than at low redshift. The location of composite and AGN galaxies also strongly depends on the ISM conditions adopted for the narrow-line region, in particular the inner gas metallicity. The exploratory approach adopted by Kewley et al. (2013) does not allow them to draw any conclusion on the relative influence of these different physical quantities on the cosmic evolution of [O III]/H β and [N II]/H α.

Our simulation results are consistent with the finding of Kewley et al. (2013) that ISM conditions can strongly affect the cosmic evolution of [O III]/H β at fixed stellar mass. Our results further suggest that (specific) SFR – which controls the ionization parameter – is the primary driver of this evolution, with potential minor contributions by hydrogen density and the hardness of stellar ionizing radiation.

## 6 SUMMARY

In the previous sections, we have investigated theoretically the physical origin of the observed cosmic evolution of optical emission-line ratios in galaxies, employing for the first time a self-consistent modelling approach.

Specifically, we compute synthetic [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α ratios for galaxies in a cosmological framework, by coupling – in post-processing – newly developed spectral-evolution models, based on photoionization calculations, with a set of 20 high-resolution cosmological zoom-in simulations of massive galaxies. The latter are performed with the code SPHGal, a modified version of GADGET3, including sophisticated prescriptions for star formation, chemical enrichment (Aumer et al. 2013), stellar feedback (Núñez et al. 2017), BH growth and AGN feedback (Choi et al. 2017).

We include nebular emission from young stars (Gutkin et al. 2016), AGN (Feltre et al. 2016) and post-AGB stars (Section 2.2.3). We adopt direct predictions from our simulations for

the redshift evolution of global and central interstellar metallicity, C/O abundance ratio, SFR, BH accretion rate, global and central average gas densities, and the age and metallicity of post-AGB stellar populations. Based on these, we select SF, AGN and post-AGB nebular-emission models for each galaxy and its most massive progenitor at any redshift. By default, we adopt fixed dust-to-metal mass ratio, ionized-gas hydrogen/electron density and power-law index of AGN ionizing radiation.

We can summarize our main results as follows:

(i) The synthetic [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α emission-line ratios predicted by our simulations are in excellent agreement with observations of both star-forming and active SDSS galaxies in the local universe.

(ii) Towards higher redshifts, at fixed galaxy stellar mass, [O III]/H β is predicted to increase and [N II]/H α, [S II]/H α and [O I]/H α to decrease. These evolutionary trends are consistent with observations by Yabe et al. (2012) and Steidel et al. (2014).

(iii) The physical origin of the cosmic evolution [O III]/H β, [N II]/H α, [S II]/H α and [O I]/H α is a complex mix of different evolving ISM and ionizing-radiation properties governing the nebular emission from young stars, AGN and post-AGB stars.

(iv) When considering the entire sample of simulated galaxies and their main progenitors, interstellar metallicity appears to be a main driver of the cosmic evolution of optical-line ratios, along with (specific) SFR, which controls the ionization parameter. This dominant role of metallicity arises primarily from the intrinsic stellar-mass evolution of simulated galaxies, combined with the correlation between mass and metallicity. Instead, at fixed stellar mass, interstellar metallicity evolves only weakly with redshift and has a negligible influence on the evolution of optical-line ratios.

(v) At fixed stellar mass, the drop in [O III]/H β from high to low redshift in our simulations is driven primarily by that in (specific) SFR, via the ionization parameter. Nebular emission from the growing population of post-AGB stars can also play a minor role in the cosmic evolution of [O III]/H β, but not that from accreting BHs.

(vi) AGN feedback appears to play a key role in the predicted cosmic evolution of [O III]/H β, as test simulations not including AGN feedback exhibit much flatter star formation histories and hardly any redshift evolution of [O III]/H β.

(vii) At fixed stellar mass, the rise in [N II]/H α from high to low redshift follows primarily from the decline in SFR (and ionization parameter), which reduces the probability of multiply ionizing nitrogen at the expense of $N^+$. At redshift $z < 1$, a drop in average central gas density and rise in BHAR/SFR ratio make the contribution by AGN emission contribute more significantly to the cosmic evolution of [N II]/H α. Nebular emission from post-AGB stellar populations hardly affects this evolution.

(viii) At fixed stellar mass, the rise in [S II]/H α and [O I]/H α from high to intermediate redshift ($z \sim 1.5$) also follows from the drop in SFR (and ionization parameter). The continued rise at lower redshift is driven in roughly equal parts by nebular emission from AGN and post-AGB stellar populations.

(ix) Applying observational flux limits to our sample of simulated galaxies and their main progenitors indicates that evolution effects are likely to dominate over flux-selection effects in determining the cosmic evolution of optical-line ratios (although this result might change for a cosmologically representative sample). The ability with our approach to draw such a conclusion is particularly noteworthy, given the difficulty in disentangling these competing effects in emission-line studies of distant galaxies (Juneau et al. 2014).
(x) We have checked that the dust-to-metal mass ratio, ionized-gas hydrogen/electron density and power-law index of AGN ionizing radiation, which are fixed by default in the nebular-emission models, have only a minor (or even negligible) influence on the cosmic evolution of $[O\text{iii}]/H\beta$, $[N\text{ii}]/H\alpha$, $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$. Based on our investigation, we can speculate that adopting a harder ionizing radiation, higher ionized-gas density or higher dust-to-metal mass ratio [which implies higher $(N/O)_{\text{gas}}$ in our self-consistent modelling of metals and their depletion on to dust grains] in high-redshift galaxies relative to local ones may strengthen the redshift dependence of $[O\text{iii}]/H\beta$, but not account for the bulk of the evolution.

The theoretical results presented in this paper provide useful insight into the physical origin of observed cosmic evolution of optical-line ratios. Nevertheless, it is important to keep in mind the sparse statistics of our sample of 20 simulated massive galaxies and their main progenitors, which is likely to affect, for example, the predicted fraction of AGN-dominated galaxies at any given cosmic epoch, and thus, the AGN contribution to the redshift dependence of $[N\text{ii}]/H\alpha$, $[S\text{ii}]/H\alpha$ and $[O\text{i}]/H\alpha$ (but not $[O\text{iii}]/H\beta$). In addition, our sample does not include galaxies with present-day low masses, for which conclusions may be different. The contribution by radiative shocks to the cosmic evolution of optical-line ratios, which we neglected in this study, must also be quantified in detail in future work. This paper is the first in a series. In follow-up studies, we plan to investigate ultraviolet-line diagnostics to help characterize the nature of ionizing radiation in very distant galaxies observed through near-infrared spectroscopy. We also plan to explore the contribution by different ionizing sources to nebular emission in different regions of a galaxy, producing spatially resolved emission-line maps to improve the interpretation of modern integral-field spectroscopic observations in terms of galaxy physical parameters.

## ACKNOWLEDGEMENTS

We thank the referee, A. Inoue, for carefully reading and providing helpful comments on our manuscript. We also thank Stéphanie Juneau and Emma Curtis-Lake for helpful advice and Dan Stark, Bodo Ziegler, Christian Maier and the NEOGAL team for fruitful discussions. MH, SC and AF acknowledge financial support from the European Research Council (ERC) via an Advanced Grant under grant agreement no. 321323-NEOGAL. AF acknowledges support from the ERC via an Advanced Grant under grant agreement no. 339659-MUSICOS. TN acknowledges support from the DFG priority program 1573 'Physics of the interstellar medium' from the DFG Cluster of Excellence 'Origin and structure of the Universe'. RSS is grateful for the generous support of the Downsborough family, and acknowledges support from the Simons Foundation through a Simons Investigator grant.

## REFERENCES

* Agertz O. et al., 2007, MNRAS, 380, 963

* Allen M. G., Groves B. A., Dopita M. A., Sutherland R. S., Kewley L. J., 2008, ApJS, 178, 20

* Andrews B. H., Martini P., 2013, ApJ, 765, 140

* Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013, Earth Planets Space, 65, 213

* Aumer M., White S. D. M., Naab T., Scannapieco C., 2013, MNRAS, 434, 3142

* Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5

* Belfiore F. et al., 2016, MNRAS, 461, 3111

* Bian F., Fan X., McGreer I., Cai Z., Jiang L., 2017, ApJ, 837, L12

* Bondi H., 1952, MNRAS, 112, 195

* Brinchmann J., Pettini M., Charlot S., 2008, MNRAS, 385, 769

* Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

* Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchibergmann T., 2000, ApJ, 533, 682

* Carton D. et al., 2017, MNRAS, 468, 2140

* Chabrier G., 2003, PASP, 115, 763

* Charlot S., Longhetti M., 2001, MNRAS, 323, 887

* Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125

* Choi E., Ostriker J. P., Naab T., Oser L., Moster B. P., 2015, MNRAS, 449, 4105

* Choi E., Ostriker J. P., Naab T., Somerville R. S., Hirschmann M., Núñez A., Hu C.-Y., Oser L., 2017, ApJ, 844, 31

* Cowie L. L., Barger A. J., Songaila A., 2016, ApJ, 817, 57

* Cullen L., Dehnen W., 2010, MNRAS, 408, 669

* Cullen F., Cirasuolo M., Kewley L. J., McLure R. J., Dunlop J. S., Bowler R. A. A., 2016, MNRAS, 460, 3002

* Daddi E. et al., 2007, ApJ, 670, 156

* de Barros S. et al., 2016, A&A, 585, A51

* De Cia A., Ledoux C., Mattsson L., Petitjean P., Srianand R., Gavignaud I., Jenkins E. B., 2016, A&A, 596, A97

* Dehnen W., Aly H., 2012, MNRAS, 425, 1068

* Dopita M. A., 2003, PASP, 20, 31

* Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

* Durier F., Dalla Vecchia C., 2012, MNRAS, 419, 465

* Elbaz D. et al., 2007, A&A, 468, 33

* Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354

* Ferland G. J. et al., 2013, RMxAA, 49, 137

* Genel S. et al., 2014, MNRAS, 445, 175

* Genzel R. et al., 2011, ApJ, 733, 101

* Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76

* Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757

* Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays Modelling the UV/X-ray Cosmic Background with CUBA

* Hainline K. N., Shapley A. E., Kornei K. A., Pettini M., Buckley-Geer E., Allam S. S., Tucker D. L., 2009, ApJ, 701, 52

* Hainline K. N., Hickox R. C., Greene J. E., Myers A. D., Zakamska N. L., Liu G., Liu X., 2014, ApJ, 787, 65

* Hasinger G., 2008, A&A, 490, 905

* Hayashi M. et al., 2015, PASJ, 67, 80

* Hirschmann M. et al., 2013, MNRAS, 436, 2929

* Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304

* Hirschmann M., De Lucia G., Fontanot F., 2016, MNRAS, 461, 1760

* Hopkins P. F., 2013, MNRAS, 428, 2840

* Hopkins P. F., Cox T. J., Hernquist L., Narayanan D., Hayward C. C., Murray N., 2013, MNRAS, 430, 1901

* Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, MNRAS, 443, 1173

* Inoue A. K., 2003, PASJ, 55, 901

* Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439

* Izotov Y. I., Thuan T. X., 1999, ApJ, 511, 639

* Izotov Y. I., Thuan T. X., Guseva N. G., 2017, MNRAS, 471, 548

* Jaskot A. E., Oey M. S., 2013, ApJ, 766, 91

* Juneau S. et al., 2014, ApJ, 788, 88

* Karakas A. I., 2010, MNRAS, 403, 1413

* Kashino D. et al., 2017, ApJ, 835, 88

* Kauffmann G. et al., 2003, MNRAS, 346, 1055

* Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

* Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, ApJ, 556, 121

* Kewley L. J., Dopita M. A., Leitherer C., Davé R., Yuan T., Allen M., Groves B., Sutherland R., 2013, ApJ, 774, 100

* Kimura M. et al., 2010, PASJ, 62, 1135
Kobulnicky H. A., Kennicutt R. C., Jr, Pizagno J. L., 1999, ApJ, 514, 544

Kornei K. A., Shapley A. E., Martin C. L., Coil A. L., Lotz J. M., Schiminovich D., Bundy K., Noeske K. G., 2012, ApJ, 758, 135

Lehnert M. D., Nesvadba N. P. H., Le Tiran L., Di Matteo P., van Driel W., Douglas L. S., Chemin L., Bournaud F., 2009, ApJ, 699, 1660

Li C., Kauffmann G., Heckman T. M., White S. D. M., Jing Y. P., 2008, MNRAS, 385, 1915

Magorrian J. et al., 1998, AJ, 115, 2285

Maier C., Ziegler B. L., Lilly S. J., Contini T., Pérez-Montero E., Lamareille F., Bolzonella M., Le Floc'h E., 2015, A&A, 577, A14

Maiolino R. et al., 2008, A&A, 488, 463

Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

Masters D. et al., 2014, ApJ, 785, 153

Masters D., Faisst A., Capak P., 2016, ApJ, 828, 18

McLean I. S. et al., 2010, in Ground-based and Airborne Instrumentation for Astronomy III Vol. 7735 of Proc. SPIE, Design and development of MOSFIRE: the Multi-Object Spectrometer for Infrared Exploration at the Keck Observatory. pp. 77351E-77351E-12.

Morisset C. et al., 2016, A&A, 594, A37

Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903

Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

Nagao T., Maiolino R., Marconi A., 2006, A&A, 459, 85

Newman S. F. et al., 2012, ApJ, 761, 43

Núñez A., Ostriker J. P., Naab T., Oser L., Hu C.-Y., Choi E., 2017, ApJ, 836, 204

Orsi Á., Padilla N., Groves B., Cora S., Tecce T., Gargiulo I., Ruiz A., 2014, MNRAS, 443, 799

Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, 2nd edn., Vol. 35. Univ. Science Books, p. 40

Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642

Panuzzo P., Bressan A., Granato G. L., Silva L., Danese L., 2003, A&A, 409, 99

Pettini M., Pagel B. E. J., 2004, MNRAS, 348, L59

Popping G., Somerville R. S., Galametz M., 2017, MNRAS, 471, 3152

Read J. I., Hayfield T., 2012, MNRAS, 422, 3037

Rémy-Ruyer A. et al., 2014, A&A, 563, A31

Rich J. A., Dopita M. A., Kewley L. J., Rupke D. S. N., 2010, ApJ, 721, 505

Rich J. A., Kewley L. J., Dopita M. A., 2011, ApJ, 734, 87

Rigby J. R., Wuyts E., Gladders M. D., Sharon K., Becker G. D., 2011, ApJ, 732, 59

Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743

Saitoh T. R., Makino J., 2009, ApJ, 697, L99

Saitoh T. R., Makino J., 2013, ApJ, 768, 44

Sanders R. L. et al., 2016, ApJ, 816, 23

Sazonov S. Y., Ostriker J. P., Ciotti L., Sunyaev R. A., 2005, MNRAS, 358, 168

Schaye J. et al., 2015, MNRAS, 446, 521

Shapley A. E., Steidel C. C., Erb D. K., Reddy N. A., Adelberger K. L., Pettini M., Barmby P., Huang J., 2005, ApJ, 626, 698

Shapley A. E. et al., 2015, ApJ, 801, 88

Shapley A. E., Steidel C. C., Strom A. L., Bogosavljević M., Reddy N. A., Siana B., Mostardi R. E., Rudie G. C., 2016, ApJ, 826, L24

Sharp R. G., Bland-Hawthorn J., 2010, ApJ, 711, 818

Shimizu I., Inoue A. K., Okamoto T., Yoshida N., 2016, MNRAS, 461, 3563

Singh R. et al., 2013, A&A, 558, A43

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Soto K. T., Martin C. L., Prescott M. K. M., Armus L., 2012, ApJ, 757, 86

Spergel D. N. et al., 2003, ApJS, 148, 175

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Stasińska G., 1980, A&A, 85, 359

Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289

Steidel C. C. et al., 2014, ApJ, 795, 165

Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R. F., 2016, ApJ, 826, 159

Steinborn L. K., Dolag K., Hirschmann M., Prieto M. A., Remus R.-S., 2015, MNRAS, 448, 1504

Steinborn L. K., Dolag K., Comerford J. M., Hirschmann M., Remus R.-S., Teklu A. F., 2016, MNRAS, 458, 1013

Stinson G. S., Brook C., Macciò A. V., Wadsley J., Quinn T. R., Couchman H. M. P., 2013, MNRAS, 428, 129

Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS

Strom A. L., Steidel C. C., Rudie G. C., Trainor R. F., Pettini M., Reddy N. A., 2017, ApJ, 836, 164

Vanzella E. et al., 2016, ApJ, 821, L27

Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295

Vincenzo F., Belfiore F., Maiolino R., Matteucci F., Ventura P., 2016, MNRAS, 458, 3466

Weinberger R. et al., 2017, MNRAS, 465, 3291

Weistrop D., Nelson C. H., Angione R., Bachilla R., Hancock M., Kaiser M. E., 2012, AJ, 143, 98

Wiseman P., Schady P., Bolmer J., Krühler T., Yates R. M., Greiner J., Fynbo J. P. U., 2017, A&A, 599, A24

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Wright S. A., Larkin J. E., Graham J. R., Ma C.-P., 2010, ApJ, 711, 1291

Yabe K. et al., 2012, PASJ, 64, 60

Yabe K. et al., 2014, MNRAS, 437, 3647

Yabe K., Ohta K., Akiyama M., Iwamura F., Tamura N., Yuma S., Dalton G., Lewis I., 2015, ApJ, 798, 45

## APPENDIX A: THE ROLE OF AGN FEEDBACK IN THE COSMIC EVOLUTION OF [O III]/H β

Throughout this study, we have investigated predictions from cosmological zoom-in simulations including a prescription for AGN feedback. Previous work (e.g. Choi et al. 2017) has shown that AGN feedback (and in particular our specific implementation of this process) can strongly affect galaxy properties, such as SFR and central and global gas densities, which control nebular emission from stars and AGN. Thus, we may expect a priori AGN feedback to be imprinted in the nebular emission from galaxies. Our aim in this appendix is to test this hypothesis, focusing on the evolutionary trends of optical-line ratios.

To achieve this, we performed a new suite of zoom-in simulations of 20 massive haloes based on the same initial conditions as described in Section 2.1.2, but switching off BH growth and AGN feedback in our simulation code SPHGal. By design, this set of re-simulated galaxies does not include any AGN contribution to nebular emission. As [N II]/H α, [S II]/H α and [O I]/H α are known to be directly affected by the presence of nebular emission from an AGN, unlike [O III]/H β (Fig. 8), it is most meaningful here to examine the impact of AGN feedback on the cosmic evolution of [O III]/H β in SF-dominated galaxies.

Fig. A1 shows the analogue of Fig. 3, i.e. the evolution of the average [O III]/H β of SF galaxies in bins of [N II]/H α (left-hand panel), [S II]/H α (middle panel) and [O I]/H α (right-hand panel), for the new simulations without AGN feedback. Compared to Fig. 3 (based on simulations including AGN feedback), the drop in [O III]/H β from high to low redshift is significantly reduced at given [S II]/H α and [O I]/H α. At fixed [N II]/H α, we even find a slightly reversed trend of higher [O III]/H β in local than distant galaxies.

To understand the origin for the hardly evolving SF branch in these optical diagnostic diagrams, we plot in Fig. A2 the analogue of the first two rows of Fig. 6, i.e. the redshift evolution of the different physical quantities used to select SF nebular-emission models, for the 20 massive galaxies and their main high-redshift progenitors
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Same as Fig. 3, but for the simulation set of galaxies *without* AGN feedback, as described in Appendix A.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure A2.** Same as the first two rows of Fig. 6, but for the simulation set of galaxies *without* AGN feedback, as described in Appendix A.</figcaption></figure>

simulated *without* AGN feedback. In contrast to the predictions including AGN feedback, where the drop in (specific) SFR has been identified as the main driver for the cosmic evolution of [O III]/H β, Fig. A2 shows much flatter histories of star formation, and hence, ionizing-photon production ($U_{sim,*}$), whether galaxies are selected by mass or not (thick black solid and dashed lines). In fact, at $z=0$, all massive galaxies in Fig. A2 are still highly star-forming, with $\log(\text{sSFR/yr}^{-1}) > -10$. This is because of the lack of mechanical and radiative AGN feedback, which, when present, can very efficiently heat and expel cold, star-forming gas from massive galaxies, thereby reducing late *in situ* star formation (see e.g. Choi et al. 2015, 2017). As a result, without AGN feedback, [O III]/H β in local massive galaxies is predicted to be typically as high as in distant galaxies.

From this analysis, we conclude that the drop in [O III]/H β from high to low redshift in our sample of simulated massive galaxies is ultimately caused by AGN feedback being the main responsible factor for the strong decrease in (specific) SFR, and thus, in the SF ionization parameter.

This paper has been typeset from a TeX/LATeX file prepared by the author.