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

TL;DR

This study models nebular emission in massive galaxies across cosmic time, revealing how galaxy properties and stellar populations influence the evolution of optical emission-line ratios, aligning well with observational data.

Contribution

It introduces a self-consistent modeling approach coupling nebular emission with cosmological galaxy simulations, elucidating drivers of emission-line ratio evolution.

Findings

Simulated line ratios match local galaxy observations.

[OIII]/Hβ increases with redshift at fixed stellar mass.

Star formation history is the main driver of line ratio evolution.

Abstract

Galaxies occupy different regions of the [OIII]$\lambda5007$/H$\beta$-versus-[NII]$\lambda6584$/H$\alpha$ 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 cosmic evolution of the optical-line ratios [OIII]$\lambda5007$/H$\beta$, [NII]$\lambda6584$/H$\alpha$, [SII]$\lambda\lambda6717,6731$/H$\alpha$ and [OI]$\lambda6300$/H$\alpha$. The line ratios of simulated galaxies agree well with observations of both star-forming and active local SDSS galaxies. Towards higher redshifts, at fixed galaxy stellar mass, the average [OIII]/H$\beta$ increases and [NII]/H$\alpha$, [SII]/H$\alpha$ and [OI]/H$\alpha$ 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 [OIII]/H$\beta$ and [NII]/H$\alpha$. For [SII]/H$\alpha$ and [OI]/H$\alpha$, this applies only to redshifts above $z=1.5$, the evolution at lower redshift being driven in roughly equal parts by nebular emission from AGN and post-AGB stars. 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.

Tables

Table 1. 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 ( α = − 1.7 𝛼 1.7 \alpha=-1.7 ), fixed dust-to-metal mass ratio ( ξ d = 0.3 subscript 𝜉 d 0.3 \xi_{\mathrm{d}}=0.3 ), and fixed ionized-gas density ( n H , ⋆ = 10 2 subscript 𝑛 H ⋆ superscript 10 2 n_{\mathrm{H},\star}=10^{2}\, cm -3 , n H , ∙ = 10 3 subscript 𝑛 H ∙ superscript 10 3 n_{\mathrm{H},\bullet}=10^{3}\, cm -3 and n H , ⋄ = 10 subscript 𝑛 H ⋄ 10 n_{\mathrm{H},\diamond}=10\, 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.

Parameter space	SF models
(Gutkin et al. 2016)	AGN models
(Feltre et al. 2016)	PAGB models
(this work)
Ionizing spectrum
(matched/fixed)	10 Myr-old stellar population with const SFR (fixed),
stellar metallicity same as that of gas (matched)	UV slope $\alpha$ = $-1.2$, $-1.4$, $-1.7$, $-2.0$
(fixed)	3, 5, 7, 9 Gyr-old stellar populations
(matched)
$Z_{\diamond ,\mathrm{stars}}$ = 0.008, 0.014, 0.017, 0.02 (matched)
Interstellar metallicity
$Z$ (matched)	$Z_{\star }=$
0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.008, 0.014, 0.017, 0.02, 0.03	$Z_{\bullet }=$
0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.008, 0.014, 0.017, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07	$Z_{\diamond }=$
0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.008, 0.014, 0.017, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07
Ionization parameter $\log U$, function of the average gas density
(matched)	$\log U_{\star }=-0.65,-1.15,-1.65,-2.15,$
$-2.65,-3.15,-3.65$	$\log U_{\bullet }=-0.65,-1.15,-1.65,-2.15,$
$-2.65,-3.15,-3.65,-4.65$	$\log U_{\diamond }=-2.15,-2.65,-3.15,-3.65,$
$-4.15,-4.65$
Dust/metal mass ratio
${\xi }_d$ (fixed)	0.1, 0.3, 0.5	0.1, 0.3, 0.5	0.1, 0.3, 0.5
Ionized-gas density $\log \left(n_H/{\mathrm{cm}}^3\right)$
(fixed)	$\log \left(n_{H,\star }\right)$=2.0, 3.0, 4.0	$\log \left(n_{H,\bullet }\right)$= 2.0, 3.0, 4.0	$\log \left(n_{H,\diamond }\right)$= 1.0, 2.0, 3.0
C/O abundance ratio in solar units
(matched)	(C/O)⋆/(C/O)${}_{\odot }{}^{}=0.1,0.2,0.27,0.38,0.52,0.72,1.0$	(C/O)∙/(C/O)${}_{\odot }{}^{}=0.1,0.2,0.27,0.38,0.52,0.72,1.0$	(C/O)⋄/(C/O)${}_{\odot }{}^{}=1.0$
Model normalization
(matched)	Star formation rate
SFR	AGN luminosity
$L_{\mathrm{AGN}}$	Mass of evolved stars
$M_{\diamond ,\mathrm{stellar}}$

Equations

$$U_{\mathrm{sim},\star }=\frac{3{\alpha }_B^{2/3}}{4c}{\left(\frac{3Q_{\mathrm{sim},\star }{ϵ}^2{n}_{\mathrm{H},\star }}{4\pi }\right)}^{1/3}.$$

$$ϵ=\frac{n_{\mathrm{gas},\mathrm{glob}}}{n_{\mathrm{H},\star }},$$

Full text

††pagerange: Synthetic nebular emission from massive galaxies I: origin of the cosmic evolution of optical emission-line ratios–LABEL:lastpage††pubyear: 2002

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 et al., 1999; Kauffmann et al., 2003; Nagao et al., 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$\alpha$, H$\beta$, [O i]$\lambda 6300$, [O ii]$\lambda\lambda 3726,3729$, [O iii]$\lambda 5007$, [N ii]$\lambda 6584$ and [S ii]$\lambda\lambda 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 et al. (1981, hereafter BPT) and Veilleux & Osterbrock (1987), relate the [O iii]/H$\beta$ ratio to the [N ii]/H$\alpha$, [S ii]/H$\alpha$ and [O i]/H$\alpha$ 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>1$ have systematically larger [O iii]/H$\beta$ ratio, at fixed [N ii]/H$\alpha$ 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 et al., 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 et al., 2016). Instead, Steidel et al. (2014, 2016, see also ) 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$\beta$ versus [N ii]/H$\alpha$ BPT diagram, as well as the potential influence of an AGN.111Kewley 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. Kewley et al. (2013) find that the SF sequence can be shifted to higher [O iii]/H$\beta$ 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 the present 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., 2016; 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 onto dust grains (Gutkin et al. 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 et al. (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 the present study: can we account for the observed evolution of optical emission-line ratios, in particular the systematically larger [O iii]/H$\beta$ ratio of high-redshift galaxies at fixed [N ii]/H$\alpha$ 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}\hbox{M{}{\odot}}\,h^{-1} and 2.7\times 10^{13}\hbox{M{}{\odot}}\,h^{-1} ($H_{0}=100\,h\,\mathrm{km\,s}^{-1}$). We performed these simulations with a modified version of the highly parallel, smoothed particle hydrodynamics (SPH) code GADGET3 (Springel et al., 2005), SPHGal (Hu et al., 2014, see also Choi et al. 2016; 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. (2016), 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 C4 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 photo-ionization background. Specifically, star formation and chemical evolution is modelled as described in Aumer et al. (2013), assuming chemical enrichment via type-Ia and type-II supernovae (SNe) 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_{\mathrm{th}}=n_{0}(T_{\mathrm{gas}}/T_{0})^{3}(M_{0}/M_{\mathrm{gas}})^{2}$, where $T_{\mathrm{gas}}$ and $M_{\mathrm{gas}}$ are the temperature and mass of the gas particle, and $n_{0}=2\,\mathrm{cm}^{-3}$ and $T_{0}=30000\,$K (see Section 2.1.2). Gas particles with densities above $n_{\mathrm{th}}$ are Jeans unstable. Their star formation rate is calculated as $d\rho_{*}/dt=\eta\rho_{\mathrm{gas}}/t_{\mathrm{dyn}}$, where $\rho_{*}$, $\rho_{\mathrm{gas}}$ and $t_{\mathrm{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_{\mathrm{out,SN}}=4500\,$km s*-1*).222This is a simplified version of the full, 3-phase blast-wave model adopted in Núñez et al. (2017) and Choi et al. (2016). 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_{\mathrm{out,AGB}}=10\,$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}\hbox{M{}{\odot}} being placed at the density minimum of any dark-matter halo whose mass exceeds 10^{11}\hbox{M{}{\odot}}.333These 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). 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 onto 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 Naab & Ostriker 2016; Somerville & Davé 2015). Instead, we rely on a more physically motivated approach including both mechanical and radiative feedback (Ostriker et al., 2010; Choi et al., 2016). 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 km s*-1* perpendicular to the gaseous disk. 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, 2016) for more details about AGN feedback modelling. We note that our set of 20 zoom-in simulations do not include any metallicity-dependent heating prescription. This is justified by the fact that, as shown by Choi et al. (2016), such refinements are not found to have any significant impact on basic properties of massive galaxies.

Choi et al. (2016) 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 co-moving periodic box length $L=100\ \mathrm{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_{\mathrm{m}}=0.26$, $\Omega_{\Lambda}=0.74$ 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_{\mathrm{DM}}=2\times 10^{8}\hbox{M{}_{\odot}}\,h^{-1} and a fixed co-moving gravitational softening length of $2.52\ h^{-1}\mathrm{kpc}$. We refer the reader to the original papers of Oser et al. (2010, 2012) for more details about the simulation setup.

From this simulation, Oser et al. (2010) selected 39 haloes with masses in the range $7\times 10^{11}$–2.7\times 10^{13}\hbox{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 as well as gas particles at higher resolution ($\Omega_{\mathrm{b}}=0.044,\Omega_{\mathrm{dm}}=0.216$). The new dark matter particles have a mass m_{\mathrm{dm}}=2.5\times 10^{7}\hbox{M{}{\odot}}\,h^{-1}, i.e., 8 times smaller than the original ones, while the gas particles have a mass m_{\mathrm{gas}}=4.2\times 10^{6}\hbox{M{}{\odot}}\,h^{-1}, equal to that of star particles. The co-moving gravitational softening length of the dark matter particles is $890\ h^{-1}\mathrm{pc}$, and that of the gas and star particles $400\ h^{-1}\mathrm{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}\hbox{M{}{\odot}} and 3\times 10^{13}\hbox{M{}{\odot}}, and associated central galaxy masses (computed as the stellar mass within a tenth of the virial radius) between 3\times 10^{10}\hbox{M{}{\odot}} and 3\times 10^{11}\hbox{M{}_{\odot}} ($h=0.72$).