Bridging Optical and Far-Infrared Emission-Line Diagrams of Galaxies from Local to the Epoch of Reionization: Characteristic High [O III] 88 $\mathrm{\mu m}$/SFR at $z>6$

We present photoionization modeling of galaxy populations at $z\sim0$, $2$, and $>6$ to bridge optical and far-infrared (FIR) emission-line diagrams. We collect galaxies with measurements of optical and/or FIR ([O III] 88 $\mathrm{\mu m}$ and [C II] 158 $\mathrm{\mu m}$) emission line fluxes and plot them on the [O III]$\lambda5007/\mathrm{H\beta}$--[N II]$\lambda6585/\mathrm{H\alpha}$ (BPT) and L([O III]88)/SFR--L([C II]158)/SFR diagrams, where SFR is the star-formation rate and L([O III]88) and L([C II]158) are the FIR line luminosities. We aim to explain the galaxy distributions on the two diagrams with photoionization models that employ three nebular parameters: the ionization parameter $U$, hydrogen density $n_\text{H}$, and gaseous metallicity $Z_\text{gas}$. Our models successfully reproduce the nebular parameters of local galaxies, and then predict the distributions of the $z\sim0$, $2$, and $>6$ galaxies on the diagrams. The predicted distributions illustrate the redshift evolution on all the diagrams; e.g., [O III]$/\mathrm{H\beta}$ and [O III]88/[C II]158 ratios continuously decrease from $z>6$ to $0$. Specifically, the $z>6$ galaxies exhibit $\sim\!0.5$ dex higher $U$ than low-redshift galaxies at a given $Z_\text{gas}$ and show predicted flat distributions on the BPT diagram at $\log{\mathrm{[O III]/H\beta}} = 0.5$-$0.8$. We find that some of the $z>6$ galaxies exhibit high L([O III]88)/SFR ratios. To explain these high ratios, our photoionization models require a low stellar-to-gaseous metallicity ratio or bursty/increasing star-formation history at $z>6$. The James Webb Space Telescope will test the predictions and scenarios for the $z>6$ galaxies proposed by our photoionization modeling.

1. INTRODUCTION Emission-line ratios reflect physical properties of the interstellar medium (ISM) of galaxies.In the local universe, optical to far-infrared (FIR) spectroscopy has provided various emission-line measurements to probe essential ISM properties, including the electron temperature, electron density, ionization state, and elemental abundances.However, direct comparisons of the ISM properties at different redshifts are challenging due to limited sensitivity and wavelength coverage of observational instruments.
Comparisons with local galaxies have unveiled the stellar and ISM properties of galaxies at z ∼ 2 and z > 6; however, galaxies in the two high-redshift ranges cannot be directly compared with each other because no emission lines are commonly observed for the two populations.One of the key tools to overcome this difficulty is photoionization models.Photoionization models simplify ISM structures of galaxies under some assumptions and predict emissionline intensity ratios from an input ionizing spectrum and nebular physical parameters.They are widely adopted to galaxy population at each redshift, including local galaxies like dwarf galaxies (Cormier et al. 2019) and (ultra-)luminous infrared galaxies (U/LIRGs; e.g., Nagao et al. 2011;Inami et al. 2013;Pereira-Santaella et al. 2017) and z ∼ 2 galaxies (e.g., Steidel et al. 2014;Sanders et al. 2016;Trainor et al. 2016;Strom et al. 2018).For z > 6 galaxies, Harikane et al. (2020b)  line luminosities, respectively, and SFR is the star-formation rate.They concluded that high-redshift galaxies have high ionization parameter, low PDR covering fraction, or both.By applying a similar model to Harikane et al., Sugahara et al. (2021) constrained the N/O ratio of B14-65666 at z = 7.15 as a function of the metallicity.
In this paper we attempt to bridge local to z > 6 galaxy populations by modeling their distributions on diagrams of various emission-line ratios.Our attempt derives average ISM properties of these galaxy populations and infers any evolutionary trends as a function of redshift.This paper is organized as follows.Section 2 describes samples of galaxy populations at z ∼ 0, 2, and > 6 that have optical and/or FIR emission-line measurements.Section 3 explains setups of photoionization models.Section 4 shows results of photoionization modeling of galaxy populations in emission-line diagrams.Section 5 discusses implications of line ratios, parameter dependence of photoionization models, and caveats.Finally Section 6 summarizes our findings.Throughout this paper, we use the Chabrier (2003) initial mass function (IMF) in the mass range of 0.1-100 M and convert SFR from Salpeter (1955) and Kroupa (2001) IMFs by multipying factors of 0.63 and 0.94(= 0.63/0.67),respectively (Madau & Dickinson 2014).In this paper the two diagrams are referred to as the BPT and FIR diagrams, respectively.We collected galaxies with emissionline measurements that can be plotted on the BPT and FIR diagrams.Section 2.1 constructs a sample of local galaxies with both optical and FIR line measurements.Section 2.2 constructs a statistical sample of galaxies at z ∼ 0 and 2 with only optical line measurements.Section 2.3 constructs a sample of z > 6 galaxies observed with ALMA.

OptFIR sample: local dwarfs and U/LIRGs
The first sample consists of local galaxies with both optical (Hβ λ4861, [O III] λ5007, Hα λ6563, and [N II] λ6585) and FIR ([O III] 88 and [C II] 158 ) line measurements, being referred to as the optFIR sample.We firstly describe FIR observations, with the PACS instrument (Poglitsch et al. 2010) on board Herschel, and then explain optical line measurements.This sample is mainly divided into two galaxy populations: local dwarfs and local U/LIRGs.
The local dwarfs were taken from the Dwarf Galaxy Survey (DGS, PI: Madden;Madden et al. 2013;Cormier et al. 2015), which is targeting 50 dwarf galaxies.We used 36 galaxies for which Cormier et al. (2015)  Here, values from De Looze et al. have priority if galaxies are listed in both references.The stellar mass spans from log M * /M = 6.5 to 10.5 with a mean of 8.6 and a standard deviation of 1.0 dex (Madden et al. 2013).
The local U/LIRGs were composed of 61 galaxies taken from the Great Observatories All-Sky LIRG Survey (GOALS; Armus et al. 2009;Díaz-Santos et al. 2013) and 14 galaxies taken from the Survey of Far-infrared Lines The red circles show the z > 6 galaxies; the filled circles show the measurements mainly taken from Harikane et al. (2020b), while the hashed circles from Carniani et al. (2020).The open red circles are taken from Harikane et al. (2020b), but their SFR is measured with the SED fitting.The gray triangles show the local dwarfs (Cormier et al. 2019) and the orange symbols show the local U/LIRGs including the GOALS galaxies (square, Díaz-Santos et al. 2017) and the SHINING galaxies (cross, Herrera-Camus et al. 2018a).The solid polygons depict the regions representing the distributions of the galaxy populations, which are used in our photoionization modeling.Left: The black and blue contours illustrate the distributions of the z ∼ 0 (SDSS) and z ∼ 2 (Strom et al. 2017;Shivaei et al. 2018) galaxies, respectively, which include 68, 95, and 99.5% of the z ∼ 0 galaxies and 68 and 90% of the z ∼ 2 galaxies.The dashed gray lines depict the criteria between star-forming galaxies and AGNs (Kauffmann et al. 2003;Kewley et al. 2001a).Right: The horizontal and vertical error bars indicate the [C II]158 and [O III]88 measurement errors, respectively, and the diagonal error bars does the SFR errors.For local galaxies, the black edges mean the data points with optical measurements, while the gray edges mean ones without them.
in Nearby Galaxies (SHINING, PI: Strum;Herrera-Camus et al. 2018a,b).The selected galaxies have both [C II] 158 and [O III] 88 line flux measurements.We excluded from the sample active galactic nuclei (AGNs) classified in the literature (Rich et al. 2015;Herrera-Camus et al. 2018a), from the bolometric AGN fractional contribution of α bol AGN ≥ 0.5 (Díaz-Santos et al. 2017), and on the BPT diagram (Kewley et al. 2001a).For the GOALS galaxies, the line luminosities were converted from "Best line flux" given by Díaz-Santos et al. ( 2017) using the luminosity distance in Armus et al. (2009).SFR was taken from Howell et al. (2010), who computed it from GALEX FUV and IRAS L TIR (Kennicutt 1998).For the SHINING galaxies, the line luminosities were converted from "integrated line fluxes" given by Herrera-Camus et al. (2018a) using given redshifts.SFR was computed from L TIR (Murphy et al. 2011), where L TIR = 1.75LFIR (Herrera-Camus et al. 2018a) and L FIR is the FIR luminosity at 42.5-122.5µm (Helou et al. 1988).The stellar mass is log M * /M 10.6-11.6 with a mean of 11.1 for the GOALS galaxies (Howell et al. 2010), which is similar to the mass of the SHINING galaxies (Herrera-Camus et al. 2018b).We note that the U/LIRGs defined here include galaxies that do not satisfy the definition of LIRGs, L IR > 10 11 L , but we categorized them as U/LIRGs in this work for simplicity.Optical emission-line ratios, [O III]/Hβ and [N II]/Hα, are taken from the literature (Veilleux et al. 1999;Kewley et al. 2001b;Moustakas & Kennicutt 2006;Moustakas et al. 2010;Rich et al. 2015;De Vis et al. 2017;Perna et al. 2021) 1 .When a galaxy is listed in multiple studies, the measurements from the integrated spectroscopic observations (Moustakas & Kennicutt 2006;Moustakas et al. 2010;Rich et al. 2015;Perna et al. 2021) have priority, considering a consistency with a large field of view of FIR observations.In this way, optical measurements were taken for 32 out of the DGS galaxies, 29 out of the GOALS galaxies, and 9 out of the SHINING galaxies.
Figure 1 shows the distributions of the optFIR sample on the BPT and FIR diagrams.The dwarfs exhibit higher [O III]/Hβ and L([O III] 88 )/SFR values than the U/LIRGs, suggesting higher ionization states and lower metallicities of the dwarfs.The low L([C II] 158 )/SFR and L([O III] 88 )/SFR ratios in the U/LIRGs are known as "line deficit" (Díaz-Santos et al. 2017;Herrera-Camus et al. 2018a,b), lower line-to-L FIR ratios at higher L FIR or higher SFR.

Optical sample: z ∼ 0 and 2 galaxies
The second sample consists of galaxies at z ∼ 0 and 2 that have measurements of the optical [O III]/Hβ and [N II]/Hα ratios but not the FIR observations.We refer to this sample as the optical sample.
The z ∼ 0 galaxies were drawn from the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data Release 7 (Abazajian et al. 2009) main galaxy sample (Strauss et al. 2002).The emission-line fluxes are given in the MPA-JHU catalog 2 .We selected ∼20,000 star-forming or starburst galaxies detected in Hβ, [O III], [N II], and Hα lines with S/N > 5. Composites of AGNs and starbursts were removed based on the criterion of Kauffmann et al. (2003) on the BPT diagram.The mean and standard deviation of the stellar mass taken from the same catalog are log M * /M = 10.0 and 0.7 dex, respectively (Kauffmann et al. 2003).
The z ∼ 2 galaxies consisted of Keck Baryon Structure Survey observed with MOSFIRE (KBSS-MOSFIRE; Steidel et al. 2014) and the MOSFIRE Deep Evolution Field survey (MOSDEF; Kriek et al. 2015).We took the emission-line ratios of 226 galaxies from Strom et al. (2017) for the KBSS-MOSFIRE survey and those of 223 galaxies from Shivaei et al. (2018) for the MOSDEF survey.Upper limits of the line ratios were not included in the sample.We also did not include objects as AGN candidates whose [O III]/Hβ values are > 0.2 dex higher than the criterion of Kewley et al. (2001a).The stellar mass for the KBSS-MOSFIRE galaxies spans from log M * /M 8.6 to 11.4 with a mean of 10.0 (Steidel et al. 2014) and the stellar mass for the MOSDEF galaxies similarly spans from 9.0 to 11.5 (MOSDEF; Kriek et al. 2015).
The contours in the left panel of Figure 1 illustrate the distributions of the z ∼ 0 and 2 galaxies on the BPT diagram.The black contours include 68, 95, and 99.5% of the z ∼ 0 galaxies, while the blue contours include 68 and 90% of the z ∼ 2 galaxies.The figure shows a clear offset of the distributions between the z ∼ 0 and 2 galaxies, indicating the evolution of ISM properties as described in Section 1.
The z > 6 galaxies are provided by two compilations, Harikane et al. (2020b) and Carniani et al. (2020).Harikane et al. (2020b)  nosities and SFR of nine UV-selected galaxies (Lyman-break galaxies and Lyman-alpha emitters) and three submillimetergalaxies (SMGs) from the literature.Carniani et al. (2020) reanalyzed ALMA data of the nine UV-selected galaxies with uv-tapering to report [C II] 158 detections for the galaxies in which previous studies did not detect.Carniani et al. also presented SFR UV+IR , SFR computed from the sum of the UV and IR luminosity, L UV+IR (Kennicutt & Evans 2012).In this work, we used the [C II] 158 and [O III] 88 line luminosities and SFR from the both compilations.For J1211-0118, J0235-0532, and J0217-0208 (Harikane et al. 2020b), we used the SFR values updated by Ono et al. in prep.Moreover, we added z7 GSD 3811 (Binggeli et al. 2021) to our sample.Although deriving the stellar mass is difficult for the z > 6 galaxies, it was inferred from the spectral-energydistribution (SED) fitting and empirical relations of UV photometry in the literature.The inferred stellar mass spans from log M * /M ∼ 9 to 10.5 (Inoue et al. 2016;Laporte et al. 2017;Marrone et al. 2018;Hashimoto et al. 2018Hashimoto et al. , 2019;;Tamura et al. 2019;Harikane et al. 2020b;Binggeli et al. 2021).
The sample finally consists of 13 galaxies with 22 data points.In the right panel of Figure 1, the red circles show the z > 6 galaxies on the FIR diagram.Regarding the galaxies that have different measurements from the two compilations, the filled red circles show measurements for Harikane et al. (2020b) and the hashed red circles for Carniani et al. (2020).The open red circles depict the galaxies whose SFR is computed from SED fitting, SFR SED .We mainly analyze galaxies with SFR UV+IR in this work.For convenience, Figure A1 corresponds the names of the z > 6 galaxies to the data points on the FIR diagrams.
The L([O III] 88 )/SFR ratios of the z > 6 galaxies are comparable to or higher than the local dwarfs, despite their high SFR.These high L([O III] 88 )/SFR ratios at a given SFR are consistent with previous results of higher line-to-L FIR ratios at z > 1 than local line deficit at a given L FIR (e.g., Herrera-Camus et al. 2018a;Sugahara et al. 2021)

FIDUCIAL PHOTOIONIZATION MODELS
We used CLOUDY version 17.02 (Ferland et al. 1998(Ferland et al. , 2017) ) to construct photoionization models.One of the major assumptions in our analysis is that a single nebular parameter set can reproduce the optical and FIR emission-line ratios simultaneously.Therefore, we ignored contributions from the secondary component like diffuse ionized gas (DIG).
Following assumptions are adopted as fiducial models, which partly share the concepts with Nagao et al. (2011Nagao et al. ( , 2012)), Harikane et al. (2020b) and Sugahara et al. (2021).The adopted assumptions are summarized in Table 1.A nebular structure is a plane-parallel geometry under the constant The colors of the four grids indicate log nH/cm −3 = 0.5, 1.5, 2.5, and 3.5 from purple (dark) to yellow (light).The symbols are the same as Figure 1, but we removed the symbols with gray edges, removed the error bars, and lightened colors.The gray and cyan contours in the left panel show the 95 and 68 percentile contours for the z ∼ 0 and 2 galaxies, respectively.The model grids can explain most of the data points well with the three parameters, U , Zgas, and nH.
pressure.Elemental and dust grain abundances are the H II region abundances, but for helium and nitrogen.A helium abundance reflects a combination of Big Bang nucleosynthesis and stellar yields, given by He/H = 0.0737 + 0.0293Z gas /Z (1) (Groves et al. 2004), where Z gas is the gaseous metallicity and Z is the solar metallicity of 0.02.A nitrogen abundance would reflect a combination of primary and secondary nucleosynthesis as a function of metallicity (Pilyugin et al. 2012;Andrews & Martini 2013).Here, all over the redshift, N/O ratios keep constant at low Z gas and then follow a local relation of extragalactic H II regions (Pilyugin et al. 2014): where Z ≡ log Z gas /Z , the solar metallicity is assumed to be 12 + log(O/H) = 8.69 (Asplund et al. 2009), and thus the criterion Z = −0.59corresponds to 12 + log(O/H) = 8.1.
Equation 2 gives N/O values comparable to those of the KBSS-MOSFIRE galaxies at z ∼ 2 (Strom et al. 2018) and that of B14-65666 at z = 7.15 (Sugahara et al. 2021).A carbon abundance is as default in the fiducial models.Our calculations extend to the PDR regions, from which a [C II] 158 emission line mainly arises (Russell et al. 1980;Tielens & Hollenbach 1985).Therefore, the fiducial models require parameters important for the PDR regions: the cosmic radioto-X-ray background as a function of redshift (Ostriker & Ikeuchi 1983;Ikeuchi & Ostriker 1986;Vedel et al. 1994), ionization by cosmic rays background, the magnetic field of 30 µG, and the PDR covering fraction (C PDR ) of unity.As strengths of the cosmic rays background and magnetic fields are uncertain for high-redshift galaxies and even for local galaxies, a default value in CLOUDY (Indriolo et al. 2007) and a typical value in nearby spiral galaxies (Beck 2015) are applied, respectively.The code stops the calculations at the V-band dust extinction of 100 mag (Abel et al. 2005).
Input ionizing spectra varies according to galaxy populations.We used two simple stellar population (SSP) models: STARBURST99 (SB99; Leitherer et al. 1999Leitherer et al. , 2014) ) with the non-rotating standard Geneva tracks and the Binary Population and Spectral Synthesis code (BPASS, Eldridge et al. 2017) version 2.2.1 (Stanway & Eldridge 2018).Both models aim to reproduce young starbursts, but the BPASS models exhibit harder ionizing spectra than the SB99 models owing to massive star binaries and thus produce higher [O III]/Hβ and L([O III] 88 )/SFR ratios than the SB99 models.The BPASS models are applicable to high-redshift galaxies (e.g., Steidel et al. 2016), which agrees with an expectation that young low-metallicity galaxies, typically seen at high-redshift, host many massive binaries (Stanway 2020).For the local U/LIRGs (optFIR sample) and z ∼ 0 galaxies (optical sample), we applied a SB99 300 Myr constant starformation (cSF) model.For the local dwarfs (optFIR sample), which would be analogs of high-redshift galaxies, we applied a BPASS 300 Myr cSF model.These local dwarfs, U/LIRGs, and z ∼ 0 galaxies are assumed to have stel-  b The PDR covering fraction is not an input parameter of CLOUDY (see Section 5.2).
lar metallicity, Z * , equal to Z gas .On the other hand, FUV and optical spectra of galaxies at z ∼ 2 indicate an importance of low Z * /Z gas ratios, which would reflect the abundance pattern of core-collapse supernovae yields, as well as massive star binaries (Steidel et al. 2016;Trainor et al. 2016).Accordingly, we used BPASS 300 Myr cSF model with Z * /Z gas = 0.2 for the z ∼ 2 galaxies (optical sample).For the z > 6 galaxies (FIR sample), we assumed the identical ionizing spectra as for the z ∼ 2 galaxies.This assumption is supported by low Z * /Z gas ratios for galaxies at z > 2 (Cullen et al. 2019;Harikane et al. 2020a;Kashino et al. 2022).Section 5.1 will revisit our choices of the input ionizing spectra.The fiducial models have three free nebular parameters: the ionization parameter U , the hydrogen density n H , and the gas metallicity Z gas .The ionization parameter at the illuminated surface, U , is a dimensionless parameter that expresses relative amount of ionizing photons to gas density, defined as U = ṅγ /n H c, where ṅγ is the number flux of the ionizing photons and c is the speed of light.U ranged from log U = −4.0 to −0.5 at 0.5 steps.The hydrogen density at the illuminated surface ranged from log n H /cm −3 = 0 to 4.0 at 0.5 steps.Finally, the gas metallicity ranged Z gas /Z = 0.05, 0.2, 0.4, 1.0, and 2.0, where the solar metallicity value is Z = 0.02.For each nebular parameter set of (U , n H , Z gas ), our models computed Figure 2 illustrates a part of our fiducial model grids.The illustrated model uses the SB99 300 Myr cSF ionizing spectrum and shows grids at log U = −1.5, −2.5, and −3.5; Z gas /Z = 0.2, 0.4, and 1.0; and log n H /cm −3 = 0.5, 1.5, 2.5, and 3.5.Our model grids can reproduce most observational data points in both of the diagrams.While the model grids are relatively insensitive to n H in the BPT diagram, a wide range of n H is necessarily to explain the distributions in the FIR diagram.This is because the FIR fine-structure lines are more sensitive to the electron (hydrogen) densities than the optical lines due to their lower critical densities.Some previous studies fixed hydrogen densities at a typical values in their photoionization models (e.g., Strom et al. 2018;Sanders et al. 2020); this assumption would be reasonable for the BPT diagram, but not appropriate for the analysis of the FIR diagram.The more detailed model dependencies on the BPT and FIR diagrams are discussed in Kewley et al. (2013) and Harikane et al. (2020b), respectively.

Modeling galaxy distributions
We search model solutions of the nebular parameters (U , n H , and Z gas ) that can reproduce the distributions of galaxies in both of the diagrams.As seen in Section 2, galaxy populations with similar characteristics are distributed in groups on the diagrams.We aim to model not each galaxy but distribu-tions of galaxy populations with the fiducial photoionization models.
In our modeling, the distributions of the galaxy populations are represented by regions defined as follows.We defined the regions by hand for the optFIR (the local dwarfs and U/LIRGs) and FIR (the z > 6 galaxies) samples.These regions include most of the data points, but do not include some outliers.The regions for the optical sample are defined from the number-density contours.We used contours including 95% and 68% of the z ∼ 0 and 2 galaxies, respectively.We chose the relatively low-percentage contour for the z ∼ 2 galaxies because the number of data points are still small and their measurement errors are large, compared with the z ∼ 0 galaxies.These defined regions are shown in Figure 1 with the solid lines.
The searched nebular-parameter space is −4.0 < log U < −0.5, 0.0 < log n H /cm −3 < 4.0, and 0.1 < Z gas /Z < 2.0.The model grids of U , n H , and Z gas are interpolated into 20 steps to give finer grids.The parameter space of each sample is additionally restricted according to the literature.The hydrogen density is conservatively restricted to be 0.5 < log n H < 2.5 for the local dwarfs (Cormier et al. 2019;Spinoglio et al. 2022), 0.5 < log n H < 3.5 for the U/LIRGs (Graciá-Carpio et al. 2011;Inami et al. 2013;Herrera-Camus et al. 2016;Zhao et al. 2016), 0.5 < log n H < 2.5 for the z ∼ 0 galaxies, and 1.5 < log n H < 3.0 for the z ∼ 2 galaxies (Masters et al. 2014;Sanders et al. 2016).For the z > 6 galaxies, we applied the same constraint as for the z ∼ 2 galaxies because there are few measurements of n H for the z > 6 galaxies.The ranges of U and Z gas were restricted for the optical and FIR samples.As shown in Section 4.2, the nebular parameters for the optical sample are poorly constrained only from the optical measurements on the BPT diagram.To complement the lack of the observational constraints on the FIR diagram, we restricted U and Z gas for the optical sample to follow U -Z gas relations observed at z 2 (e.g., Pérez-Montero 2014).Whereas a local U -Z gas relation seems to hold at z ∼ 2 at Z gas /Z < 1/3 (Sanders et al. 2020), U values at z 1 seem to be 0.2-0.5 dex higher than local at Z gas /Z ∼ 1 (Nakajima & Ouchi 2014;Kashino et al. 2017).Therefore, we analytically expressed a U -Z gas relation of Pérez-Montero (2014) and employed it for the z ∼ 0 galaxies, We then modified the relation for the z ∼ 2 galaxies, where the upper limit is 0.3 dex higher than Equation 3 at Z = 0.Even if the increase of the upper limit becomes 0.2-0.5 dex, our results do not change.We used the same lower limit as Equation 3 because the increase in the lower limit makes it difficult to model the z ∼ 2 galaxies exhibiting relatively low [O III]/Hβ ratios.These modifications are consistent with a recent observational study of the U -Z gas relation at z ∼ 1-2 (Papovich et al. 2022).For the FIR sample, the z > 6 galaxies were assumed to have sub-solar metallicities, because high-redshift galaxies likely have low metallicities according to the redshift evolution of the massmetallicity relation (e.g., Zahid et al. 2013).and sub-solar metallicity (Z gas /Z < 1.0), as shown in the right panel.The range of the obtained nebular parameters is consistent with those measured for the DGS galaxies (Cormier et al. 2015(Cormier et al. , 2019)), although some of the dwarfs are measured to have very low metallicities of Z gas /Z < 0.1 (Madden et al. 2013).We note that there are model solutions around (log Z gas /Z , log U ) ∼ (0.2, −1.3), but we ignore these solutions because these are far from the main solutions.
The bottom panels show the results for the local U/LIRGs.The U/LIRGs exhibit high [N II]/Hα and low [O III]/Hβ in the BPT diagrams and low L([O III] 88 )/SFR in the FIR diagrams, indicating that these galaxies have nebular properties different from the dwarfs.In the right panel, the positive correlation between U and Z gas is tight despite a wide range of n H ; log U increases from −3.5 to −2.0 as log Z gas /Z increases from −0.25 to 0.3.This positive U -Z gas correlation is opposite to the negative correlations of the local dwarfs, H II regions, and galaxies at z ∼ 0 (Pérez-Montero 2014; Kashino & Inoue 2019), but these negative U -Z gas relations smoothly connected to the distributions of the U/LIRGs.Specifically, the nebular parameter space for the U/LIRGs do not share with that for the dwarfs.Our results of U and Z gas are roughly consistent with previous results using MIR and FIR lines (−3.2 log U −2.0 and 0.7 Z gas /Z 2.0; Graciá-Carpio et al. 2011;Inami et al. 2013;Díaz-Santos et al. 2017;Pereira-Santaella et al. 2017).However, there is a caution that emission from the U/LIRGs may be still contaminated by AGNs.As our models assume only ionizing spectrum of star formation, we would be careful to interpret the results for the U/LIRGs.
Comparisons with previous studies support that our fiducial photoionization models can reproduce the galaxy nebular parameters from the BPT and FIR diagrams.In the next sections, we apply these models to the optical and FIR samples to predict galaxy distributions on the FIR and BPT diagrams, respectively.4.2.FIR diagram at z ∼ 0-2 for the optical sample The top panels of Figure 4 show the results for the z ∼ 0 galaxies.We searched the nebular parameters that satisfy the region in the BPT diagram (left panel) and the restrictions of U and Z gas (right panel).The range of the obtained nebular parameters are −3.5 < log U < −2.5 and −0.45 < log Z/Z < 0.1.The middle panel shows a distribution of the z ∼ 0 galaxies in the FIR diagram predicted from the obtained nebular parameters.To check the accuracy of our predictions, we plotted local star-forming galaxies observed with ISO (Brauher et al. 2008) with the cross symbols.Their SFR was taken from De Looze et al. ( 2014); their optical measurements were taken from the same literature as for the optFIR sample; and AGNs, U/LIRGs, and resolved objects were not included.The distribution of the ISO galaxies in the FIR diagram agrees well with the prediction for the z ∼ 0 galaxies, while four ISO galaxies located on the region of the local U/LIRGs.This good agreement supports results and predictions of our photoionization modeling.
The bottom panels show the results for the z ∼ 2 galaxies.The obtained U and Z gas values are similar to those for the z ∼ 0 galaxies, but slightly lower Z gas values are dom-inant.The predicted distribution in the FIR diagram shows a weak positive correlation between L([O III] 88 )/SFR and L([C II] 158 )/SFR, similar to the local dwarfs and U/LIRGs.We note that, if there exist galaxies at z ∼ 2 having high hydrogen density of log n H /cm −3 > 3.0 (Lehnert et al. 2009;Hainline et al. 2009;Bian et al. 2010), they would be located at lower L([O III] 88 )/SFR and L([C II] 158 )/SFR than the predicted distributions.Compared with the z ∼ 0 galaxies, the z ∼ 2 galaxies are predicted to exhibit higher For both of the z ∼ 0 and 2 galaxies, the gray points in the U -Z gas diagrams (right panels) show the non-selected nebular parameters that explain the regions in the BPT diagram but do not satisfy the restrictions of the U and Z gas .These non-selected nebular parameters demonstrate that it is difficult to constrain the nebular parameters only from the BPT diagram, due to the parameter degeneracy between (low U , low Z gas , high n H ) and (high U , high Z gas , low n H ) (see Figure 2).
One may concern why almost the same U and Z gas values for the z ∼ 0 and 2 galaxies can reproduce the different dis- .Solutions of the photoionization modeling for the optical sample: the z ∼ 0 galaxies (Top) and the z ∼ 2 galaxies (Bottom).We searched the nebular parameters that satisfy the regions in the BPT diagrams (Left) and the constraints on the nebular parameters and plot the obtained solutions in the FIR (Middle) and U -Zgas (Right) diagrams.The cross symbols in the top panels indicate the ISO galaxies (Brauher et al. 2008).Middle and Right: The gray data points show the model parameters that satisfy the regions in the BPT diagram but do not the constraints on the nebular parameters (solid lines in the right panels).The obtained regions are illustrated with the dashed polygons.Other symbols are the same as Figure 3.
tributions in the BPT diagrams.These differences actually originate from the different assumptions in our photoionization models: the input ionizing spectra and hydrogen density.In Section 3, we assumed ionizing spectra of BPASS models with Z * /Z gas = 0.2 for the z ∼ 2 galaxies.These massive star binaries and low stellar metallicities make the spectrum harder than the ionizing spectra of SB99 models, which was used for the z ∼ 0 galaxies.The hard spectra and high hydrogen (i.e., electron) densities increase both [O III]/Hβ and [N II]/Hα in the BPT diagram (Kewley et al. 2013).In our photoionization models, thus, we explain the redshift evolution from z ∼ 0 to 2 with assumed harder ionizing spectra and higher hydrogen density at z ∼ 2, as well as slightly higher U and lower Z gas .

Nebular parameters and BPT diagram at z > 6 for the FIR sample
Figure 5 shows the results for the FIR sample (i.e., the z > 6 galaxies).We searched the nebular parameters that satisfy the region in the FIR diagram for the z > 6 galaxies.
The right panel of Figure 5 shows that U and Z gas have a negative correlation similar to those seen in galaxies at z ∼ 0 (e.g., Pérez-Montero 2014; Kashino & Inoue 2019).Interestingly, this negative correlation is obtained from only the region on the FIR diagram.The ionization parameter at low Z gas is comparable to those of Lyman-alpha emitters at z ∼ 2 (Nakajima & Ouchi 2014;Trainor et al. 2016;Kojima et al. 2017).The metallicity spans all the allowed range of −1.0 < log Z gas /Z < 0.0, implying that our analysis cannot find characteristic Z gas for the z > 6 galaxies.To clarify the metallicity dependence of the z > 6 galaxies, we divided the region into the three Z gas range as illustrated in Figure 6.The FIR diagram (bottom right panel) shows that an allowed range of Z gas depends on L([O III] 88 )/SFR; namely, Z gas is constrained to be high for galaxies with high L([O III] 88 )/SFR whereas not constrained for galaxies with low L([O III] 88 )/SFR.This dependence of the allowed Z gas range is consistent with predictions of the analytical models in Yang & Lidz (2020)   the predicted log [N II]/Hα varies from < −2.5 to −0.3.This predicted flat distribution is different from those of other low-redshift galaxy populations that follow the curved locus of star-forming galaxies (i.e., [O III]/Hβ decreasing with [N II]/Hα).Upcoming NIR observations with the James Webb Space Telescope (JWST) will test our predictions for the BPT diagram.Moreover, [N II]/Hα significantly depends on Z gas as shown in Figure 6.This dependence suggests a possibility that the [N II]/Hα ratio is still a metallicity tracer in the JWST era (N2 index; Pettini & Pagel 2004;Nagao et al. 2006), although the different nebular parameters at z > 6 may cause systematics in metallicity measurements (Bian et al. 2018).
We would like to stress that our results are consistent with conclusions of Harikane et al. (2020b), who used similar photoionization models.They concluded that high U , low C PDR , or both are necessarily to explain the distributions in the FIR diagram for galaxies at z > 6, especially for galaxies with high [O III] 88 /[C II] 158 ratios.Our results similarly predict high ionization parameters of log U −2.5 to −1.5 for the z > 6 galaxies at C PDR = 1.0.A major difference is that the region for the photoionization modeling in our work does not include data points with the highest [O III] 88 /[C II] 158 ratios, SXDF-NB1006-2 and MACS1148-JD1 (open red circles), because these measurements are based on SFR SED .Due to this difference, higher ionization parameters like log U = −1 are unnecessary in our modeling.

Bridging galaxy populations
Figure 7 summarizes the regions of the galaxy populations, observed and predicted from the modeling.At z ∼ 0, in all the panels, the distributions of the galaxy populations continuously shift from the local U/LIRGs through the normal z ∼ 0 galaxies to the local dwarfs.Although this trend has been already observed in the BPT diagram, the predicted distributions of the z ∼ 0 galaxies bridge a bi-modality between the dwarfs and U/LIRGs in the FIR and U -Z gas diagrams.The gas metallicity monotonically increases from the dwarfs to the U/LIRGs and the ionization parameter gradually changes along with metallicity.Combining the dwarfs and z ∼ 0 galaxies (that follow Equation 3) presents a negative U -Z gas correlation, which supports that the nebular parameters for the dwarfs, obtained from our photoionization modeling, agree with the local U -Z gas relation (Pérez-Montero 2014).In contrast, the U/LIRGs show a positive diagrams.The regions show the distributions of the dwarfs (gray), U/LIRGs (orange), z ∼ 0 (black), z ∼ 2 (blue), and z > 6 (red) galaxies.These regions are either drawn by hand on the basis of the galaxy distributions or obtained from the fiducial photoionization modeling.The gray dashed lines depict the criteria between star-forming galaxies and AGNs (Kauffmann et al. 2003;Kewley et al. 2001a).
correlation that smoothly connects with the z ∼ 0 galaxies at log Z gas /Z ∼ 0.1.
The z ∼ 0, 2, and > 6 galaxies show the redshift evolution of the galaxy distributions in all the diagrams.In the BPT diagram, the z ∼ 2 galaxies exhibit higher [O III]/Hβ ratios than the z ∼ 0 galaxies, and the z > 6 galaxies will exhibit higher [O III]/Hβ ratios than the z ∼ 2 galaxies on average.The z > 6 galaxies exhibit lower [N II]/Hα ratios than the z ∼ 2 galaxies; however, as shown in Figure 6, [N II]/Hα ratios strongly depend on Z gas and similar Z gas (log Z gas /Z = −0.6 to 0) will predict similar [N II]/Hα ratios at all the redshifts.Thus, higher-redshift galaxies exhibit higher [O III]/Hβ ratios at given Z gas (or [N II]/Hα).In the FIR diagram, higher-redshift galaxies show higher [O III] 88 /[C II] 158 ratios on average, which is consistent with previous findings at z ∼ 0 and > 6 (Inoue et al. 2016;Hashimoto et al. 2019;Harikane et al. 2020b).The evolution in the BPT and FIR diagrams would reflect the redshift evolution of the nebular parameters and the hardness of ionizing spectra.In the U -Z gas diagram, the z > 6 galaxies exhibit ∼ 0.5 dex (up to 1 dex) higher U than the z ∼ 0 and 2 galaxies at a given Z gas , in good agreement with the results of Harikane et al. (2020b).This higher U and assumed harder ionizing spectra at higher redshift give rise to the higher [O III]/Hβ and [O III] 88 /[C II] 158 ratios in our photoionization models.
The stellar mass is important in discussions of the redshift evolution.The average stellar mass of the z > 6 galaxies may be several times less than those of the z ∼ 0 and 2 galaxies (Section 2), but an increase in U caused by the low stellar mass would be less than 0.25 dex at given Z gas and n H according to local scaling relations (Kashino & Inoue 2019).Therefore, the predicted increases in U for the z > 6 galaxies are larger than expected from the difference of the stellar mass.
Specifically, the z > 6 galaxies and local dwarfs share a large part of their distributions on the diagrams.This similar-ity may support an argument that, regarding ISM properties, nearby dwarf galaxies are local analogues of high-redshift galaxies, despite differences of their masses.The ionizing spectra we assumed for the local dwarfs and z ∼ 2 galaxies are the BPASS 300 Myr cSF models with Z * /Z gas = 1 and 0.2, respectively; these hard spectra help to reproduce high [O III]/Hβ ratios observed in these galaxy populations.Both of the galaxy populations include galaxies with log [O III]/Hβ > 0.7 in their regions on the BPT diagram.The SB99 cSF models are often used for modeling local galaxies (e.g., Inami et al. 2013), including the DGS galaxies (Cormier et al. 2019), and they are actually suitable to explain low [O III]/Hβ ratios in the z ∼ 0 galaxies and U/LIRGs in this study.However, the SB99 models are unable to reproduce the high [O III]/Hβ ratios in the local dwarfs because the ionization parameters are constrained by the region on the FIR diagram.For this reason, the local dwarfs require harder input ionizing spectra than the SB99 models and prefer the BPASS models.The requirement becomes more stringent for the z ∼ 2 galaxies, which exhibit higher [N II]/Hα ratios than the local dwarfs.The BPASS cSF models at Z * /Z gas = 1 can explain only a half of the region of the z ∼ 2 galaxies on the BPT diagram with satisfying the U -Z gas relation.To explain the remaining galaxies having high [O III]/Hβ, harder ionizing spectra due to low Z * /Z gas are necessarily, which is consistent with results in previous emission-line studies (Steidel et al. 2016;Trainor et al. 2016;Shapley et al. 2019).We note that in our models the N/O ratio follows a single relation of Equation 2, although observed N/O ratios have a large dispersion (Pilyugin et al. 2012).Including the N/O dispersion in the models may help to explain high [O III]/Hβ ratios (Strom et al. 2018;Curti et al. 2022) by moving the model grids to high [N II]/Hα direction (see Figure 2).
A possible solution is low Z * /Z gas , which is adopted in our fiducial photoionization models.Low Z * /Z gas increases ionizing photon flux illuminating metal-enriched gas per SFR as well as hardness of ionizing spectra.These result in high L Hβ /SFR and high [O III] 88 /Hβ, respectively, and hence high L([O III] 88 )/SFR ratios (see Appendix B).Indeed, the fiducial models assume Z * /Z gas = 0.2 and successfully explain the z > 6 galaxies with the high L([O III] 88 )/SFR ratios within line measurement errors.Although low Z * /Z gas at high redshift was inferred from UV-to-optical emissionline studies (Cullen et al. 2019;Harikane et al. 2020a;Kashino et al. 2022), our photoionization models highlight the necessity of low Z * /Z gas ratios at z > 6 from FIR L([O III] 88 )/SFR ratios.
Bursty or increasing star-formation history is an alternative solution to explain the high L([O III] 88 )/SFR ratios.For a part of the z > 6 galaxies, short star-formation ages of 5 Myr are inferred from the SED fittings (Inoue et al. 2016;Hashimoto et al. 2018;Tamura et al. 2019;Hashimoto et al. 2019).Such a bursty or increasing star-formation history leads to underestimating SFR UV+IR and thus overestimating L([O III] 88 )/SFR UV+IR , as illustrated in Figure A1 where all the SFR SED are higher than SFR UV+IR estimated by Carniani et al. (2020).Moreover, this solution is in line with theories that bursty star formation can explain high [O III] 88 /[C II] 158 ratios of the z > 6 galaxies (Arata et al. 2020;Vallini et al. 2021).The bursty star-formation history is compatible with low Z * /Z gas ; therefore, the high L([O III] 88 )/SFR ratios support both of these two properties regarding the ionizing spectrum for the z > 6 galaxies.
Finally, we mention possible explanations of high L([O III] 88 )/SFR ratios other than low Z * /Z gas and bursty star-formation history.
Photoionization models with log n H /cm −3 < 1.0 can also explain the high L([O III] 88 )/SFR ratios, as shown in Yang & Lidz (2020).However, this models require high U (log U > −2.0) and Z gas /Z ∼ 1 as well as low n H , which are totally different from the nebular parameters of the z ∼ 0 and 2 galaxies and would contrast with predictions that high-redshift galaxies exhibit higher electron densities than galaxies at z ∼ 0 (e.g., Sanders et al. 2016;Kaasinen et al. 2017).We note that highly-ionized diffuse gas (having high U and low n H ) may contribute [O III] 88 emission as observed in the local universe (Kawada et al. 2011;Lebouteiller et al. 2012;Polles et al. 2019).High dust temperature is another possibility to be considered.The dust temperature is highly uncertain for the z > 6 galaxies because only one or two ALMA measurements are available in the most cases.Higher dust temperature yields higher L FIR and higher SFR UV+IR at given dust continuum fluxes.If all the z > 6 galaxies with the high L([O III] 88 )/SFR have dust temperature of 60 K (Harikane et al. 2020b, Ono et al. in prep.),L([O III] 88 )/SFR becomes small enough to be explained by photoionization models.
One may expect that a top-heavy IMF can explain the high L([O III] 88 )/SFR ratios of the z > 6 galaxies.To test this possibility, we adopted BPASS models using double powerlaw IMF with a massive slope of −2.0 at 0.5-100 M to construct a photoionization model.This top-heavy model failed to reproduce the high L([O III] 88 )/SFR ratios because predicted [O III] 88 /Hβ is similar to the Chabrier IMF models.

Parameter dependence of photoionization models
Although we assumed several parameters in the fiducial photoionization models (Section 3), the assumed parameters for the z > 6 galaxies are the same as for the z ∼ 2 galaxies because most of them have not been constrained observationally.In this section, we change parameters assumed for the z > 6 galaxies in the fiducial models to assess uncertainties of our model predictions.This section discusses the BPT diagram, but all the BPT, FIR, and U -Z gas diagrams for the changed models are shown in Appendix C.
SSP model and stellar age-The fiducial photoionization models adopt BPASS 300 Myr cSF model.We changed the input ionizing spectrum to SB99 300 Myr cSF and BPASS 10 Myr cSF models.The panel (a) of Figure 8 compares predicted BPT diagrams for the z > 6 galaxies.We note that the panel (a) only shows a range of log[N II]/Hα > −1.8 for the SB99 model because the minimum Z * in SB99 of 0.05Z prevents modeling at Z gas /Z < 0.25 using the models with Z * /Z gas = 0.2.Although there are slight differences at high [N II]/Hα (i.e, high metallicity), all the models predict the same [O III]/Hβ values at low [N II]/Hα (low metallicity).The differences of SSP models and stellar ages do not change the predicted BPT diagram for the z > 6 galaxies.
Stellar-to-gaseous metallicity ratio-The fiducial models for the z > 6 galaxies assume Z * /Z gas = 0.2 by following previous studies on galaxies at z 2. We changed Z * /Z gas from 0.2 to 1.0.In the panel (b) of Figure 8, the Z * /Z gas ratios do not affect predicted distributions at low [N II]/Hα (i.e., low metallicity), while they slightly affect the distributions at log [N II]/Hα > −1.5 (high metallicity).These results are similar to the case when SSP models and stellar ages change.
C  The dashed gray lines depict the criteria between star-forming galaxies and AGNs (Kauffmann et al. 2003;Kewley et al. 2001a Cosmic ray and magnetic field-The background cosmicray intensity flux and the magnetic field amplitudes in the ISM is quite unclear at high redshift.As an experiment, we changed the two parameters to 1-dex higher and lower than the fiducial values.The panel (d) of Figure 8 illustrates that the two parameters change the predicted distributions by < 0.2 dex in [O III]/Hβ at most, which is much lower than the changes generated by C/O and C PDR .Therefore, the effects of the cosmic rays and magnetic fields are small despite their unclear amplitudes at high redshifts.
A prediction of our fiducial models, log [O III]/Hβ 0.7 for the z > 6 galaxies, holds even when many assumed parameters change, while low C/O and C PDR values make a dispersion of [O III]/Hβ larger.These predicted BPT diagrams will be tested by the upcoming JWST observations.JWST GTO and GO cycle 1 programs plan to perform follow-up observations of the z > 6 galaxies with the NIRSPEC integrated field units.In the NIRSPEC observations, the dispersion of [O III]/Hβ may be helpful to distinguish whether the ISM of the z > 6 galaxies have low C/O and C PDR .The MIRI MIR observations targeting Hα and [N II] emission lines will also contribute to estimate the metallicity of the galaxies.

Caveats
Aperture differences in optical and FIR observations-The optical and FIR line fluxes used in this work were measured with different aperture sizes.Optical spectroscopy usually uses slits or fibers, which can observe only a part of a nearby galaxy, whereas FIR instruments have large apertures to observe a whole galaxy.These aperture differences may cause systematics in estimates of the nebular parameters from optical and FIR observations, even though our models assumed that the same nebular parameters can simultaneously explain galaxy distributions on the BPT and FIR diagrams.We attempted to reduce the aperture differences by using the integrated spectroscopic observations as much as possible.In the near future, the z > 6 galaxies will be observed with JWST integral field units (IFU) that covers the entire of them within its field of view.Therefore, observing local galaxies with optical IFU will be important for a fair comparison with the z > 6 galaxies.
Geometrical variation of dust attenuation-This work ignores effects of the dust attenuation because FIR lines are not attenuated by dust and the BPT diagram takes ratios of lines with close wavelengths.However, the geometrical variation of dust attenuation may bias optical and FIR radiation sources within galaxies.Although we can observe FIR emission lines from a whole galaxy, the optical line fluxes may be dominated by emission from less dusty regions.Thus, this spatial variation of attenuation may break our assumption that we can reproduce both optical and FIR emission-line ratios with the same nebular parameters.Some studies report spatial offsets between UV and FIR line/continuum detections at z > 5 (Carniani et al. 2018;Bowler et al. 2022).These effects should be considered in the future analyses.
Contributions from DIG and AGN-Our models ignore the DIG, which contributes low-ionization emission lines such as [N II] and [C II] lines in local galaxies (e.g., Martin 1997;Kaufman et al. 2006).The DIG increases the [N II]/Hα ratio by a few times 0.1 dex and systematically biases metallicity estimates (Zhang et al. 2017), but its effect to [N II]/Hα (and to the BPT diagram in this work) is relatively smaller than those to [S II]/Hα and [O II]/Hβ (Sanders et al. 2017).At z ∼ 2, and possibly at higher redshift, DIG contributions to low ionization lines would be negligible because of high SFR surface densities of galaxies at the redshift (Sanders et al. 2017;Shapley et al. 2019).Our models also ignore the hidden AGN contributions.Although we cannot reject a possibility that AGNs affects the line ratios of the z > 6 galaxies, currently there are no clear signatures of AGNs in the z > 6 galaxies.
Geometry of the ionized and photodissociated regions-Most photoionization models, including ours, assume idealized nebular geometries such as spherical or plane-parallel geometry.However, the actual H II regions and galaxies have complex gaseous geometry, stellar distributions, and temperature/density structures.This discrepancy between models and observations may complicate interpretations of observed emission-line ratios of galaxies.More sophisticated photoionization modeling will be necessarily in the era of multiwavelength, IFU observations (e.g., Jin et al. 2022).U/LIRGs), the optical sample (z ∼ 0 and 2 galaxies), and the FIR sample (z > 6 galaxies).The constructed photoionization models have three free nebular parameters of the ionization parameter U , hydrogen density n H , and gaseous metallicity Z gas .We have defined the regions representing distributions of the galaxy populations on the BPT and FIR diagrams and have searched nebular parameters that satisfy the regions.
Our fiducial photoionization models successfully reproduce the nebular parameters of the local dwarfs and U/LIRGs consistent with results of previous studies.Then we have applied the photoionization models to the optical and FIR samples.For the z ∼ 0 galaxies, the predicted distribution on the FIR diagram is consistent with ISO observations.The z ∼ 2 galaxies are predicted to exhibit higher [O III] 88 /[C II] 158 ratios than the z ∼ 0 galaxies on average.For the z > 6 galaxies, the predicted nebular parameters have a negative U -Z gas correlation and U and Z gas values are comparable to those of Lyman-alpha emitters at z ∼ 2. The distribution on the BPT diagram for the z > 6 galaxies is relatively flat; log [O III]/Hβ = 0.5-0.8 while [N II]/Hα strongly depends on Z gas .
Comparing the galaxy populations illustrates continuous distributions from low-to high-mass galaxies at z ∼ 0 and continuous shifts from z ∼ 0 to z > 6, on all the diagrams (Figure 7).The z ∼ 0 galaxies bridge the distributions of the local dwarfs and U/LIRGs on the FIR and U -Z gas diagrams, as well as on the BPT diagram.These diagrams show continuous transitions of the stellar and ISM properties from dwarfs through normal star-forming galaxies to U/LIRGs.The z > 6 galaxies have higher U than the z ∼ 0 and 2 galaxies at given Z gas .Thanks to high U and assumed hard ionizing spectrum, the z > 6 galaxies have higher [O III]/Hβ on the BPT diagram and higher [O III] 88 /[C II] 158 on the FIR diagram than the z ∼ 0 and 2 galaxies.These continuous shifts on the diagrams demonstrate the redshift evolution of the stellar and ISM properties of galaxies from the epoch of reionization to present.In addition, the z > 6 galaxies share a large part of the distributions with the local dwarfs, indicating a similarity of ISM properties between them.
We find that some of the z > 6 galaxies exhibit high L([O III] 88 )/SFR ratios, which cannot be reproduced by BPASS 300 Myr constant star-formation models with Z * /Z gas = 1 and 1.5 < n H /cm −3 < 3.0.To reproduce such high L([O III] 88 )/SFR ratios, our photoionization models require: 1) low stellar-to-gaseous metallicity ratios of Z * /Z gas = 0.2, consistent with results of UV-to-optical emission line studies; or 2) bursty or increasing star-formation history, in agreement with the SED-fitting results of some z > 6 galaxies.We mentioned other possibilities of low hydrogen density and high dust temperature.Similarly, high [O III]/Hβ ratios of the local dwarfs and z ∼ 2 galaxies require BPASS models with Z * /Z gas = 1 and 0.2, respectively, rather than SB99 non-rotating models.
Upcoming JWST plans to detect rest-frame optical emission lines of the z > 6 galaxies.These new NIR observations can test our model predictions, including the distribu-tions on the BPT diagrams and the nebular parameters of the z > 6 galaxies, and improve our understanding of the early universe.A. THE z > 6 GALAXIES USED IN THIS WORK Figure A1 shows names of the z > 6 galaxies for convenience of the reader.As nine UV-selected galaxies have two measurements given by Harikane et al. (2020b) and Carniani et al. (2020), their two data points are connected with the dashed lines.where x denotes the tracer of SFR, L x the tracer luminosity, C x the conversion factor, t the age.Although most of the references cited in this work (De Looze et al. 2014;Madden et al. 2013;Howell et al. 2010;Herrera-Camus et al. 2018a;Carniani et al. 2020;Harikane et al. 2020b) derived SFR using formulae of Kennicutt (1998), Hao et al. (2011)  and Wilkins et al. (2019).This section revisits the SFR conversion factors to make them useful for this work.
The production rate of ionizing photons is proportional to the starformation rate, (B5) The bottom panel of Figure B1 shows the dependence of C FUV on the parameters.Equation B4 is available to convert line-to-Hβ ratios output by CLOUDY to line-to-SFR ratios; However, this SFR is inconsistent with those estimated from observations used in this work because C Q(H 0 ) depends on SSP, Z * , and t.In the observations, SFR is often derived from the UV (and FIR) luminosity, following where C KE12 FUV is the conversion factor in Kennicutt & Evans (2012) which connect the Hβ luminosity in CLOUDY to SFR that can be compared with the observations.Figure B2 illustrates dependence of C obs Hβ on stellar age and metallicity.Hβ luminosity at given SFR increases as the stellar age and metallicity decrease.Regarding the z > 6 galaxies, while the observed SFR UV+IR can be compared with SFR obs of Equation B7, the observed SFR SED can be compared with SFR of Equation B4 because the SED fitting already took stellar age and metallicity into account.

C. RESULTS OF NON-FIDUCIAL PHOTOIONIZATION MODELING
In Section 5.2, we changed assumptions of the photoionization models from the fiducial ones.Figure C1 shows panels that are the same as Figure 5 but for photoionization models with different assumptions.

Figure 1 .
Figure 1.BPT diagram (Left) and FIR diagram(Right).The red circles show the z > 6 galaxies; the filled circles show the measurements mainly taken fromHarikane et al. (2020b), while the hashed circles fromCarniani et al. (2020).The open red circles are taken fromHarikane et al. (2020b), but their SFR is measured with the SED fitting.The gray triangles show the local dwarfs(Cormier et al. 2019) and the orange symbols show the local U/LIRGs including the GOALS galaxies (square, Díaz-Santos et al. 2017) and the SHINING galaxies (cross,Herrera- Camus et al. 2018a).The solid polygons depict the regions representing the distributions of the galaxy populations, which are used in our photoionization modeling.Left: The black and blue contours illustrate the distributions of the z ∼ 0 (SDSS) and z ∼ 2(Strom et al. 2017;Shivaei et al. 2018) galaxies, respectively, which include 68, 95, and 99.5% of the z ∼ 0 galaxies and 68 and 90% of the z ∼ 2 galaxies.The dashed gray lines depict the criteria between star-forming galaxies and AGNs(Kauffmann et al. 2003;Kewley et al. 2001a).Right: The horizontal and vertical error bars indicate the [C II]158 and [O III]88 measurement errors, respectively, and the diagonal error bars does the SFR errors.For local galaxies, the black edges mean the data points with optical measurements, while the gray edges mean ones without them.

Figure 2 .
Figure2.A part of fiducial photoionization model grids on the BPT (Left) and FIR (Right) diagrams.The input ionizing spectrum is the SB99 300 Myr constant star-formation model.Each model grid indicates log U = −1.5, −2.5, and −3.5 and Zgas/Z = 0.2, 0.4, and 1.0.The colors of the four grids indicate log nH/cm −3 = 0.5, 1.5, 2.5, and 3.5 from purple (dark) to yellow (light).The symbols are the same as Figure1, but we removed the symbols with gray edges, removed the error bars, and lightened colors.The gray and cyan contours in the left panel show the 95 and 68 percentile contours for the z ∼ 0 and 2 galaxies, respectively.The model grids can explain most of the data points well with the three parameters, U , Zgas, and nH.
[O III]λ5007, [N II]λ6548, Hα, [O III] 88 , and [C II] 158 line intensities relative to Hβ line.As the FIR diagram requires SFR in denominators of the two axes, we computed model line-to-SFR ratios by converting Hβ line intensity to SFR with Equation B7.Section B discusses the derivation of the conversion factor in detail.
Figure3shows the obtained nebular-parameter solutions for the optFIR sample that satisfy the regions in the BPT (left) and FIR (middle) diagrams.The right panels show the obtained nebular parameters (U , Z gas , n H ). The optFIR sample helps to check whether our photoionization model reproduces proper nebular parameters of local galaxies from the BPT and FIR diagrams.The top panels show the results for the local dwarfs.The dwarfs exhibit low [N II]/Hα and high [O III]/Hβ in the BPT diagram and high L([O III] 88 )/SFR in the FIR diagram, suggesting high ionization and low metallicity environments.In fact, our photoionization model gives −3.0 < log U < −1.5 and sub-solar metallicity (Z gas /Z < 1.0), as shown in the right panel.The range of the obtained nebular parameters is consistent with those measured for the DGS galaxies(Cormier et al. 2015(Cormier et al. , 2019)), although some of the dwarfs are measured to have very low metallicities of Z gas /Z < 0.1(Madden et al. 2013).We note that there are model solutions around (log Z gas /Z , log U ) ∼ (0.2, −1.3), but we ignore these solutions because these are far from the main solutions.The bottom panels show the results for the local U/LIRGs.The U/LIRGs exhibit high [N II]/Hα and low [O III]/Hβ in the BPT diagrams and low L([O III] 88 )/SFR in the FIR diagrams, indicating that these galaxies have nebular properties different from the dwarfs.In the right panel, the positive correlation between U and Z gas is tight despite a wide range of n H ; log U increases from −3.5 to −2.0 as log Z gas /Z increases from −0.25 to 0.3.This positive U -Z gas correlation is opposite to the negative correlations of the local dwarfs, H II regions, and galaxies at z ∼ 0 (Pérez-Montero 2014; Kashino & Inoue 2019), but these negative U -Z gas relations smoothly connected to the distributions of the U/LIRGs.Specifically, the nebular parameter space for the U/LIRGs do not share with that for the dwarfs.Our results of U and Z gas are roughly consistent with previous results using MIR and FIR lines (−3.2 log U −2.0 and 0.7 Z gas /Z 2.0;Graciá-Carpio et al. 2011;Inami et al. 2013; Díaz-Santos  et al. 2017;Pereira-Santaella et al. 2017).However, there is a caution that emission from the U/LIRGs may be still contaminated by AGNs.As our models assume only ionizing spectrum of star formation, we would be careful to interpret the results for the U/LIRGs.Comparisons with previous studies support that our fiducial photoionization models can reproduce the galaxy nebular parameters from the BPT and FIR diagrams.In the next sections, we apply these models to the optical and FIR samples to predict galaxy distributions on the FIR and BPT diagrams, respectively.

Figure 3 .
Figure 3. Solutions of the photoionization modeling for the optFIR sample: the dwarfs (Top) and the U/LIRGs (Bottom).We searched the nebular parameters that satisfy the regions in the BPT (Left) and FIR (Middle) diagrams and plot the obtained parameters in the U -Zgas diagram (Right).Left and Middle: The obtained solutions are plotted with being color-coded by nH.The symbols and contours are the same as Figure 2, but the symbols of the target galaxies are highlighted.Right: The colored circles show the obtained solutions, of which outer and inner colors show the maximum and minimum nH, respectively, at a given U and Zgas.The gray data points are model solutions that satisfy one of the regions in the BPT or FIR diagrams, but not the other.The obtained nebular-parameter regions are illustrated with the dashed polygons.We note that model solutions for the dwarfs around (log Zgas/Z , log U ) ∼ (0.2, −1.3) are ignored because these are far from the main solutions

Figure 4
Figure4.Solutions of the photoionization modeling for the optical sample: the z ∼ 0 galaxies (Top) and the z ∼ 2 galaxies (Bottom).We searched the nebular parameters that satisfy the regions in the BPT diagrams (Left) and the constraints on the nebular parameters and plot the obtained solutions in the FIR (Middle) and U -Zgas (Right) diagrams.The cross symbols in the top panels indicate the ISO galaxies(Brauher et al. 2008).Middle and Right: The gray data points show the model parameters that satisfy the regions in the BPT diagram but do not the constraints on the nebular parameters (solid lines in the right panels).The obtained regions are illustrated with the dashed polygons.Other symbols are the same as Figure3.
. As proposed inMoriwaki et al. (2018),Jones et al. (2020), andYang &Lidz (2020), [O III] 52 µm and [O III] λ5007 lines are important to resolve a degeneracy in U , n H , and Z gas at z > 6.The predicted distribution of the z > 6 galaxies in the BPT diagram is shown in the left panel of Figure5.The predicted log [O III]/Hβ is almost constant at 0.5-0.8whereas

Figure 5 .
Figure 5. Solutions of the photoionization modeling for the z > 6 galaxies (the FIR sample).We searched the nebular parameters that satisfy the regions in the FIR diagram (Middle) with Z < 1.0Z and plot the obtained solutions in the BPT (Left) and U -Zgas (Right) diagrams.Left and Right: The gray data points show model parameters that satisfy the region in the FIR diagram, but at Zgas/Z ≥ 1.The obtained regions are illustrated with the dashed polygons.The symbols are the same as Figure 3, but the symbols of the z > 6 galaxies are highlighted.

Figure 7 .
Figure7.Graphical summary of the regions where the galaxy populations are distributed on the BPT (Left), FIR (Middle), and U -Zgas (Right) diagrams.The regions show the distributions of the dwarfs (gray), U/LIRGs (orange), z ∼ 0 (black), z ∼ 2 (blue), and z > 6 (red) galaxies.These regions are either drawn by hand on the basis of the galaxy distributions or obtained from the fiducial photoionization modeling.The gray dashed lines depict the criteria between star-forming galaxies and AGNs(Kauffmann et al. 2003;Kewley et al. 2001a).
Implications of high [O III]/Hβ and L([O III] 88 )/SFR In the fiducial photoionization models, we have assumed the ionizing spectra depending on the galaxy populations (Section 3).This section discusses the necessity of the input ionizing spectra to represent line ratios, especially for [O III]/Hβ and L([O III] 88 )/SFR.

Figure 8 .
Figure 8. Photoionization model solutions on the BPT diagrams obtained for the z > 6 galaxies under the different assumptions (i.e., non-fiducial models).The obtained solutions are displayed with a line and bars, representing means and widths of the distribution at several [N II]/Hα values, respectively.The fiducial model, which is the same as the left panel of Figure 5, is shown in black in all the panels.(a) Different SSP models and ages of the input ionization spectra: BPASS 300 Myr, SB99 300 Myr, and BPASS 10 Myr constant star-formation models.The SB99 model only shows log[N II]/Hα > −1.8 (see text).(b) Different stellar-to-gaseous metallicity ratios: Z * /Zgas = 1.0 and 0.2.(c) Different C/O abundance ratio and PDR covering fraction: R C/O CPDR = 1.0, 0.3, 0.1, and 0.03.(d) Different cosmic-ray background intensity and strength of magnetic field: ×10, ×1, and ×0.1, and log B/G =−3.5, −4.5, and −5.5.The dashed gray lines depict the criteria between star-forming galaxies and AGNs(Kauffmann et al. 2003;Kewley et al. 2001a).Only low C/O and CPDR increase a dispersion of [O III]/Hβ while other assumptions do not affect the model predictions.
photoionization modeling to distributions of galaxy populations from z ∼ 0 to z > 6 on the BPT ([O III]λ5007/Hβ-[N II]λ6585/Hα) and FIR (L([O III] 88 )/SFR-L([C II] 158 )/SFR) diagrams.The galaxy samples are divided into three according to their available measurements: the optFIR sample (local dwarfs and Figure A1.Same as the FIR diagram in Figure1, but names of the z > 6 galaxies are denoted near the data points.The black dashed lines connects two measurements for the same objects listed inHarikane et al. (2020b) andCarniani et al. (2020).
B. SFR CONVERSION FACTORS DEPENDENT ON SSP MODELS, METALLICITY AND STELLAR AGERelations between SFR and its tracers depend on the SSP models, star-formation history, metallicity, and IMF.Even under constant star-formation history as assumed in this work, the relation depends on stellar population age; namely, L x /SFR = C x (SSP, Z * , t, IMF), (B1)

Figure B1 .
Figure B1.SFR conversion factors as a function of the stellar metallicity and stellar age.The solid and dashed lines show the BPASS and SB99 SSP models.The line colors depend on metallicity.Top: conversion factor for Q(H 0 ).Bottom: conversion factor for LFUV.The gray dotted line indicates the conversion factor in Kennicutt & Evans (2012).

Figure B2 .
Figure B2.SFR conversion factor for L Hβ (Equation B7) as a function of the stellar metallicity and stellar age.The solid and dashed lines show the BPASS and SB99 SSP models.The gray dotted line indicates the conversion factor in Kennicutt & Evans (2012), which is obtained by dividing C KE12Hα by the Hβ/Hα ratio of 2.87 under Case B recombination.This factor C obs Hβ converts Hβ intensity in models to SFR that can be compared with observations.

Figure C1 .
Figure C1.Same as Figure 5, but for photoionization models with assumptions that differs from the fiducial ones.The changed assumptions are written in the top left of the left panels.

Table 1 .
Assumed parameters of the fiducial photoionization models.
/O abundance ratio and PDR covering fraction-These two parameters change [C II] 158 line intensity at a given metallicity.The C/O ratio directly scales [O III] 88 /[C II] 158 ratio.We define R C/O ≡ (C/O)/(C/O) ISM , where ).Only low C/O and CPDR increase a dispersion of [O III]/Hβ while other assumptions do not affect the model predictions.C PDR to 0.3, 0.1, and 0.03 by reducing [C II] intensity at those fractions.In the panel (c) of Figure 8, the predicted [O III]/Hβ value decreases as R C/O C PDR decreases.Figure C1 in the appendix shows that the decrease in [O III]/Hβ is attributed to a decrease in U at a given Z gas .Although the fiducial models explain high [O III] 88 /[C II] 158 ratios of the z > 6 galaxies by increasing U (Section 4.4), low C PDR and R C/O can explain the high [O III] 88 /[C II] 158 ratios with keeping U low (Harikane et al. 2020b).However, Figure C1 also shows that it is more difficult for low R C/O C PDR models to explain galaxies with high L([O III] 88 )/SFR than for the fiducial models.For this reason, not all the z > 6 galaxies would exhibit low C/O and C PDR .If R C/O C PDR values vary among the z > 6 galaxies, [O III]/Hβ values have a dispersion of more than 0.5 dex, larger than a prediction of the fiducial models.