Nebular Continuum and Line Emission in Stellar Population Synthesis Models

To calculate photoionization models, we use version 13.03 of Cloudy,⁴ last described by Ferland et al. (2013). Cloudy simulates physical conditions within a gas cloud to predict the thermal, ionization, and chemical structure of the cloud, and the resultant spectrum of the diffuse emission. Users must describe the physical properties of the gas cloud and provide an external source of radiation to photoionize the cloud. Each model must specify (1) the geometry and (2) the chemical content of the gas cloud, and (3) the spectrum and (4) the intensity of the ionizing radiation source. Our choices for these parameters are as follows.

2.1.1. Assumed Geometry

2.1.2. Gas Chemical Content

We adopt the gas phase abundances specified by Dopita et al. (2000), which are based on the solar abundances from Anders and Grevesse (1989). We assume that the gas phase metallicity scales with the metallicity of the stellar population, given that the metallicity of the most massive stars should be identical to the metallicity of the gas cloud from which the stars formed. Metal abundances are solar-scaled, with the exception of nitrogen, which has a known secondary nucleosynthetic contribution. For the scaling of nitrogen with metallicity, we follow the piecewise relationship between nitrogen and oxygen specified by Dopita et al. (2000).

Elements like carbon and nitrogen are observed to be heavily depleted onto dust grains in H ii regions. This alters the chemical composition of the nebula, but also the thermal properties of the nebula, since these elements are important gas coolants. To account for this, the relative abundances used in this work include the effect of depletion onto dust grains derived from observations, as defined in Dopita et al. (2000). The applied depletion factors are constant factors applied to specific elements and do not scale with the metallicity of the model. The depletion factors are applied regardless of whether the model nebula includes dust grains, due to their important effect on the temperature structure of the cloud. A complete description of the elemental abundances and depletion factors used in this work can be found in Table 1. In Table 2 we compare the abundances used in this work with those used in other nebular emission models.

Table 1.  Solar-metallicity (${{\rm{Z}}}_{\odot }$) and Depletion Factors (D) Adopted for Each Element

Element	${\mathrm{log}}_{10}Z/{Z}_{\odot }$	log (D)
H	0	0
He	−1.01	0
C	−3.44	−0.30
N	−3.95	−0.22
O	−3.07	−0.22
Ne	−3.91	0
Mg	−4.42	−0.70
Si	−4.45	−1.0
S	−4.79	0
Ar	−5.44	0
Ca	−5.64	−2.52
Fe	−4.33	−2.0

Note. Solar abundances are from Anders & Grevesse (1989) and depletion factors are from Dopita et al. (2000).

Download table as:  ASCIITypeset image

Table 2.  Abundances of Important Gas Coolants in Various Nebular Emission Models

Abundance Set	Solar Abundance Set	log (C/H)	log (D)	log (N/H)	log (D)	log (O/H)	log (D)
Cloudy $\langle$ Orion Nebula $\rangle$	Anders & Grevesse (1989)	−3.52	⋯	−4.15	⋯	−3.40	⋯
Dopita et al. (2013)	Grevesse et al. (2010)	−3.87	(−0.30)	−4.65	(−0.05)	−3.38	(−0.07)
Dopita et al. (2000)	Anders & Grevesse (1989)	−3.74	(−0.30)	−4.17	(−0.22)	−3.29	(−0.22)
Charlot & Longhetti (2001)	Grevesse & Noels (1993)	−3.45	⋯	−4.30	(−0.27)	−3.18	(−0.05)
Levesque et al. (2010)	Anders & Grevesse (1989)	−3.70	⋯	−4.22	⋯	−3.29	⋯

Note. Values in this table reflect absolute abundance at solar metallicity and include the indicated depletion factors.

Download table as:  ASCIITypeset image

We produce two separate nebular models: one that includes grains within the nebula and one that does not. The majority of the models discussed hereafter do not include dust grains within the nebula. For the dusty models, we use a dust grain model with a size distribution and abundances pattern appropriate for the ISM of the Milky Way. The grain prescription includes both a graphite and silicate component and generally reproduces the observed extinction properties for a ratio of extinction per reddening of ${R}_{{\rm{v}}}\equiv {A}_{{\rm{v}}}/E(B-V)=3.1$. We note that in real galaxies, ${R}_{{\rm{v}}}$ will not necessarily equal the canonical Milky Way value, and depletion patterns may differ from the ones adopted in this work, which could in turn significantly alter the physical properties of the model H ii region.

2.1.3. Ionizing Spectra

2.1.4. Ionizing Spectrum Intensity

For photoionization models, we care primarily about the intensity of the ionizing spectrum blueward of $912\,{\rm{\mathring{\rm A} }}$, the ionization threshold for hydrogen. We quantify the intensity of each ionizing spectrum with ${Q}_{{\rm{H}}}$, the total number of photons emitted per second that are capable of ionizing hydrogen:

$$\begin{eqnarray}&&{Q}_{{\rm{H}}}\equiv {\int }_{{\nu }_{0}}^{\infty }\displaystyle \frac{{f}_{\nu }}{h\nu }d\nu =\displaystyle \frac{1}{{hc}}{\int }_{0}^{{\lambda }_{0}}\lambda {f}_{\lambda }d\lambda .\end{eqnarray} \tag{ 1 }$$

${Q}_{{\rm{H}}}$ can be calculated directly from the SSP ionizing spectra from FSPS, which assume the formation of one solar mass of new stars. Over the full range of SSP ages and metallicities used in this work, ${Q}_{{\rm{H}}}$ calculated from the $1\,{M}_{\odot }$ spectra varies from 10⁴³–10⁴⁷ (s⁻¹). However, $1\,{M}_{\odot }$ SSPs do not provide enough ionizing photons to reproduce the ionizing photon rates observed in massive H ii regions. A single O7V star with ${M}_{{\rm{i}}}\sim 25\,{M}_{\odot }$ produces ${Q}_{{\rm{H}}}\sim {10}^{48.75}$, while a young star cluster with $M\sim 3\times {10}^{4}\,{M}_{\odot }$ produces ${Q}_{{\rm{H}}}\sim {10}^{51}$. Thus to achieve values of ${Q}_{{\rm{H}}}$ consistent with observed massive H ii regions, we must increase the intensity of the ionizing spectra by increasing the total stellar mass of the SSPs. We first describe how we specify the intensity of the ionizing spectra within Cloudy, and then describe the relation to SSP stellar mass.

When running Cloudy, we set the intensity of the ionizing spectrum using the ionization parameter, ${ \mathcal U }$, a dimensionless quantity that gives the ratio of hydrogen ionizing photons to total hydrogen density,

$$\begin{eqnarray}&&{ \mathcal U }\equiv \displaystyle \frac{{Q}_{{\rm{H}}}}{4\pi {R}^{2}\cdot {n}_{{\rm{H}}}\cdot c},\end{eqnarray} \tag{ 2 }$$

where R is radius of the ionized region, ${n}_{{\rm{H}}}$ is the number density of hydrogen (${\mathrm{cm}}^{-3}$), and c is the speed of light. The ionization parameter used in this work, ${ \mathcal U }$, differs from the ionization parameter q ($\mathrm{cm}\cdot {{\rm{s}}}^{-1}$) used by Levesque et al. (2010) and Dopita et al. (2013) by a factor of c, the speed of light: $q={ \mathcal U }\cdot c$.

In its derivation, ${ \mathcal U }$ is computed at $R={R}_{{\rm{S}}}$, the Strömgren radius, the location where ionization and recombination rates are balanced in thermal equilibrium, which can only be calculated after a photoionization model is computed. Cloudy defines ${ \mathcal U }$ at $R={R}_{\mathrm{inner}}$, the distance from the ionizing source to the illuminated inner face of the cloud. To avoid confusion, we define ${{ \mathcal U }}_{0}$ as the ionization parameter calculated at the inner radius of the gas cloud. The distinction does not matter for a thin spherical shell, the geometry assumed in this work (for details, see Charlot & Longhetti 2001).

${{ \mathcal U }}_{0}$ conveniently folds in both the intensity of the ionizing source (${Q}_{{\rm{H}}}$) and the geometry of the gas cloud (${R}_{\mathrm{inner}}$, ${n}_{{\rm{H}}}$), which allows us to make a simplification to reduce the dimensionality of our model grid. For a fixed EUV shape and metallicity, any combination of ${Q}_{{\rm{H}}}$, ${R}_{\mathrm{inner}}$, and ${n}_{{\rm{H}}}$ that produces the same value of ${{ \mathcal U }}_{0}$ will produce the same nebular spectrum, a simplification that holds at typical H ii region densities and sizes (${n}_{{\rm{H}}}\sim 10\mbox{--}1000$, ${R}_{\mathrm{inner}}\sim 0.1\mbox{--}10$ pc).

Each model is run through Cloudy at seven different ionization parameters, $-4\leqslant {\mathrm{log}}_{10}{{ \mathcal U }}_{0}\leqslant -1$, in steps of 0.5, a range consistent with ionization parameters observed in local starburst galaxies (Rigby & Rieke 2004). To produce the various values of ${{ \mathcal U }}_{0}$, we choose to fix ${R}_{\mathrm{inner}}$ and ${n}_{{\rm{H}}}$ and only vary ${Q}_{{\rm{H}}}$. Other groups have taken different approaches: Moy et al. (2001) fix ${R}_{\mathrm{inner}}$, ${n}_{{\rm{H}}}$, and ${Q}_{{\rm{H}}}$, and produce different ionization parameters by varying the filling factor of the surrounding gas cloud.⁶ Levesque et al. (2010) fix the total mass of the instantaneous bursts at ${10}^{6}\,{M}_{\odot }$ and vary the inner radius of the cloud to produce different values of ${{ \mathcal U }}_{0}$.

We note that running models at different values of ${{ \mathcal U }}_{0}$ and fixed ${R}_{\mathrm{inner}}$ and ${n}_{{\rm{H}}}$ by varying ${Q}_{{\rm{H}}}$ implicitly varies the effective mass of the SSPs. Older SSPs produce fewer ionizing photons and must be more massive than a younger SSP to achieve the same value of ${{ \mathcal U }}_{0}$. To distinguish the ${Q}_{{\rm{H}}}$ input to Cloudy from the one measured from the $1\,{M}_{\odot }$ ionizing spectrum, we define ${\hat{Q}}_{{\rm{H}}}$ as the ionizing photon rate per solar mass: ${\hat{Q}}_{{\rm{H}}}\equiv {Q}_{{\rm{H}}}/{M}_{\odot }$. To account for the different intensities required to produce the desired range in ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}$, we simply normalize the output from Cloudy by the total stellar mass, ${\hat{Q}}_{{\rm{H}}}$/${Q}_{{\rm{H}}}$.

The range in total stellar mass required to produce ionization parameters in the range of $-4\leqslant {\mathrm{log}}_{10}{{ \mathcal U }}_{0}\leqslant -1$ varies from 10¹ to ${10}^{8}\,{M}_{\odot }$, depending on the age and metallicity of the SSP. For a 1 Myr solar-metallicity SSP, $-4\leqslant {\mathrm{log}}_{10}{{ \mathcal U }}_{0}\leqslant -1$ implies variation in stellar mass from 10¹ to ${10}^{4}\,{M}_{\odot }$ at the assumed ${R}_{\mathrm{inner}}={10}^{19}$ cm and ${n}_{{\rm{H}}}=100$ ${\mathrm{cm}}^{-3}$ of our model. At 7 Myr, the oldest SSP used in our model, this requires larger masses, 10⁵ to ${10}^{8}\,{M}_{\odot }$.

2.1.5. Other Model Specifications

The nebular emission model samples the following values for SSP age, SSP metallicity, and ionization parameter, ${{ \mathcal U }}_{0}$:

• ${\mathrm{log}}_{10}{{\boldsymbol{ \mathcal U }}}_{0}$: −4.0, −3.5, −3.0, −2.5, −2.0, −1.5, −1.0
• ${\mathrm{log}}_{10}Z/{Z}_{\odot }$: −2.0, −1.5, −1.0, −0.6, −0.4, −0.3, −0.2, −0.1, 0.0, 0.1, 0.2
• Age: 0.5, 1, 2, 3, 4, 5, 6, 7, 10 million years (Myr)

The nebular model has one "free" parameter, ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}$. For each SSP of age t and metallicity Z, we run photoionization models at each ionization parameter, ${{ \mathcal U }}_{0}$. We normalize the line and continuum emission by ${Q}_{{\rm{H}}}$ input to Cloudy. The normalized line and continuum emission are recorded into separate look-up tables.

For a given SSP $(t,Z)$ and specified ${{ \mathcal U }}_{0}$, FSPS returns the associated line and continuum emission associated with that grid-point from the look-up table. This maintains the model self-consistency, such that the nebular emission is added to the same spectrum that was used to ionize the gas cloud.

The normalized line and continuum emission are rescaled by ${Q}_{{\rm{H}}}$ calculated from the SSP's spectrum; FSPS then removes the ionizing photons from the SED to enforce energy balance. FSPS includes a parameter, frac_obrun, which allows some fraction of the ionizing luminosity to escape from the H ii region. However, in this work we assume an ionizing photon escape fraction of zero.

The nebular model is implemented for SSPs with age $t\lt {t}_{\mathrm{esc}}$, a parameter within FSPS that specifies the time a given SSP is surrounded by its birth cloud. For complex stellar populations, the total nebular contribution is the sum of emission from all SSPs that contribute to the SFH with age $t\lt {t}_{\mathrm{esc}}$. In practice, the emission from H ii regions surrounding star clusters is a relatively short-lived phenomenon ($\lesssim {10}^{7}$ years), with most of the contribution coming from SSPs 3 Myr and younger, though deviations from this are discussed in Section 6.1. In general, we neglect the contribution from old planetary nebula and hot intermediate-aged stars. In Section 6.2, however, we assess the importance of the contribution from post-AGB stars. We do not consider the contribution from AGN.

For giant H ii regions, the ionization is provided by groups of rapidly evolving massive stars. The ionizing spectrum will thus change with time, with stars of different masses dominating the spectrum at different times. For H ii regions, we are primarily concerned with the evolution of flux at wavelengths shorter than $912\,{\rm{\mathring{\rm A} }}$, as these are the photons with energies high enough to ionize hydrogen in the surrounding gas cloud.

Figure 1 shows the relative flux contribution of stars in different mass ranges to the ionizing spectrum at solar metallicity as a function of time. In all panels, the radiation capable of ionizing hydrogen is produced by the most massive stars: those with initial masses $\gtrsim 10\,{M}_{\odot }$. At 1 Myr, stars with initial mass $\gtrsim 35\,{M}_{\odot }$ dominate the ionizing spectrum. These stars are extremely short-lived, and by 3 Myr the spectrum is dominated by stars with initial masses $\gtrsim 25\,{M}_{\odot }$. O-type stars ($\gt 25\,{M}_{\odot }$) have evolved off of the main sequence by 5 Myr, and stars with masses $10\mbox{--}25\,{M}_{\odot }$ (B-type stars) dominate the spectrum from $6\mbox{--}20\,\mathrm{Myr}$. By 10 Myr there are not enough stars left with sufficiently high temperatures to produce significant amounts of ionizing radiation.

Zoom In Zoom Out Reset image size

Figure 1 shows the evolution in the intensity of ionizing photon production with age. We quantify this in Figure 2, where we show typical values of ${\hat{Q}}_{{\rm{H}}}$, the production rate of photons that are capable of ionizing hydrogen ($\lambda \geqslant 912\,{\rm{\mathring{\rm A} }}$) as a function of the age and metallicity of the stellar population. As expected, the youngest stellar populations produce the most ionizing photons and have the highest values of ${\hat{Q}}_{{\rm{H}}}$. As the population ages, cooler stars dominate the SED, producing less light at higher energies and decreasing the overall ionizing photon rate.

Zoom In Zoom Out Reset image size

The SSP metallicity has a second-order effect on the ionizing photon rate, attributed to (1) metallicity-dependent changes in the stellar atmospheres and (2) metallicity-dependent changes in stellar evolution. First, metals in stellar atmospheres absorb heavily, diminishing the UV flux for high-metallicity SSPs, lowering ${\hat{Q}}_{{\rm{H}}}$. Second, low-metallicity stellar populations have longer main sequence lifetimes, and can thus produce photons capable of ionizing hydrogen for longer periods.

Not only is the absolute number of ionizing photons important, but the exact distribution of ionizing photon energies is important as well. Photoionization ejects an electron with kinetic energy proportional to the energy of the ionizing photon, and thus the spectrum of initial electron velocities in the nebula reflects the spectrum of ionizing photons. The more high-energy photons are present, the "harder" a spectrum is, which ultimately affects the temperature and ionization structure of the H ii region. To quantify the hardness of the SSPs in our models, we calculate the slope of the extreme-ultraviolet (EUV) portion of the SED, as measured between the ionization threshold for helium (HeI, 24.6 eV or 505 ${\rm{\mathring{\rm A} }}$) and hydrogen (13.6 eV or 912 ${\rm{\mathring{\rm A} }}$). A large slope implies relatively few high-energy photons (a "soft" ionizing spectrum) and a smaller, flatter slope implies relatively more high-energy photons (a "hard" ionizing spectrum).

We show the time evolution of the EUV slopes in Figure 3. To the first order, the slope of the EUV spectrum is a function of SSP age. The youngest populations have spectra dominated by hot stars, which emit relatively more high-energy photons. Older populations have spectra dominated by cooler stars, which emit relatively fewer high-energy photons. Thus, as the population ages, the slope gradually steepens, or "softens." The abrupt change in spectral slope at $\mathrm{log}t\sim 6.5$ (∼3 Myr) roughly corresponds to the lifetime of an O-type star; once the O-type stars have evolved off of the main sequence, the slope suddenly becomes much "softer."

Zoom In Zoom Out Reset image size

The EUV slope is also a function of metallicity, since metallicity affects stellar evolutionary timescales and mass-loss rates. Low-metallicity populations are generally hotter, which hardens the EUV spectrum. Low-metallicity populations have weaker line-driven winds and experience less mass-loss, affecting main sequence lifetimes. The low-metallicity models produce a more gradual softening of the EUV spectrum, maintaining hard ionizing spectra for 1–2 Myr longer than the most metal-rich population.

Flux in the EUV portion of the spectrum is primarily produced by stars with $M\gt 10\,{M}_{\odot }$, and is thus closely tied to massive star evolution. As such, the reliability of the photoionization models is tied to the reliability of evolutionary models of massive stars. Although the ionizing spectrum tends to become both softer and less intense with time, there are noticeable departures from these trends at $\sim 5\,\mathrm{Myr}$. These fluctuations are due to Wolf–Rayet (WR) stars, which are massive stars (initial mass $\gtrsim 30\,{M}_{\odot };$ highlighted as $\sim 35\mbox{--}75\,{M}_{\odot }$ in Figure 1) that have evolved off of the main sequence. WR stars are stellar cores exposed from extreme mass loss, and their hot temperatures significantly harden the ionizing spectrum of the stellar population, seen clearly at $\mathrm{log}t\sim 6.7$ in Figure 3. While the effect of the WR phase appears huge in Figure 3, it occurs on such short timescales that it will likely not be detectable compared with the predictions of classical models of massive stars.

We note, however, that models of WR evolution are extremely uncertain, and that the exact masses and evolutionary pathways are not well-constrained. The extreme non-LTE nature of their atmospheres make WR stars challenging to model, and variations in mass loss and the strength of stellar winds make it difficult to predict ionizing fluxes (see review by Crowther 2007, and references within). This limitation likely imparts significant uncertainties into the evolution of the ionizing spectrum, particularly at late ($\sim 5$ Myr) times.

The Padova+Geneva isochrones used in this work do not include the effects of stellar rotation or binarity, both of which have important effects on the ionizing spectrum and main sequence lifetimes of massive stars (e.g., Eldridge & Stanway 2012; Levesque et al. 2012; Stanway et al. 2016). We discuss this issue in detail in Section 6, and in Section 6.1 we include an analysis of the ionizing spectra produced by the MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016), which include stellar rotation.

Ionizing spectra from single bursts may be a good approximation for a single massive H ii region, but real galactic stellar systems can show more complexity. Star clusters do not form instantaneously, and may be better modeled by a population with a range of ages spanning a few million years. Likewise, galaxies are subject to prolonged bursts of star formation involving many clusters, and their integrated spectra may be better represented with even more extended, complex SFHs.

We illustrate the effect of extended star formation by comparing the instantaneous burst models to models generated with a constant star formation rate (CSFR) of 1 ${M}_{\odot }$ per year. In this analysis, the constant SFR spectra are generated with FSPS and then fed as input to Cloudy. To calculate the total emission spectrum for complex or extended SFHs, FSPS adds the nebular emission from contributing SSPs from the nebular look-up tables generated by instantaneous bursts.

In the top panel of Figure 4 we compare the time evolution of ${\hat{Q}}_{{\rm{H}}}$ for instantaneous bursts and constant SFR populations. For instantaneous burst models, ${\hat{Q}}_{{\rm{H}}}$ decreases steadily with age as the massive star population evolves off the main sequence. In models with constant SFR, however, the rate of stars forming and the rate of stars dying eventually reach an equilibrium, after which there is little evolution in the ionizing spectrum. The constant SFR models reach an equilibrium around 4 Myr, after which the ionizing photon rate is essentially constant.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 4 shows the time evolution of the EUV slope. The slope of the constant-SF model reaches equilibrium around 6 Myr, with a slope that is roughly equivalent to that of a 2 Myr instantaneous burst. We note that these identified equilibrium ages correspond to a specific IMF and set of evolutionary tracks, and may be different for stellar populations generated with different properties.

Cloudy calculates the full radiative transfer through the gas cloud, so each individual H ii region model has internal structure, with radial variations in ionization state and temperature, which in turn affect the location within the nebula where various emission lines are produced. Before discussing the emission properties of the model H ii regions, we first explore their internal structure to better understand the physical conditions driving the global spectrum.

4.1.1. Model H ii Region Temperatures

The emergent emission spectrum is sensitive to the kinetic temperature of the free electrons, ${T}_{e}$, since collisions between electrons and metal ions are responsible for producing some of the most prominent emission lines observed in H ii region spectra. In this section, we describe how the equilibrium temperature of model H ii regions is affected by the input ionizing spectrum and gas phase metallicity. Our model links the metallicity-dependent changes in the ionizing spectrum with the coolants in the gas cloud, which drive the bulk trends in cloud temperature and ionization structure.

The equilibrium temperature of an H ii region is set by the balance of heating and cooling processes. Photoionization of hydrogen is the dominant source of heating in an H ii region. The net heating from photoionization is determined by (1) the photoionization rate and (2) the average energy of the liberated electrons, both of which depend on the intensity and shape of the ionizing radiation field by means of the number of incident ionizing photons and the energy injected per photoionization. The volumetric heating rate from photoionization, ${{\rm{\Gamma }}}_{\mathrm{ion}}$, is thus given by the photoionization rate weighted by the energy of the freed electron,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ion}}={n}_{{\rm{H}}}{\int }_{{\nu }_{0}}^{\infty }\displaystyle \frac{4\pi {J}_{\nu }}{h\nu }(h\nu -h{\nu }_{0}){\sigma }_{\nu }({\rm{H}})d\nu ,\end{eqnarray} \tag{ 3 }$$

where ${n}_{{\rm{H}}}$ is the neutral hydrogen density, ${\sigma }_{\nu }({\rm{H}})$ is the photoionization cross section of neutral hydrogen, ${\nu }_{0}$ the ionization energy, and ${J}_{\nu }$ is the mean specific intensity of the radiation field. This form is similar to integrating a monochromatic version of ${Q}_{{\rm{H}}}$ weighted by the energy of each photoelectron. Thus, for a given ${Q}_{{\rm{H}}}$, if there are relatively more high-energy photons (i.e., if the ionizing spectrum is harder), each photoionization will inject more kinetic energy, increasing the heating rate.

Cloud cooling is radiative, through a combination of line and continuum emission. We represent the total cooling rate, ${{\rm{\Lambda }}}_{\mathrm{rad}}$, as a sum of each of the major cooling processes:

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{rad}}={{\rm{\Lambda }}}_{\mathrm{ce}}+{{\rm{\Lambda }}}_{\mathrm{fb}}+{{\rm{\Lambda }}}_{\mathrm{ff}},\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Lambda }}}_{\mathrm{ce}}$ is the cooling from collisionally excited metal ions, ${{\rm{\Lambda }}}_{\mathrm{fb}}$ is the cooling from free-bound emission, and ${{\rm{\Lambda }}}_{\mathrm{ff}}$ is the cooling from free–free emission.⁸ In Equation (4) the cooling rates are ordered by their contribution to the total cooling rate; for H ii regions, the dominant cooling process is emission from collisionally excited metal ions.

In Figure 5 we show the fractional contribution of collisionally excited metal ions to the total cloud cooling as a function of model age and metallicity. At near-solar and super-solar metallicities, the cooling is dominated by forbidden and fine-structure transitions from collisionally excited metal ions, which can provide up to 95% of the total cooling. However, for metal-poor nebulae, the contribution from collisionally excited metal ions is less than 1% of the total cooling rate. In this regime, hydrogen provides the bulk of the cooling emission via free-bound line and continuum emission.

Zoom In Zoom Out Reset image size

The model H ii regions are in thermal equilibrium, so the net energy gained and lost by the nebula is equal at every age and metallicity. While the net heating and cooling rates in each cloud are always equal, the efficiency of the various cooling processes is highly variable, resulting in different model equilibrium temperatures. Metal line cooling is a particularly efficient coolant and produces much lower equilibrium temperatures than cooling from recombination emission.

In the top panel of Figure 6 we show the volume-averaged equilibrium temperatures of model H ii regions as a function of metallicity and ionization parameter for a 2 Myr model. ${T}_{e}$ varies from ∼4000–20,000 K, with the lowest cloud temperatures found in the most metal-rich models. Metal line cooling is the dominant cooling process for nebulae with metallicities $-1.0\lt \mathrm{logZ}/{{\rm{Z}}}_{\odot }\lt 0.1$, and the cooling efficiency is a strong function of the gas cloud metallicity. Scaling up the gas phase abundances increases the cooling efficiency, producing lower equilibrium temperatures, as expected from Figure 5.

Zoom In Zoom Out Reset image size

Below ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-1.0$, the shift in the dominant coolant produces a secondary dependence on the ionization parameter. In these metal-poor gas clouds, hydrogen is the only available coolant, with most of the cloud cooling emission produced in recombination lines and continuum emission. The strength of the recombination emission is strongly dependent on the number of incident ionizing photons; thus below ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-1.0$, ${T}_{e}$ depends on both metallicity and ionization parameter.

The bottom panel of Figure 6 shows the time evolution of the model H ii region temperatures at several different metallicities and ionization parameters. As the ionizing spectra soften with age, the equilibrium temperatures decrease. However, variations in metallicity and ionization parameter drive much larger variations in equilibrium temperature.

4.1.2. Ionization Structure and Line Emissivity of the Model H ii Regions

The ionization state of the cloud is critical for determining a number of processes within the nebula: the rate of radiative cooling, the rate at which the cloud absorbs photons from stars, and the chemical processes that can proceed within the cloud. The ionization structure sets the accessibility of critical emission lines; for example, [O iii] emission can only be produced if oxygen is doubly ionized. It is therefore necessary to understand what drives the ionization structure of H ii regions in order to understand the global emission line spectrum. Photoionization is by far the most important ionization process, and the frequency dependence of the ionization cross section means that the spectral shape will thus play an important role in determining the ionization structure within the model H ii region.

In Figure 7 we show the ionization structure and line emissivity of oxygen in a model H ii region at 1, 2, and 3 Myr at a fixed ionization parameter. At all ages, the high ionization species are found in the inner region of the nebula, while lower ionization species are found in the outer region of the nebula. Spectral shape regulates the spatial extent of high ionization species and the prevalence of partial-ionization zones, which in turn controls where emission from collisionally excited transitions is produced.

Zoom In Zoom Out Reset image size

While oxygen emission is strong at all ages, the age-dependent softening of the ionizing spectrum changes the location within the cloud where each emission line is produced. In the 1 Myr model, O++ is appreciably present in $\sim 60 \%$ of the cloud, with [O iii] emission produced at nearly all radii. In the 3 Myr model, however, O++ is only present in the innermost 25% of the cloud, and the emission from [O iii] is entirely contained within the doubly ionized oxygen zone. While the [O i], [O ii], and [O iii] lines are produced in overlapping physical regions in the 1 Myr model, the lines are produced in physically distinct regions of the cloud in the 3 Myr model.

Ideally we would like temperature and density diagnostics for the same ionic state, which would probe the same physical region within the nebula. [O iii] is a commonly used temperature diagnostic, but probes the temperature in the O++ zone, near the inner edge of the nebula. [O ii] is the only optically accessible oxygen density diagnostic, which measures the density in the the O+ zone, in the outer region of the nebula. The region probed in either case may not be representative of the nebula as a whole. This has implications for the interpretation of line ratios, as lines produced at different radii probe physically distinct regions of the nebula.

The strength of emission lines is not set solely by the total abundance of an ionic species in the nebula (i.e., the ionization structure shown in Figure 7), but is modulated by the energetic accessibility of energy levels as well. Collisional excitation involves an interaction between an ion and a free electron, and the rate of collisions will thus depend on the kinetic energy of the free electrons in the nebula.⁹

This is especially important in the context of our nebular model, which links the gas phase abundances to the metallicity of the ionizing population. In Section 4.1.1, we demonstrated that ${T}_{e}$ is a strong function of metallicity through the efficiency of metal line cooling. Both the hardness of the ionizing spectrum and the timescales associated with stellar evolution are metallicity-dependent, which has important effects on the thermodynamic properties of the model H ii regions.

In Figure 8 we show the oxygen ionization structure and line emissivity of a 3 Myr model at constant ionization parameter and different metallicities. In Figure 7 the age-induced softening of the ionizing spectrum changed the spatial extent of different ionization zones. Here we see a similar effect from the metallicity-induced changes in the ionizing spectrum, where the harder ionizing spectra in the low-metallicity model extends the O++ ionization zone.

Zoom In Zoom Out Reset image size

4.1.3. Hydrogen Recombination Lines

Hydrogen recombination lines are typically the strongest lines in H ii region spectra, and are widely used diagnostics of H ii region conditions when combined with metal lines. However, the hydrogen lines have a fundamentally distinct emission mechanism than the collisionally excited metal lines discussed in Section 4.1.2. Recombination lines are produced by the radiative recombination of a free electron with a proton into an excited state, followed by a radiative cascade to lower levels. They therefore have a different dependence on the properties of the nebula.

Analytic prescriptions for nebular emission generally assume that the number of ionizing photons and the number of recombination photons are equal, which yields a simple equation that converts ${Q}_{{\rm{H}}}$ to an Hα luminosity. It is thus unsurprising that the total power in the recombination lines closely tracks ${Q}_{{\rm{H}}}$, as measured from the ionizing spectrum.

However, the conversion between ${Q}_{{\rm{H}}}$ and ${L}_{{\rm{H}}\alpha }$ is not one to one, because the recombination coefficient, ${\alpha }_{{\rm{B}}}$, is a function of temperature. To get the correct luminosity of Hα and Hβ, it is essential to factor in both the number of ionizing photons and the temperature of the nebula. This is explicitly self-consistent in our nebular model, which adds the nebular emission to the same SSP (t, Z) as the one input to Cloudy, preserving the temperature sensitivity of the recombination lines as a function of age and metallicity. We note that the PopStar models do include a metallicity-dependent ${T}_{e}$ in their recombination coefficient, but not an age-dependent ${T}_{e}$.

In Figure 9 we show the Hα luminosity of the model H ii regions as a function of age, metallicity, and ionization parameter. We have normalized the Hα luminosities by ${Q}_{{\rm{H}}}$ to demonstrate the magnitude of the variations in Hα luminosity driven by changes in temperature. At young ages, metallicity changes produces 10%–20% variations in Hα line strength. The softening of the ionizing spectrum with age can change the recombination line strengths by 40%–50% Both of these will impact Hα-based SFR indicators that do not account for the temperature dependence of recombination lines.

Zoom In Zoom Out Reset image size

4.1.4. Nebular Continuum

The nebular emission model has two components: line emission and continuum emission from ionized gas. The implementation within FSPS allows the user to choose whether to include one or both in the output spectrum. For young, metal-poor populations, the continuum can contribute a significant portion of the total flux at optical and IR wavelengths, with a relative contribution that can be comparable with the stellar flux (e.g., Reines et al. 2010).

For H ii regions, the most important continuum emission processes are free-bound and free–free transitions: free-bound emission is responsible for most of the nebular continuum emission at optical and NIR wavelengths, while the free–free emission is more important at longer wavelengths.

The free-bound continuum is produced when a free electron recombines into an excited level of hydrogen, followed by the radiative cascade that produces recombination lines. The resultant continuum spectrum has a sharp "edge" at the ionization energy, followed by continuous emission to higher energies. As a hydrogen recombination process, the free-bound continuum is sensitive to model ionization parameter, which will set the overall intensity of the continuum emission. However, the general shape of the continuum is determined by the distribution of electron velocities and the recombination cross section, and will thus be sensitive to the temperature of the H ii region as well.

The free–free continuum, which is the result of a free electron scattering off of an ion or proton, produces a roughly power-law ($\propto {\nu }^{-2}$) distribution of photon energies and is also sensitive to the temperature of the H ii region. The two-photon continuum, while smaller in magnitude, can be important in the UV. The two-photon continuum is the result of a bound-bound process, where the excited 2s state of hydrogen decays to the 1s state by the simultaneous emission of two photons. The energy of the two photons totals to $h{\nu }_{\mathrm{Ly}\alpha }$, producing a bump in the UV portion of the spectrum.

In the top panel of Figure 10 we show the nebular continuum spectrum for a 1 and 3 Myr solar metallicity model. The overall intensity of the nebular continuum spectrum is set by the model ionization parameter, since recombination emission depends on the number of incident ionizing photons. At fixed metallicity and age, the continuum spectrum is nearly identical modulo a scaling factor, ${\hat{Q}}_{{\rm{H}}}$/${Q}_{{\rm{H}}}$, which scales the intensity of the model, as described in Section 2.1.4. Several spectral features are easily discernible: bound-free transitions produce the characteristic sawtooth edges across the spectrum; the two-photon continuum is responsible for the bump at $\sim 1500\,{\rm{\mathring{\rm A} }}$.

Zoom In Zoom Out Reset image size

While ${{ \mathcal U }}_{0}$ sets the continuum normalization, the relative height and steepness of the recombination edges are sensitive to the nebular temperature and thus the metallicity of the model. The temperature sensitivity of the recombination continuum is well-known, and the strength of the edge features have long been used as temperature diagnostics (e.g., Peimbert 1967).

In the bottom panel of Figure 10 we show the nebular continuum for models with different metallicities, again normalized by ${\hat{Q}}_{{\rm{H}}}$/${Q}_{{\rm{H}}}$. Low metallicity models have relatively more continuum emission. The emission edges are much broader, and there is less difference in amplitude among recombination edges. This behavior is primarily due to temperature changes driven by the changing cooling efficiencies.

Nebular continuum emission is clearly strongest at high ionization parameters and low metallicity. The harder ionizing spectra associated with young stellar populations further enhances this. In Figure 11 we show the total SED for 1 and 5 Myr SSPs at high and low metallicities and ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}=-1.0$. The 1 Myr model has much stronger line and continuum emission, which is strongest in the low-metallicity model.

Zoom In Zoom Out Reset image size

Reines et al. (2010) found that the nebular continuum contributed significantly ($\sim 40 \%$) to the NUV-NIR broadband flux of young ($\lt 3$ Myr) clusters, especially near the Balmer break. In Figure 12 we show the fractional contribution to the total flux as a function of model age for several wavelength ranges from the EUV to the NIR at fixed ionization parameter (${\mathrm{log}}_{10}{{ \mathcal U }}_{0}=-1.0$), and at two metallicities (solar and ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-1$). Our emission models agree with Reines et al. (2010), with continuum and line emission contributing at least 30% of the total flux at young ages in every wavelength range considered.

Zoom In Zoom Out Reset image size

The largest contribution from line and continuum emission relative to the stellar emission is seen in the optical and NIR. Nebular line and continuum emission contributes 95% of the total flux from $1.0\mbox{--}5.3$ μm for several million years, with continuum emission producing 70% of the total flux.

We show the fractional contribution to the broadband flux for a ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-1$ model in the bottom panel of Figure 12. The broad behavior in the low metallicity model is nearly identical to the solar metallicity model. We might expect line emission to contribute a larger fraction of the total flux in the low-metallicity model, since oxygen emission is generally stronger at sub-solar metallicities; this effect is only a few percent. The most noticeable difference between the high and low metallicity models is the time dependence of the nebular emission. In the low metallicity model, nebular emission contributes a significant portion of the total flux for longer, by a few Myr. Lower metallicity stellar populations have extended main sequence lifetimes, which in turn extends the timescale for nebular emission. This is an important feature of our nebular model, which links the stellar and nebular abundances, and ultimately allows for a more accurate accounting of the total flux from galaxies.

Figure 13 shows the same fractional contributions, but for models assuming a constant SFR. As the young stars evolve off of the main sequence and a substantial population of older stars is built up, the stellar contribution increases. Nebular emission remains a significant fraction of the total light at 10 Myr and beyond, and is especially important at low metallicity. Thus, for galaxies with a burst of recent star formation or star formation extended over several million years, it is essential to include the contribution from nebular line and continuum emission to accurately predict the underlying properties of the stellar population from broadband photometry.

Zoom In Zoom Out Reset image size

With an understanding of the physical properties of the model H ii regions, we can better assess the processes driving variations in nebular emission. In this section, we characterize the emission line properties of the model H ii regions and their dependence on the model parameters.

4.2.1. Broad Trends in Emission Line Strength

Figure 14 shows global trends in line strength for the strongest optical emission lines: Hβλ4861, [O iii]λ5007, Hαλ6563, [N ii]λ6584, and [S ii]λ6713. To showcase the results of the implementation within FSPS, the emission lines plotted in Figure 14 were generated directly from python-fsps by generating SSPs with and without emission lines and subtracting the stellar-only spectrum.¹⁰

Zoom In Zoom Out Reset image size

As discussed in Section 2.2, when FSPS adds in nebular emission, the emission lines are scaled by the number of ionizing photons in the incident SSP. To first order, the model age sets the overall intensity of the emission lines through this ${Q}_{{\rm{H}}}$ normalization. The 1 Myr SSP is much brighter and produces many more ionizing photons than the 4 Myr SSP and has stronger emission lines. Despite the normalization by ${Q}_{{\rm{H}}}$, emission lines sensitive to the overall ionization parameter still vary with ${{ \mathcal U }}_{0}$. For example, [S ii] is a well-known ionization-sensitive line, and is strongest at low ${{ \mathcal U }}_{0}$.

For models at constant age and ionization parameter, changing the model metallicity has two effects: it changes (1) the shape of the ionizing spectrum and (2) the gas phase abundances. The first of these effects is reflected in the metallicity-dependent variability in Hα and Hβ line strengths. The metal-poor models produce hotter H ii regions, which in turn produces stronger recombination line emission. The second of these effects is demonstrated in the [N ii] and [S ii] line strengths, which are strongest in the highest metallicity models, which reflects the decreasing metal abundance.

The [O iii] emission is less straightforward to describe due to the well-known double-valued and nonlinear relationship between oxygen line strength and metallicity (Pilyugin & Thuan 2005; Kewley & Ellison 2008). For the selection of metallicities plotted in Figure 14, [O iii] emission is strongest at ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-0.3$. At the lowest metallicities, oxygen is not abundant enough to produce significant emission; at the highest metallicities, the temperature is too low to collisionally excite a substantial population of O++ ions to make the [O iii]λ5007 transition.

4.2.2. Emission Line Ratios

In 4.2.1, we established broad trends in metal and recombination line strengths and their dependence on model parameters. In this section, we will build upon that intuition to understand the behavior of various line ratios used to probe the physical state of the nebula.

Most line ratios involve the comparison of a metal line to a Balmer line. The total power emitted in collisionally excited metal lines is proportional to the total cooling in metal lines, and thus sensitive to metallicity and the ionization parameter (Section 4.1.1). As discussed in Section 4.1.3, the total power emitted in the hydrogen recombination lines is primarily driven by ${Q}_{{\rm{H}}}$ and should thus trace the total luminosity of the ionizing source.

In practice, however, we never have access to reliable measurements of the total power in metal line emission and instead measure the fluxes of a small number of individual lines from a particular species. For example, in the left panel of Figure 15, we show the the ratio of [O iii]λ5007/Hβ as a function of metallicity and ionization parameter. The [O iii] line strengths do not scale directly from the total metal line emission, and the resultant ratio between [O iii] and Hβ does not monotonically scale with metallicity. Instead, the line ratio is double-valued with metallicity.

Zoom In Zoom Out Reset image size

Some of this behavior is due to changes in the ionization state of oxygen from O+ to O++. At high metallicity, temperatures are too low to collisionally excite an appreciable population of O++, decreasing the [O iii] emission. The middle panel of Figure 15 compares the line strengths of oxygen from two different ionization states, [O iii]λ5007 and [O ii]λ3727. While the [O iii]λ5007 emission decreases, the [O ii] emission stays relatively strong. The ratio of just these two lines proves to be an excellent probe of the overall excitation of the H ii region.

In the right panel of Figure 15, we show the [N ii]λ6584/Hα line ratio as a function of model metallicity and ionization parameter. The [N ii] line strength is not double-valued with metallicity at the temperatures found in our model H ii regions.

We test our predicted emission line luminosities against those measured from spectroscopic observations of both H ii regions and star-forming galaxies, to show that the FSPS nebular model can reproduce observed emission from both single H ii regions and complex populations of multiple H ii regions. Even though we do not know the ${{ \mathcal U }}_{0}$, Z, and t values of the observed objects independently, theoretical grids should at least be able to span the range of observed line ratios using reasonable physical parameters. The H ii region comparison spectra from van Zee et al. (1998) consist of observations of massive H ii regions in nearby galaxies, and are a standard comparison set for nebular emission models (Dopita et al. 2000; Kewley et al. 2006; Levesque et al. 2010; Dopita et al. 2013). For star-forming galaxies, SDSS galaxy spectra are another common comparison set used in optical line ratio diagnostic diagrams (Dopita et al. 2000; Kewley et al. 2006; Levesque et al. 2010). The SDSS spectra typically cover large areas of a galaxy and are likely to contain emission from multiple H ii regions.

The van Zee et al. (1998) catalog reports emission line intensities with contributions from both doublet lines for [N ii] (6548, 6584) and [O iii] (4959, 5007). Standard diagnostic diagrams only use the stronger of the two, [N ii]λ6584 and [O iii]λ5007; we remove the contribution from the weaker line using the theoretical intensity ratio between the two lines (e.g., ${I}_{5007}=2.88\ {I}_{4959}$). This results in a decrease in the [N ii]λ6584/Hα and [O iii]λ5007/Hβ ratio of $\sim {\mathrm{log}}_{10}(3/4)$ or ∼0.1 dex.

The star-forming galaxy sample is derived from galaxy spectra from the Sloan Digital Sky Survey Data Release 7 (SDSS DR7; York et al. 2000; Abazajian et al. 2009) and emission line fluxes measured from the publicly available SDSS DR7 MPA/JHU catalog4 (Kauffmann et al. 2003a; Brinchmann et al. 2004; Salim et al. 2007). We use the emission line sample presented in Telford et al. (2016), briefly summarized here. The sample includes ∼135,000 galaxies with redshifts between 0.07 and 0.30. Galaxies are required to have S/N of 25 in the Hα line, 5 in the Hβ line, and 3 in the [S ii] lines. Emission line fluxes are corrected for dust extinction using the Balmer decrement and the Cardelli et al. (1989) extinction law, assuming ${R}_{{\rm{V}}}=3.1$ and an intrinsic Balmer decrement of 2.86. AGNs are removed from the sample according to the empirical BPT diagram classification of Kauffmann et al. (2003b).

In Figures 16 and 17, we show how each of the model parameters affects the location of a model H ii region in the BPT diagram. In Figure 16 we vary ${{ \mathcal U }}_{0}$ and ${\mathrm{log}}_{10}Z/{Z}_{\odot }$ for a 1 and 2 Myr model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

${{ \mathcal U }}_{0}$ variations: Increasing model ${{ \mathcal U }}_{0}$ moves the model along the star-forming sequence to higher values of [O iii]λ5007/Hβ and lower values of [N ii]λ6584/Hα. Nitrogen and oxygen have similar ionization potentials, so comparing a doubly ionized population, [O iii], with a singly ionized population, [N ii], probes the ionization state of the gas cloud. Thus increasing ${{ \mathcal U }}_{0}$ will enhance the doubly ionized population at the expense of the singly ionized population.

Abundance variations: Decreasing the gas phase abundances moves the model away from the star-forming sequence, to lower values of [N ii]λ6584/Hα and [O iii]λ5007/Hβ. However, decreasing the gas phase metallicity also results in a higher equilibrium temperature. Initially, the changing line ratios reflect the change in ionization state, with enhanced [O iii] emission and decreased [N ii] emission. Decreasing the gas phase abundances below ${\mathrm{log}}_{10}Z/{Z}_{\odot }=-0.6$ results in a decrease in both [O iii]λ5007/Hβ and [N ii]λ6584/Hα. For high-metallicity models, the temperature is too low to collisionally excite O++ ions to make the [O iii]λ5007 transition, and so [O iii]λ5007/Hβ decreases at roughly constant [N ii]λ6584/Hα.

Age variations: We show the age evolution of the ${\mathrm{log}}_{10}Z/{Z}_{\odot }$–${\mathrm{log}}_{10}{{ \mathcal U }}_{0}$ grid for instantaneous bursts in the top panel of Figure 17. The grid is well-matched to the star-forming sequence for the 1 Myr models. By 3 Myr, however, very few models can produce line ratios consistent with the star-forming locus or H ii regions. Those that do have very high ionization parameters (${\mathrm{log}}_{10}{{ \mathcal U }}_{0}\gt 2$) and a very narrow range of metallicities (${\mathrm{log}}_{10}Z/{Z}_{\odot }\sim -0.6$). The WR phase at ∼5 Myr significantly hardens the ionizing spectrum.

The bottom panel of Figure 17 shows the time evolution of the ${\mathrm{log}}_{10}Z/{Z}_{\odot }$–${\mathrm{log}}_{10}{{ \mathcal U }}_{0}$ grid for constant SFR models. At the youngest ages, the ionizing spectra produced by the instantaneous burst models and the constant star formation models are very similar, and produce very similar line ratios. As discussed in Section 3.5, the CSFH models reach an equilibrium around ∼4 Myr. Thus the CSFH models continue to produce line ratios that match the observed star-forming locus to later ages.

We compare the constant SFR model grid with the model grid presented in Dopita et al. (2013; hereafter referred to as D13), which uses Starburst-99 ionizing spectra and the Mappings-III photoionization code. We do not run the photoionization models ourselves; rather we use the line ratios from the D13 paper directly. To make a cleaner comparison, we run our constant SFR model through Cloudy at the same ionization parameters and metallicities are used in the D13 model, and match the gas phase abundances to those used in D13 (see Table 1). We attempt to match the density of the D13 models at ${n}_{{\rm{H}}}=300$. However, as discussed previously, Mappings-III is a pressure-based photoionization code, and ${n}_{{\rm{H}}}=300$ is just the initial average density in the Dopita et al. (2013) model. The D13 paper tested several different electron velocity distributions; we compare our model to their model with $\kappa =\infty$, which corresponds to the standard Maxwell–Boltzmann distribution assumed in Cloudy.

We show the BPT diagram for the two constant SFR model grids in Figure 18. Despite the different approaches, both models show substantial overlap and are able to reproduce most of the observed SF-sequence. While the two grids show similar trends with ionization parameter and metallicity, the overall model coverage is different. The FSPS model shows a wider spread in the [N ii]λ6584/Hα and [O iii]λ5007/Hβ line ratios, while the D13 grid extends to more extreme values in [O iii]λ5007/Hβ. This difference is likely due to the fact that the D13 grids extend to higher metallicities than the FSPS grid, up to ${\mathrm{log}}_{10}Z/{Z}_{\odot }=0.7$, or $5{{\rm{Z}}}_{\odot }$. Our model ties the gas phase abundances with the metallicity of the stellar population, and the maximum metallicity in the Pavoda+Geneva isochrones is ${\mathrm{log}}_{10}Z/{Z}_{\odot }=0.2$.

Zoom In Zoom Out Reset image size

While the BPT diagram is the most widely used diagnostic diagram, there are many combinations of line ratios used to probe physical properties of H ii regions and star-forming galaxies. Theoretical grids of photoionization models are able to reproduce observed BPT diagram line ratios quite easily, but these same models can fail to reproduce observed line ratios in other diagnostic diagrams (e.g., Telford et al. 2016). In this section we assess several diagnostic diagrams and showcase our model's ability to reproduce observed line ratios.

N2O2: The [N ii] λ6548/[O ii] λ3727 (N2O2) line ratio correlates with metallicity and is widely used as an abundance indicator (Veilleux & Osterbrock 1987; Dopita et al. 2000; Levesque et al. 2010). The line ratio is especially useful as a diagnostic when paired with an ionization-sensitive line ratio like log[O iii] λ5007/[O ii] λ3727 (O3O2), discussed in Section 4.2.2. In Figure 19 we show the N2O2 and O3O2 line ratios for a 1 Myr instantaneous burst and a 4 Myr constant SFR model. The model grids show good overall agreement with the data, and can reproduce most of the observed H ii region line ratios. The more extreme N2O2 line ratio values may be produced by gas phase metallicities not reached by our model grid. The turnover in ionization parameter at the highest metallicity in our model, ${\mathrm{log}}_{10}Z/{Z}_{\odot }=0.2$, suggests that the ionizing spectra at high metallicity are not hard enough.

Zoom In Zoom Out Reset image size

Ne3O2: Levesque & Richardson (2014) presented the utility of the [Ne iii] λ3869/[O ii] λ3727 (Ne3O2) line ratio as an ionization parameter diagnostic. The O3O2 line ratio is a well-known probe of excitation, but the wavelength separation of the lines makes the diagnostic quite sensitive to reddening and requires observations that cover a wide wavelength range. The Ne3O2 remedies both of these issues, and shows a tight correlation with O3O2. However, Levesque & Richardson (2014) found considerable offset ($\sim 0.6$ dex) between the models and the observations of Ne3O2 and O3O2 line ratios, suggesting an underlying deficiency in the predicted emission line fluxes, which they attributed to an insufficiently hard ionizing spectrum.

In Figure 20, we show the Ne3O2 and O3O2 line ratios produced by our ${{ \mathcal U }}_{0}$ − Z model grid for a 1 Myr instantaneous burst and a 4 Myr constant SFR model compared with the massive extragalactic H ii regions from van Zee et al. (1998). Our grids show considerable improvement from the models used in Levesque & Richardson (2014), reducing the offset between model and data to ∼0.2 dex. Both the instantaneous burst and constant SFR models are able to reproduce 75% of the observed H ii region line ratios.

Zoom In Zoom Out Reset image size

For completeness, in Figures 21 and 22 we show model grids for two additional common diagnostic diagrams: R23 versus O3O2 and [S ii]/Hα versus [O iii]/Hβ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recent work has suggested that the redshift evolution in the location of the star-forming sequence in the standard BPT diagram can be attributed to harder ionizing spectra (e.g., Steidel et al. 2014). The intensity and hardness of the EUV portion of the spectrum dictates the overall ionization structure and temperature of the nebula, which in turn affects the overall location of the model grid in various diagnostic diagrams presented in Section 5.3. SPS models have a number of knobs to turn that can alter the amount of EUV flux significantly, through stellar atmospheres (winds, opacities) and stellar evolution (mass loss, rotation, binarity).

We do not compare different stellar atmospheric models, since this is not currently an easily flexible aspect of FSPS. We do, however, consider different treatments of stellar evolution, by varying the evolutionary tracks used to produce the ionizing spectra input to Cloudy. We first analyze bulk differences in the ionizing spectra generated with different isochrones, and then discuss the effect on the resultant photoionization models. Within FSPS, there are four isochrone sets to choose from: Padova (Bertelli et al. 1994; Girardi et al. 2000; Marigo et al. 2008) + Geneva (Schaller et al. 1992; Meynet & Maeder 2000); BaSTI (Pietrinferni et al. 2004; Cordier et al. 2007); PARSEC (Bressan et al. 2012); and MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016). Here we focus on the comparison between the Padova+Geneva isochrones used in Section 4, which reflect the "industry standard" in this context, and the MIST isochrones (Choi et al. 2016; Dotter 2016), the newest stellar evolution calculations included in FSPS.

The Padova+Geneva model uses 2007 Padova isochrones with high mass-loss rate Geneva isochrones adopted for $M\gt 70\,{M}_{\odot }$. For the comparisons made in this work, both the Padova+Geneva models and the MIST models use a Kroupa IMF and identical stellar atmospheres. O-star spectra are from WMBasic (Pauldrach et al. 2001; updates from J. J. Eldridge et al. 2017, in preparation); Wolf–Rayet stellar spectra are from CMFGEN (Hillier & Lanz 2001); and post-AGB stellar isochrones are from Vassiliadis & Wood (1994), with spectra from Rauch (2003).

The MIST models include stellar rotation, which affects stellar lifetimes, luminosities, and effective temperatures through rotational mixing and mass loss (see Choi et al. 2016 for further details). Rotating stars tend to have harder ionizing spectra and higher luminosities. In Figure 23 we show the time and metallicity evolution of ${\hat{Q}}_{{\rm{H}}}$ and the EUV slope for the MIST evolutionary tracks, identical to the ones shown for the Padova+Geneva models in Figure 3. The left panel of Figure 23 shows the evolution of ${\hat{Q}}_{{\rm{H}}}$, the ionizing photon rate, which shows the same qualitative behavior as the Padova+Geneva models: ${\hat{Q}}_{{\rm{H}}}$ gradually decreases with time and the lower metallicity models produce higher values of ${\hat{Q}}_{{\rm{H}}}$. Quantitatively, however, the MIST values of ${\hat{Q}}_{{\rm{H}}}$ are larger by a factor of 2–3 and maintain larger values for longer. For MIST, the lowest metallicity spectra produce elevated values of ${\hat{Q}}_{{\rm{H}}}$ until nearly 10 Myr.

Zoom In Zoom Out Reset image size

The right panel of Figure 23 shows the time evolution of the EUV slope. At the earliest ages, the MIST and Padova+Geneva ionizing spectra have similar slopes, but the evolution with time is significantly different. The Padova+Geneva spectra begin to soften around 3 Myr, while the MIST models stay relatively hard until 5 Myr. This is likely due to the extended main sequence lifetimes afforded by rotational mixing.

The prolonged high values of ${Q}_{{\rm{H}}}$ seen in the MIST models imply that the MIST SSPs could sustain ionizing radiation for a longer period of time compared with the Padova+Geneva models. We show the time evolution of the BPT diagram in Figure 24. While the two models are qualitatively similar at the youngest ages, they produce very different line ratios at later ages. The Padova+Geneva grid begins to fall away from the star-forming locus on the BPT diagram at 2 Myr; at 3 Myr the Padova+Geneva models fail to reproduce H ii region line ratios, consistent with observed star-forming regions. Conversely, the MIST models can match the star-forming locus until at least 4 Myr, due to the increased number of ionizing photons and harder ionizing spectra.

Zoom In Zoom Out Reset image size

The effect of WR stars at 5 Myr on the Padova+Geneva models in Figure 24 is striking. WR stars produce an extremely hard ionizing spectrum, which rejuvenates an ionizing spectrum that would otherwise be too soft to produce emission line ratios consistent with observed star-forming galaxies. However, the short timescales associated with the WR phase imply that this effect would be minimal on galaxy scales.

6.1.1. Alternative Line Ratios

He2: [He ii] emission is produced in H ii regions but can also be produced in the atmospheres of WR stars. Steidel et al. (2016) discussed the discrepancy between predicted [He ii] λ1640/Hβ (He2) line ratios and those observed for a sample of massive star-forming galaxies near $z\sim 2$. Predictions for the nebular [He ii] emission from Starburst-99 models were too weak to match observed [He ii] line ratios. Only the BPASS models (Eldridge 2012), which include stellar and nebular [He ii] emission, could match the observed [He ii] flux. Thus Steidel et al. (2016) deduced that photospheric emission must also contribute to the [He ii] flux.

In Figure 25 we show the He2 line ratios predicted by our model grids. The Padova+Geneva models are shown in gray, and are unable to produce enough [He ii] emission to explain the observed He2 line ratios from Steidel et al. (2016). However, the MIST models produce significant nebular [He ii] emission from 3–5 Myr, which can fully account for the observed He2 line ratio without the need to also include stellar emission. The harder ionizing spectra produced by the MIST models show clear improvement from the Padova+Geneva models for emission lines associated with high excitation species.

Zoom In Zoom Out Reset image size

For completeness, in Figure 26 we show MIST model grids for several additional common diagnostic diagrams, including R23 and [S ii]/Hα.

Zoom In Zoom Out Reset image size

The region of the BPT diagram between the star-forming sequence and the AGN sequence is occupied by objects classified as "low ionization emission regions" (LIERs; Belfiore et al. 2016).¹¹ LIER-like emission is characterized by strong low-ionization forbidden lines (e.g., [N ii], [S ii], [O ii], [O i]) relative to recombination lines, and was originally discovered in the nuclear regions of galaxies and attributed to AGN-related activity (Kauffmann et al. 2003b; Kewley et al. 2006; Ho 2008). While some cases are certainly still driven by the presence of a weak AGN, the discovery of spatially extended (∼kpc scales) LIER-like emission has led to work suggesting that hot, evolved stars could be responsible for the ionizing radiation in other cases (Singh et al. 2013; Belfiore et al. 2016). The leading candidate for the ionizing source is the population of post-AGB stars (Binette et al. 1994; Sarzi et al. 2010; Yan & Blanton 2012). Extreme horizontal branch stars have also been suggested candidates, but considering their effect is beyond the scope of this paper.

Post-AGB stars are stars that have left the AGB, evolving horizontally across the HR diagram toward very hot temperatures (∼10⁵ K) before cooling down to form white dwarfs, with a fraction of the post-AGB population forming planetary nebulae. The exposed cores of post-AGB stars are hot enough to ionize hydrogen, and thus could produce a radiation field capable of ionizing the surrounding ISM, provided that there are enough post-AGB stars.

Post-AGB stars are not capable of matching the ionizing flux produced by a single O-star; their strength lies in numbers. The progenitors of these stars have initial masses from 1 to 5 ${M}_{\odot }$, and are much more common than O-type stars. For early-type galaxies, where the bulk of the stellar population is old, AGB and post-AGB stars can contribute a significant portion of the total galaxy luminosity. Yan & Blanton (2012) measured the ionization parameter and gas density for galaxies with LIER-like emission and compared them with the typical numbers of ionizing photons produced by post-AGB stars. They deduced that the ionization parameter for post-AGB stars falls short of the required value by a factor of 10, implying that the gas geometry must be quite close to the stars themselves.

The geometry between the stars and gas is one of the major uncertainties associated with interpreting LIER-like emission. We first determine if populations including post-AGB stars are capable of producing line ratios consistent with LIER-like emission. The geometry previously assumed in this work (${\mathrm{log}}_{10}{R}_{\mathrm{inner}}=19;$ ${n}_{{\rm{H}}}=100$) is appropriate for massive H ii regions found in star-forming galaxies, where the gas is associated with the natal gas cloud but may not accurately describe the gas geometry in a scenario where old stars provide the ionizing radiation.

To test the sensitivity to geometry, we generate ionizing spectra for older SSPs ($\gt 1$ Gyr) that include post-AGB stars to use as input to Cloudy, and run photoionization models at different values of ${R}_{\mathrm{inner}}$ and ${n}_{{\rm{H}}}$. In Figure 27 we show the BPT diagram line ratios for the post-AGB star ionizing spectra at several different ionization parameters, with ${n}_{{\rm{H}}}$ varied from 10 to 1000 and the inner radius varied from 10¹⁸ to 10²⁰ cm (0.3–30 pc).

Zoom In Zoom Out Reset image size

The post-AGB models produce line ratios consistent with LIER-like emission, well outside of the "pure star-forming" region of the BPT diagram, as identified by Kauffmann et al. (2003a). The ionization parameters required to produce observable line ratios are ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}\sim -5$ to ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}\sim -3$. At ${R}_{\mathrm{inner}}={10}^{19}$ and ${n}_{{\rm{H}}}=100$ cm⁻³, this implies a total initial stellar mass of 10⁶–${10}^{8}\,{M}_{\odot }$. The required stellar mass would be higher for models where the inner radius is further from the ionizing source. We note that the ionizing spectrum was based on a stellar population at solar metallicity; low-metallicity post-AGB stars are hotter and would likely enhance the ionizing radiation.

The line ratios in Figure 27 show little sensitivity to the star-gas geometry, a result of our simplified model in which the gas exists in a plane-parallel shell surrounding a central point source of ionizing radiation. If the gas is produced by the AGB stars themselves or has a spatial distribution that differs from the distribution of stars, the geometry will differ substantially from the simplified shell used in this work. In future work we plan to test the effects of model geometry in more detail.

Another major uncertainty associated with the interpretation of LIER-like emission lies in our incomplete understanding of the late stages of stellar evolution. In addition to the poor constraints on the age distribution and lifetimes of post-AGB stars, it is uncertain what fraction of post-AGB stars contribute to the large-scale photoionizing radiation field. A recent census of the old, UV-bright stellar population in M31 was unable to reproduce predicted numbers of post-AGB stars (Brown et al. 2008).

To understand the sensitivity of LIER-like emission to the underlying stellar model, we generate ionizing spectra with varying contribution from post-AGB stars using the pagb parameter in FSPS. This parameter specifies the weight given to the post-AGB star phase, where pagb = 0.0 turns off post-AGB stars and pagb = 1.0 implies that the Vassiliadis & Wood (1994) tracks are implemented as is, the default behavior in FSPS. In the top panel of Figure 28 we show the ionizing spectrum produced by a 3 Gyr SSP with pagb set to 1, 10⁻¹, 10⁻², and 10⁻³. Scaling down the implemented post-AGB stars scales down the emergent EUV radiation in the spectrum.

Zoom In Zoom Out Reset image size

We show the resultant BPT diagram for the photoionization models that vary in both ${{ \mathcal U }}_{0}$ and pagb in the bottom panel of Figure 28. We find that the luminosity of the post-AGB stars could be reduced by a factor of two and still produce LIER-like emission, provided there are enough stars.

Belfiore et al. (2016) found that the [S ii]/Hα line ratio provided a clean separation between LIER-like emission and Seyfert-like emission. In Figure 29 we show the [S ii] emission produced by our post-AGB models. The line ratios produced by our post-AGB star models show the elevated [S ii]/Hα ratios observed in LIER galaxies, and confirm the utility of [S ii] as a means of identifying low-ionization emission.

Zoom In Zoom Out Reset image size

FSPS includes two separate nebular emission models: one that is dust-free and one that includes the effects of dust within the gas cloud itself. We note that neither the dusty nor the dust-free models include the effect of reddening or extinction. The effect of including grains is primarily on the thermodynamic properties of the cloud, and extinction and reddening should be applied to the emergent spectrum.

Users should pair the dusty nebular grids with the dust model used in Conroy et al. (2009) from Charlot & Fall (2000), already included within FSPS. This dust model has two separate dust attenuation prescriptions: a component that affects the entire stellar population and a component that applies additional extinction to the young stellar populations to mimic the effect of young clusters embedded in a birth cloud. The wavelength dependent extinction, ${\hat{\tau }}_{\lambda }$, is given by

$$\begin{eqnarray*}{\hat{\tau }}_{\lambda }(t)\equiv \left\{\begin{array}{ll}{\hat{\tau }}_{1}{(\lambda /5500)}^{-0.7} & t\leqslant {t}_{\mathrm{esc}}\\ {\hat{\tau }}_{2}{(\lambda /5500)}^{-0.7} & t\gt {t}_{\mathrm{esc}},\end{array}\right.\end{eqnarray*}$$

with ${\hat{\tau }}_{1}=1.0$ for the galaxy-wide attenuation and ${\hat{\tau }}_{2}=0.3$ for the young populations. ${t}_{\mathrm{esc}}={10}^{7}$ years, roughly matched to the timescale of birth cloud disruption (Conroy et al. 2009).

The combination of the dusty nebular model with the Charlot & Fall (2000) dust model assumes a very simple gas-dust geometry for the H ii regions, where the ionizing source is embedded within a gas cloud, with a layer of "natal" dust exterior to the ionized region. We do not include the effects of dust attenuation or extinction on the ionizing spectrum input to Cloudy, and the nebula experiences the full, unattenuated ionizing photon flux. The emergent spectrum is subsequently attenuated and reddened before leaving the galaxy. This is an unrealistic take on the dust geometry, since the dust and gas is likely to be clumpy and unevenly distributed around the stellar population.

The main difference between the dusty and dust-free models is in their thermodynamic properties. Dust changes the temperature structure of the cloud, since it provides alternative mechanisms for heating and cooling the gas. If cloud temperatures are low enough to allow the survival of dust grains, grain photoionization and thermionic emission can be an efficient heating mechanism ($\lesssim 10 \%$ of the total heating).

In Figure 30 we show the change in temperature produced by model H ii regions that include dust grains when compared with the grain-free model as a function of metallicity and ionization parameter for instantaneous burst models at 1 Myr. The models that include grains produce slightly hotter H ii regions; at solar metallicity and ${\mathrm{log}}_{10}{{ \mathcal U }}_{0}=-1$, the grain model produces a volume-averaged electron temperature that is 25% hotter than the grain-free model. This effect is only important at high metallicity, where the nebula is cool enough such that dust grains survive and contribute additional heating. In the high temperatures associated with the low-metallicity models, dust grains are destroyed and the inclusion of grains has little effect on the thermal properties of the H ii region.

Zoom In Zoom Out Reset image size

The inclusion of dust grains is most important for highly ionized gaseous regions at high metallicity, where they have a non-negligible effect on emission line strengths. In Figure 31 we show the BPT diagram for the model H ii regions with and without dust grains. Both grids do a reasonable job of covering the range of observed emission line strengths. However, for models at high metallicity and ionization parameter, the higher equilibrium temperatures produced in the dusty models produce enhanced [O iii]λ5007/Hβ at the same [N ii]λ6584/Hα. The line ratios differ by less than 10% at their most extreme values, and still fail to reproduce the star-forming locus beyond ∼3 Myr. A more careful analysis is required to determine which model provides the best fit to the data; we leave this for future studies.

Zoom In Zoom Out Reset image size