Ionization State of Inter-Stellar Medium in Galaxies: Evolution, SFR-M_star-Z Dependence, and Ionizing Photon Escape

We present a systematic study for ionization state of inter-stellar medium in galaxies at z=0-3 with ~140,000 SDSS galaxies and 108 intermediate to high redshift galaxies from the literature, using an ionization-parameter sensitive line ratio of [OIII]5007/[OII]3727 and photoionization models. We confirm that z=2-3 galaxies show an [OIII]/[OII] ratio significantly higher than a typical star-forming galaxy of SDSS by a factor of>~10, and the photoionization models reveal that these high-z galaxies have an ionization parameter of log(qion/cm s^{-1})~7.6-9.0, a factor of ~4-10 higher than local galaxies. For galaxies at any redshift, we identify a correlation between the [OIII]/[OII] ratio and galaxy global properties of star-formation rate (SFR), stellar mass (M_star), and metallicity (Z). We extend the fundamental metallicity relation (FMR) and develop the fundamental ionization relation (FIR), a four-dimensional relation of ionization parameter, SFR, M_star, and Z. The intermediate and high-z galaxies up to z~3 follow the FIR defined with the local galaxies, in contrast with the FMR whose possible evolution from z~2 to 3 is reported. We find that the FMR evolution of z~2-3 appears, if one omits ionization parameter differences, and that the FMR evolution does not exist for an average metallicity solution of z~3 galaxies with a high-ionization parameter. Interestingly, all of two local Lyman-continuum emitting galaxies (LyC leakers) have a high [OIII]/[OII] ratio, indicating a positive correlation between [OIII]/[OII] and ionizing photon escape fraction (fesc), which is successfully explained by our photoionization models. Because [OIII]/[OII] ratios of z~2-3 galaxies, especially Ly-alpha emitters (LAEs), are comparable to, or higher than, those of the local LyC leakers, these high-z galaxies are candidates of Lyman-continuum emitting objects. (abridged)


INTRODUCTION
For galaxy formation studies, it is essential to understand physical conditions of inter-stellar medium (ISM) of galax-⋆ E-mail: kimihiko.nakajima@nao.ac.jp ies. Physical properties of hot ISM are characterized by gas-phase metallicity and ionization parameter. A gas-phase metallicity is a record of galaxy's star formation history and gas infall+outflow, and is commonly defined by oxygen abundance, 12 + log(O/H). Metallicity estimates can be made with metal lines divided by a hydrogen recombi-nation line, such as ([O iii]λλ5007, 4959+[O ii]λ3727)/Hβ (R23-index; e.g., Pagel et al. 1979;Kewley & Dopita 2002) and [N ii]λ6584/Hα ratio (N 2-index; e.g., Storchi-Bergmann et al. 1994;Kewley & Dopita 2002;Pettini & Pagel 2004). An ionization parameter, qion, is sensitive to the degree of excitation of an H ii-region. It is defined as the ratio of mean ionizing photon flux to mean atom density, and its dimensionless form, U = qion/c, is also used. Ionization parameters are estimated with emission lines of different ionization stages of the same element, such as [O iii]λ5007/[O ii]λ3727 ratio (e.g., Kewley & Dopita 2002).
Metallicities of galaxies have been well studied in the local and high-z universe in the diagram of stellar-mass vs. metallicity (M⋆-Z) relation. In the local universe, a tight M⋆-Z relation is recognized, where a metallicity increases with stellar mass (e.g., Tremonti et al. 2004). At high-z, similar M⋆-Z relations are reported at z = 0.5 − 3.5 (e.g. Savaglio et al. 2005;Erb et al. 2006;Maiolino et al. 2008;Mannucci et al. 2009;Zahid et al. 2011;Yabe et al. 2012). These studies find the evolution of M⋆-Z relation; a metallicity decreases towards high-z at a given M⋆. Recently, Mannucci et al. (2010) and Lara-López et al. (2010) identify that the M⋆-Z relation depends on star formation rate (SFR), which is initially suggested by Ellison et al. (2008), and propose the fundamental metallicity relation (FMR) that forms a single plane in the stellarmass, metallicity, and SFR space. The FMR holds over z ≃ 0 − 2.5 with no evolution. This would indicate that there are physical processes commonly found in local and high-z galaxies. Many follow-up studies test the universality of the FMR over the wide range of stellar mass, SFR, and redshifts (e.g, Mannucci et al. 2011;Richard et al. 2011;Nakajima et al. 2012;Niino 2012;Yabe et al. 2012;Cresci et al. 2012;Wuyts et al. 2012;Christensen et al. 2012a;Belli et al. 2013;Lara-López et al. 2013). Although galaxies at z 2.5 follow the FMR, galaxies at z 3 appear to fall below the FMR (Mannucci et al. 2010). Mannucci et al. (2010) claim that this disagreement suggests that there are physical mechanisms different from those of low-z galaxies at z ∼ 3. However, this is still an open question.
In contrast to metallicities, ionization parameters of galaxies are poorly understood. Recently, Nakajima et al. (2013) use a diagram of [O iii]/[O ii] ratio and R23-index to study ionization parameters and metallicities for galaxies at z = 0−3, focusing on the properties of Lyα emitters (LAEs) at z ∼ 2. Nakajima et al. (2013) report that a typical ionization parameter increases towards high-z, and suggest the evolution of ionization parameter (see also, Lilly et al. 2003;Hainline et al. 2009;Richard et al. 2011). However, until recently, statistical studies for ionization parameters do not cover (i) physical relations between ionization parameter and galaxy global properties, such as stellar mass and SFR, and (ii) evolution of the physical relations with ionization parameter. Dopita et al. (2006a) present that in the local universe an ionization parameter seems to gradually decrease with M⋆, but do not address the issue about evolution of this trend (see also Brinchmann et al. 2008).
There is another diagnostics of ionization parameter, using a plot of [N ii]/Hα and [O iii]/Hβ ratios that is referred to as the BPT diagram (Baldwin et al. 1981). Lo-cal star-forming galaxies form a tight sequence on the BPT diagram (e.g., Kauffmann et al. 2003a). In contrast, high-z galaxies depart from the local star-forming sequence towards higher [O iii]/Hβ ratio (e.g., Shapley et al. 2005;Erb et al. 2006;Liu et al. 2008;Hainline et al. 2009;Finkelstein et al. 2009;Yabe et al. 2012). Similarly, recent spectroscopic surveys report that high-z galaxies on average have [O iii]/Hβ ratios much higher than local similar mass galaxies (Cullen et al. 2014;Holden et al. 2014). One interpretation of that departure is evolution of ionization parameter (e.g., Brinchmann et al. 2008;Kewley et al. 2013a,b). Kewley et al. (2013a) have compared high-z galaxy line ratios and theoretical predictions to investigate evolution of the star-forming sequence on the BPT diagram, and claimed that high-z galaxies would have a larger ionization parameter, a higher electron density, and/or a harder ionizing radiation field (see also Brinchmann et al. 2008;Kewley et al. 2013b). Shirazi et al. (2013) have also argued a possibility that the typical gas density of H ii-regions would be much higher in high-z galaxies.
Ionization parameters are also important for accurate estimates of metallicity from nebular emission lines. Metallicity estimation methods based on local empirical relations (e.g., Maiolino et al. 2008) are often applied to high-z galaxies. These methods implicitly assume ionization parameters of typical local galaxies. Since various studies suggest that ionization parameter evolves with redshift, the estimated metallicities from local empirical relations would be biased.
Recently, a relationship between ionization parameter and escape fraction of ionizing photons (fesc) and its importance for cosmic reionization are discussed. Nakajima et al. (2013) have pointed out a possibility that fesc positively correlates with ionization parameter, since an optically-thin nebula could show a high [O iii]/[O ii] ratio (e.g., Giammanco et al. 2005;Brinchmann et al. 2008;Kewley et al. 2013b). If this correlation of ionization parameter and fesc is established, ionization parameters allow us to investigate ionizing photon escapes that cannot be directly observed for galaxies at the epoch of reionization (z > 6; Fan et al. 2006;Dunkley et al. 2009). Studies of galaxies at z ∼ 6 − 7 have identified a tension between the ionizing photon production rate and cosmic reionization (e.g., Ouchi et al. 2009;Robertson et al. 2010). A moderately high average fesc of 0.2 would be required for galaxies in the early universe to keep the inter-galactic medium (IGM) ionized. Nakajima et al. (2013) have suggested that LAEs exhibit higher ionization parameters than other highz galaxy populations such as Lyman-break galaxies (LBGs). Because observations reveal that a fraction of LAEs to star-forming galaxies of LBGs increases with redshift (e.g., Ouchi et al. 2008;Vanzella et al. 2009;Stark et al. 2010Stark et al. , 2011, more LAEs with a high ionization parameter, i.e., a high fesc, would ease the tension of the deficit of ionizing photons at the early universe. However, the claim of Nakajima et al. (2013) is made only with two LAEs, and it is unclear whether this is a general feature of LAEs. One needs to (i) investigate the trend that LAEs typically show an ionization parameter higher than LBGs, and (ii) understand the physical reason of the correlation between ionization parameter and fesc with theoretical models. Although predicting fesc is challenging by using simulations or theoretical models due to the large uncertainties in the assumptions about the ISM and the limited resolutions (e.g., Rahmati et al. 2013), we investigate the correlation with photoionization models assuming a simple, homogeneous gas cloud as the first step.
This paper is organized as follows. We describe our samples of galaxies in Section 2, and test the main sample of SDSS in Section 3. In Section 4, we show all of our galaxies in the [O iii]/[O ii] vs. R23-index diagram, and investigate galaxy properties in the diagram with the aid of photoionization models. We then clarify the dependencies of ionization parameter on galaxy global properties of stellar mass, SFR, and metallicities. We introduce a new relation named fundamental ionization relation (FIR) that is described in the four-dimensional parameter space of [O iii]/[O ii] ratio, R23-index, stellar mass, and SFR. We discuss these observational results in Section 5, addressing the issues of (i) the evolution of ionization parameter, (ii) the FMR evolution at z 3, (iii) the relations between stellar mass,and SFR,and (iv) the correlation between ionization parameter and fesc. Finally, we summarize our results in Section 6. Throughout this paper, we estimate stellar masses and SFRs based on a Chabrier (2003) initial mass function (IMF), and assume an Asplund et al. (2009) solar metallicity, where log(Z/Z⊙) = 12 + log(O/H) −8.69. We adopt a concordance ΛCDM cosmology with (Ωm, ΩΛ, H0) = (0.3, 0.7, 70 km s −1 Mpc −1 ).

SDSS galaxies
We make an SDSS sample from the SDSS Data Release 7 MPA/JHU catalog 1 including emission line properties and stellar masses (Kauffmann et al. 2003b;Brinchmann et al. 2004;Salim et al. 2007). We select galaxies with the criteria similar to those adopted by Mannucci et al. (2010); (1) 0.07 < z < 0.30, (2) a signal-to-noise ratio (S/N) of Hα > 25, (3) AV < 2.5, (4) Hα/Hβ > 2.5, (5) active galactic nuclei (AGN) are removed based on the BPT diagram (Baldwin et al. 1981) with the classification of Kauffmann et al. (2003a) 2 , (6) metallicities estimated from the N 2-and R23-index (Maiolino et al. 2008) agree within 0.25 dex. We adopt the average value of the two metallicities 3 , and use total stellar masses listed in the updated MPA/JHU stellar mass catalog 4 . All the emission line fluxes are corrected for extinction using Hα/Hβ and the extinction 1 http://www.mpa-garching.mpg.de/SDSS/DR7/ 2 Although Brinchmann et al. (2004) recommend that requiring S/N > 3 in all four BPT diagram lines is useful, we do not adopt the criteria because the requirement could bias the sample regarding metallicity. Indeed, the criteria of S/N ratio of Hα > 25 is strong enough to remove lots of galaxies with low-S/N emission lines in the sample. By imposing S/N ratios of all four BPT diagram lines > 3, the sample size decreases to ∼ 120, 000. Nevertheless, we have checked that our main results do not depend significantly on the criteria. In order to fairly compare with the FMR of Mannucci et al. (2010), we do not impose the requirement of S/N ratios for the BPT diagram. 3 These metallicities are only used for the consistency check of the FMR of Mannucci et al. (2010) (Section 3). 4 http://home.strw.leidenuniv.nl/~jarle/SDSS/ curve given by Cardelli et al. (1989). SFRs are estimated from the Hα luminosities by using the Kennicutt (1998)'s relation, and (total) stellar masses from the broadband photometry. All of the SFRs and stellar masses are corrected to the Chabrier (2003) IMF. We remove the same objects in the catalog that are found by the MPA/JHU group. We thus obtain 136, 067 SDSS galaxies.
In this sample, we highlight two populations of galaxies shown below.
• Green Pea galaxies (GPs; Cardamone et al. 2009) that are local star-forming galaxies with a very strong [O iii] emission line whose EW is ∼ 1000Å. They are found in the SDSS spectroscopic sample of galaxies with green colors in the SDSS gri composite images, where the strong [O iii] emission fall in the r band. Because the strong [O iii] emission has to be redshifted to the bandpass of r, the redshift of GPs ranges from 0.112 to 0.360. We use 80 GPs given by Cardamone et al. (2009) that are free from AGN contamination. We cross-match the 80 GPs to our SDSS galaxies, and identify 51 GPs. To include GPs at z > 0.30, we perform the same cross-matching to SDSS galaxies that do not meet the redshift criterion of (1), and obtain additional 15 GPs at z > 0.3. We use a total of 66 = (51 + 15) GPs in our analysis.
• Lyman Break Analogs (LBAs; Heckman et al. 2005;Hoopes et al. 2007;Basu-Zych et al. 2007;Overzier et al. 2008Overzier et al. , 2009) that are super-compact UV-luminous local galaxies whose characteristics are similar to those of highz LBGs. The UV imaging survey performed by the Galaxy Evolution Explorer (GALEX) is used to select star-forming galaxies at 0.1 < z < 0.3 with a high luminosity of LFUV > 10 10.3 L⊙ and a high far UV surface brightness of IFUV > 10 9 L⊙ kpc −2 . We compile catalogs of Hoopes et al. (2007), Basu-Zych et al. (2007), and Overzier et al. (2009), and obtain 51 LBAs. Cross-matching them to our SDSS galaxies and removing AGN candidates, we make an LBA sample composed of 37 LBAs.
We remove sources in our GP and LBA samples from our SDSS galaxies, and define a sample of SDSS galaxies with no GPs and LBAs as SDSS sample.

LyC leakers
There are two local star-forming galaxies emitting Lymancontinuum (LyC) radiation, Haro 11 (Bergvall et al. 2006;Leitet et al. 2011) and Tol 1247-232 (Leitet et al. 2013). They are initially selected as blue compact galaxies or H ii galaxies, and found to leak LyC radiation from the Far Ultraviolet Spectroscopic Explorer (FUSE) data. In order to study the relation between ionization state and ionizing photon escape, we investigate the two LyC leaking galaxies. Hereafter we refer to these two galaxies as LyC leakers. Emission line fluxes for Haro 11 and Tol 1247-232 are obtained from Bergvall &Östlin (2002) and Terlevich et al. (1993), respectively, and stellar masses from Leitet et al. (2013). Their SFRs are derived in the same manner as those of our SDSS galaxies (Section 2.1.1). Neither of them is contaminated by AGN (Leitet et al. 2013). We note that Haro 11 can be classified as an LBA (Heckman et al. 2011). The mass-metallicity relations of galaxies from our SDSS sample. Each line denotes a SFR subsample from our SDSS sample whose SFR is presented at the bottom of the panel. Center: The SFR-metallicity relation given by our SDSS sample. The lines represent the relations for subsamples with M⋆ values shown in the legend. Right: Same as the left panel, but for the µ 0.32 -metallicity relation, a two-dimensional view of the FMR. The color code is the same as the one in the left panel. The black solid curve is the quadratic best-fit of our SDSS sample, while the black dotted curve is the original FMR obtained by Mannucci et al. (2010) and Mannucci et al. (2011). The black solid and dotted lines are almost identical. Thus, our SDSS sample reproduces the original FMR.

Intermediate and High Redshift galaxies
We make samples of intermediate-and high-z galaxies that have firm estimates of stellar mass and spectroscopic measurements of nebular emission lines (hydrogen Balmer lines, . We particularly include strongly lensed galaxies, BX/BM galaxies, LBGs, and LAEs. • Intermediate redshift galaxy sample (z = 0.5 − 1.5) is composed of 24 galaxies from the Gemini Deep Deep Survey (GDDS) and 15 galaxies from the Canada-France Redshift Survey (CFRS) (Savaglio et al. 2005;Maier et al. 2005). Moreover, we add 26 low-mass galaxies recently reported by Henry et al. (2013), and 4 strongly lensed galaxies (Christensen et al. 2012a;Queyrel et al. 2009) to the sample. In total, the intermediate-z sample contains 69 galaxies.
We use a total of 108 = (69 + 39) intermediate-and high-z galaxies in our analysis.
We divide the high-z galaxy sample into two subsamples, LAE and LBG samples, to investigate the relationship between ionization state and Lyα photon escape. For high-z galaxies with a Lyα flux measurement, we classify galaxies with EWrest (Lyα) > 20Å as LAE; Lynx arc (Fosbury et al. 2003), BX418 (Erb et al. 2010), Abell 1689 31.1 arc, and SMACS J2031 arc (Christensen et al. 2012a,b). We include the two narrow-band selected LAEs, CDFS-3865 and COSMOS-30679 (Nakajima et al. 2013), to the LAE sample. For high-z galaxies with EWrest (Lyα) < 20Å or no Lyα flux measurement, we regard these galaxies as LBGs for simplification, given the observational fact that continuum-selected galaxies on average have an EW of Lyα ∼ 0Å (e.g., Shapley et al. 2003;Kornei et al. 2010). Indeed, most of the z ∼ 3 galaxies in the sample of Maiolino et al. (2008) are originally spectroscopically-confirmed based on rest-frame UV absorption lines, not Lyα emission line (see also Vanzella et al. 2006).
We derive physical quantities for these intermediateand high-z galaxies with the same procedures as those for our SDSS galaxies. Emission line fluxes are corrected for extinction with Hα/Hβ and the extinction curve given by Calzetti et al. (2000). For galaxies with no Hα or Hβ measurement, we correct for extinction estimated from spectral energy distribution (SED) fitting, assuming that the amount of extinction is the same in continuum and emission lines. For galaxies neither with Hα/Hβ nor dust extinction from SED fitting 5 , we adopt dust extinction inferred from stellar mass based on the empirical relation (Gilbank et al. 2010 (Richard et al. 2011). We have confirmed that this assumption does not affect our results significantly. We estimate intrinsic SFRs and stellar masses, using lens magnification factors for the strongly lensed galaxies. We place the lower (upper) limits on the SFR and stellar mass of SGAS1527 (SGAS1226) that only has the upper (lower) limit of the magnification factor (Wuyts et al. 2012). We adopt the lower limit of stellar mass given by Villar-Martín et al. (2004) for Lynx arc (Fosbury et al. 2003). For high-z galaxies with no [O ii] detections (Fosbury et al. 2003;Erb et al. 2010;Richard et al. 2011 We note that all the properties for the local, intermediate-z, and high-z galaxies used in this paper are and R23-index. The gray dots are the SDSS galaxies. The red and blue diamonds represent z = 2 − 3 LAEs and LBGs, respectively. Specifically, the diamonds with a circle denote galaxies at z > 3. The gray triangles are z ∼ 1 star-forming galaxies. The green and purple circles show GPs and LBAs, respectively, and the orange pentagons are LyC leakers in the local universe. The black solid curves present photoionization model tracks in the range of 12 + log(O/H)= 7.5 − 9.5 for q ion = 3 × 10 8 , 8 × 10 7 , and 2 × 10 7 cm s −1 (Kewley & Dopita 2002). The numbers with the small open circles on the model tracks denote metallicities in 12 + log(O/H). The dotted lines connect the photoionization model curves with the same metallicity value. derived in a consistent manner, and that they can be fairly compared to each other.

TESTING FMR WITH OUR SDSS SAMPLE
Before we present our analysis results, we test whether our SDSS sample reproduces the FMR (Mannucci et al. 2010). For this purpose, we derive metallicities based on the Maiolino et al.'s (2008) empirical relations, as adopted by Mannucci et al. (2010) (Section 2.1.1). It should be noted that our main discussion is based on metallicities derived with ionization parameter from the photoionization models of Kewley & Dopita (2002) (e.g., Section 4.3). In Figure  1, three panels from left to right show the M⋆-metallicity, SFR-metallicity, and µ0.32-metallicity plots with our SDSS sample. The variable µα is defined as a combination of M⋆ and SFR, where α is a free parameter. Mannucci et al. (2010) obtain α = 0.32 that minimizes the scatter of their SDSS galaxies in the µα-metallicity plane. The colored lines in each panel represent a mean metallicity of our SDSS galaxies for subsamples defined with SFR (left and right panels) and M⋆ (center panel) in a 0.15 dex range. We display subsamples containing > 50 galaxies.
In the right panel of Figure 1, the solid black curve indicates the quadratic best-fit of our SDSS sample, while the dotted curve presents the original FMR (Mannucci et al. 2010(Mannucci et al. , 2011. The original FMR is expressed by a combination of a linear and a 4th-order polynomial functions (Mannucci et al. 2011). From this plot, the FMR of our SDSS sample is almost identical to that of Mannucci et al. (2010Mannucci et al. ( , 2011. Therefore, we conclude that our SDSS sample successfully reproduces the original FMR of Mannucci et al. (2010) vs. R23-index diagram is useful to investigate ionization parameter and metallicity of galaxies (see also Kewley & Dopita 2002). Figure 2 is the diagram The underlying gray dots are the SDSS galaxies. Figure  2 indicates that galaxies with a higher R23-index have a higher [O iii]/[O ii] ratio. The red and blue diamonds denote z ∼ 2 − 3 LAEs and LBGs, respectively. Compared to the SDSS galaxies, these high-z galaxies clearly present an [O iii]/[O ii] ratio higher than most of galaxies by a factor of 10. The majority of z ∼ 1 galaxies (the gray triangles) are placed at the middle between the SDSS galaxies and the z ∼ 2 − 3 LBGs. LBGs at z 3 (blue diamonds with a circle) have an [O iii]/[O ii] ratio higher than those at z 3 on average. In Figure 2, LAEs have the This trend is first claimed by Nakajima et al. (2013), but we confirm this trend with our LAE and LBG samples significantly larger than those used in Nakajima et al. (2013). The two local galaxies of GPs (green open circles) and LBAs (purple circles) are placed near the high-z galaxies. These galaxies present the highest ratio, a) spectral shape of ionizing source, b) gas temperature, c) metallicity, d) ionization parameter, and e) gas density. First, if the incident ionizing source includes an AGN whose spectrum is hard, high-ionization lines are stronger than low-ionization lines. However, since our samples are free from a significant AGN contamination (Section 2.1.1), the ionizing source difference does not significantly change the [O iii]/[O ii] ratio. Second, gas temperature changes the emissivity of [O iii] and [O ii] lines. Gas temperature is determined by the metal abundance, because coolings are dominated by collisionally excited metal ions such as O + , O 2+ , N + etc. We simply regard metallicity as a control parameter of temperature, and investigate the temperature dependence on the basis of gas metallicity. Thus, there are three key parameters, metallicity, ionization parameter, and gas density, that affect the In Figure 2, we plot the photoionization models given by Kewley & Dopita (2002) with the black curves. The models predict line ratios at a given metallicity and ionization parameter with a fixed gas pressure 8 of P/k = 10 5 cm −3 K assuming isobaric conditions. The models in Figure  ( 8 For electron temperatures of 10 4 K. This pressure corresponds to a density of order 10 cm −3 . (3) R23-index increases in a low-metallicity regime (low-Z branch) from 12 + log(O/H) = 7.5 to ∼ 8.5, and decreases in a high-metallicity regime (high-Z branch) from 12 + log(O/H) ∼ 8.5 towards a higher metallicity.
The trend of (1) is explained by the boost of high-energy ionizing photons as the ionization parameter increases. Such high-energy photons produce more abundant [O iii]. The trend of (2) is caused by a decrease of gas temperature in an H ii-region. As gas metallicity increases, gas temperature drops due to cooling of metal lines. For gas clouds with high metallicity, the far-infrared fine-structure cooling of O 2+ is more efficient than those in optical wavelength (e.g., Charlot & Longhetti 2001;Dopita et al. 2000 ratio depends strongly on ionization parameter and moderately on metallicity. R23-index is mainly governed by metallicity, but also depends, albeit weakly, on ionization parameter. In addition to metallicity and ionization parameter, emission line ratios would be changed by gas density. In fact, the recent study of Shirazi et al. (2013) argues that a high gas density would be an origin of the large [O iii]/[O ii] ratio of high-z galaxies. This argument is based on the following ideas. The ionization parameter is defined by where Q H 0 is the number of hydrogen ionizing photons per second, Rs is the Strömgren radius, and nH is the total hydrogen density. The hydrogen ionizing photon production rate of Q H 0 is balanced by a recombination rate under the assumption of ionization equilibrium, where αB is a coefficient of the total hydrogen recombination to the n > 1 levels, and ǫ is the volume filling factor of the ionized gas. The combination of Equations (2) and (3) yields (see also Charlot & Longhetti 2001;Shirazi et al. 2013). By the comparison with low and high redshift galaxies with similar M⋆ and specific SFR (and thus with similar metallicity under the assumption of the FMR being in place), Shirazi et al. (2013) have found that the high-z galaxies show an [O iii]/[O ii] ratio higher than low-z counterparts by ∼ 0.5 dex, indicating that high-z galaxies have a higher qion than the low-z counterparts by ∼ 0.3 dex. Based on Equation (4), Shirazi et al. (2013) have argued that nH increases by a factor of ∼ 8 to high-z, assuming that high-z galaxies and their low-z counterparts have similar Q H 0 and ǫ values. Indeed, Shirazi et al. (2013) have confirmed that high-z galaxies tend to have a high electron density, which is measured directly from emission lines. Because ionization parameter is determined by Q H 0 , nH, and ǫ (Equation 4), ionization parameter is a fundamental parameter of photoionization models, which allows the potential evolutions of nH, Q H 0 , and ǫ. The assumption of no evolution of IMF would also be an open question (e.g., Davé 2008;Narayanan & Davé 2013 where x = log(R23) and the subscript low (high) corresponds to a metallicity value in the low-Z (high-Z) branch of metallicity.  (5)-(7). We repeat the process 500 times, and de-  (Table 1). fine one sigma errors of metallicity and ionization parameter by the 68 percentile of the distribution.
The estimated values are summarized in Table 1 and shown in Figure 3. We list the ranges of values in both the high-Z and low-Z branches for the intermediate and high-z galaxy samples but for the local galaxy samples. For the SDSS sample, only the values in the high-Z branch are given, since the average metallicity for the low-Z branch, 12 + log(O/H) = 7.7 − 8.1, is significantly lower than that estimated with R23-index and N 2-index ( Figure 1). We have checked that the SDSS high-Z branch metallicity of 12 + log(O/H) = 8.92 − 9.12 (Table 1) is consistent with our estimate in Section 3, as well as those estimated in previous studies (e.g., Tremonti et al. 2004;Mannucci et al. 2010). For GPs, LBAs, and LyC leakers, we calculate the metallicities and ionization parameters for individual objects, whose two-branch solutions are distinguished by other metallicity indicators of gas temperature and/or N 2-index (see also Section 5.2). Their 68 percentile of the distributions are listed in Table 1 and shown in Figure 3.
Although only weak constraints are placed on metallicity due to the degeneracy of the two branch solutions, the difference of ionization parameters between galaxy samples is clearly found in Figure 3. The SDSS galaxies have the average value of log(qion/cm s −1 ) ∼ 7.3 (see also e.g., Dopita et al. 2006a). In contrast, high-z galaxies have an ionization parameter of log(qion/cm s −1 ) ∼ 7.6 − 9.0. High-z LBGs have an ionization parameter higher than the SDSS galaxies by a factor of ∼ 4, and high-z LAEs show the highest ionization parameter among the galaxy samples, which is about an order of magnitude higher than the SDSS galaxies. In the local galaxies, GPs and the LyC leakers exhibit ionization parameters much higher than the SDSS galaxies. Similarly, the ionization parameter of LBAs is high. These extreme local populations of GPs, LyC leakers, and LBAs have an ionization parameter as high as LBGs and LAEs at high redshifts. Moreover, the ionization parameter of GPs is comparable to that of LAEs. In this sense, GPs could be local counterparts of high-z LAEs. Figure 2 suggests that there is a variation of ionization parameter in the SDSS sample. In fact, the SDSS galaxy distribution sequence departs from the model curve of ∼ 2 × 10 7 cm s −1 . In the large R23-index regime, we recognize the slope of the SDSS galaxy sequence steeper than that of the photoionization model predictions. This trend indicates that galaxies with a low metallicity of 12 + log(O/H) ∼ 8.5 have a high ionization parameter (e.g., Dopita et al. 2006a;Nagao et al. 2006). If the M⋆-Z relation is in place (e.g., Tremonti et al. 2004), this would imply that low-M⋆ galaxies have a high ionization parameter. We examine the possible dependencies of M⋆ and other galaxy global properties on ionization parameters in Sections 4.4 and 4.5. Similarly, the important characteristics can be found in the small R23-index regime. The sequence of SDSS galaxies presents an upturn from high to low R23-index values. This trend is probably caused by the presence of a hidden AGN activity in a galaxy. In this regime, the SDSS galaxies departing from the model curve would be also dominated by objects with a low SFR and a high M⋆ (Section 4.5; see also Dopita et al. 2006a).

Dependence of [O iii]/[O ii] on M⋆ and SFR
In Section 4.3, we have estimated the typical ionization parameters for our galaxy samples from the comparisons of the photoionization models with the measurements of [O iii]/[O ii] and R23-index, and compared them ( Figure 3). However, these comparisons do not focus on possible dependencies on galaxy global properties such as M⋆ and SFR. Figure 2 presents     Before moving to the next section, we comment here on the variations of global properties for the galaxies shown in Figure 4. First, GPs are the least massive and the most actively star-forming galaxies in the local universe. Their low µ0.32 values are probably due to their low metallicities (e.g., Amorín et al. 2010). LBAs are the next extreme population. The LyC leakers have SFRs as high as GPs and LBAs. Highz galaxies have a M⋆ comparable with the SDSS galaxies, but a SFR higher than these local galaxies. This tendency is consistent with the evolution of the star-formation main sequence (e.g., Daddi et al. 2007). The SFRs and M⋆ of highz galaxies are little higher than those of GPs and LBAs, while the sSFR of high-z galaxies are comparable with those of GPs and LBAs. LAEs exhibit remarkable properties on average, having the least M⋆ and the highest SFR among galaxies at any redshifts.

Fundamental Ionization Relation
As shown in Figure 2, Figure 2 would be explained by the differences of the galaxy global properties.
In order to investigate the dependencies of M⋆, SFR, sSFR, and µ0.32 on the [O iii]/[O ii] vs. R23-index plane, we make subsamples of our SDSS sample. We choose one quantity from the four galaxy global properties. The chosen quantity is referred to as a '3rd quantity' hereafter. A bin size of a subsample is 0.15 dex in the 3rd quantity. For a given subsample, we calculate a median [O iii]/[O ii] ratio of galaxies in a range of ∆logR23 = 0.1. We use subsamples with > 50 galaxies for our analysis. The left panels of Figure 5 show the subsamples of M⋆, SFR, sSFR, and µ0.32 on the [O iii]/[O ii] vs. R23-index plane. Figure 5

indicates that the relation between [O iii]/[O ii] and R23-index depends on the galaxy global properties. A high [O iii]/[O ii
] ratio can be found in less massive, more efficiently star-forming, and probably more metal-poor galaxies (Dopita et al. 2006a). We introduce a parameter, ζ β (Θ), which combines R23-index and one of the 3rd quantities: where β is a free parameter and Θ = log(M⋆/10 10 M⊙), log(SFR/M⊙ yr −1 ), log(sSFR/Gyr −1 ), and (µ0.32 −10.0) for a 3rd quantity of M⋆, SFR, sSFR, and µ0.32, respectively. We where a0-a2 are free parameters. The best-fit parameters are given in Table 2. We show the best-fit quadratures in the right panels of Figure 5. These plots suggest that the  3rd quantity can reduce the systematic dispersions. However, it is unclear whether the new parameter of β reduces the dispersion by a physical reason or a simple statistical reason. Here we use µα defined in Equation (1) as a 3rd quantity and search for an optimal combination of M⋆ and SFR by varying parameters of both α and β 9 . Note that the parameter α denotes the best combination of M⋆ and SFR for metallicity (Mannucci et al. 2010), while the parameter β defines the best combination of these quantities for both metallicity and ionization parameter. Based on Equation (8) with Θ = (µα −10.0), we calculate χ 2 of the free parameters, α and β. Here we assume that 1σ [O iii]/[O ii] dispersions of our subsamples are originated from measurement errors. These 1σ dispersions are used to obtain the χ 2 values. Figure 6 shows χ 2 values of α and β. The minimum χ 2 is provided by the parameters of (α, β) = (0.29, 0.25). The best-fit parameters are α = 0.29 +0.04 −0.26 and β = 0.25 +0.06 −0.12 . Interestingly, the best-fit value of α = 0.29 +0.04 −0.26 is consistent with the one of the FMR (α = 0.32). Thus, µ0.32 (µα 9 Although µα does not represent a function of SFR with no dependence of M⋆, we have found in our analysis that a ζ function with Θ = Θ(SFR) provides the best-fit value of β ≃ 0 (Table 2), which indicates that a SFR alone does not reduce the scatter of the SDSS galaxy distribution.  Figure 5. In this plot, we compare the relation of the SDSS galaxies with all of our galaxy samples, high-z galaxies, GPs, LBA, and the LyC leakers.
for α = 0.32) is the best 3rd quantity to describe the relation. Moreover, we have newly constrained the β parameter that is tightly related to ionization parameter. The best-fit β rules out β = 0 at the > 95 % confidence level. The results of our analysis suggest that ionization parameters given with the [O iii]/[O ii] vs. R23-index diagram is important to characterize galaxies, and that ionization parameters are closely related to the galaxy global properties of M⋆ and SFR as well as metallicity.
We compare our various galaxy samples with the bestfit function of the SDSS galaxies in Figure 7. This figure is the same as the bottom right panel of Figure 5, but with high-z galaxies and local extreme populations, GPs, LBAs, and the LyC leakers. Although the best-fit function is derived with the SDSS sample alone, high-z galaxies and local extreme populations follow the same relation or fall on the extrapolation of the function. It should be noted that no significant differences are found between galaxies at z ∼ 2 and 3.

The relation of best-fit function between [O iii]/[O ii]
and ζ β (Θ) (Figure 7) is thus thought to be a fundamental relation consisting of ionization parameter, metallicity, M⋆, and SFR over cosmic time. Hereafter we refer to the relation as fundamental ionization relation (FIR).

Evolution of Ionization State
We have confirmed the trend that high-z galaxies have ionization parameters significantly higher than local galaxies. The SDSS galaxies typically show qion ∼ 2 × 10 7 cm s −1 , while the high-z LBGs exhibit an ionization parameter higher than the SDSS galaxies by a factor of ∼ 4 (Table 1, Figure 3). This evolution of ionization parameter is consistent with the trend that high-z galaxies depart from the local star-forming sequence in the BPT diagram towards higher [O iii]/Hβ ratio (e.g., Brinchmann et al. 2008;Kewley et al. 2013a,b), and with the recent spectroscopic surveys of high-z galaxies (Cullen et al. 2014;Holden et al. 2014). The evolution of ionization parameter can be closely related to the evolution of SFR and M⋆. As shown in Figure 5, a high [O iii]/[O ii] ratio, i.e., a high ionization parameter, is found in low M⋆, low µ0.32, and high sSFR galaxies. Since highz galaxies tend to have a higher sSFR (or smaller µ0.32) than local galaxies (Figure 4; cf. the evolution of the starformation main sequence; Daddi et al. 2007), they are supposed to show an ionization parameter typically higher than local galaxies. Note that ionization parameter is, by definition, the ratio of ionizing photon and hydrogen atom densities. Since the amount of ionizing photons and hydrogen atoms would be positively correlated with SFR and stellar mass, respectively, the evolution of ionization parameter would be explained by the increase of sSFR, i.e., the high SFR and low mass for high-z galaxies.
The local extreme populations of GPs, LyC leakers, and LBAs have an ionization parameter as high as z ∼ 2 − 3 LBGs and LAEs. Moreover, the ionization parameter of GPs is comparable to that of LAEs. GP's low metallicities (see also Amorín et al. 2010) and SFR (Izotov et al. 2011) are also analogous to those observed in LAEs. In these senses,  (Maiolino et al. 2008). The black solid curves denote photoionization models (Kewley & Dopita 2002). The labels on the curves present values of log(q ion ).
GPs could be local counterparts of z ∼ 2−3 LAEs. However, GPs are quite rare objects in the local universe, occupying only ∼ 0.06 % of the SDSS sample (Nakajima et al. 2013;Cardamone et al. 2009). The number density of GPs is approximately two orders of magnitude smaller than that of z ∼ 2 LAEs. Interestingly, Heckman et al. (2005) report a similar difference of abundances between LBAs and z ∼ 3 LBGs. These results indicate that galaxies with a high ionization parameter would emerge in a high redshift universe, and that these galaxies are more dominant in high-z than the local universe.
Recently, some studies have found galaxies with very strong nebular emission lines that are called extreme emission line galaxies (EELGs; e.g., van der Wel et al. 2011;Atek et al. 2011 ( 3 − 100). They argue that USELs typically have ionization parameters as high as qion 10 8 cm s −1 . Obviously, ionization parameter is a key quantity to characterize these extreme populations of galaxies.

Revisiting the Issue of FMR Evolution from
z ∼ 2 to 3 Mannucci et al. (2010) have suggested that the FMR can describe properties of galaxies with no redshift evolution up to z ∼ 2.5, and claimed that galaxies at z 3 fall significantly below the FMR by ∼ 0.6 dex. The analysis of Mannucci et al. (2010) indicates the evolution of FMR from z ∼ 2 to 3, but it is not clear why the FMR evolves only between z ∼ 2 and z ∼ 3. In contrast to the FMR evolution, we have argued in Section 4.5 that no significant differences are found between galaxies at z ∼ 2 and 3 in the FIR (Figure 7). We address the question why the FMR shows the evolution but no evolution is found in the FIR. Mannucci et al. (2010) use metallicities of z 3 galaxies given by Maiolino et al. (2008) and Mannucci et al. (2009) that are estimated with the local empirical relations of Maiolino et al. (2008). The local empirical relations of Maiolino et al. (2008) are obtained from two samples, (i) 259 local galaxies with 12 + log(O/H) < 8.4 whose metallicities are estimated by the electron temperature Te method, and (ii) 22, 482 star-forming SDSS galaxies with 12 + log(O/H) > 8.2 given with the photoionization model (Kewley & Dopita 2002; see also Nagao et al. 2006). Figure  8 presents R23-index and [O iii]/[O ii] ratio of the two samples. The red dashed curves denote the local empirical relations determined with these two samples (Maiolino et al. 2008). Figure 8 compares the photoionization models (black solid curves) with the local empirical relations, and indicates that the local empirical relations implicitly assume an ionization parameter. Given the fact that high-z galaxies have an ionization parameter higher than local galaxies (Section 5.1), the local empirical relations assuming a relatively low ionization parameter would provide biased metallicity estimates for high-z galaxies. At z 2.5, the [N ii]λ6584/Hα ratio can be used to estimate metallicities and/or to discriminate one from two metallicity solutions given by R23-index. However, at z 2.5 [N ii]λ6584 and Hα lines cannot be detected from the ground. The [O iii]/[O ii] ratio is needed to determine one from two metallicity solutions of R23-index. The bottom panel of Figure 8 shows that the photoionization model predicts a lower metallicity if a lower ionization parameter is assumed. Thus, one may underestimate metallicities of z 3 galaxies with the local empirical relations, since the ionization parameter of z 3 galaxies is significantly higher than those of local galaxies (Section 5.1).
We evaluate the bias of z ∼ 3 galaxy metallicity estimates, which is originated from the local empirical relation. We estimate metallicities of z ∼ 3 galaxies with the Kobulnicky & Kewley (2004) method in the same manner as Section 4.3, allowing the ionization parameter evolution. We compile a total of 15 z ∼ 3 galaxies from AMAZE (Maiolino et al. 2008) and LSD (Mannucci et al. 2009) projects. First, for a consistency check, we estimate metallicities of our 15 z ∼ 3 galaxies with the local empirical relation used by Mannucci et al. (2010). We find that the average metallicity of the z ∼ 3 galaxies is 12 + log(O/H) ∼ 8.1 (left panel of Figure 9), which is consistent with that of Mannucci   . z ∼ 3 galaxies in the FMR plot. Left: The black dashed curve and the gray shaded area indicate the FMR and its typical dispersion, respectively (Mannucci et al. 2010(Mannucci et al. , 2011. The black open and filled circles represent galaxies at z = 1-2 and z ∼ 3, respectively, that are compiled by Mannucci et al. (2010). The green circles denote z ∼ 3 galaxies (Maiolino et al. 2008;Mannucci et al. 2009) whose metallicities are estimated with the local empirical relations in the same manner as Mannucci et al. (2010). These green circles are consistent with the average value shown in Mannucci et al. (2010). Right: Same as the left panel, but for metallicity estimates using the Kobulnicky & Kewley (2004) method that is free from the bias of ionization parameters. The red and blue circles denote metallicities of the z ∼ 3 galaxies in the cases of high-and low-Z branches, respectively. Four z ∼ 3 galaxies are omitted from the plot, since their R23-index values are out of the range where the Kobulnicky & Kewley (2004)  obtain average metallicities of our 11 z ∼ 3 galaxies for the low-Z and high-Z branches that are 12 + log(O/H) = 8.20 and 8.65, with ionization parameters of qion ∼ 6.6 × 10 7 and 1.1 × 10 8 cm s −1 , respectively. Typically ∼ 0.2 − 0.3 dex errors are associated with the estimates of metallicity and ionization parameter for each of the galaxies. Note that four out of fifteen galaxies are omitted from these metallicity estimates, since their R23-index values fall out of the range covered by Equations (6) and (7). The right panel of Figure 9 compares the metallicity estimate of Mannucci et al. (2010) with ours allowing the ionization parameter evolution. The high-Z branch metallicities appear to follow the FMR (black dashed curve; Mannucci et al. 2010Mannucci et al. , 2011. In contrast, the low-Z branch metallicities fall below the FMR, which is consistent with the result of Mannucci et al. (2010).
Because there is no useful line ratio, such as [N ii]λ6584/Hα, that determines one from two metallicity solutions of low and high-Z branches, we cannot conclude which metallicities are correct and whether z ∼ 3 galaxies follow the FMR. However, we find that the high ionization parameter of ∼ 1.1 × 10 8 cm s −1 for the high-Z branch solution agrees with the our findings of ionization parameter evolution from z = 0 to 3 (Section 5.1). Moreover, there are no physical reasons for the sudden departure from the FMR from z = 2 to 3, e.g., the smooth evolution of cosmic SFR/stellar mass densities at z 2 (Bouwens et al. 2011;Stark et al. 2013). Theoretical studies do not find the departure from the FMR at z ∼ 3. Our metallicity estimates of high-Z branch of z ∼ 3 galaxies are comparable with those predicted by Davé   simple analytic model. These pieces of evidence imply that the z ∼ 3 galaxy metallicities of high-Z branch with the high-ionization parameter are correct, and that z ∼ 3 galaxies follow the FMR.
As local counterparts of high-z galaxies, we also examine LBAs in the FMR plot in Figure 10. Their metallicities are estimated in the same manner as those for z ∼ 3 galaxies. However, their two metallicity solutions can be distinguished, since other metallicity indicators of gas temperature and/or [N ii]λ6584/Hα ratio are available for the local galaxies. The filled circles display metallicities that agree better with those estimated by the other metallicity indicators. Figure 10 suggests that LBAs follow the FMR if their metallicities are estimated with ionization parameter. We have also checked that the FMR is valid for GPs, whose metallicities are estimated in the same manner as those for LBAs. These findings confirm the validity of the FMR over the wide ranges of stellar mass and SFR in the local universe. Another notable trend in Figure 10 is that the LBAs metallicities of high-Z branch with high-ionization parameter are generally correct. This is what we argue for the z ∼ 3 galaxies. Since LBAs are local UV-selected galaxies whose characteristics are similar to those for LBGs (e.g., Overzier et al. 2009), their trends in the FMR plot support our argument that z ∼ 3 LBGs have metallicities of high-Z branch with high-ionization parameter and follow the FMR.
For more conclusive discussion of z ∼ 3 galaxies, one needs deeper NIR spectroscopy to provide an ionization parameter indicator of [Ne iii]λ3869/[O ii]λ3727 (e.g., Nagao et al. 2006;Richardson et al. 2013) or IR (λ > 3µm) spectroscopy from the space to obtain [N ii]λ6584/Hα that determines one from two metallicity solutions. In particular, the [Ne iii]/[O ii] ratio is observable from the ground up to z ∼ 5. We have confirmed the utility of the [Ne iii]/[O ii] ratio for LAEs, whose ionization parameters are very high (Fosbury et al. 2003;Christensen et al. 2012b).

Fundamental Ionization Relation
As demonstrated in Section 5.2, ionization parameter is important for metallicity estimates. Although the FMR reproduces properties of galaxies with no significant evolution, input metallicities to the FMR would be sometimes biased due to the evolution of ionization parameters. In this sense, the FIR is advantageous because the FIR allows the difference of ionization parameter. The FIR consisting of the fiducial µ0.32 parameter requires non-zero β, where β = 0 is ruled out at the > 95 % confidence level (Section 4.5). This indicates that the four-dimensional relation of FIR between ionization parameter, metallicity, SFR, and M⋆ is useful to characterize gas properties of galaxies.

Relationship between Ionization State and Ionizing Photon Escape
We have found in Figure 2   To examine the dependence of [O iii]/[O ii] ratio on fesc, we perform photoionization model calculations with cloudy (version 13.02; Ferland et al. 1998Ferland et al. , 2013. We assume a constant-density homogeneous inter-stellar gas cloud with a spherically closed geometry. A volume filling factor of the ionized gas is assumed to be unity for simplification. We include dust physics and the depletion factors of the various elements from the gaseous phase in the same manner as the analysis of Dopita et al. (2006b). In non-solar metallicities, we assume that both the dust model and the depletion factors are unchanged, but the dust abundance is assumed to scale linearly with the gas metallicity, as adopted in Nagao et al. (2011). All elements except for nitrogen, carbon, and helium are taken to be primary nucleosynthesis elements. For nitrogen, we use a form given by López-Sánchez et al. (2012) to take account of its secondary nucleosynthesis component in a high metallicity range. For carbon and helium, we use forms in Dopita et al. (2006b). Photoionization models are calculated with the parameters of Z = (0.05, 0.1, 0.2, 0.5, 1.0) Z⊙, log(qion/cm s −1 ) = (7.5, 8.0, 8.5), and nH = 10 2 cm −3 10 . The incident ionizing radiation fields are generated by Star- 10 We have checked that changing the density to n H = 10 cm −3 and 10 3 cm −3 does not alter our results. Since both [O iii] and [O ii] are collisionally excited lines, their intensities are similarly enhanced if the density is smaller than their critical den-!"" !" !""

Ionization-bounded nebula Density-bounded nebula
(a) (b) Figure 12. Schematic illustrations of two H ii regions. A photoionized H ii region is presented with blue, and the outer H i region is shown with gray. Yellow stars are central ionizing sources.
(a) An ionization-bounded nebula whose radius is determined by the ionization equilibrium. The cyan shows the Strömgren radius. (b) A density-bounded nebula whose radius is determined by the distribution of gas cloud. The surrounding H i cloud is small enough that the central sources can ionized it completely. The nebula thus cannot form a complete Strömgren sphere.
burst99 (Leitherer et al. 1999). We use instantaneous burst models at zero age with a Kroupa (2001) IMF and lower and upper mass limits of 0.1 and 120 M⊙, respectively. Stellar metallicities are matched to the gas-phase ones. Under these assumptions and conditions, we evaluate fesc with the number of ionizing photons escaping from a thin cloud with a low column density. We adopt a neutral hydrogen column density, NHI, as the stopping criterion of our calculations that ranges from NHI = 10 15 to 10 20 cm −2 . The escape fraction is defined as a ratio of transmitted ionizing photons (λ < 912Å) to input ionizing photons. Figure 11 shows (1) Ionization-bounded nebula whose radius is determined by the ionization equilibrium between the ionizing photon production rate and the recombination rate (Equation 3; Figure 12a). The radius of the ionized region is called the Strömgren radius.
(2) Density-bounded nebula whose radius is determined by the distribution of gas cloud (Figure 12b). The densitybounded nebula does not expend all of ionizing photons, but leak ionizing photons that are not used for the ionization of the nebula.
The ionizing photon escape takes place in the densitybounded nebula under the assumption of homogeneous gas cloud. A density-bounded nebula with a less NHI emits more ionizing photons, and has a higher fesc (Giammanco et al. 2005 1990; Oey & Kennicutt 1997;Pellegrini et al. 2012). Thus, the O + zone of density-bounded nebula is smaller than that of an ionization-bounded nebula, but the size of O 2+ zone is similar in the density-bounded and ionization-bounded nebulae. and fesc. We note that we consider a very simple case, with the assumption that H iiregions are homogeneous density-bounded nebulae. In reality, conditions of H ii-regions with ionizing photon escape could be more complex, e.g., ionization-bounded nebulae with highly ionized low density holes (e.g., Zackrisson et al. 2013; see also the end of this section). However, we present in this paper the simple case to investigate the possibility that the [O iii]/[O ii] ratio could be affected by the condition of H ii-regions where ionizing photon escapes take place. It is beyond the scope of this paper whether or not the assumption of the density-bounded nebulae is reasonable.
The assumption remains to be tested by future observations.
In Figure 11 In Figure 7, we have found that most of galaxies at z = 0 − 3 follow the FIR within the uncertainties. However, some galaxies clearly depart from the FIR at a few sigma level. The two z ∼ 0 LyC leakers have an [O iii]/[O ii] ratio higher than the FIR by a factor of ∼ 2 − 3 at high significance. The enhancement of [O iii]/[O ii] ratio may be caused by an ionizing photon escape. If it is true, some GPs and LAEs lying above the FIR would be good candidates of ionizing photon emitting objects. In Figure 7, LAEs particularly have a very high [O iii]/[O ii] ratio, some of which are departing from the FIR more significantly than LBGs over the sizes of error bars. This physical characteristics would be originated from a very large fesc given by a low NHI of LAE. The low NHI of LAE is also supported by the gasdynamics study of Hashimoto et al. (2013) who claim that LAEs typically have a NHI lower than LBGs from the analysis of velocity offsets between a Lyα line and low-ionization interstellar absorption lines with respect to the systemic velocity. This tendency is confirmed with a larger number of LAEs by Shibuya et al. (2014a). In this companion paper, Shibuya et al. (2014b) examine LAE structures, suggesting that LAEs with a large Lyα EW tend to have a small ellipticity. This is consistent with the theoretical results that Lyα photons can more easily escape from face-on disks having a small ellipticity, due to a low NHI (e.g., Verhamme (Feldmeier et al. 2013). Since most of the observed Lyα emission in the diffuse halos is probably originated from the galaxy H ii-regions and scattered to the line of sight by extended H i gas in the galaxy's circum-galactic medium (CGM), LAEs could exhibit Lyα halos weaker than LBGs due to their smaller amounts of H i gas in the CGM. This scenario is consistent with our argument of low NHI of LAE. Interestingly, that interpretation also consistently explains the results of direct LyC detections of z ∼ 3 galaxies that reveal a high fesc for LAEs (∼ 0.1 − 0.3) than LBGs (∼ 0.05) (Iwata et al. 2009;Nestor et al. 2011Nestor et al. , 2013. Thus, an origin of the high fesc in LAEs would be their optically thin H i gas surrounding star-forming regions. A similar suggestion of the high fesc is inferred for GPs (Jaskot & Oey 2013 LAEs with a high fesc would play a key role in supplying ionizing photons for cosmic reionization in the early universe. Previous z ∼ 6 − 7 galaxy surveys have suggested that the universe could not be totally ionized by galaxies alone at z ∼ 6 − 7 if one assumes similar properties of galaxies at z ∼ 3 that include fesc ∼ 0.05 (e.g., Ouchi et al. 2009;Robertson et al. 2010). Because the fraction of LAEs in star-forming galaxies increases with redshift (at least up to z ∼ 6; e.g., Ouchi et al. 2008;Vanzella et al. 2009;Stark et al. 2010Stark et al. , 2011Hayes et al. 2011), the average fesc for z ∼ 6 galaxies would be higher than that for z ∼ 3 galaxies. The higher fesc for z ∼ 6 galaxies may resolve the problem of ionizing photon deficit in cosmic reionization. Note that the number of observed LAEs suddenly drops from z ∼ 6 to 7. This could be due to the Lyα damping absorption of partly neutral IGM (Vanzella et al. 2011;Pentericci et al. 2011;Ono et al. 2012;Schenker et al. 2012) or optically thick systems in the intervening IGM (Bolton & Haehnelt 2013). Alternatively, the decreasing fraction of LAEs could be explained by the intrinsic decrease of galaxies with strong Lyα emission at z > 6 ( Dijkstra et al. 2014). If we assume the former scenarios, galaxies with an intrinsically bright Lyα line would not decrease from z ∼ 6 to 7. However, this is still an open question. A similar implication is provided by Schenker et al. (2013) on the basis of a small velocity offset of Lyα line with respect to the systemic velocity for z = 3−4 galaxies. We also note that this picture would be consistent with the observational fact that gas fraction of star-forming galaxies increases with redshift, or more properly, with sSFR (e.g., Tacconi et al. 2010Tacconi et al. , 2013Santini et al. 2014). If the molecular gas exists in dense cores of clumps embedded in density-bounded nebulae ionized by young stars, the galaxy could have a high gas fraction and high fesc. A dense molecular gas clump could reside inside the H ii-region, since the dense gas could be self-shielded from ionizing radiation. A good example is found in the Orion nebula (see, e.g., Genzel & Stutzki 1989 for a review). Such clumpy geometry of dense molecular gas could explain the high gas fraction, active star-formation, and high fesc in high-z LAEs.
In the discussion of high fesc for LAEs, there is one caveat that Lyα emission would become weaker in a high fesc, because less ionizing photons are spent for making ionized gas that emit Lyα. Therefore, LAEs may not be a population of high fesc. The dashed curve in Figure 13 shows the relation between EW(Lyα) and fesc predicted by our cloudy models. Here we assume qion = 10 8 cm s −1 and Z = 0.2 Z⊙, typical values observed in GPs (Amorín et al. 2010) and z ∼ 2 LAEs (Nakajima et al. 2013). The dashed curve presents that the EW(Lyα) decreases sharply if its fesc are over ≃ 0.8. However, in fesc ≃ 0.0 − 0.8, EW(Lyα) remains unchanged within a factor of 2. We have checked that the overall trends are the same, adopting other values of qion and Z. For example, EW(Lyα) increases (decreases) by a factor of ∼ 2 − 3 for Z = 0.05 Z⊙ (1 Z⊙). LAEs are usually defined as galaxies with EW(Lyα) > 20Å. If an escape fraction is very high, fesc 0.9, a galaxy cannot be observed as an LAE. However, LAEs, or galaxies with EW(Lyα) > 20Å, can exist in the range of fesc 0.8. Since a very high fesc of 0.8 is not required for cosmic reionization but only fesc ≃ 0.2 , the weakening of Lyα in density-bounded nebulae does not rule out the possibility that major ionizing sources are LAEs. Figure 13 also compares our model and the recent observational results for z ∼ 3 LAEs and LBGs whose fesc are measured from LyC fluxes (Nestor et al. 2013) 11 . EW(Lyα) values of our model are regarded as upper limits, because EW(Lyα) decreases by dust extinction and hydrogen scattering. The observational results fall below our model, and are consistent with our model. The observations indicate a positive correlation between fesc and EW(Lyα) in fesc 0.5, while our model shows a nearly constant EW(Lyα) in this regime. This difference would infer that LAEs have a less dust extinction, a more homogeneous geometry for isotropic Lyα scattering, a less metallicity, and more ionizing photons than LBGs.
We note that our calculations assume a densitybounded nebula with a spherically closed geometry, and predict dependencies of NHI on the emergent spectra. A volume filling factor of the ionized gas is assumed to be unity. Thus, we have presented this specific case of EW(Lyα) as a function of fesc. However, the EW(Lyα) and fesc involve various complex physical effects. EW(Lyα) depends on gas and dust geometries (e.g., Neufeld 1991;Scarlata et al. 2009;cf. Duval et al. 2014) and galactic scale outflow (e.g., Kunth et al. 1998;Verhamme et al. 2006). fesc anti-correlates with a covering fraction of neutral gas. Recent spectroscopic studies suggest a low covering fraction of neutral gas in high-z star-forming galaxies (e.g., Shapley et al. 2003;Steidel et al. 2010;Jones et al. 2013). There exist the other complex physics of EW(Lyα) and fesc. Nevertheless, our analysis reveals the relation of EW(Lyα) and fesc for the case that H ii-regions are homogeneous density-bounded nebulae.

SUMMARY
We investigate ionization state of hot ISM in galaxies at z = 0 − 3 with ∼ 140, 000 SDSS galaxies and 108 intermediate and high redshift galaxies, using an ionization parameter and metallicity sensitive line ratios of [O iii] • We confirm that z ∼ 2 − 3 galaxies show an [O iii]/[O ii] ratio significantly higher than a typical star-forming galaxy of SDSS by a factor of 10 ( Figure 2). The photoionization models reveal that these high-z galaxies have an ionization parameter of log(qion/cm s −1 ) ∼ 7.6−9.0, a factor of ∼ 4−10 higher than local galaxies (Figure 3).
• We find a high ionization parameter in local extreme populations of green pea galaxies (GPs), Lyman-continuum emitting galaxies (LyC leakers), and Lyman break analogs (LBAs), and the ionization parameters of these galaxies are comparable with the ones of z = 2 − 3 galaxies (Figure 3). These local extreme populations are quite rare, and their number densities are nearly two orders of magnitude smaller than those of z = 2 − 3 galaxies.
• We identify a correlation between the [O iii]/[O ii] ratio and galaxy global properties of star-formation rate (SFR), stellar mass (M⋆), and metallicity ( Figure 4). A high [O iii]/[O ii] ratio is found in less massive, efficiently starforming, and metal-poor galaxies.
• We develop a four-dimensional relation of ionization parameter, SFR, M⋆, and metallicity ( Figure 5), extending the fundamental metallicity relation (FMR; Mannucci et al. 2010;Lara-López et al. 2010) with ionization parameter. The intermediate and high-z galaxies up to z ≃ 3 as well as the local extreme populations follow the same relation defined with the SDSS galaxies (Figure 7). The relation is thus referred to as the fundamental ionization relation (FIR).
• We particularly find no significant differences between galaxies at z ∼ 2 and 3 in the FIR (Figure 7). This is in contrast with the FMR whose possible evolution from z ∼ 2 to 3 is reported (Mannucci et al. 2010). We suggest that the FMR evolution can arise, if one omits ionization parameter differences and adopts local empirical metallicity relations for high-z galaxies. We indicate that the FMR evolution does not exist for one out of two average metallicity solutions of z ∼ 3 galaxies with a high ionization parameter of log(qion/cm s −1 ) 8 (right panel of Figure 9).
• We find the z ∼ 0 LyC leakers exhibit [O iii]/[O ii] ratios significantly higher than typical low-z galaxies, although there are only two LyC leakers. This trend suggests a positive correlation between [O iii]/[O ii] and ionizing photon escape fraction (fesc). Our photoionization models of cloudy reproduce the trend (Figure 11). Since [O iii]/[O ii] ratios for z = 2 − 3 galaxies, especially Lyα emitters (LAEs), are comparable to, or higher than, those of the low-z LyC leakers, these high-z galaxies are candidates of high fesc objects.
• Our cloudy photoionization models reveal that EW(Lyα) remains almost unchanged if its fesc is ≃ 0 − 0.8 (Figure 13), indicating that galaxies are identified as LAEs of a Lyα emitting population if their fesc values are below ≃ 0.8. Faint galaxies with intrinsically-bright Lyα having a high fesc could significantly contribute to the ionizing photon production for cosmic reionization, since recent observations report that a fraction of LAEs in star-forming galaxies increases with redshift at the epoch with no IGM damping absorption (e.g., Ouchi et al. 2008;Vanzella et al. 2009;Stark et al. 2010Stark et al. , 2011. If it is true, this would ease the tension of ionizing photon deficit problem (e.g., Ouchi et al. 2009;Robertson et al. 2010).