Reconciling escape fractions and observed line emission in Lyman continuum leaking galaxies

Finding and elucidating the properties of Lyman continuum (LyC) emitting galaxies is an important step for our understanding of cosmic reionization. Although the z - 0.3-0.4 LyC emitters found recently show strong optical emission lines, no consistent quantitative photoionization model taking into account the escape of ionizing photons and inhomogenous interstellar medium (ISM) geometry of these galaxies has yet been constructed. We hence construct one- and two-zone photoionization models accounting for the observed LyC escape, which we compare to the observed emission line measurements. We find that one-zone density-bounded photoionization models cannot reproduce the emission lines of the LyC leakers since they systematically underpredict the lines of species of low ionization potential, as [OI] and [SII]. Introducing a two-zone model, with differing ionization parameter and a variable covering fraction and where one of the zones is density-bounded, we show that the observed emission line ratios and escape fractions of the LyC emitters are well reproduced. The [OI] excess, which is observed in some LyC leakers, can be naturally explained in this model, e.g., by emission from low ionization and low filling-factor gas. LyC emitters with high escape fraction (fesc>38%) are deficient both in [OI]6300A and in [SII]6716,6731A. We also confirm that a [SII] deficiency can be used to select LyC emitter candidates. Finally, we find indications for a possible dichotomy in term of escape mechanisms for LyC photons between galaxies with relatively low (fesc<10%) and higher escape fractions. We conclude that two-zone photoionization models are needed to and capable of explaining the observed emission line properties of z - 0.3-0.4 LyC emitters. These models provide a first step towards the use of optical emission lines and their ratios as quantitative diagnostics of LyC escape from galaxies.


Introduction
Significant progress has been made in recent years on searches for and studies of star-forming galaxies, which emit copious amount of Lyman continuum (LyC) radiation and whose ionizing radiation can escape the interstellar medium (ISM) of these galaxies to ionize their surroundings, and thus contribute to intergalactic ionizing radiation. If sufficient population of such galaxies exists in the early Universe, they could explain cosmic reionization, which is known to have started shortly after the Big Bang and ended ∼ 1 Gyr later, at z ∼ 6.
After intense efforts to search for LyC emitters or LyCleaking galaxies ("leakers" in short), several groups have identified such galaxies from low redshift up to z ∼ 4 (e.g., Leitherer et al. 2016;Izotov et al. 2016aIzotov et al. ,b, 2018aVanzella et al. 2016Vanzella et al. , 2018Vanzella et al. , 2020de Barros et al. 2016;Steidel et al. 2018;Fletcher et al. 2019;Bian et al. 2017;Wang et al. 2019, etc.) with LyC escape fractions, f esc , ranging from a few percent up to ∼ 72 % in few cases (Izotov et al. 2018a;Vanzella et al. 2020). Although different methods and criteria have been used to search for these galaxies, some of the distinct observational features shared by many LyC emitters are strong Lyα emission Article number, page 1 of 16 arXiv:2009.09882v1 [astro-ph.GA] 21 Sep 2020 A&A proofs: manuscript no. leakers_EL_0 with a double-peaked Lyα line profiles (see Verhamme et al. 2015Verhamme et al. , 2017Vanzella et al. 2020) and intense optical emission lines with a high [O iii]/[O ii] 1 ratio (Izotov et al. 2018b;de Barros et al. 2016;Vanzella et al. 2020). In fact, after Jaskot & Oey (2013) and Nakajima & Ouchi (2014) suggested that high values of [O iii]/[O ii] could pinpoint density-bounded, i.e., LyC emitting galaxies, this simple prediction was put to test by Izotov and collaborators, who selected compact star-forming galaxies with intense emission lines and [O iii]/[O ii]> 5, and who were thus able to find 11 LyC emitters at z ∼ 0.3 − 0.4 with 100% success rate using the COS spectrograph on-board the Hubble Space Telescope (HST). The galaxies found in this way have properties comparable to the numerous galaxies at z > 6. They are strong producers of ionizing photons, show significant LyC escape ( f esc ∼ 2 − 72 %), and are therefore good analogs of the sources of cosmic reionization (cf. Schaerer et al. 2016 Izotov et al. 2018b;Naidu et al. 2018;Bassett et al. 2019;Nakajima et al. 2020, and discussions therein). Similarly, Wang et al. (2019) recently selected some galaxies for being [S ii]-deficient, which could indicate a density-bounded ISM, and were able to detect LyC emission from two low-z galaxies with fairly low excitation ([O iii]/[O ii] ∼ 1.3 − 1.5). However, they could not find a correlation between the [S ii]-deficiency and the LyC escape fraction of the known leakers. These findings are not entirely surprising, since it is well known that the [O iii]/[O ii] ratio depends also on several other factors, including, in particular, ionization parameter and metallicity. Whether density-bounded models are applicable to known LyC emitters remains therefore to be examined quantitatively.
On the other hand, the detailed analysis of UV absorption lines (including the Lyman series lines of hydrogen and other low ionization lines) of the z ∼ 0.3−0.4 leakers of Izotov and collaborators undertaken by Gazagnes et al. (2018) and Chisholm et al. (2018), has shown the presence of high column density gas, optically thick to LyC radiation, co-existing with lines of sight, which allow the escape of UV ionizing and non-ionizing radiation from these galaxies. This implies that the ISM of these galaxies cannot be described by simple one-zone photoionization models, and that a more complex ISM structure should be accounted for to consistently describe these objects. Therefore, testing one-zone photoionization models, tailored to the properties (metallicities etc.) of the LyC emitters, and exploring more complex structures to consistently describe their emission line spectra are needed. This is one of the main goals of our study.
Specifically, the first basic question posed is if and how well simple, i.e., single, 1D photoionization models can reproduce the observed emission line properties of the LyC emitters. Comparing grids of photoionization models with observational diagrams of ∼ 800 low to moderate metallicity star-forming galaxies, Stasińska et al. (2015) concluded that the strong lines indicate no strong LyC leakage (i.e., f esc < ∼ 10%) in galaxies with high [O iii]/[O ii], since the observed emission line ratios are not compatible with density-bounded H ii regions and, therefore, with a significant ionizing photon escape. The same authors also stressed the importance of considering high and low excitation lines. Furthermore, they already conjectured that some leakage might occur because of a nebular covering factor smaller than unity, and that [O i] emission could come from optically thick clumps embedded in a density-bounded medium (Stasińska et al. 2015). Related to this, Plat et al. (2019)

have recently noted that low-z LyC emitters show a higher [O i]/[O iii] ratio than their ionization-bounded photoionization models, in contrast to expectations from density-bounded models, where [O i]/[O iii]
should be reduced. To solve this discrepancy they invoke a contribution of hard radiation from AGN or radiative shocks. However, no consistent model, allowing for the observed leakage and escape fraction, including the geometrical constraints derived from the UV lines (cf. above), and explaining the intensities of the major emission lines has yet been done.
To progress further on these issues, we here compute both one-and two-zone photoionization models with metallicities tailored to those of the observed LyC emitters, to quantitatively examine how such models compare to the various detailed observational constraints, including measured LyC escape fractions. First, we re-examine the most important observed emission line ratios of the z ∼ 0.3 − 0.4 leakers, and compare them to a large sample of star-forming galaxies of comparable metallicity, serving as a meaningful reference sample, which is presumably dominated by non-leaking galaxies. We then compare the data to standard one-zone photoionization models to see the successes and failures of these models. Finally, and since one-zone models cannot properly reproduce the data, we are led to construct two-zone models which show a lot of interesting features and provide useful insight on the ISM and emission line properties of the LyC emitters.
If successful, consistent models should also allow us to better understand the origin of the different optical emission lines, their dependence on geometrical factors, leaking, and others. Hopefully, they could lead to the development of a new diagnostic of LyC leakage using rest-frame optical emission lines, applicable to other observations, especially of high-z galaxies, for which many spectra will be obtained in the near future with the James Webb Space Telescope (JWST).
The need for multi-zone (or multi-component) models is not only of importance for LyC emitters. Indeed, numerous papers have pointed out the necessity to account for variations in density, ionization parameter or geometry to simultaneously explain the observed strong emission line ratios of galaxies including different elements and ionization stages (see Stasińska & Leitherer 1996;Stasińska et al. 2006). Similarly, detailed studies of nearby galaxies with a large set of observations have clearly shown the limitations of one-zone models and the need for more complex and realistic models (e.g., Stasińska & Schaerer 1999;Cormier et al. 2012Cormier et al. , 2019Lebouteiller et al. 2017). We here explore some of these aspects with a focus on LyC emitting galaxies.
This paper is organized as follows. In Sect. 2 we present our observational samples. Sect. 3 shows the main optical emission line diagrams used to compare to "classical" single-component photoionization models (both density-and ionization-bounded). Then we construct two-component photoionization models including LyC leaking emission and compare them with the observations (Sect. 4). Our results are discussed in Sect. 5. Section 6 summarizes and concludes the paper. Throughout this paper, we assume a solar chemical composition following Asplund et al. (2009) where log(Z/Z ) = 12 + log(O/H) − 8.69. HST and whose properties have been amply discussed in several papers (Izotov et al. 2016a(Izotov et al. ,b, 2018aSchaerer et al. 2016;Verhamme et al. 2017;Gazagnes et al. 2018;Chisholm et al. 2018;Schaerer et al. 2018). The metallicities of these galaxies are in the range 12 + log(O/H) = 7.65 − 8.16, with a median oxygen abundance of 12 + log(O/H) of 7.91.

Optical emission line properties of LyC emitters
These sources represent currently the best sample of the local LyC emitters, for which a wide range of observations is available, including LyC, UV, and optical spectroscopy. Furthermore, Gazagnes et al. (2018) and Gazagnes et al. (2020) have undertaken a detailed analysis of their UV absorption lines (including Lyman series lines and other low ionization absorption lines), which provides important geometrical constraints on the ISM of these galaxies. We here analyze the optical emission line properties of these LyC emitters.
Before discussing the observed emission line properties of the LyC emitters, we briefly summarize the observational data used, including the main galaxies of interest (11 LyC emitters) and an appropriate comparison sample.

Lyman continuum emitters
The emission line measurements from the SDSS spectra of the LyC emitters are reported in the discovery papers (Izotov et al. 2016a(Izotov et al. ,b, 2018a. We have recently obtained new VLT spectra of five of the eleven leakers (J0925+1403, J0901+2119, J1011+1947, J1154+2443 and J1442-0209) with the XShooter spectrograph. We use the improved emission line measurements from these spectra, and in particular also the [S iii] λ9068 measurements, which are present in the near-IR spectra of these galaxies. A more detailed description of these observations and results are presented in Guseva et al. (2020). All emission line ratios reported here refer to extinction-corrected ratios.

Comparison sample
For comparison we use a sample of compact star-forming galaxies from the Sloan Digital Sky Survey (SDSS) Data Release 14 (Abolfathi et al. 2018) analyzed in earlier publications , and for which the [O iii] λ4363 line is detected with an accuracy better than 4σ, allowing direct abundance determinations using the T e -method. In practice, from the list of 5607 galaxies with direct metallicity measurements, we retain 4571 galaxies in the metallicity range 12 + log(O/H) = 7.6 − 8.2, for comparison with the LyC emitters discussed here.

Optical emission line properties of LyC emitters
The emission line properties of the LyC emitters and the comparison sample are shown in Figs. 1, 2, and 4. By selection all the galaxies studied here, including the LyC emitters, have emission lines ratios compatible with those of star-forming galaxies, i.e., they fall in the range of stellar photoionization sources in the classical diagrams by Veilleux & Osterbrock (1987) used for spectral classification (including the [N ii] BPT diagram introduced by Baldwin et al. 1981). Also, by selection, the sources studied here have a relatively low metallicity, which implies a fairly high excitation, reflected, e.g., by a high [O iii]/ Hβ ratio. At a finer level, however, the LyC emitters show some interesting distinctive properties, which we now discuss.

A [S ii]-deficiency in LyC emitters
In Fig. 1   To quantify this comparison we derive an SDSS "ridge-line" which marks the median location of normal star forming galaxies in the same range of O/H (metallicity) as our LyC emitters. It was obtained by binning the data in logarithmic bins of 0.1 dex along the [O iii]/Hβ axis with each bin having at least 10 galaxies. Each bin was then fitted by a gaussian distribution of a center x 0 and a standard deviation σ. The ridge-line defined by the x 0 was then fitted with an eight-order polynomial, where y =[O iii]λ5007/Hβ, and the coefficients α n are reported in Table 1.   Guseva et al. 2020). Therefore, the mean trends described by this ridge-line should be quite independent of the exact S/O abundance. As shown in Fig. 1 finding already pointed out by Wang et al. (2019). Several reasons could be invoked to explain such a behavior: S/H abundance variations, a decrease of [S ii] due to the loss of the outer, neutral regions related to the escape of LyC radiation, a lower contribution of diffuse ionized gas (DIG), which is in general known to "boost" [S ii] emission (cf. Sanders et al. 2020), or possibly other explanations. The [S ii]-deficit and its link with the observed LyC escape fraction will be discussed below more quantitatively (Sect. 4.4).

An excess of [OI]/[OIII]
Another peculiarity of the LyC emitters is found in the diagram comparing the three ionization stages of oxygen traced by the optical line ratios [ Fig. 2. This empirical finding clearly shows the presence of neutral oxygen in the LyC leakers, which already indicates that the ISM in these galaxies can probably not be described by simple density-bounded models, and whose presence and intensity needs to be explained quantitatively by models.

A shift in low ionization line properties
In Fig. 3, we show the displacement relative to the ridge-lines in O3Hβ-O1Hα and O3Hβ-S2Hα diagrams (see columns 2 and 3 in Table 1). In this plane, we see a clear shift from the SDSS comparison sample of presumably non-leaking galaxies. Interestingly, the three most leaking galaxies (the 3 reddest circles, f esc > ∼ 38%) in our sample stand in the lower-left corner with low ∆[S ii] and ∆[O i] (< -0.2 dex). On the other-hand, most of the small leakers in our sample (yellow circles with f esc < ∼ 10%) exhibit moderate deficiency in . We note, however, that three galaxies in our sample with large deficit in [S ii] (∆[S ii] < ∼ -0.3) have intermediate (J1011+1947, f esc =11,4%) to low escape fractions (J1248+4259, J1333+6246, labeled on the plot). One of them (J1333+6246) strikingly stands apart with a non-detection of [S ii] line but the largest [O i]-excess in the sample. A more detailed analysis of these shifts in [O i] and [S ii], including a comparison to models, is given in Sect. 4.4.

Other emission line peculiarities of LyC emitters
In the [N ii] BPT diagram (O3Hβ-N2Hα, Figs. 4 and 8) it is possible to note a slight offset of most of the leakers with respect to the average of the comparison sample. Probably, related to this, Guseva et al. (2020) have noted that the N/O abundance of LyC emitters is again higher than the average of star-forming galaxies at the same O/H abundance. The origin of this behavior is currently unknown.
In another diagram discussed later, Finally, Izotov et al. (2017) have proposed that LyC emitting galaxies show peculiar line ratios of weak He i lines. As shown by Guseva et al. (2020), our Xshooter observations confirm the trend suggested by this study.
To gain further insight on the observed behavior of the emission line ratios we now proceed to comparisons with photoionization models, first with simple and then more complex geometries.

Comparison to classical photoionization models
The first natural step in a comparison between the observed emission line properties of LyC emitters is to examine how simple 1D photoionization models allowing for varying degrees of ionizing photon escape behave and compare with the data.

The models: BOND models from the 3MdB database
To do so, we have used Cloudy photoionization models assembled in the BOND (Bayesian Oxygen and Nitrogen abundance Determinations) grid of models (Vale Asari et al. 2016), which is available on the Mexican Million Models database (3MdB, Morisset et al. 2015). This grid was designed to be used with the BOND code, a bayesian code, which simultaneously derives oxygen and nitrogen abundances in giant H ii regions. One of the specificity of the BOND code is to get rid of any assumption on relations between O/H and N/O. Hence, the grid spans a wide range in O/H, N/O and ionization parameter U, and covers different starburst ages and nebular geometries (spherical bubble or thin shell). The spectral energy distribution of the ionizing radiation is obtained from the population synthesis code PopStar (Mollá et al. 2009) for a Chabrier (Chabrier 2003) stellar initial mass function between 0.15 and 100 M . Varying starburst ages describe variations in the ionizing radiation field hardness, which arise due to the evolution of the stellar populations ionizing the H ii regions. A subset of the BOND model grid has been used by Stasińska et al. (2015), who studied also the effect of LyC leakage. In the current study, we use models from the 3MdB_17 database 2 (ref='BOND_2'), which is an updated version of 3MdB using Cloudy v17.02 (Ferland et al. 2017) and covering a somewhat refined grid.
For our study, we used a sub-grid of BOND models, with parameters appropriate for the LyC emitting galaxies and the comparison sample. Specifically, the metallicity range is limited to 12 + log(O/H)=7.6-8.2 (with a step of 0.2), covering the observed metallicities of the LyC emitters. The ionization parameter log U varies between -4 and -1 (with a step of 0.25), the age between 1 and 3 Myr (in steps of 1 Myr). The models are spherical filled bubbles with fixed density of 100 cm −3 , an abundance ratios of N/O= 10 −1 and S/O=10 −1.66 . The choice of a spherical filled bubble geometry influences the ionization structure of the modeled H ii regions: based on Stasińska et al. (2015), we expect thin shell models to produce fainter emission from low ionization potential species (e.g., . Concretely, this would induce a lateral shift to the left-hand side in the O32-O13 and O32-S23 diagnostics plots ( Fig. 5 and 9) that we will examine in Section 3.2.
To describe density-bounded models, we used a parameter H β f rac which cuts the nebula at different radii corresponding to a range of 10-100% of the total Hβ luminosity. 3 This parameter corresponds to the fraction of Hβ luminosity emitted by the entire modeled galaxy relative to the Hβ luminosity that would be emitted by the same model if it was ionization-bounded. Since the Hβ emissivity is roughly constant throughout the H ii region, this quantifies the ratio of the radius of the truncated, densitybounded model relative to the Strömgren radius of its ionizationbounded analog. We directly translate H β f rac into the corresponding LyC escape fraction by defining f esc (Hβ) = 1 -H β f rac that we use as a proxy to estimate how much leakage is allowed by a given density-bounded model. We refer to the models used here as simple, one-zone models.

The spectral diagnostics
We now compare the one-zone photoionization models to the observations, using the classical diagrams introduced by Veilleux & Osterbrock (1987) for spectral classification, the abovementioned O32-O13, and other diagrams. Figure 4 shows the LyC emitters, the comparison sample and the one-zone models in the classical diagrams from Veilleux & Osterbrock (1987), including [N ii] BPT diagram. As before, we also indicate the ridge-line, describing the location of the average SDSS galaxies in the corresponding plots (cf. Sect. 2.2). . Classical diagrams with BOND models and observations. Photoionization models are represented by the colored symbols: the colorbar represents the variation of the escape fraction and the size of symbol is proportional to log(U). The shape either represents the variation of metallicity (first row: for the age of 1 Myr) or stellar age (second row: for the metallicity of 8.2). The observed data points (gray dots) show the SDSS galaxy sample in the restricted metallicity range, 12 + log(O/H) = 7.6 − 8.2. The plain blue lines show the ridge-lines of the entire sample. Dashed lines correspond to the AGN delimitations from Kewley et al. (2001).

Classical diagrams
For illustration, dashed lines indicate the theoretical delimitations between star-forming galaxies (below) and AGN (above), according to the models of Kewley et al. (2001).
Overall, the ionization-bounded ( f esc = 0) photoionization models cover reasonably well the observed line ratios of the entire galaxy sample including the LyC emitters. The O3Hβ-N2Hα diagram is better covered by these models than the O3Hβ-S2Hα and O3Hβ-O1Hα diagrams. For example, the ridge-line in the [S ii] diagram is located outside the parameter space examined here; discussing the origin of these offsets is beyond the present scope. The position of the LyC leakers indicates high ionization parameters and young stellar ages of the ionizing populations; the latter is in good agreement with the young ages found from the observed UV spectral features and from spectral modeling of these galaxies (see, e.g., Izotov et al. 2019).
The second most important conclusion from Fig. 4 is that the domain of the emission line ratios observed in the LyC emitters is only covered by ionization-bounded photoionization models, i.e., models with f esc = 0 or clearly only for low values of f esc < 0.1, which is in contradiction with the observed escape fractions reaching up to ∼ 70% for some of the leakers. As already shown by Stasińska et al. (2015), increasing the escape fraction in these simple one-zone, 1D photoionization models, decreases strongly the [S ii] and [O i] emission, which results from trimming the outer parts of the H ii region. This shows that overall the observed LyC emitters (or at least those with a significant LyC escape fraction) cannot be reproduced by simple density-bounded models, as these systematically underpredict lines of species with low ionization potentials, such as [O i] and [S ii].

Diagnostics of the low ionization regions
In Fig. 5 we examine complementary and probably more powerful emission line diagnostics to probe the ionization structure and hence different zones of the ionized ISM in the LyC emitters. By focusing on different ionization stages of the same element, one also avoids dependences on elemental abundance ratios and metallicity, at least to first order.
These O32-O13 and O32-S23 diagnostics highlight two main points that will motivate the introduction of more complex models in the section 4: (1) Although the oxygen line ratios (O32-O13) predicted by simple ionization-bounded models reproduce quite well the SDSS sequence with ionizing spectra of young populations (∼ 1 − 2 Myr), the models are not able to account for galaxies with a certain [O i] λ6300 excess, i.e., a higher-than- ratio, which is observed for the LyC emitters and other sources (see Fig. 5

left). For [S ii]/[S iii],
shown on the right panel of this figure, the deficiency of the models is even more striking, as they fail to reproduce the bulk of the observations. The two discrepancies might arise since the models do not take into account a diffuse ionized gas (DIG) component, which can significantly contribute to the emission of [O i] and [S ii] (Sanders et al. 2017(Sanders et al. , 2020. Alternatively, these discrepancies could also be attributed to the presence of shocks, which are not included in our models (see, e.g., Plat et al. 2019), due to problems with atomic data for sulfur (Izotov et al. 2006;Kewley & Dopita 2002;Kewley et al. 2019), or others. Sulfur lines ratio [S ii]/[S iii] is also sensitive to more subtle effects including stellar winds, outflows or photoevaporation of photodissociation regions (PDR) in galactic dense gaseous environments (e.g., Westmoquette et al. 2013;McLeod et al. 2015).
(2) The density-bounded models allowing for LyC photon escape do not at all reproduce the observed oxygen and sulfur line ratios of the leakers, as already noted above. On the contrary, while leakers populate the upper part of the SDSS sequence, the density-bounded models produce very low emission in [O i] and [S ii], which drives them in the opposite direction to the bottom-left corner of our diagnostic plots. Figure 6  In short, although density-bounded models would explain the Lyman continuum escaping photons, they clearly do not match the spectral signatures arising from the outer part of the H ii region and which are observed in LyC emitters. On the other hand, while ionization-bounded models can fairly well describe the overall properties of most galaxies (the comparison sample), except for sulfur line ratios, they underpredict the emission from low ionization species (e.g., the [O i] line) in the LyC emitters, and most importantly, they do not account the escape of LyC photons. These failures of simple, one-zone photoionization models, clearly indicate that more complex models are required to describe the observed emission line properties of the LyC emitters. We have, therefore, examined composite photoionization models based on the combination of two distinct zones with different properties.

Motivation and method
As just demonstrated (Sect. 3.2), single-component (one-zone) photoionization models fail to simultaneously reproduce the observed line ratios and escape fractions of LyC leakers. To answer our problem (1) (insufficient emission in low ionization lines), we introduce a component of low ionization (low U) that would mimic a DIG contribution in star-forming galaxies. The relative contribution of this component is described by a factor ω that varies from one galaxy to another. In the case of LyC emitters, although they are expected to be relatively poor in DIG (see Oey et al. 2007), this component could effectively represent the effect of low ionization channels carved into the ISM through which some LyC photons escape, and which may contain gas with a low volume filling factor, such as clumps or filaments. As an answer to problem (2) (inconsistency between the escape fraction and the presence of low ionization lines), we include in our two-zone model the possibility for LyC photon escape from one region, while the second remains ionization-bounded.
A two-zone, "picket-fence" type of model, is also suggested by other observations. Indeed, Gazagnes et al. (2018), Chisholm et al. (2018), and Gazagnes et al. (2020) have analyzed the UV absorption lines of hydrogen and low ionization metal lines of the LyC emitters studied here. Their work clearly shows the presence of absorption lines formed in the neutral gas, with column densities too high for the escape of LyC photons. Furthermore, the depths of the saturated absorption lines show that channels of low column density must be present in these galaxies, allowing LyC escape. For our purpose, their model is basically equivalent to our two-zone model. As we will show below, such an idealized two-zone model can also capture and describe the observed optical emission lines of the LyC emitters.
In practice, we have considered two different scenarios that are represented in Fig. 7. Each scenario is based on a linear combination of two models from the BOND database (see Sect. 3.1), one with a high ionization parameter weighted by a factor ω, and one with a low U weighted by a factor (1 − ω). In each scenario, we allow only one component to be leaking.
The high U component corresponds to the main body of the galaxy which is highly ionized and homogeneous (filling factor ∼ 1). We briefly discuss here the motivation for introducing a low ionization component. The BOND models are parameterized by an input volume-averaged ionization parameter that, for a given geometry, is proportional to the ionization parameter of the corresponding ionization-bounded model, at the Strömgren radius (see Equation 4 in Stasińska et al. 2015). For ionizationbounded models, the ionization parameter at the Strömgren radius is defined by the following equations: where Q is the number of ionizing photons, n the density of hydrogen, α(T e ) the recombination coefficient of hydrogen, and the filling-factor for an inhomogeneous medium (0 < ∼ < ∼ 1). Equation 3 leads to: where A is almost constant (A ≈ 1.3 × 10 −21 cm s 1/3 for T e = 10 4 K), which means that low density or clumpy medium ( 1) will both correspond to a model with a low input ionization parameter. We note that the density-bounded models created by imposing a cut in radius to the full ionization-bounded model do not necessarily have a mean ionization parameter that matches their input ionization parameter. The two components do not have to share exactly the same Strömgren radius but one can ensure that the physical sizes are compatible by adjusting the filling factor of one component to match the size of the other. We now further describe the scenarios in Fig. 7.
Density-bounded channels (DBC): ionization-bounded galaxy with density-bounded channels. The galaxy is ionizationbounded but LyC photons can escape through channels, which are matter-bounded. Such holes could be carved in the ISM by supernovae explosions or could result from hydrodynamical fragmentation effects induced by turbulence or stellar feedback (e.g., Kakiichi & Gronke 2019;Hogarth et al. 2020). We expect such mechanisms to efficiently clear out channels, either producing low density channels or almost emptied holes with small clumps of matter remaining within. In both cases we model such channels by a density-bounded component with low input ionization parameter (see Equation 4). Since the LyC emitters studied here are clearly dominated by very young stellar populations with ages < ∼ 3 − 4 Myr, too young for significant feedback from supernovae (Izotov et al. 2018b), we suspect radiation pressure and related feedback effects to be the main driver for the creation of these channels or holes, although the explanation may be more complex (cf. Hogarth et al. 2020).
Density-bounded galaxy (DBG): density-bounded galaxy with ionization-bounded clumps. The core of the galaxy is completely ionized by massive young stars, which produce a large number of ionizing photons exceeding what can be absorbed by the available nebular gas. This is represented by a high U, density-bounded component, from which photons can directly escape at the edge. We add a low ionization component to represent denser, ionization-bounded filaments or small clumps, which can naturally emit the [O i] and [S ii] emission lacking in one-zone models (see section 3) and where the observed UV low ionization lines are formed. Both components have the same chemical composition and the same ionizing spectrum, emitted by a central source.
For a first approach, we combined two arbitrary photoionization models with different ionization parameters: one with a high  . Schematic illustration of our two-zone photoionization models for LyC emitters, combining an ionization-and density-bounded regions, chosen to have different ionization parameters ("high" and "low" U). In the "density-bounded channels" (DBC) scenario, ionizing photons escape from the low U region, in the "density-bounded galaxy" (DBG) scenario, they escape from the high U region with f esc =0 for the non-leaking region and f esc between 0-80% for the leaking region.
ionization parameter U h and one with a low ionization parameter U l . The high and low ionization parameters were chosen "by hand" so that the non-leaking combined models reproduce well the upper part of the SDSS sequence on our diagnostic plots (see Sect. 4.2) and the tight SDSS sequence in the O32-O13 and O32-S23 plots (Fig. 5 and 9). We used U h = 10 −1 and U l = 10 −3.5 . We briefly comment on this choice and an optimizing procedure below (Sect. 5.1). We computed both scenarios for an identical stellar population of 1 Myr, a metallicity 12 + log(O/H) = 8, and N/O=0.1, as above. Since our combined model is meant to represent a galaxy ionized by the same stellar population, we have to make sure that the photoionization models describing the two components also correspond to the same number of ionizing photons Q h and Q l , respectively. This is simply achieved by scaling the predicted line luminosities of the U l model with the factor s = Q h /Q l . Then, the two components are combined with relative weights ω for the U h component, and (1-ω) for the U l component. Any ratio of emission lines A and B of the combined model then becomes For each scenario, we allowed the LyC escape fraction of the leaking component to vary from 0 to 80% in steps of 10%, provided by the 3MdB_17 database. The total f esc of the twozone model, is computed consistently from the relative weight of zones ω, by defining Since only one component is leaking in each scenario, Eq. 6 simply becomes f c esc = (1−ω) f l esc in the DBC scenario and f c esc = ω f h esc in the DBG scenario. In the absence of leakage ( f esc = 0), the two scenarios are trivially identical.
The resulting emission line ratios predicted from our twomodels with varying relative weights and LyC escape fractions are shown in Figs. 8 and 9, allowing thus a comparison with the one-zone models (Figs. 4 and 5). Non-leaking models are represented by the yellow symbols having 0% of escaping photons. In this case the two scenarios (DBC and DBG) are identical, by construction, and their symbols (yellow triangles on the first row, respectively yellow diamonds on the second row) are thus superposed. As mentioned before, we chose the ionization parameter of each component such that the combined non-leaking models cover well the upper part of the SDSS sequence.

[O i] and [S ii] λλ6716,6731 emission in normal galaxies explained
As we can see from Fig. 8, the non-leaking two-zone models follow a bended or sometimes straight "mixing line" from very high to moderate excitation (in [O iii]λ5007/Hβ) and covering a wide range in the [N ii]/Hα, [S ii]/Hα, and [O i]/Hα ratios. Despite our somewhat arbitrary choice of models describing the individual zones, we already see that observed line ratios of the bulk of the SDSS galaxy sample (which primarily hosts galaxies with very low or zero f esc ) can be covered with the two-zone models combining two different ionization parameters. In particular, there seems to be no difficulty to explain the relatively high observed intensities of [O i] λ6300, and even the comparatively high values of [S ii] λλ6716,6731/[S iii] λ9068 in normal galaxies, which cannot be reproduced by classical one-zone models (compare Fig. 5 and Fig. 9). The inclusion of a second ISM com-  Fig. 4) showing the predictions from the two-zone photoionization models. The upper and lower pannels represent DBC and DBG models, respectively. ponent with a lower ionization parameter seems thus both essential and sufficient to solve the above-mentioned discrepancies of classical models when compared to normal galaxies.

LyC emitters explained by two-zone models
In contrast to the one-zone models, the combined two-zone models are able to cover the entire range of emission line ratios of the LyC emitters, as shown by Figs. 8 and 9. Most importantly, the models can -at least qualitatively -explain the observed [O i] excess and the [S ii]-deficiency (see Fig. 1 and Fig. 2), which characterizes the leakers, and which cannot be understood by one-zone models (see Sects. 2.2 and 3.2). Furthermore, the twozone models account for the escape of LyC photons, and are thus consistent with this important observational property. Below we will further explore the consistency of different line ratios, and to which extent the models also correctly reproduce the observed LyC escape fractions.
Let us now examine how the predictions from the two-zone models in classical diagrams vary with LyC escape fraction. In the O3Hβ-N2Hα diagram (Fig. 8 left), most of the models are superposed quite independently of their f esc value. This is due to the fact that [N ii] originates from the highly excited region in the inner part of the H ii region, which remains unaffected when trimming the outskirts of the galaxy. Only the most densitybounded models from the DBG scenario, where the highly ionized component dominates almost all the galaxy (ω = 99%) and is leaking, affect the [N ii]/Hα ratio, when f esc exceeds 10%. In contrast, the [S ii] and [O i] line ratios (shown in the middle and right panels) are very sensitive to LyC leakage. We note that these emissions primarily arise from the low ionization component, i.e., the one leaking in the DBC scenario (upper pannel). This explains why these line ratios are barely affected in the DBG scenario (lower pannel), even for high escape fractions. The same trend is visible on the O32-O13 plot (Fig. 9). The SDSS sequence seems bounded by f esc ∼ 0 − 10% for the DBC scenario (first row). Such values are sufficient to explain the observed escape fraction of some of the leakers, represented by yellow circles, but it fails to reproduce higher leakage of ionizing photons, represented by the orange to brown circles. Meanwhile, the DBG scenario can produce escape fractions up to 80% without decreasing strongly the emission of [O i] and [S ii]. On the other-hand, the O32-S23 plan seems to be better covered by models from the DBC scenario, while the DBG scenario tends to overestimate the [S ii]/[S iii] ratio. For the five leaking galaxies with [S iii] measurements, the observed escape fractions seem consistent with predictions from the DBC scenario, even for the most leaking galaxy J1154+2443 (red circle, f esc = 46%). Indeed, we note that this galaxy stands on the left-hand side of the median line, as expected for the highest escape fractions produced by the DBC scenario. It falls between two model points with escape fractions of 30-40%, which is roughly consistent with the observed value. We also note that models with a weight ω < ∼ 0.7 fail to reproduce the leakers studied here. This is easily understood if we keep in mind that these galaxies show a high excitation, as indicated by their high O32 ratio: decreasing the contribution ω of the U high region below a certain threshold produces two-zone models that are dominated by the low ionization (U low ) component, and hence are unable to reproduce the observed line ratios of the leakers.

Links with f esc
In Fig. 10 we now revisit the predicted O32 and O31 line ratios in models with LyC escape presented earlier for one-zone models (see Fig. 6). As already mentioned and clearly shown in this Figure, the two-zone models solve the [O i] problem of simple models, and they cover the entire range of observed line ratios of the LyC emitters. A comparison of the two panels also shows that, at least approximately and without adjustment of the components of the two-zone models, the models are capable of consistently reproducing the two oxygen line ratios and the observed escape fraction. This Figure also shows that LyC emitters with very high f esc > ∼ 20% can only be reproduced by our DBG scenario (where the high U component is leaking), whereas both scenarios are possible for the galaxies with lower escape fractions. Figure 11 (left) illustrates the evolution of f esc as a function of ∆[S ii]. This plot suggests that a sufficiently strong [S ii]deficiency (e.g., ∆[S ii] < ∼ −0.2) effectively picks out LyC leakers, although not all of them. Indeed, several LyC emitters are also located within the normal range of the comparison sample (presumably dominated by non-leakers), as also pointed out by Wang et al. (2019). Using ∆[S ii] to determine the escape fraction is, however, not obvious. Compared to the dispersion of our comparison sample represented by the gray-shaded areas, 6 out of 11 leakers show a deficit in [S ii] larger than 1σ (5 out of 11 with deficit larger than 2σ). If we include the 2 leaking galaxies from Wang et al. (2019), this amounts to 8 out of 13 galaxies with clear [S ii]-deficit at 1σ (6 out of 13 at 2σ). The three strongest leakers (J1256+4509, J1154+2443 and J1243+4646), which have f esc > ∼ 30% all exhibit large deficiency in [S ii] of at least 0.2 dex from the ridge-line (> 1σ). J1011+1947, with an intermediate escape fraction of 11.4% also shows a large [S ii]deficit. However, the most [S ii]-deficient galaxy in our sample (J1248+4259) and J1333+6246, which is undetected in [S ii], have escape fractions of only a few percents. This could mean that the observed f esc for these galaxies is underestimated, for example, due to that LyC radiation escapes only along other, unobserved lines of sight. Interestingly, in the analysis of the UV absorption lines of more than 30 low-z galaxies from Gazagnes et al. (2020), J1248+4259 stands out as an outlier with an exceptionally strong Lyα emission (high equivalent width) despite showing a large covering fraction in the H i absorption lines. Demonstrating if any of these sources indeed emits more ionizing photons in other directions seems difficult. In any case, our analysis confirms that the method proposed by Wang et al. (2019) could be a promising way to identify new leaking galaxies, although the determination of f esc from the [S ii]-deficit remains challenging.
From Fig. 11 (left) we notice once more that only the DBG scenario is able to reproduce the more extreme leakage above 20%. The DBC scenario is able to cover the space, where reside most of the galaxies with small f esc , including the [S ii]-deficient galaxies selected by Wang et al. (2019) (green color along the contours of circles). We note, however, that the two galaxies hav-ing the stronger [S ii] emission (∆[S ii] close to 0) stand outside of the contours of the DBC scenario. This indicates that leakage powered exclusively by channels of low ionization is unable to produce, at the same time, strong [S ii] emission and escaping LyC photons. This comes from the fact that in this scenario leakage is allowed exclusively by trimming the low ionization component, which always results in decreasing [S ii] emission. This could either suggest that density-bounded LyC could still be the preferred mechanism even for rather small ( < ∼ 10%) escape fractions, or that a more complex geometry needs to be considered to allow leakage without systematically decreasing [S ii].
In Fig. 11 (right) we show ∆[O i], calculated from the O3Hβ-O1Hα plots and defined earlier (see Fig. 4 which traces the photo-dissociated interface between ionized and atomic gas or cold dense clumps in the H ii region might be linked to the disappearance of the remaining Lyman continuum, leading to lower escape fractions. We notice, however, that two galaxies (J1248+4259 and J0901+2119) that are [O i]-deficient exhibit very small escape fractions ( f esc ∼ 2%). This indicates that other effects are likely to cause the drop of Lyman continuum even in [O i] deficient environment. However, the overall trend shown by the observations in Fig. 11 suggests that selecting [O i]-deficient galaxies could also be an effective way to target galaxies with high leakage and might correlate better with f esc than [S ii]-deficiency. On the other hand, the [S ii] lines are generally stronger than [O i], making the former easier to measure.
In Fig. 11 (right) the contours from the DBC scenario include only the two galaxies with the lowest escape fractions. We conclude that ∆[S ii] and ∆[O i] could be interesting proxies to estimate the escape fraction but only in the case of ionization-bounded leakage (DBC scenario). In the case of density-bounded leakage (DBG scenario), likely producing the highest escape fractions, [O i] and [S ii] seem uncorrelated with f esc , and other diagnostics, sensitive to density boundaries, are to be preferred (see, e.g., Fig. 5 and 9). It seems that of the ionizing cluster, its evolutionary stage, the gas content of each galaxy, its star formation rate and metallicity, all of which could produce other combinations for a LyC leaking cluster (e.g., Wang et al. 2019).

Discussion
We now briefly discuss further improvements, potential caveats, implications and future prospects resulting from the present work.

Can the LyC escape fraction be detected and measured from the optical emission lines ?
In this study, we used the observed LyC escape fractions and optical emission lines to build simple, but consistent two-zone photoionization models of the LyC emitters, and to discriminate between different leakage scenarios. An obvious further step would be to use these two-component models to infer the expected es-cape fraction from the emission line observations only. Such indirect methods are badly needed in order to optimally exploit upcoming JWST observations of rest-frame optical emission lines for galaxies during the epoch of reionization. A first automated routine was implemented in an effort to extract best-fit values for the parameters of the two-component models. Although very promising, this method requires a robust implementation and a careful choice of the tracers to use that we postpone to future work. We present here the first trends obtained from this study and discuss possible caveats of this method.
We ran a χ 2 -minimization routine on a sub-grid of our combined models, which includes the ones presented in Sect. 4. We defined a 4-dimensional grid with the following free-parameters: U high and U low (varied from 10 −4 to 10 −1 ), the weight of the component of high ionization ω, and the escape fraction of the leaking component f holes esc (both varied between 0 and 1). As observational constraints we used six line ratios: the four ratios used in the classical diagrams (Fig. 8), and two ratios from the O32-O13 diagrams (Fig. 9). We did not include upper limits. As our Article number, page 13 of 16 A&A proofs: manuscript no. leakers_EL_0 prescription for f combined esc is fully determined by the parameters of each model (see Eq. 6), we could also determine the corresponding value of the total LyC escape fraction for the χ 2 minimum. We also examined the weight of the ionization-bounded component ω in the DBC scenario (respectively (1 − ω) in the DBG scenario) to compare it to the UV covering fractions deduced from the Lyman series absorption lines Chisholm et al. 2018;Schaerer et al. 2018). Hereafter we summarize some interesting trends that will motivate further examination: -The models minimizing χ 2 always favor two distinct values for the ionization parameters U high and U low . Inspection of their marginalized probability distribution functions shows that typically U high ∼ 10 −1 − 10 −2 , and U low ∼ 10 −3 − 10 −4 . We interpret this as a confirmation of our conclusion in Sect. 4 that two-component models do a better job at reproducing the observed line ratios than single component models. -Regardless of the considered scenario, the weight of the high-ionization component is always much higher than the low-ionization one (ω ∼ 0.8). This means that no matter which component is leaking, the high-ionization component always dominates the total emission, in the same proportion of about 80% (respectively 20% for the low-ionization component). This suggests that the evolution from one scenario to the other could result from a continuous evolution (cf Kakiichi & Gronke 2019) leaving the global proportions of the galaxy's structure unchanged. However, we note that the two scenarios produce very different covering fractions. We note that for galaxies with a low escape fraction ( f esc < ∼ 10%), the ω obtained from the best-fitting model has a value comparable to the covering fractions, C f , obtained from absorption lines analysis Chisholm et al. 2018;Gazagnes et al. 2020). This is consistent with what we expect from the DBC scenario, where the covering fraction should be close to ω but cannot be reconciled with the DBG scenario, where C f ∼ (1 − ω) and which is a better fit for some of these galaxies. On the other hand, for strong leakers ( f esc > ∼ 10%), we find C f < ω. This strengthens the conclusion in Sect. 4.4 that the DBC scenario is unable to reproduce strong leakage, as it predicts too high a covering fraction, and hence too low escape fractions. However, we note that the DBC scenario sometimes better reproduces the six line ratios we considered, even for strong leakers. -Finally, the escape fraction obtained for the best-fit models allows us to distinguish the four galaxies with the strongest leakage in our sample. In fact, the best-fit models predict f esc ≈ 0 for all the weak leakers but one, and nonzero escape fractions ( f esc ∼ 0.6 − 38%) for the galaxies with high escape fractions. We note that the highest predicted escape fraction corresponds to J1243+4646, the strongest leaker in our sample, although the predicted value of f esc is a factor ∼ 2 below the observed value.
In short, this exercise supports our conclusion in Sect. 4.4 that the geometry of weak leakers resembles that of the DBC scenario, and strong leakers more the DBG scenario. However, we find that a robust distinction between those scenarios, based only on emission lines is not easy. Furthermore, accurate predictions of f esc from the optical emission line ratios do not appear to be straightforward. This will probably require a careful choice of the relevant line fluxes or ratios to be used in such an analysis, a proper treatment of upper limits and uncertainties, a robust statistical method to overcome the caveats of a simple χ 2 minimization, the exploration of SEDs allowing also for different ages, and other improvements, which are beyond the scope of the present publication.

Low ionization line diagnostics: progress and caveats
In Sect. 3, we discussed the inadequacy of simple models to reproduce simultaneously the high and low ionization lines, as revealed, for example, in the O32-O13 and O32-S23 diagnostics (Fig. 5). Specifically, we emphasized the offset of BOND models in O32-S23, which systematically underestimate the [S ii]/[S iii] ratio. This offset of standard photoionization models has been shown in previous work (e.g., Pérez-Montero & Vílchez 2009;Kehrig et al. 2006), in which such models fail to explain the hardening in ionizing radiation shown by most low mass, low metallicity H ii galaxies in the O32-S23 plane. Problems with the sulfur line strengths have often been reported in the past (see, e.g., Garnett 1989;Kewley et al. 2019;Mingozzi et al. 2020, and references therein). The shortcomings of several Cloudy and Mappings models to reproduce the [S ii]/[S iii] ratio have, for example, recently been shown by Mingozzi et al. (2020). Considering that this only seems to affect sulfur lines, these authors suggest that this is due to limitations in stellar atmosphere modeling and/or the atomic data for sulfur. Based on the observation that the discrepancy holds for single H ii regions (from Kreckel et al. 2019), they exclude the role of diffuse ionized gas as a cause of this issue. However, in our study, we see in Fig. 9 that adding a low ionization component, which provides additional emission of [S ii] without producing [S iii], helps to solve this discrepancy. We conclude that for both our reference sample of star-forming galaxies and the LyC emitting galaxies, the addition of such a low ionization component is sufficient to reproduce the observed [S ii]/[S iii] ratios. Furthermore, this component also contributes to low ionization emission of oxygen, simultaneously reproducing the observed line strength of [O i] λ6300 and lines tracing the other ionization stages of oxygen. The exact nature, spatial distribution, and other properties of such a component remain, however, to be understood.
While this success of our two-zone models is notable and certainly calls for further studies and applications, it must be noted that other physical processes can also play a role for the low ionization lines. For example, it is well known that shocks in the ISM or the presence of sources with hard ionizing spectra (far-UV and X-ray emission) lead to enhanced emission of highly ionized species (e.g., [Ne v], [Fe v], [O iv], He ii, and others) and also of low ionization lines (e.g., Schmitt 1998;Stasińska et al. 2015;Plat et al. 2019). The effect of a contribution from shock models with different magnetic field strengths (tuned to reproduce a typical He ii/Hβ ratio) on the O32-O13 diagram is shown in Stasińska et al. (2015). Clearly, some model combinations of shock with both density-or ionization-bounded H ii regions can be found, which reproduce the three ionization stages of oxygen. How these models would affect the sulfur lines, and if they could also be consistently reproduced remains to be seen. Given different degeneracies and uncertainties it seems currently difficult to exclude a contribution from shocks (see Stasińska et al. 2015). Similarly, a weak AGN component (contribution ∼ 3% of the hydrogen ionizing photons) cannot be excluded based on the observed emission line diagrams for the galaxies of our reference sample (Stasińska et al. 2015). Again, the sulfur line ratios used here may provide an additional handle on these questions.
Finally, the low ionization lines are also affected by diffuse ionized gas (DIG) in galaxies, as recently discussed by Mingozzi et al. (2020), Shapley et al. (2019), Sanders et al. (2020) and others. In the LyC emitters studied here, the contribution of the DIG would presumably be low, if we adopt the observed scaling of the DIG contribution with surface density of star-formation (e.g., Oey et al. 2007;Shapley et al. 2019), since the LyC emitters are compact by selection, and have a very high SFR surface density (Izotov et al. 2018b).
Whatever the exact origin of emission of the low ionization lines (such as [O i] λ6300 and [S ii] λλ6716,6731) is in LyC emitting galaxies, the UV observations of these galaxies clearly reveal the presence of both low H i column density (N HI ) lines of sight towards the observer allowing the escape of LyC radiation, and lines of sight with high N HI , where other neutral and low ionization metal absorption lines are also formed Chisholm et al. 2018;Gazagnes et al. 2020). Hence, their ISM contains at least two "zones", i.e., regions with different conditions, even if we do not know on which spatial scales this exactly happens. Furthermore, it is very plausible that properties, which influence the ionization parameter (e.g. density, clumpiness, filling factor) differ between these regions. This justifies our two-zone approach using two different ionization parameters, U high and U low , for the LyC leakers. If, and to what extent such an approach is also motivated for "normal" galaxies (most of which are presumably not LyC emitters) remains to be seen. In any case given the ubiquitous complexity of the ISM in resolved star-forming galaxies it is probably not surprising that simple, 1-D photoionization models cannot reproduce all the observables.

Future improvements
Some of the methodology used and some of the questions addressed here certainly call for other applications and for future improvements. Indeed, although photoionization models are frequently used to predict and analyze emission lines in starforming galaxies, they are generally 1-D, spherical or plane parallel models computed with Cloudy or similar codes, which then allow for exploration in a wide range of parameter space, in terms of ionization parameters, radiation field, abundances etc. Such models, which typically assume acsingle source of ionizing radiation and a very simple geometry, when compared to observations of integrated galaxy spectra, are obvious oversimplifications of the true complexity and diversity of the radiation field, ISM, and geometry.
Relatively few more complex models have so far been used for comparisons with observations. Stasińska et al. (2015) have, for example, explored the combination of stellar photoionization and other emission sources (AGN, shocks and others) with variable but arbitrary contributions, and confronted them to a large sample of star-forming galaxies from the SDSS. Cormier et al. (2012Cormier et al. ( , 2019 have used one and two-component models to consistently analyze emission lines originating in the ionized and neutral ISM, as well as in adjacent photo-dissociation regions (PDR). These have been compared to detailed multi-wavelength observations of 39 nearby dwarf galaxies. A very complex tailored model including up to five components, allowing also for ionization-and density-bounded parts, has been constructed for the well known metal-poor galaxy I Zw 18 by Lebouteiller et al. (2017). Clearly other models can and need to be envisaged, although sufficient independent observational constraints are necessary to avoid degeneracies in a multi-dimensional parameter space.
Another approach is to use high resolution hydrodynamic simulations of galaxies, ideally including radiation transfer and a proper treatment of the ISM, to predict the resulting emission lines and compare them to observations. Such models, with fairly different degrees of sophistication have, for example, been presented by Shimizu et al. (2016), Hirschmann et al. (2017, Katz et al. (2019), Wilkins et al. (2020), and by Pellegrini et al. (2020) for smaller scales. They are able to capture several observational trends and allow useful predictions, in particular, for high-redshift galaxies. However, direct and detailed comparisons to observations are not trivial to undertake or to interpret. For example, Katz et al. (2020) predicts optical and IR emission lines of LyC leakers and non-leakers, which show a good agreement for the [O iii] λλ4959,5007, [O ii] λλ3726, 3728, and Hβ emission lines when compared to observations. In contrast, their predictions for [S ii] λλ6716,6731/Hα and [O iii]/Hβ are quite different than those observed in low-z LyC emitters and in z ∼ 2 − 3 galaxies. Understanding the origin and importance of such a discrepancy is, however, not straightforward. Different, complementary modeling approaches should help to develop more robust diagnostics and thus a better insight on the properties of the ISM, radiation field, geometry and other properties of star-forming galaxies and LyC emitters.
In parallel, new observations, such as those from the HST Low-z LyC survey (PI A. Jaskot), will provide for the first time a statistical sample of more than sixty z ∼ 0.2 − 0.4 galaxies with direct LyC observations. This sample of confirmed LyC emitters and also non-LyC emitters should represent a useful base to test several of the trends observed here, and an ideal sample for detailed comparisons with emission line models from different approaches. The first comparison with simple two-zone models and a small sample presented here is very encouraging, but incomplete. Nevertheless, there is hope that optical emission lines can to some extent be used to identify LyC emitters, and ideally also to allow a determination (even approximate) of the escape fraction of ionizing photons from these galaxies. Clearly, future improvements on these lines will be very welcome and of great interest, especially in view of the upcoming launch of the JWST, which will routinely observe rest-frame optical emission lines out to the highest redshifts.

Conclusion
We have analyzed the main optical emission lines of low-z LyC emitting galaxies (leakers) discovered recently with the HST (Izotov et al. 2016a(Izotov et al. ,b, 2018aWang et al. 2019)  We have then compared the observations to simple one-zone and two-zone photoionization models accounting the escape of ionizing radiation to construct, for the first time, consistent models which reconcile their observed emission lines and the measured LyC escape fractions. The main results from this analysis can be summarized as follows: -Classical one-zone photoionization models underpredict the low-excitation line emission of [O i] and [S ii] arising from the outskirts of H ii regions, strongly for density-bounded models with LyC escape, and also somewhat for ionizationbounded models (see Sect. 3). -Two-zone models, combining regions with a high and low ionization parameter (U high , U low ), where one of which is density-bounded, have allowed us to resolve these discrepancies and to reproduce the observed line ratios of the LyC emitters. Furthermore, the models show LyC escape fractions comparable to those measured directly from the UV LyC observations. -From the two-zone models we also found a possible dichotomy between galaxies with different LyC escape fractions. Indeed, the models indicate that galaxies with f esc < 10% can be explained by the presence of low columndensity, low filling factor regions or channels with a low ionization parameter, through which the ionizing photons escape. The observed [O i] emission can originate from these filaments/channels, in what we called the DBC scenario (Fig. 7). For much higher f esc , our photoionization models show that LyC leakage must occur from the region with high ionization parameter (DBG scenario) in order to explain the observed emission line ratios.
Our two-zone models are inspired by and compatible with the "picket-fence" type of model, which successfully reproduces the observed UV absorption lines of hydrogen and low ionization metal lines of the LyC emitters, as demonstrated by the studies of Gazagnes et al. (2018Gazagnes et al. ( , 2020. Our finding of two different scenarios, favored by LyC emitters with low/high escape fraction, is also compatible with the recent analysis of Gazagnes et al. (2020) combining information from the Lyα line profiles and absorption lines. Our results are also consistent with the numerical simulations of Kakiichi & Gronke (2019), which suggest that the mechanism powering LyC escape may evolve in a continuous way from the DBC scenario, producing low escape fractions, to the DBG scenario, which appears as a necessary configuration to reach high levels of leakage.
Finally, we have discussed the possibility of improved analysis tools to predict escape fractions based on optical emission lines and using two-zone photoionization models (Sect. 5). We conclude that robust predictions would require multi-phase and multi-sectors modeling, including some density-bounded components to successfully capture the complexity of the ISM structure that leads to photons leakage. We also suggest that lowexcitation tracers from the outer part of H ii regions such as [O i] and [S ii] lines should be included in the list of relevant lines to study the morphology and dynamical state of leaking and nonleaking galaxies.