Lisa J Kewley

Michael A Dopita

Claus Leitherer Stsci

Romeel Davé

Tiantian Yuan

Mark Allen

Brent Groves

Ralph Sutherland

THEORETICAL EVOLUTION OF OPTICAL STRONG LINES ACROSS COSMIC TIME

Keywords: galaxies, starburst-galaxies, abundances-galaxies, fundamental parameters

We use the chemical evolution predictions of cosmological hydrodynamic simulations with our latest theoretical stellar population synthesis, photoionization and shock models to predict the strong line evolution of ensembles of galaxies from z = 3 to the present day. In this paper, we focus on the brightest optical emission-line ratios, [N II]/Hα and [O III]/Hβ. We use the optical diagnostic Baldwin-Phillips-Terlevich (BPT) diagram as a tool for investigating the spectral properties of ensembles of active galaxies. We use four redshift windows chosen to exploit new near-infrared multi-object spectrographs. We predict how the BPT diagram will appear in these four redshift windows given different sets of assumptions. We show that the position of star-forming galaxies on the BPT diagram traces the ISM conditions and radiation field in galaxies at a given redshift. Galaxies containing AGN form a mixing sequence with purely star-forming galaxies. This mixing sequence may change dramatically with cosmic time, due to the metallicity sensitivity of the optical emission-lines. Furthermore, the position of the mixing sequence may probe metallicity gradients in galaxies as a function of redshift, depending on the size of the AGN narrow line region. We apply our latest slow shock models for gas shocked by galactic-scale winds. We show that at high redshift, galactic wind shocks are clearly separated from AGN in line ratio space. Instead, shocks from galactic winds mimic high metallicity starburst galaxies. We discuss our models in the context of future large near-infrared spectroscopic surveys.

INTRODUCTION

Understanding how galaxies formed and evolved is one of the primary drivers of modern astronomical research. Measuring the fundamental properties of galaxies as a function of cosmic time is now possible by combining kewley@mso.anu.edu.au spectroscopy from the world's largest telescopes with state-of-the-art theoretical simulations.

The collisionally excited emission-lines provide key diagnostics of the gas-phase chemical abundance, the ionization state of the gas, the dust extinction, and the ionizing power source of the galaxy. [START_REF] Asplund | Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis[END_REF] showed that a diagram of [N II]/Hα and [O III]/Hβ (com-monly referred to as the BPT diagram) can be used to classify galaxies dominated by AGN from those dominated by star formation. Star-forming galaxies form a tight abundance sequence on the BPT diagram (Dopita & Evans 1986;Dopita et al. 2000). The location of the star-forming abundance sequence on the [N II]/Hα versus [O III]/Hβ diagram at a given redshift probes (1) the spread in global metallicity of the observed starforming galaxy population, (2) the stellar ionizing radiation field of the star-forming population, and (3) the conditions of the interstellar medium (ISM) surrounding the star-forming regions. Therefore, the [N II]/Hα versus [O III]/Hβ diagram can be used as a tool for investigating the metallicity and ionizing radiation field in starforming galaxies as a function of cosmic time, independent of the large systematic errors that plague the chemical abundance scale (Kewley & Ellison 2008;Bresolin et al. 2009;Kudritzki et al. 2012).

Galaxies containing an Active Galactic Nucleus (AGN) are easily separated from purely star-forming galaxies using the BPT diagram. The extreme ultra-violet (EUV) "hard" radiation field from the accretion disk of an AGN ionizes the [O III] and [N II] lines, producing larger [O III]/Hβ and [N II]/Hα line ratios than usually seen in star-forming galaxies. [START_REF] Osterbrock | [END_REF] and Veilleux & Osterbrock (1987, hereafter VO87) derived the first semi-empirical AGN classification schemes for the BPT diagram based on a combination of observations and and photoionization models. These diagnostics were refined by Kewley et al. (2001b) using stellar population synthesis, photoionization and shock models. Power-law AGN models show that position and spread of the AGN region on the BPT diagram traces the metallicity in the extended narrow-line region of AGN (Groves et al. 2004b).

Clearly, the gas-phase metallicity has a critical influence on the location of both star-forming galaxies and AGN on the BPT diagram. Both theory and observations indicate that the mean metallicity of star-forming galaxies rises with cosmic time (e.g., Nagamine et al. 2001;De Lucia et al. 2004;Kobulnicky & Kewley 2004;Kobayashi et al. 2007;Maiolino et al. 2008;Davé et al. 2011a;Yuan et al. 2012;Zahid et al. 2012). The metallicity evolution of galaxies will therefore change the position of galaxies on the BPT diagram as a function of redshift. Because optical classification of starburst and AGN is based on either empirical fits to local galaxies, or theoretical models developed for local galaxies, current BPT classification methods may not be applicable at high redshift.

Testing BPT classification methods at intermediate or high redshift has been difficult in the past. At z > 0.4, the [N II] and Hα lines are redshifted into the nearinfrared. With single-slit near-infrared spectroscopy, [N II]/Hα and [O III]/Hβ ratios have now been observed for small numbers of individual galaxies at high redshift (Teplitz et al. 2000;Finkelstein et al. 2009;Hainline et al. 2009;Bian et al. 2010;Rigby et al. 2011;Yabe et al. 2012). Stacked [N II]/Hα and [O III]/Hβ ratios of large numbers of galaxies have also been measured (Erb et al. 2006), and the first large NIR spectroscopic surveys are now being conducted (Trump et al. 2013, e.g.,).

The majority of galaxies at z > 0.4 show an offset towards larger [N II]/Hα and [O III]/Hβ ratios compared with local galaxies. This offset may be caused by a higher (2×) ionization parameter in high redshift galaxies. A larger ionization parameter may be produced by high nebular electron densities, a higher rate of star formation, a top-heavy initial mass-function (IMF), a high volume filling factor, and a large escape fraction of UV photons (Brinchmann et al. 2008b). Lehnert et al. (2009) suggest that the offset results from high gas densities and pressures that are similar to the most intense nearby SF regions locally, but spread over scales of 10-20 kpc in high redshift galaxies. Higher electron densities and ionization parameters have been measured for several high redshift galaxies (Hainline et al. 2009;Bian et al. 2010;Liu et al. 2008), but a large electron density by itself cannot explain the offset in all cases (see Rigby et al. 2011). Groves et al. (2006) use photoionization models to show that the emission-line ratios in some high redshift galaxies could be explained by a combination of starburst and AGN activity. A similar conclusion was reached by Trump et al. (2011) after combining HST spectra with Chandra X-ray data for a sample of galaxies at z ∼ 2. Concurrent starformation and AGN activity has also been found in some high redshift lensed galaxies (Wright et al. 2010).

The astronomical community is now on the cusp of obtaining the BPT diagnostic emission lines for large samples of z > 0.4 galaxies for the first time, thanks to new near-infrared multi-object spectrographs. New near infrared multi-object spectrographs include MOSFIRE on Keck (McLean et al. 2010), FMOS on Subaru (Kimura et al. 2010), MMIRS on Magellan [START_REF] Mcleod | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF]), and FLAMINGOS II on Gemini [START_REF] Dopita | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF]). These near-infrared multi-object spectrographs will allow the ionizing source in galaxies to be analyzed as a function of redshift in statistically significant numbers of galaxies, for the first time.

To interpret these new observations, a theoretical understanding of how the key BPT features may change at high redshift is needed. Such an understanding is essential for separating samples of star-forming galaxies from AGN at high redshift, and for tracking fundamental physical properties of the active galaxy population with redshift. This theoretical understanding is now achievable through modern cosmological hydrodynamic simulations, stellar evolution and photoionization models.

In this paper, we combine the predictions of cosmological hydrodynamic simulations models with stellar evolution, photoionization and shock models to predict how the BPT diagram may change with redshift. In Section 2, we describe the theoretical chemical evolution predictions used. Our simulations for the star-forming galaxy abundance sequence are given in Section 3. In Section 4, we model how the the line ratios of galaxies containing AGN may evolve with redshift. Section 5 gives theoretical predictions of the position of the BPT diagram within specific redshift windows. We investigate contamination from shocks, and we discuss the limitations of this work. Our conclusions are presented in Section 6. Throughout this paper, we adopt the flat Λ-dominated cosmology as measured by the 7 year WMAP experiment (h = 0.72, Ω m = 0.29; Komatsu et al. 2011).

CHEMICAL EVOLUTION OF ACTIVE GALAXIES

We utilize the chemical evolution estimates from Davé et al. (2011b,a). The Davé et al. (2011b) models use the GADGET-2 N-body + Smoothed Particle Hydrodynamic code (Springel et al. 2005) in a ΛCDM cosmology. This code incorporates gas cooling and heating processes, including the effects of metal line cooling (Oppenheimer & Davé 2006). Density-driven star formation is calculated using a Schmidt law (Schmidt 1959). Chemical enrichment from Type II supernovae, Type I supernovae, and AGB stars is included. Dave et al. use a Monte Carlo approach to model galactic outflows, where the mass loss due to outflows is related to the star formation rate and a variable mass loading factor.

We use the predicted evolution in the gas-phase chemical abundance for galaxies with stellar mass M * > 10 9 M ⊙ across 0 < z < 3. We parameterize the relative change in chemical abundances, ∆(log O/H), from z = 0 to an arbitrary redshift z by a 3rd order polynomial:

∆(log O/H) = -0.0013 -0.2287z + 0.0627z 2 -0.0070z 3 .
(1) This equation applies to star-forming galaxies with M * > 10 9 M ⊙ (chosen to avoid Malmquist bias). The 1σ error about equation 1 is ±0.1 dex at z = 0, falling to ±0.05 dex at z=3. Equation 1 assumes that there is no mass dependence in the chemical evolution for M * > 10 9 M ⊙ and over 0 < z < 3. Theory suggests that more massive galaxies undergo more rapid enrichment than less massive galaxies due to additional contributions from wind recycling (Davé et al. 2011a). However, we show in Yuan et al. (2012) that equation 1 fits the slope of the current metallicity history of galaxies for M * > 10 9 M ⊙ to within the observational errors. In future, when the chemical enrichment history of galaxies as a function of stellar mass is understood, a mass term could be included in equation 1.

THE STAR-FORMING ABUNDANCE SEQUENCE

We combine stellar evolutionary synthesis models with our MAPPINGS III photoionization models to generate theoretical limits to the expected ionizing radiation field as a function of redshift. We have tested this combination of stellar population synthesis and photoionization models extensively for local star forming galaxies (e.g., Kewley et al. 2001a;Levesque et al. 2010). Zero-age and 1 Myr old models are able to reproduce the [N II]/Hα and [O III]/Hβ line ratios in the majority of star-forming galaxies. However, the [S II]/Hα and [O I]/Hα line ratios require a harder ionizing radiation field than is available in current stellar population synthesis models (Levesque et al. 2010). This situation may be resolved when the effects of stellar rotation are incorporated into the stellar evolutionary tracks used by the population synthesis models. Initial investigations into the effect of stellar rotation on the ionizing radiation field at solar metallicity is promising (Levesque et al. 2012). Until the full set of stellar tracks with rotation become available, we limit our analysis of BPT evolution to the [N II]/Hα and [O III]/Hβ emission-line ratios which can already be reproduced using current stellar population synthesis models.

The Local Star-forming Abundance Sequence

The Starburst99 (SB99) models that we use are described in detail in Levesque et al. (2010) and Nicholls et al. (2012). Briefly, we apply a Salpeter IMF (Salpeter 1955) with an upper mass limit of 100 M ⊙ . The choice of IMF makes negligible difference on the optical emission-line ratios used in this analysis. We use the Pauldrach/Hillier model atmospheres, which employ the WMBASIC wind models of Pauldrach et al. (2001) for younger ages when O stars dominate the luminosity (< 3 Myr), and the CMFGEN atmospheres from Hillier & Miller (1998) for later ages when Wolf-Rayet (W-R) stars dominate. These stellar atmosphere models include the effects of metal opacities. We use the Geneva group "high" mass-loss evolutionary tracks [START_REF] Meynet | The Feedback of Chemical Evolution on the Stellar Content of Galaxies[END_REF]). These tracks include enhanced mass-loss rates that are applicable to low-luminosity W-R stars and can reproduce the blue-to-red supergiant ratios observed in the Magellanic Clouds (Schaller et al. 1992;[START_REF] Meynet | The Feedback of Chemical Evolution on the Stellar Content of Galaxies[END_REF]). Starburst99 generates a synthetic FUV spectrum using isochrone synthesis (e.g., Charlot & Bruzual 1991), in which isochrones are fitted to the evolutionary tracks across different masses rather than discretely assigning stellar mass bins to specific tracks. We use the zero-age instantaneous burst models because these models provide the best fit to the SDSS starforming galaxy sequence at z ∼ 0 [START_REF] Dopita | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF]. To match the nebular metallicities of our photoionization code, we interpolate between the STARBURST99 model grids as a function of metallicity.

We use our Mappings IV photoionization code (Binette et al. 1985;Sutherland & Dopita 1993;[START_REF] Dopita | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF] to model the interstellar medium surrounding the SB99 ionizing radiation field. We assume the solar abundance set of [START_REF] Asplund | Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis[END_REF]. This abundance set includes revised solar abundances for key elements, including oxygen and carbon. As in Levesque et al. (2010), the α-element abundances are assumed to scale linearly with metallicity, with the exception of Helium and Nitrogen. For helium, we include the stellar yield in addition to the primordial abundance from Pagel et al. (1992). For nitrogen, we assume primary and secondary nucleosynthetic components as measured by Mouhcine & Lançon (2002) and Kennicutt et al. (2003). The resulting N/H ratio is parameterized in Groves et al. (2004a). Metals are depleted out of gas phase and onto dust grains. The dust depletion factors are given in Groves et al. (2006), and are based on Kimura et al. (2003), who examined the metal absorption along several lines of sight within the local interstellar cloud.

Mappings uses either a plane parallel or spherical geometry with the ionization parameter defined at the initial edge of the nebula. For a spherical geometry, an effective ionization parameter q can be defined that takes into account the spherical divergence of radiation at the Stromgren radius R s (Strömgren 1939):

q ef f = Q H 0 4πR s 2 n H (2)
where Q H 0 is the flux of ionizing photons above the Lyman limit. If the thermal gas only occupies a fraction of the available volume, the ionization parameter can be defined in terms of a volume filling factor.

The ionization parameter q has units of velocity (cm/s)

and can be thought of as the maximum velocity ionization front that an ionizing radiation field is able to drive through a nebula. This dimensional ionization parameter is related to the dimensionless ionization parameter U through the identity U ≡ q/c. The dimensional ionization parameter is typically -3.2 < log U < -2.9 for local H II regions (Dopita et al. 2000) and star-forming galaxies (Moustakas 2006;Moustakas et al. 2010). In practice, all models with a similar effective ionization parameter produce very similar spectra, assuming all other parameters are held constant (e.g., Dopita et al. 2000).

To minimize small uncertainties produced by particular geometries, we calculate spherical models in which q is determined at the inner radius. The average ionization parameter is lower than this initial value, and is dependent on the ionization parameter of the initial radius. Models were run with pressure P/k = 10 5.5 cm -3 K, where k is the Boltzmann constant. In an ionized nebula, electron temperatures are ∼ 10 4 K, yielding a density of 10 -30 cm -3 , typical of giant extragalactic H II regions. Detailed photoionization, excitation, and recombination are calculated at increments (step size 0.03) throughout the nebula. The model completes when the hydrogen gas is fully recombined. A full description of the models, including geometry, is given in López-Sánchez et al. (2012).

Unlike the Kewley & Dopita (2002) models, our current models include a sophisticated treatment of dust, including the effects of absorption, grain charging, radiation pressure, and photoelectric heating of the small grains (Groves et al. 2004a). The latest version of Mappings incorporates a Kappa temperature distribution that is a more realistic representation of the electron temperature distribution in a turbulent ISM than a Stefan-Boltzmann distribution (Nicholls et al. 2012).

Figure 1 (left panel) shows our model grid in comparison to the local star-forming galaxy sequence from the Sloan Digital Sky Survey (SDSS) from Kewley et al. (2006b). The SDSS galaxies were selected within the redshift range of 0.04 < z < 0.1 to minimize aperture effects and Malmquist bias (Kewley et al. 2004).

We fit an equation of the form log[OIII]/Hβ = a log([NII]/Hα)+b +c to our model grid, where a, b, c are constants. This polynomial form was chosen for consistency with previous fits to theoretical model grids, including the hard ionizing radiation field from the Pègase stellar population synthesis models (Kewley et al. 2001a;Kauffmann et al. 2003). According to our models, the mean position of a galaxy along the local star-forming sequence is 2002) metallicity scale. We note that the following analysis depends only on the relative change in metallicity along the local starburst abundance sequence (equation 3). Different metallicity diagnostics have systematic offsets in their absolute abundance scale but relative metallicity measurements are conserved to within ±0.03 dex on average (Kewley & Ellison 2008). Thanks to this conservation of relative metallicities, the relative change in position along equation 3 is independent of the metallicity calibration used.

Factors that affect the Star-Forming Abundance

Sequence Individual galaxies that lie on the BPT diagram at high redshift are likely to evolve off the BPT diagram before z = 0. Thus, an observed change in the starforming abundance sequence will track the change in intrinsic galaxy properties across the star forming galaxy population at different epochs.

The position of our theoretical star-forming galaxy abundance sequence is determined by: (1) the shape of the ionizing radiation field, (2) the geometrical distribution of gas with respect to the ionizing sources, (3) the metallicity range, and (4) the electron density (pressure) of the gas. We discuss the effect of changing each of these quantities below.

Shape of the Ionizing Radiation Field

The stellar ionizing radiation field may change with redshift as a result of a change in the fraction of ionizing photons produced by the young stellar population. In a pure star forming galaxy, the hardness of the ionizing radiation field is related to the slope of the IMF, the age of the stellar population and the metallicity of the galaxy. A stellar population with a shallow initial mass function produces a hard ionizing radiation field, but there is no solid evidence for a change in IMF with redshift (see Bastian et al. 2010;Greggio & Renzini 2012).

The stellar population age is directly related to the shape of the ionizing radiation field. Hard ionizing radiation fields can be produced at ∼ 3 -5 Myr when the stellar population may be dominated by Wolf-Rayet stars (e.g., Schaerer 1996;Kehrig et al. 2008). Broad He II λ1640 emission has been observed in stacked spectra of Lyman Break Galaxies (Shapley et al. 2003) and in some individual high redshift gravitationally lensed galaxies (Cabanac et al. 2008;Dessauges-Zavadsky et al. 2010). The broad He II feature has been attributed to a significant contribution from O and Wolf-Rayet stars to the ionizing EUV radiation field at low metallicity (Brinchmann et al. 2008a). We note that radiative shocks can produce a narrow, nebular He II feature (Dopita et al. 2011;Lagos et al. 2012) which may be blended with the broad component produced in stellar atmospheres of luminous stars. Such blending could be difficult to distinguish at high redshift.

A hard ionizing radiation field has been linked with low metallicity in star forming galaxies (e.g., Campbell et al. 1986;Galliano et al. 2005 (2006b). Left: The colored curves show our new theoretical stellar population synthesis and photoionization model grid for star-forming galaxies based on a κ electron temperature distribution. Right: The red solid curve shows the mean star-forming sequence for local galaxies.

The shape of the red solid curve is defined by our theoretical photoionization models, while the position is defined by the best-fit to the SDSS galaxies. The ±0.1 dex curves (dashed lines) represent our model errors and contain 91% of the SDSS star-forming galaxies.

2010; Levesque et al. 2010

). There are several potential reasons for the correlation between metallicity and the hardness of the ionizing radiation field:

(1) High energy photons produced by metal-rich stars are absorbed by metals in the stellar atmosphere, known as metal blanketing. The preferential absorption of high energy photons yields a softer ionizing radiation field (e.g., González Delgado et al. 2005).

(2) The Hayashi track shifts to hotter effective temperatures at low metallicities, enabling metal-rich massive stars to maintain a higher effective temperature compared with metal-poor stars of similar spectral types (Elias et al. 1985;Levesque et al. 2006).

(3) Low metallicities correspond to lower mass loss rates, allowing low metallicity stars to remain on the main sequence for longer timescales [START_REF] Meynet | The Feedback of Chemical Evolution on the Stellar Content of Galaxies[END_REF]Maeder & Conti 1994).

(4) In isolated stars, rotational mixing causes heavy mass loss. This mass loss produces bluer colors in the red supergiant phase, lowering the mass limit required for a star to enter the Wolf-Rayet phase (Levesque et al. 2012). Thus, a population of rotating massive stars will contain a larger fraction of hot, massive stars to contribute ionizing photons to the stellar radiation field. This rotational hardening is a function of metallicity, with more significant hardening in metal-poor environments [START_REF] Leitherer | IAU Symposium[END_REF].

(5) In binary stars, efficient mass-transfer can spin-up the rotation of the companion star, causing similar mixing effects as in rotating isolated stars (de Mink et al. 2009;[START_REF] Eldridge | [END_REF]. Whether a rapidly rotating star can spin down depends on stellar winds, which are weaker at low metallicities due to metal opacity. Rotation and binarity are not yet included in stellar population synthesis models for a range of metallicities.

A hard ionizing radiation field can also be produced by contamination from an active galactic nucleus or radiative shocks. Slow shocks (100-200 km/s) associated with galactic-scale winds have been observed both locally (Rich et al. 2010(Rich et al. , 2012) ) and at high redshift (Yuan et al. 2012). We address these two possibilities in Section 4 and Section 5.1.

It is unclear whether high redshift star-forming galaxies have a harder radiation field than local galaxies at the same metallicity. If a harder ionizing radiation field exists in high redshift galaxies, both the [O III]/Hβ and [N II]/Hα line ratios will be affected. A larger fraction of photons with energies above the ionization potentials of [N II] and [O III] will raise the [N II]/Hα and [O III]/Hβ line ratios because there are more photons available to ionize nitrogen and oxygen. The [O III] emission-line is substantially more sensitive to the hardness of the EUV ionizing radiation field than [N II] because the difference in ionization potentials is large (35eV c.f. 14.5eV). In Figure 2, we illustrate the effect of raising the hardness of the ionizing radiation field in galaxies along the local abundance sequence (orange dashed line). The harder radiation field moves the entire abundance sequence above and to the right on the BPT diagram. A harder ionizing radiation field could account for part (or all) of the large [N II]/Hα and [O III]/Hβ ratios seen at high redshift.

-2.0 -1.5 -1.0 -0.5 0.0 0.5

LOG ([NII]/Hα

Geometrical distribution of the gas

The geometrical distribution of the gas around the ionizing source can be changed by shocked stellar winds. Stellar winds clear ionized gas from the interior of an H II region. Because highly ionized species, such as [O III], form preferentially at the inner radii of H II regions, shocked stellar winds will lower the effective ionization parameter of the nebula (Yeh & Matzner 2012).

Our theoretical photoionization models assume a uniform medium, where the geometrical distribution of the gas is approximated by a volume filling factor. In reality, H II regions may be clumpy and porous. Ultracompact H II regions in the Milky Way contain a porous interstellar medium (Kurtz et al. 1999;Kim & Koo 2001). Radio observations indicate that young H II regions in low metallicity galaxies also have a clumpy and porous dis-tribution of gas (Johnson et al. 2009). A clumpy, porous medium allows some ionizing photons to escape the nebula without being absorbed by the interstellar gas. In this scenario, the effective ionization parameter may be lowered, depending on the escape fraction and the optical thickness of the porous medium.

Our models are calculated for radiation-bounded H II regions. In the radiation-bounded scenario, the model finishes when hydrogen is completely recombined. If H II regions are density bounded (i.e. the H II region density is sufficiently low that the stars can ionize the entire nebula), the [O III] ionization zone is likely to be largely unaffected, but the [N II] excitation and Hydrogen recombination zones may be shortened. Therefore, the [O III]/Hβ ratio may be larger, while the [N II]/Hα ratio will be similar to or smaller than observed in a radiationbounded nebula. Detailed ionization parameter mapping of the H II regions in nearby galaxies can constrain whether the H II regions are radiation bounded or density bounded. A mixture of radiation-bounded and density-bounded H II regions have been observed in the local group (Pellegrini et al. 2012). Nakajima et al. (2012) suggest that Lymanα emitters at high redshift contain density-bounded H II regions. It is unclear whether density-bounded nebulae are common in normal star-forming galaxies, either locally or at high redshift.

The geometrical distribution issues described above can be considered in terms of an effective ionization parameter. We discuss the effect of changing the ionization parameter on the [N II]/Hα and [O III]/Hβ line ratios in Section 3.2.5.

The Metallicity Range of Galaxies

The spread in metallicity across a galaxy sample determines the length of the star-forming abundance sequence. Galaxy samples that span a small range of metallicities occupy only a portion of the star-forming abundance sequence. For example, interacting or merging galaxies typically have lower central metallicities due to largescale gas infall (Kewley et al. 2006a(Kewley et al. , 2010;;Rupke et al. 2010b;Ellison et al. 2010;Rich et al. 2012;Scudder et al. 2012). On the other hand, low metallicity galaxies, such as blue compact dwarfs, occupy the highest positions (lowest metallicities) on the star-forming abundance sequence (Levesque et al. 2010).

The metallicity range of high redshift galaxies is unknown. Samples selected from rest-frame blue colors or the Lyman break may be missing a population of faint, low metallicity galaxies and a population of dusty, metalrich star forming galaxies. Gravitationally lensed samples probe fainter galaxy samples and a broader range of metallicities, within current instrumentation detection limits (see Figure 5 in Yuan et al. 2012, for a comparison of current instrumentation detection limits for lensed and non-lensed samples).

The Electron Density

In an isobaric density distribution, the density is defined in terms of the ratio of the mean ISM pressure, P , and mean electron temperature, T e , through n e = P Tek . For an ionized gas, the electron temperature is ∼ 10 4 K and the density is simply determined by the ISM pressure.

The SDSS abundance sequence is fit by our photoionization models with an electron density of n e = 10 -10 2 cm -3 , typical of local H II regions [START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF]. Brinchmann et al. (2008b) shows that the distance away from the star-forming abundance sequence correlates strongly with electron density; SDSS galaxies with larger electron densities (n e = 10 2 cm -3 ) lie above and to the right of the mean SDSS star-forming abundance sequence. Highly active star-forming galaxies, such as warm infrared galaxies and luminous infrared galaxies have dense ionized gas (10 2 -10 3 cm -3 ; Kewley et al. 2001b; Armus et al. 2004) and lie above and to the right of the SDSS star-forming abundance sequence (Yuan et al. 2010).

Many high redshift galaxies also have dense nebular gas (∼ 10 3 cm -3 ) (e.g. Brinchmann et al. 2008b;Hainline et al. 2009;Bian et al. 2010;Rigby et al. 2011;Wuyts et al. 2012a) and galaxies with high Hα surface brightness (Lehnert et al. 2009;Le Tiran et al. 2011). Figure 2 shows the effect of raising the electron density from n e = 10 cm -3 (solid red line) to n e = 1000 cm -3 (dotdashed green line). Both the [O III]/Hβ and [N II]/Hα ratios are larger at higher electron densities due to the increased rate of collisional excitation. We conclude that high electron densities could account for at least part of the enhanced [N II]/Hα and [O III]/Hβ ratios observed in current samples of high redshift galaxies.

Ionization Parameter

Much previous research into the differences between high redshift and local galaxy emission-line ratios has focused on the ionization parameter (Brinchmann et al. 2008b;Liu et al. 2008;Hainline et al. 2009). The ionization parameter is a useful tool for comparisons among H II regions or galaxies because it can be readily measured using the ratio of high ionization to low ionization species of the same atom (such as [O III]/[O II]). However, it is important to understand the fundamental physical properties that govern the ionization parameter. The ionization parameter is determined by the hydrogen ionizing photon flux, and the pressure of the interstellar medium.

The number of hydrogen ionizing photons per unit area can be changed by either scaling the ionizing radiation field (i.e. raising or lowering the luminosity of the stellar population) or by modifying the shape (hardness) of the ionizing radiation field. Because the ionization parameter, electron density and radiation field hardness are interrelated, isolating the cause(s) of a shift towards larger [N II]/Hα and/or [O III]/Hβ ratios at high redshift is non-trivial. Furthermore, different geometrical gas distributions can mimic a higher or lower observed ionization parameter. Several of these processes may contribute to the observed ionization parameter in starforming galaxies.

Local galaxies with high measured or inferred ionization parameters include Wolf-Rayet galaxies (Brinchmann et al. 2008b), strong Hα emitter galaxies (Shim & Chary 2012), and low metallicity galaxies (Kewley et al.

2007

). In Wolf-Rayet galaxies, the high ionization parameter is caused by the shape of the ionizing radiation field produced in Wolf-Rayet atmospheres with line blanketing (González Delgado et al. 2002). Strong Hα emitting galaxies may have high ionization parameters due to a high ISM pressure (Brinchmann et al. 2008b;Shim & Chary 2012), and/or low metallicity.

The ionization field can be scaled by raising the star formation rate of a stellar population. High redshift galaxies typically have larger star formation rates than local galaxies (Shapley et al. 2005;Hainline et al. 2009). Hainline et al. (2009) suggest that these larger star formation rates may lead to a larger reservoir of ionizing photons, leading to a higher ionization parameter.

Alternatively, a large ionization parameter may be a natural consequence of a low metallicity environment. The shape of the ionizing radiation field is a strong function of time since the most recent burst of star formation. The shape of the ionizing radiation field is also determined by the fraction of ionizing photons that are either absorbed by the stellar atmospheres and surrounding gas, and conversely, the fraction of ionizing photons that escape the nebula without being absorbed by ISM. These effects combine to produce a correlation between ionization parameter and metallicity (Dopita & Evans 1986;Dopita et al. 2000;Kewley & Dopita 2002). Dopita et al. (2006) show theoretically that the ionization parameter inversely correlates with the metallicity of an H II region, Z, according to q ∝ Z -0.8 . The reason for this correlation is twofold;

(1) At high metallicity, the stellar wind opacity is larger. This large opacity absorbs a greater fraction of the ionizing photons, reducing the ionization parameter in the surrounding H II region.

(2) At high metallicity, stellar atmospheres scatter the photons emitted from the photosphere more efficiently. This scattering efficiently converts luminous energy flux into mechanical energy flux at the base of the stellar wind. The reduced luminous energy flux available to ionize the surrounding H II region is observed as a lower ionization parameter.

Many high redshift galaxies have high ionization parameters (-2.9 < log U < -2.0 Pettini et al. 2001;Lemoine & Pelletier 2003;Maiolino et al. 2008;Hainline et al. 2009;Erb et al. 2010;Richard et al. 2011;Wuyts et al. 2012a). Nakajima et al. (2012) compare the ionization parameter and metallicity (as traced by the [O III]/[O II] vs R 23 line ratios) for a large sample of high redshift galaxies, including Lyman Break, lensed, and Lyman-α emitting galaxies. Their comparison shows that the majority of high redshift galaxies have similar ionization parameters to local galaxies at the same metallicity. Thus, the high ionization parameters at high redshift may simply be a natural consequence of a low metallicity environment.

The effect of ionization parameter on the [N II]/Hα and [O III]/Hβ ratios depends on the metallicity regime. Figure 1 (left panel) shows that at low metallicities (log([N II]/Hα)< -0.9), a larger ionization parameter will raise the [O III]/Hβ ratio, while the [N II]/Hα ratio remains constant within ±0.2 dex. At high metallicities (log([N II]/Hα)> -0.7, a larger ionization parameter will reduce the [N II]/Hα ratio while [O III]/Hβ remains roughly constant (to within ±0.1 dex). This effect is illustrated in Figure 2. For the metallicities of most high redshift galaxies (log(O/H) + 12 < 8.6; Erb et al. 2006;Wuyts et al. 2012b;Yuan et al. 2012), a high ionization parameter will raise the [O III]/Hβ ratio such that high redshift galaxies will lie above the local abundance sequence, as observed. Therefore, a high ionization parameter can account for part of the large [O III]/Hβ and [N II]/Hα ratios observed in high redshift galaxies.

The High Redshift Star-forming Abundance

Sequence Although a hard ionizing radiation field is linked with low metallicity, it is unclear whether hard ionizing radiation fields are a feature of high redshift star-forming galaxies. Figure 2 shows that the combination of a larger ionization parameter at low metallicity and a higher electron density may mimic the effects of a harder ionizing radiation field. Rest-frame UV spectroscopy of high metallicity galaxies at z > 2 is needed to determine whether the large [N II]/Hα and [O III]/Hβ ratios seen in high redshift galaxies are caused by a higher electron density and ionization parameter or the combined effects of a higher electron density, ionization parameter, and a harder EUV ionizing radiation field.

Further, selection biases may be limiting the current range of electron densities and ionization parameters observed at high redshift. The combination of surface brightness dimming and current observational sensitivities limits high redshift samples to galaxies with high emission-line surface brightness. This limit biases samples towards galaxies with high star formation rates (SFR> 4M ⊙ /yr), and more intense star formation. Local starburst galaxies are characterized by higher pressures than normal H II regions (Kewley et al. 2001b). High pressures of 10 6 < P/k < 10 7 implies gas densities of 100 -1000 cm -3 . A greater fraction of starburst galaxies in high redshift samples would produce larger [N II]/Hα and [O III]/Hβ ratios when compared with local samples of normal star-forming galaxies.

Gravitationally lensed samples can probe an order of magnitude fainter in star formation rate (SFR > 0.4M ⊙ /yr; Richard et al. 2011). Consequently, gravitationally lensed samples cover a larger range of [N II]/Hα and [O III]/Hβ ratios than non-lensed samples (Hainline et al. 2009;Richard et al. 2011;Wuyts et al. 2012a). However, current lensed samples still do not sample the full range of star formation rates seen in the global spectra of local normal star forming galaxies ( SFR > 0.01 M ⊙ /yr Kewley et al. 2002). A statistically significant sample of high redshift galaxies that covers a broad range of metallicity and stellar mass is required to determine the full range of [N II]/Hα and [O III]/Hβ ratios at high redshift. Until such samples are available, we assume a lower limit and an upper limit to the star forming abundance sequence at high redshift (z ∼ 3).

3.3.1. Lower Limit Abundance Sequence (Normal SF conditions)

Our lower limit abundance sequence is calculated assuming that galaxies at a given redshift have the same shape EUV radiation field, the same range of electron densities, and the same relationship between metallicity and ionization parameter as local star-forming galaxies. In this very conservative scenario, the star-forming abundance sequence is given by equation 3. Chemical evolution moves galaxies down along the star-forming galaxy sequence towards smaller [O III]/Hβ ratios and larger [N II]/Hα ratios.

3.3.2. Upper Limit Abundance Sequence (Extreme SF conditions)
We use PEGASE 2 models (Fioc & Rocca-Volmerange 1999) to provide the upper limit to the ionizing radiation field at z = 3. These models provide the hardest radiation field because the models use the Clegg et al. (1992) planetary nebula nucleus (PNN) atmospheres for stars with high effective temperatures (T > 50, 000K). Clearly, the ionizing spectrum of active star-forming galaxies is dominated by the emission from the atmospheres of massive stars, not planetary nebulae. However, we use the PNN atmospheres as a substitute for massive star stellar atmospheres, in the absence of the full suite of stellar atmospheres with the effects of stellar rotation. The PE-GASE 2 models use the Padova tracks (Bressan et al. 1993) and the OPAL opacities (Iglesias et al. 1992). The application of PEGASE 2 with our Mappings III code is described in detail in Kewley et al. (2001a).

If we assume that the hard radiation field sets the upper limit at z = 3 for star-forming galaxies, then the evolution of the star-forming abundance sequence between 0 < z < 3 follows log([OIII]/Hβ) = 1.1 + 0.03z + 0.61 log([NII]/Hα) + 0.08 -0.1833z

.( 5)

This upper limit is consistent with the limit that we derive from a combination of a larger ionization parameter (-2.9 < log U) and larger electron density (N e ∼ 1000 cm -3 ) for metallicities below log(O/H) + 12∼ 8.8. According to our models, the metallicity of a galaxy at redshift z on the star-forming abundance sequence given by equation 5 is

12 + log O/H = 8.97 + 0.0663z - log(O3N 2)(0.32 -0.025z) ( 6
)
where O3N 2 = [OIII]/Hβ/[NII]/Hα. In Figure 3, we show the evolution of the upper limit star-forming abundance sequence with redshift given by equation 5. Figure 3 should not be used over small redshift intervals with small samples to test the evolution of the galaxy population because the model errors (±0.1 dex) and spread in ionization parameter log(U ), (±0.1 dex) at a given redshift are larger than the predicted evolution of the sequence for an ensemble of galaxies across each redshift interval.

Abundance Sequence Shape

The metallicity range (or spread) of the star-forming population determines the length of the abundance sequence. In cosmological hydrodynamic simulations, the evolution of the lower metallicity bound is dominated by chemical enrichment of the most metal poor galaxies in the galaxy population. The upper metallicity bound may be influenced by gas inflows from the intergalactic medium such that the upper metallicity bound might fall with redshift (Nagamine et al. 2001).

Again, we assume that at least some galaxies at z ∼ 3 have reached the level of enrichment of galaxies in the local universe, and we assume that the mean metallicity of the galaxy population evolves according to equation 1.

The scatter of galaxies about the abundance sequence is dominated by the range in ionization parameters and electron densities in the galaxy population at a given redshift. How the spread in these properties changes with redshift is unknown. For simplicity, we assume that the spread in these properties about the mean is constant with redshift and is ±0.1 dex.

THE MIXING SEQUENCE

Local active galaxies form two branches on the BPT diagram. While pure star-forming galaxies lie along the abundance sequence, galaxies with a contribution from a non-thermal radiation field form a sequence that extends towards the pure AGN region of the BPT diagram (i.e. towards large [N II]/Hα and [O III]/Hβ ratios). This sequence can be produced by either a mixture of gas ionized by hot stars and gas ionized by an AGN (Groves et al. 2004b), or a mixture of gas ionized by hot stars and gas ionized by radiative shocks (Kewley et al. 2001a;Allen et al. 2008). We therefore refer to this sequence as the mixing sequence, where mixing refers to the mixture of a soft photoionizing radiation field from star-formation and a hard non-thermal radiation field from either AGN or shocks. In this section, we focus on starburst-AGN mixing. We investigate starburst-shock mixing in Section 5.1.

AGN Models

The [N II]/Hα and [O III]/Hβ line ratios of a galaxy purely ionized by an AGN are influenced by the metallicity of the surrounding narrow-line region, the shape of the AGN ionizing radiation field (characterized by a power-law), and the ionization parameter. We use the Mappings III dusty, radiation pressure dominated models of Groves et al. (2004a).

These models use the same abundance set and depletion factors as our starburst models (Section 3.1). The ionizing spectrum is based on a power-law,

F µ ∝ µ α (7)
where the frequency µ is defined over 5eV < µ < 1000eV with four values of the power-law index α (-1.2, -1.4, -1.7, -2.0).

The shape of the ionizing radiation field may change with metallicity; a smaller fraction of metals and dust may change the AGN torus structure and potentially alter the accretion disk. Unfortunately, the effect of metallicity on the AGN ionizing radiation field is unknown. For simplicity, we assume the same ionizing spectral shape for all metallicities.

In isobaric radiation-pressure dominated models, the local density varies continually throughout the models. In this case, the hydrogen density is defined as the density of the [S II] emission zone, which is located close to the ionization fronts in the Narrow Line Region (NLR). In this zone, we assume that the electron density is n H = 10 3 cm -3 .

The total radiative flux entering the narrow line region cloud is determined by the ionization parameter, U N L defined at the inner edge of the nebula. The ionization parameter range is 0.0 < log(U N L ) < -4.0 with intervals of -0.3, -0.6, and -0.1 dex. The NLR models are truncated at a column density of log[N(H I)] = 21, consistent with observations of NLR clouds (Crenshaw et al. 2003).

In Figure 4, we show the position of the dusty radiation-pressure dominated AGN models as a function of metallicity on the BPT diagram. The position of the AGN models changes substantially with metallicity. At low metallicity, the AGN models move towards lower [N II]/Hα ratios and the spread in [O III]/Hβ ratios becomes smaller. The reason for this effect is twofold; (1) at low metallicity, nitrogen shifts from a secondary nucleosynthetic element to a primary nucleosynthetic element, and (2) the rise in electron temperature offsets the fall in oxygen abundance, yielding a roughly constant mean [O III]/Hβ ratio as a function of metallicity.

Figure 4 indicates that it will be impossible to distinguish low metallicity galaxies containing AGN (i.e. metallicities log(O/H) + 12 8.4) from low-metallicity star-forming galaxies using the BPT diagram.

A sample of AGN galaxies at a given NLR metallicity may span a range of ionization parameter and power-law indices. We therefore use the full suite of ionization parameter and power-law indices to define our 100% AGN region on the BPT diagram. We have tested these photoionization models for local galaxies containing AGN; the models produce remarkable agreement with the observed position of AGN on the BPT diagram (Allen et al. 1998;Groves et al. 2004b;Kewley et al. 2006b). scale). Indeed, low metallicity AGN are extremely rare in local galaxies (Groves et al. 2006). The star forming abundance sequence and the AGN region on the BPT diagram are linked together by a "mixing sequence" which can be modelled by a combination of starburst and power-law AGN photoionization models (e.g., Groves et al. 2004b). The overlap region between the starforming abundance sequence and the AGN mixing sequence provides an additional constraint on the metallicity of local galaxies containing AGN. The local SDSS mixing sequence extends from only the most metal-rich star-forming galaxies (9.0 < log(O/H) + 12 < 9.2) to the AGN region. It is not clear whether AGN are only found in the most metal rich galaxies at higher redshift.

AGN Metallicity

The combination of metallicity gradients and surface brightness dimming may affect the observed metallicity of galaxies along the mixing sequence. The narrow-line region surrounding an AGN is typically 1-5 kpc (Bennert et al. 2006), smaller than the typical radius sampled by global spectra of local pure star-forming galaxies (∼ 8 kpc on average; Kewley et al. 2004). Therefore, in galaxies with very steep metallicity gradients, the AGN narrow-line region may be more enriched than the extended star-forming regions which probe a larger radius, on average. For example, a steep metallicity gradient (-0.15 dex kpc -1 ) can produce a 0.2 dex metallicity difference from spectra measured within circular apertures of 1 kpc and 5 kpc. On the other hand, a flat metallicity gradient would give the same metallicity for the star-forming gas and the NLR, regardless of the size of the star-forming and narrow-line regions. Flat metallicity gradients can be produced by large-scale gas flows triggered by galaxy interactions (Kewley et al. 2006a(Kewley et al. , 2010;;Rupke et al. 2010b;Rich et al. 2012). Theory indicates that metallicity gradients systematically flatten after first pericenter and flatten again during final coalescence (Rupke et al. 2010a;Torrey et al. 2012). The metallicity gradient can be steepened by a late nuclear starburst, providing that outflows have not already removed a substantial fraction of the nuclear star-forming gas (Torrey et al. in prep).

We note that local disk galaxies with typical gradients (-0.04 ± 0.09 dex kpc -1 Zaritsky et al. 1994;van Zee et al. 1998;Rupke et al. 2010a) are unlikely to produce a large difference in the NLR and star forming regions. Local early-type galaxies have even shallower gradients (Henry & Worthey 1999).

Current observations of metallicity gradients at high redshift are limited by small numbers and low angular resolution observations. Metallicity gradients in some high redshift lensed isolated star-forming galaxies appear to be significantly steeper than local galaxies (Yuan et al. 2011;Jones et al. 2013), but in other high redshift galaxies, the gradients appear to be flat or inverted (Cresci et al. 2010;Swinbank et al. 2012;Queyrel et al. 2012). The steepest gradients were obtained by observing gravitationally lensed galaxies with integral field spectrographs with laser-guided adaptive optics (Jones et al. 2010;Yuan et al. 2011;Jones et al. 2013). Observational issues such as poor angular resolution (> 0.1"F W HM ) or low S/N can lead to weak line smearing which systematically flattens or inverts metallicity gradients that are intrinsically steep (Yuan et al. 2013). In the absence of robust observational or theoretical constraints on the cosmic evolution of the metallicity gradient in galaxies, we consider two extreme cases:

• Case 1: Metal-rich NLR The AGN narrow line region at z = 3 has reached the level of enrichment seen in the narrow line region of local galaxies. This extreme case implies that galaxies at high redshift have a steeper metallicity gradient on average than observed in local galaxies. In this scenario, the metallicity gradient of galaxies would flatten with time until z = 0.

• Case 2: Metal-poor NLR The AGN narrow line region at z = 3 is at the same average metallicity as the surrounding star-forming gas. This extreme case implies a flat metallicity gradient on average at high redshift. If the average galaxy metallicity gradient is flat at high redshift, the gradient steepens with time until z = 0. . We derive starburst-AGN mixing sequences for each redshift using our AGN and starburst models. For both cases, we assume that AGN reside in the most metal-rich galaxies at any epoch, as observed locally. We assume that the fraction of the starburst sequence that intersects the AGN mixing sequence covers the top 0.2 dex in metallicity. If this assumptiondoes not hold at high redshift, then we would expect to see galaxies to the left of our mixing sequence (i.e. towards lower [N II]/Hα line ratios than spanned by the mixing sequence). We fit the mixing sequences with a 2nd or 3rd order least-square polynomial:

We interpolate between the model grids in

y = a + bx + cx 2 + dx 3 (8)
where y = log [OIII]/Hβ and x = log [NII]/Hα. The constants ad are given for each mixing line boundary at each redshift in Table 1, along with the range of [N II]/Hα and [O III]/Hβ values over which the mixing sequences are defined. The four scenarios tabulated in Table 1 are described in detail in Section 5.

THE COSMIC BPT DIAGRAM

With our theoretical predictions for the star-forming abundance sequence and the AGN mixing sequence, we can predict how the the BPT diagram might appear at different redshifts. We have two extreme cases for the star-forming galaxies and two extreme cases for the AGN NLR metallicity, which we briefly summarize below.

Star forming galaxies at high redshift (z = 3) may have ISM conditions and/or an ionizing radiation field that are the same as local galaxies (normal ISM conditions) or that are more extreme than local galaxies (extreme ISM conditions). Extreme conditions in star-forming galaxies can be produced by either a larger ionization parameter and a more dense interstellar medium, and/or a harder ionizing radiation field.

The AGN narrow line region at high redshift may have reached the level of enrichment seen in local galaxies (metal-rich), or it may evolve similarly to the starforming gas. In the latter case, the AGN narrow-line region at high redshift would be more metal-poor than local AGN narrow-line regions.

Combining our two sets of extreme cases gives four limiting scenarios for the position of galaxies on the BPT diagram at each redshift:

• Scenario 1: Normal ISM conditions, and metalrich AGN NLR at high-z.

• Scenario 2: Normal ISM conditions, and metalpoor AGN NLR at high-z.

• Scenario 3: Extreme ISM conditions, and metalrich AGN NLR at high-z.

• Scenario 4: Extreme ISM conditions, and metalpoor AGN NLR at high-z.

In Figure 5, we show how the abundance sequence and the mixing sequence are expected to evolve with redshift given our four limiting scenarios (columns 1-4 in Figure 5). The local SDSS sequence boundaries are shown for reference (blue dashed lines). We discuss the predictions given by each limiting scenario individually below.

Scenario 1: In this scenario, we have assumed that the ISM conditions (and/or the ionizing radiation field) are constant as a function of redshift, and the AGN NLR at z = 3 has already reached the level of enrichment seen in local AGN NLRs. Therefore, the only change to the position of galaxies on the BPT diagram is a broadening of the mixing sequence at the intersection with the abundance sequence. This broadening is due to the lower mean and larger range of metallicities in star-forming galaxies at higher redshift.

Scenario 2: In this picture, we have assumed that the ISM conditions (and/or the ionizing radiation field) are constant as a function of redshift, and that the AGN NLR metallicity evolves with time in the same way that pure star-forming galaxies evolve (equation 1). This implies that AGN galaxies at z = 3 would have flat metallicity gradients that steepen with time. At high redshift, the mixing sequence occupies lower [N II]/Hα ratios than local galaxies due to the metallicity sensitivity of the [N II]/Hα ratio. The 100% AGN region is broad in this scenario, spanning a large range of [N II]/Hα ratios (-1.0 <[N II]/Hα< 0). The breadth of the AGN region illustrates the effect of differing AGN power-law indices and ionization parameters at lower metallicities.

Scenario 3: Here, we have assumed that the the ionizing radiation field or the ISM conditions in star forming galaxies become more extreme at high redshift, and that the AGN NLR at z = 3 has already reached the level of enrichment seen in local AGN NLRs (i.e. a steep metallicity gradient at high-z). In this scenario, the position of the star-forming abundance sequence rises from z = 0 to z = 3 towards larger [O III]/Hβ ratios, while the AGN mixing sequence becomes shorter. At z = 3, the abundance sequence and AGN mixing sequence form almost a single horizontal sequence across the BPT diagram. If a Mixing line boundaries extending from the star-forming abundance sequence (0% AGN) to the AGN region (100% AGN) on the BPT diagram. The coefficients a, b, c, d are defined according to y = a + bx + cx 2 + dx 3 where y = log([OIII]/Hβ) and x = log([NII]/Hα). The range of x and y values (xr, yr ) over which these polynomials are valid are shown. Scenarios (1-4) correspond to the 4 limiting scenarios for our theoretical models, described in Section 5 and shown as columns in Figure 5.

this scenario exists in high redshift galaxies, classifying galaxies into star-forming and AGN using the BPT diagram would be extremely difficult for z ≥ 2.5. Scenario 4: In this scenario, we have assumed that high redshift star forming galaxies have more extreme ISM conditions and/or a harder ionizing radiation field, and that the AGN NLR metallicity evolves with time in the same way that pure-star-forming galaxies evolve (i.e. a flat metallicity gradient at high redshift). According to these predictions, the star-forming abundance sequence rises at larger redshift, while the AGN mixing sequence become substantially broader and shorter. At z = 3, it would be extremely difficult to distinguish between high metallicity star-forming galaxies and galaxies containing AGN using the BPT diagram.

Figure 5 shows that distinguishing between these four scenarios is possible at specific redshifts. Scenarios (1) and ( 2) could be easily distinguished with intermediate or high redshift samples (z ≥ 0.8). All four scenarios predict substantially different abundance and mixing sequence positions at z ≥ 2.5. Thus, observations of the rest-frame optical emission-line ratios for statistically significant samples of active galaxies at z ≥ 2.5 are likely to yield important information on the ionizing radiation field and/or the ISM conditions in star-forming galaxies, as well as the nature of metallicity gradients in galaxies containing AGN at high redshift.

The Effect of Shocks

Shocks associated with galactic winds may substantially raise the [N II]/Hα emission-line ratio observed in the global spectra of galaxies at high redshift. Radiative shocks are produced by a variety of astrophysical sources in galaxies, including starburst or AGN-driven outflows, cloud-cloud collisions from galaxy interactions, and jetcloud collisions. In Rich et al. (2010) and Rich et al. (2011), we showed that outflows in infrared luminous galaxies can drive gas clouds into the ambient gas in the outer regions of galaxies, creating shock fronts which ionize and excite the gas. This shocked gas produces strong emission lines at red wavelengths ([NII], [SII], [OI]). Similar shock excitation has recently been found in local luminous and ultraluminous infrared galaxies (Sharp & Bland-Hawthorn 2010;Soto et al. 2012;Weistrop et al. 2012). Shocks have also been observed in members of galaxy clusters (Farage et al. 2010;McDonald et al. 2012). The prevalence of shock excitation in less massive, more normal star-forming systems is unknown. Numerous kinematic and detailed spectroscopic studies indicate that galactic winds exist in a significant fraction of galaxies at z∼ 1 (Kornei et al. 2012), and that winds are prevalent in z > 2 galaxies (e.g., Pettini et al. 2000Pettini et al. , 2001;;Nesvadba et al. 2007;Steidel et al. 2010;Le Tiran et al. 2011;Genzel et al. 2011). The fraction of galactic winds that drive shock excitation is unclear. However, recent integral field spectroscopic studies have revealed the first evidence of shock excitation in high redshift gravitationally lensed galaxies (Yuan et al. 2012;Newman et al. 2012). Therefore, shock excitation may be an important contributor to the optical line ratios in high redshift galaxies.

Models of photoionized shocks were developed by Sutherland & Dopita (1993) and Dopita & Sutherland (1995, 1996), and are described in detail in Dopita et al. (2003). We apply the fast shock models of Allen et al. (2008) with shock velocities of 300-1000 kms -1 and the slow shock models described in Rich et al. (2010Rich et al. ( , 2011) ) covering velocities of 100-200 kms -1 . In slow shocks (v 200 kms -1 ) the shock front travels faster than the photoionization front emitted by the shocked gas. In this case, the optical emission-line ratios are directly dependent on the shock velocity. In fast shocks (v 200 kms -1 ) , the ionizing flux emitted by the shock front is large, producing a supersonic photoionization front that moves ahead of the shock front and preionizes the gas. This photoionization front is referred to as the photoionizing precursor. The photoionizing precursor can make a significant contribution to the optical emission-lines. The precursor is usually combined with the theoretical emission from the shocked gas. For our fast shock models, we assume a 50:50 ratio of shock to precursor emission. This ratio of shock to precursor emission produces [O III]/Hβ line ratios that are close to the maximum observed in local active galaxies (log([O III]/Hβ)∼ 1.0). We assume a magnetic field strength at pressure equipartition, and a preshock density of n ∼ 1.0 cm -3 . Allen et al. (2008) shows that changing the preshock density from n ∼ 0.01 cm -3 to n ∼ 1000 cm -3 has negligible effect on the position of the equipartition shock+precursor models on the BPT diagram.

The full range of shock parameters is explored within the BPT diagram in Allen et al. (2008). In Figure 6, we show how the position of the fast (red) and slow (blue) shocks change as a function of metallicity and shock velocity. The colored curves indicate the location of galaxy emission that is fully dominated by shocks. In reality, shocks rarely dominate a galaxy's global emission. If a galaxy contains shocks, there may also be a contribution from star-formation and/or AGN. In the case of composite starburst+shock activity, a galaxy will lie along a mix- Galaxies containing a mixture of ionizing sources, such as shocks and star-formation will lie along mixing sequences between the star-forming sequence and the 100% shock models. The position and shape of the mixing sequences will depend on the metallicity of the galaxy and the shock velocity. Examples of typical mixing sequences are given in Rich et al. (2010Rich et al. ( , 2011)).

ing sequence between the pure star-forming sequence and the 100% shock models. Examples of mixing sequences for nearby galaxies containing both star-formation and shocks are given in Rich et al. (2010Rich et al. ( , 2011)).

Figure 6 shows that the location of the shock models on the BPT diagram is a strong function of metallicity, similar to our AGN models. The [N II]/Hα ratio becomes 0.3 dex smaller between metallicities typical of local starburst galaxies (log(O/H) + 12∼ 9.0) and metallicities typical at high redshift (log(O/H) + 12∼ 8.7). Over this metallicity range, the [O III]/Hβ ratio remains roughly constant within the model errors (±0.1 dex) for shocks.

Shocks associated with galactic outflows are often spatially distributed throughout the galaxy and are therefore likely to trace regions of similar metallicity to the star forming gas. Powerful galactic winds may force metals out of galaxies, while less powerful outflows may act as "fountains", redistributing metals from the nuclear regions to larger radii (see Putman et al. 2012, for a review). Assuming that shocked regions occupy similar metallicities to star-forming regions, we calculate the position of shocks on the cosmic BPT diagram, shown in Figure 7.

The position of the 100% fast and slow shock models on Figure 7 reveals how a contribution from shock emission (black and purple curves) can mimic non-shock power sources. Galaxies with a contribution from both star formation and shocks will lie along a mixing sequence between the starburst sequence and the shock models. Although the exact location of the starburst-shock mixing sequence will depend on the metallicity and on the shock velocity, we can draw two general conclusions from Figure 7:

• Star-forming galaxies with a contribution from fast shocks can mimic a composite starburst-AGN galaxy at all redshifts between z = 0 to z = 3.0.

• Galaxies containing emission from slow shocks can masquerade as galaxies with composite starburst-AGN activity locally, but may mimic high metallicity starburst galaxies at high redshift (z > 1.5).

We recommend that galaxies containing (or likely to contain) a significant contribution from shock emission be removed from samples prior to comparison with Figure 5. Observations of shock-sensitive optical emissionlines such as [S II] λλ6717,31, [O I] λ6300, UV shock sensitive lines (Allen et al. 1998), or measurements of the velocity dispersion from integral field spectroscopy (Rich et al. 2011) can help rule out a contribution from shock excitation.

LINERs

Despite decades of study, the power source of galaxies dominated by emission from Low Ionization Narrow Emission-line Regions (LINERs) is unresolved. The power source of LINER emission depends critically on sample selection and the aperture used to observe LINER emission. Optical and X-ray studies suggest that nuclear LINERs (< 1 kpc) are predominantly inefficiently accreting AGN, analogous to the "low" state of black hole binary systems (Ho 2005;Satyapal et al. 2005;Kewley et al. 2006b;González-Martín et al. 2009;Gu & Cao 2009;Ho 2009). Extended LINER emission (> 1 kpc) may be produced by shock excitation (Heckman 1980;Lípari et al. 2004;Monreal-Ibero et al. 2006;Sharp & Bland-Hawthorn 2010;Rich et al. 2010Rich et al. , 2011)). Aged stellar populations in post-starburst galaxies can also produce line ratios typical of LINERs (Taniguchi et al. 2000;Stasińska et al. 2008;Yan & Blanton 2012;Kehrig et al. 2012).

LINERs lie along the starburst-AGN mixing sequence in the [N II]/Hα versus [O III]/Hβ diagnostic diagram. LINERs can not be distinguished from Seyfert galaxies using the [N II]/Hα versus [O III]/Hβ diagram alone. Therefore, the mixing sequence on Figure 5 may contain both Seyfert and LINER galaxies, depending on the magnitude limit, sample selection, and power source of the LINER emission. Thus, the mixing sequence at intermediate and high redshift should not be interpreted as a pure starburst-AGN mixing sequence, unless LINERs and/or shock excitation can be ruled out.

Spectral Classification at High Redshift

In local galaxies, semi-empirical and theoretical classification lines are used to successfully separate starforming galaxies from those containing an AGN or shock excitation (e.g., Veilleux & Osterbrock 1987;Kewley et al. 2001b;Kauffmann et al. 2003;Kewley et al. 2006b). These classification lines were developed for and tested on local galaxy samples. Figure 5 implies that these standard local classification lines should not be applied at high redshift if the ionizing radiation field is harder, or if the ISM conditions are more extreme (a denser environment or a larger ionization parameter). Thus, comparisons between Figure 5 and large samples of star-forming galaxies at intermediate and high redshift are critical for testing the applicability of local classification methods to higher redshift samples. Until such tests have been performed, local classification methods should not be applied to samples with z ≥ 1.5.

Discussion

Figure 5 can be used with samples of high redshift galaxies to investigate the ISM conditions in star-forming galaxies at high redshift and, potentially, to constrain metallicity gradients in starburst-AGN composite galaxies as a function of redshift. For constraining the ISM conditions in star-forming galaxies, redshifts z > 2.5 is ideal. For distinguishing between steep and flat chemical abundance gradients at high redshift, z = 1.5 is ideal. At higher redshift, if the ISM conditions in star-forming galaxies are more extreme than locally, distinguishing between steep or flat gradients (i.e. columns 3 and 4) is likely to be difficult, especially with low S/N data.

Our predictions rely on the Davé et al. (2011a) cosmological hydrodynamic simulations for the gas-phase metallicity evolution of massive star-forming galaxies. Initial tests suggest that these simulations successfully reproduce the observed metallicity history for galaxies of intermediate stellar masses (10 9 M ⊙ < M * < 10 10 M ⊙ ), but that the simulations predict a larger metallicity for the highest mass galaxies (10 10 M ⊙ < M * ) than observed for 0 < z < 3 (Maiolino et al. 2008;Yuan et al. 2012). These tests are currently limited by small samples. If The best match between a scenario and an ensemble of galaxies may change with redshift. Galaxies at a given redshift may lie within the boundaries of one scenario, while galaxies at another redshift may lie within (or closer to) the boundaries in another scenario. Such a change in location with redshift might indicate a more rapid change in properties than we have assumed. For example, if galaxies at z = 2.5 fall within the boundaries of Scenario 3 (extreme ISM conditions, steep metallicity gradient), but a statistically significant sample of galaxies at z = 0.8 fall within the boundaries of Scenario 1 (i.e. local ISM and AGN conditions), then the transition from extreme ISM to local ISM conditions are likely to have occured in the ∼ 4 Gyr between 0.8 < z < 2.5, rather than over the 11 Gyr spanned between 0 < z < 2.5.

Galaxies lying outside our sequence boundaries would provide important constraints for future models of high redshift galaxy spectra. In general, galaxies lying above the sequence boundaries at a given redshift (i.e. above the abundance sequence and the mixing sequence boundaries in Figure 5), would indicate that there are either more extreme ISM conditions (or a harder radiation field) than we have assumed in our models, or that some AGN that are more metal-poor than we have assumed.

Note that care must be taken to avoid selection effects when interpreting Figure 5 in terms of the evolution of ensembles galaxies. For example, local ultraluminous infrared galaxies systematically occupy the mixing sequence and the high metallicity (i.e. large [N II]/Hα) end of the abundance sequence compared with non-infrared selected galaxies (e.g., Yuan et al. 2010). A similar selection effect may exist at higher redshift for infraredselected samples. Locally, this offset is due to a combination of star-formation, AGN activity, and/or shock excitation (Rich et al. 2012).

Galaxies selected through strong emission-line equivalent widths may also show a bias towards particular locations on the BPT diagram. Strong Hα equivalent widths are produced by low metallicity galaxies (with little continuum from an old stellar population), and galaxies with new intense bursts of star formation. Both of these populations may have a higher ionization parameter than in the general population at the same redshift. These galaxies may lie above the abundance sequence for a volumelimited sample at a given redshift.

Magnitude-limited samples may preferentially select galaxies with the strongest optical or UV continuum, which may miss the faintest (metal-poor) galaxies and the dustiest (potentially metal-rich) galaxies at a given redshift. Samples of local galaxies or our theoretical model spectra can be used to investigate which galaxies may be missing from the BPT diagram due to selection biases.

To summarize, the predictions in Figure 5 may give important insights into both the ISM conditions and metallicity gradients in galaxies at a given redshift, providing that selection effects are well understood in comparison samples.

CONCLUSIONS

We have combined current stellar evolutionary synthesis and photoionization models with chemical evolution estimates from cosmological hydrodynamic simulations to predict how the star-forming abundance sequence and the starburst-AGN mixing sequence may appear at intermediate and high redshift on the BPT diagram. We show that:

• The position of the star-forming abundance sequence at a given redshift depends on the hardness of the ionizing radiation field and on the ISM conditions of the nebulae surrounding the active star-forming regions.

• In star-forming galaxies, a harder ionizing radiation field and/or a larger electron density moves galaxies above the normal star-forming abundance sequence on the BPT diagram. A larger ionization parameter may raise or lower the [N II]/Hα and [O III]/Hβ line ratios, depending on the metallicity of a given galaxy.

• The position of the AGN branch on the BPT diagram depends sensitively on the metallicity of the AGN narrow-line region. At low metallicities (log(O/H) + 12 8.4), the AGN region coincides with the star-forming abundance sequence. Low metallicity AGN galaxies can not be distinguished from low-metallicity star-forming galaxies using the BPT diagram alone.

• Investigations into the position of the galaxy population on the BPT diagram at a given redshift potentially provides a powerful probe of how the ISM conditions in star forming galaxies and how the metallicity gradient in AGN galaxies has changed with time.

We apply our latest fast and slow shock models to investigate how well the BPT diagram can isolate shocked galaxies at intermediate and high redshift. We conclude that:

• Galaxies dominated by fast shock emission masquerade as galaxies dominated by AGN at all redshifts.

• Galaxies dominated by slow shocks typical of galactic outflows mimic starburst-AGN composites at low redshift, and mimic high metallicity starburst galaxies at high redshift (z > 1.5).

• In galaxies where shocks are suspected, we recommend the use of high resolution integral field spectroscopy to isolate shocked regions from regions of star-formation prior to line ratio analysis.

Over the coming decade, near-infrared multi-object spectroscopy and high resolution near-infrared integral field spectroscopy will allow rest-frame optical spectroscopic analysis to be applied to large, statistically significant samples beyond the local universe, for the first time.

These surveys will allow the detailed application of theoretical population synthesis, photoionization, and shock models beyond the local universe, significantly improving our understanding of the ISM conditions in galaxies in the intermediate and high redshift universe. This work is the first in a series aimed at investigating the ISM conditions in active galaxies across cosmic time.

  right panel), we show the polynomial fit from equation 3. This polynomial represents the mean local star-forming abundance sequence with model estimated errors of ±0.1 dex in both the [N II]/Hα and [O III]/Hβ directions. The upper and lower bounds of ±0.1 dex (dashed lines) encompass 91% of the SDSS starforming galaxies. The [N II]/Hα and [O III]/Hβ line ratios along equation 3 are related to metallicity through the following equation 12 + log O/H = 8.97 -0.32x, (4) where x = log [OIII]/Hβ [N II]/Hα and metallicity is defined on the Kewley & Dopita (

Figure 1 .

 1 Figure 1. The [N II]/Hα versus [O III]/Hβ optical diagnostic diagram for the Sloan Digital Sky Survey galaxies analyzed by Kewley et al.(2006b). Left: The colored curves show our new theoretical stellar population synthesis and photoionization model grid for star-forming galaxies based on a κ electron temperature distribution. Right: The red solid curve shows the mean star-forming sequence for local galaxies. The shape of the red solid curve is defined by our theoretical photoionization models, while the position is defined by the best-fit to the SDSS galaxies. The ±0.1 dex curves (dashed lines) represent our model errors and contain 91% of the SDSS star-forming galaxies.

Figure 2 .

 2 Figure 2. An illustration of the effect of varying different galaxy parameters on the star-forming galaxy abundance sequence in the [N II]/Hα versus [O III]/Hβ diagnostic diagram. The original SDSS star-forming galaxy sequence is well-fit by the red theoretical curve. Raising the hardness of the ionizing radiation field (orange dashed line) moves the abundance sequence towards larger [N II]/Hα and [O III]/Hβ ratios. A similar effect is seen when the electron density of the gas is raised (green dot-dashed line). The relationship between ionization parameter, metallicity and the [N II]/Hα and [O III]/Hβ line ratios is more complex. At high metallicities, raising the ionization parameter causes the [N II]/Hα ratio to become smaller, while [O III]/Hβ is largely unaffected. At low metallicities, raising the ionization parameter raises the [O III]/Hβ ratio while [N II]/Hα is largely unaffected.

Figure 3 .

 3 Figure3. The BPT diagram showing the theoretical lower limit (z = 0; blue dot-dashed line) and upper limit star-forming abundance sequence as a function of redshift (z = 0.8, 1.5, 2.5, 3.0; dashed, dot-dashed, triple dot-dashed and solid lines, respectively).

Figure 4 .

 4 Figure 4. The (Groves et al. 2004a) dusty radiation pressure dominated photoionization models for AGN as a function of metallicity. Lines of constant ionization parameter (orange) and constant power-law index (red) are shown. For comparison, the SDSS sample from (Kewley et al. 2006b) is shown.

  Figure 4 to derive the range of possible [N II]/Hα and [O III]/Hβ ratios for the 100% AGN region at each redshift. The AGN model metallicity is simply the local average metallicity (log(O/H) + 12∼ 9.0) for Case 1. For Case 2, the average AGN model metallicity is determined by the redshift through equation 1. The spread of the 100% AGN region on the BPT diagram for both cases is determined by the full range of ionization parameters and power-law indices in our model grid at a given AGN model metallicity.

Figure 5 .

 5 Figure 5. The Cosmic BPT Diagram: Our theoretical predictions for the position of the star-forming abundance sequence (left curve) and the starburst-AGN mixing sequence (right curve). The primary driver for our BPT evolution model is chemical enrichment. The rows give the expected position of the two sequences as a function of redshift for four limiting model scenarios given in each column. Column (1): Normal ISM conditions, and metal-rich AGN NLR at high-z; Column (2): Normal ISM conditions, and metal-poor AGN NLR at high-z; Column (3): Extreme ISM conditions, and metal-rich AGN NLR at high-z; Column (4): Extreme ISM conditions, and metal-poor AGN NLR at high-z.

Figure 6 .

 6 Figure 6. The (Allen et al. 2008) fast shock models (red) and the (Rich et al. 2010) slow shock models (blue) as a function of metallicity. For comparison, the SDSS sample from (Kewley et al. 2006b) is shown. The fast and shock slow model positions indicatewhere galaxies that are 100% dominated by shocks are likely to lie. Galaxies containing a mixture of ionizing sources, such as shocks and star-formation will lie along mixing sequences between the star-forming sequence and the 100% shock models. The position and shape of the mixing sequences will depend on the metallicity of the galaxy and the shock velocity. Examples of typical mixing sequences are given inRich et al. (2010Rich et al. ( , 2011)).

Figure 7 .

 7 Figure7. The Cosmic BPT Diagram: Our theoretical predictions for the position of the star-forming abundance sequence (left curve) and the starburst-AGN mixing sequence (right curve) as a function of redshift (rows) for our four limiting scenarios (columns) described in Section 5. In our model, the primary drivers for BPT evolution are chemical enrichment and a change in the ISM conditions for star-forming galaxies. The(Allen et al. 2008) fast shock models (purple) and the(Rich et al. 2010) slow shock models (black) are shown for comparison.

  ;Madden et al. 2006;Hunt et al.

log(q) 8.30 8.00 7.75 7.50 7.25 7.00 6.75 6.50 log(O/H)+12 9.39 9.17 8.99 8.69 8.39 8.17

  Since the [O II] zone is shorter in a density-bounded nebula, the [O III]/[O II] ratio becomes larger. If radiation-bounded models are applied to such nebulae, then the larger [O III]/[O II] line ratio would be interpreted as a high ionization parameter.

Table 1

 1 Mixing Sequence Boundaries a

		Scenario 1	Scenario 2	Scenario 3		Scenario 4
		SF model: normal	SF model: normal	SF model: extreme	SF model: extreme
		AGN model: Metal-rich NLR AGN model: Metal-poor NLR AGN model: Metal-rich NLR AGN model: Metal-poor NLR
		z = 0		z = 0	z = 0		z = 0
		lower	upper	lower	upper	lower	upper	lower	upper
	a	0.034	0.885	0.029	0.885	0.031	0.917	0.029	0.917
	b	1.447	-0.792	1.340	-0.792	1.441	-0.491	1.332	-0.491
	c	-0.986	-6.712	-0.712	-6.712	-0.879	-6.090	-0.710	-6.090
	d	. . .	. . .	1.472	. . .	. . .	. . .	1.594
	xr [-0.45, 0.29]	[-0.45 , -0.12]	[-0.45 , 0.29]	[-0.45 , -0.12]	[-0.50 , 0.29]	[-0.44 , -0.12]	[ -0.47 , 0.29] [-0.44 , -0.12]
	yr [-0.90 , 0.38]	[ -0.20 , 0.91]	[ -0.90 , 0.38] [ -0.20 , 0.91]	[ -0.90 , 0.38] [ -0.10 , 0.91]	[ -1.05 , 0.38] [ -0.10 , 0.91]
		z = 0.8	z = 0.8	z = 0.8		z = 0.8
		lower	upper	lower	upper	lower	upper	lower	upper
	a	0.034	1.002	0.603	-13.734	0.025	1.032	0.567	-6.277
	b	1.447	0.602	1.422	-38.844	1.429	0.882	1.728	-18.664
	c	-0.986	-2.078	-1.606	-25.672	-0.693	-1.382	-2.702	-11.999
	d	. . .	. . .	5.072	. . .	. . .	. . .	6.479
	xr [ -0.45 , 0.29] [ -0.54 , -0.12] [ -0.45 , 0.30] [ -0.77 , -0.59]	[ -0.37 , 0.29] [ -0.49 , -0.12] [ -0.34 , 0.30] [ -0.77 , -0.56]
	yr [ -0.90 , 0.38] [ 0.05 , 0.91]	[ -0.90 , 1.00] [ 0.13 , 1.00]	[ -0.60 , 0.38] [ 0.26 , 0.91]	[ -0.85 , 1.00] [ 0.40 , 1.00]
		z = 1.5	z = 1.5	z = 1.5		z = 1.5
		lower	upper	lower	upper	lower	upper	lower	upper
	a	0.034	1.027	0.745	-10.585	0.047	1.022	0.696	-6.293
	b	1.447	0.902	1.365	-26.300	1.252	0.939	1.955	-15.695
	c	-0.986	-0.837	-0.233	-14.970	-0.262	0.155	-3.789	-8.433
	d	. . .	. . .	10.3905	. . .	. . .	. . .	18.5499
	xr [ -0.45 , 0.29] [ -0.61 , -0.12] [ -0.45 , 0.17] [ -0.88 , -0.67]	[ -0.35 , 0.29] [ -0.55 , -0.12] [ -0.25 , 0.17] [ -0.88 , -0.71]
	yr [ -0.90 , 0.38] [ 0.16 , 0.91]	[ -0.90 , 1.00] [ 0.25 , 1.00]	[ -0.40 , 0.38] [ 0.55 , 0.91]	[ -0.60 , 1.00] [ 0.60 , 1.00]
		z = 2.5	z = 2.5	z = 2.5		z = 2.5
		lower	upper	lower	upper	lower	upper	lower	upper
	a	0.034	1.033	0.879	-8.604	0.163	0.958	0.834	-64.678
	b	1.447	0.995	1.801	-19.285	0.722	0.415	2.479	-139.855
	c	-0.986	-0.246	4.277	-9.706	0.160	0.212	-3.271	-74.252
	d	. . .	. . .	20.374	. . .	. . .	. . .	58.818
	xr [ -0.45 , 0.29] [ -0.67 , -0.12] [ -0.45 , 0.07] [ -0.99 , -0.75]	[ -0.30 , 0.29] [ -0.70 , -0.12] [ -0.17 , 0.07] [ -1.01 , -0.99]
	yr [ -0.90 , 0.38] [ 0.25 , 0.91]	[ -0.90 , 1.00] [ 0.37 , 1.00]	[ -0.05 , 0.38] [ 0.77 , 0.91]	[ -0.20 , 1.00] [ 0.85 , 1.00]
		z = 3.0	z = 3.0	z = 3.0		z = 3.0
		lower	upper	lower	upper	lower	upper	lower	upper
	a	0.034	1.033	0.942	-8.212	0.259	0.933	0.908	2.037
	b	1.447	1.005	2.268	-17.802	0.370	0.190	2.692	1.186
	c	-0.986	-0.109	7.423	-8.621	0.208	0.103	-0.839	0.171
	d	. . .	. . .	26.057	. . .	. . .	. . .	71.964
	xr [ -0.45 , 0.29] [ -0.68 , -0.12] [ -0.45 , 0.04] [ -1.03 , -0.78]	[ -0.28 , 0.29] [ -0.80 , -0.12] [ -0.15 , 0.04] [ -1.12 , -1.03]
	yr [ -0.90 , 0.38] [ 0.30 , 0.91]	[ -0.90 , 1.00] [ 0.40 , 1.00]	[ 0.16 , 0.38]	[ 0.85 , 0.91]	[ -0.10 , 1.00] [ 0.92 , 1.00]

This research has made use of NASA's Astrophysics Data System Bibliographic Services. L. Kewley wishes to acknowledge the ANU Interstellar Plotters for stimulating discussions. L. Kewley gratefully acknowledges support by an Australian Research Council (ARC) Future Fellowship FT110101052 and NSF Early CAREER Award AST0748559. T.T. Yuan acknowledges a Soroptomist Founder Region Fellowship for Women. Dopita and Kewley would also like to thank the Australian Research Council for support under Discovery project DP0984657 and Discovery Project DP130104879.