Open Access
Volume 644, December 2020
Article Number A21
Number of page(s) 16
Section Extragalactic astronomy
Published online 24 November 2020

© L. Ramambason et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Significant progress has been made in recent years in searches for and studies of star-forming galaxies, which emit copious amounts of Lyman continuum (LyC) radiation. This ionizing radiation can escape the interstellar medium (ISM) and ionize the surroundings, thus contributing to intergalactic ionizing radiation. If a sufficient population of such galaxies exists in the early Universe, their presence could explain cosmic reionization, which is known to have started shortly after the Big Bang and to have ended ∼1 Gyr later, at z ∼ 6.

After intense efforts to search for LyC emitters or LyC-leaking galaxies (abbreviated to “leakers”), several groups have identified such galaxies from low redshift up to z ∼ 4 (e.g., Leitherer et al. 2016; Izotov et al. 2016a,b, 2018a,b; Vanzella et al. 2016, 2018, 2020; de Barros et al. 2016; Steidel et al. 2018; Fletcher et al. 2019; Bian et al. 2017; Wang et al. 2019) with LyC escape fractions, fesc, ranging from a few percent up to ∼72% in some 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 with a double-peaked Lyα line profile (see Verhamme et al. 2015, 2017; Vanzella 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). Following the suggestion by Jaskot & Oey (2013) and Nakajima & Ouchi (2014) that high values of [O III]/[O II] could pinpoint density-bounded galaxies, that is, 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 (fesc ∼ 2 − 72%), and are therefore good analogs of the sources of cosmic reionization (cf. Schaerer et al. 2016).

Despite the fact that high [O III]/[O II] was used as a basic selection criterion for these searches, it nevertheless turns out that fesc does not well correlate with the [O III]/[O II] ratio (see 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 by their [S II] deficiency, 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, these authors were not able to find a correlation between the [S II] deficiency and the LyC escape fraction of the known leakers. These findings are not entirely surprising, as it is well known that the [O III]/[O II] ratio depends also on several other factors, including in particular the ionization parameter and metallicity. Whether or not density-bounded models are applicable to known LyC emitters therefore remains 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) revealed the presence of high-column-density gas that is optically thick to LyC radiation and co-exists with lines of sight, which allows the escape of UV ionizing and nonionizing 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 is required, as well as the exploration of more complex structures to consistently describe their emission line spectra. These are two of the main goals of our study.

Specifically, we would like to know whether or not simple – that is, single – 1D photoionization models can reproduce the observed emission-line properties of the LyC emitters, and if so, how well they do this. Comparing grids of photoionization models with observational diagrams of approximately 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., fesc ≲ 10%) in galaxies with high [O III]/[O II], because 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 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) 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 these latter authors invoke a contribution of hard radiation from active galactic nuclei (AGNs) or radiative shocks. However, no consistent model allowing for the observed leakage and escape fraction, including the geometrical constraints derived from the UV lines (as mentioned above), and explaining the intensities of the major emission lines has yet been made.

To progress further on these issues, here we 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 nonleaking galaxies. We then compare the data to standard one-zone photoionization models to see the successes and failures of these models. Finally, because 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 into 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 dependency on geometrical factors, leaking, and other factors. Hopefully, these models will lead to the development of a new diagnostic of LyC leakage using rest-frame optical emission lines that is 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 in order 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 representations (e.g., Stasińska & Schaerer 1999; Cormier et al. 2012, 2019; Lebouteiller et al. 2017). Here we 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. Section 3 shows the main optical emission line diagrams used to compare to “classical” single-component photoionization models (both density- and ionization-bounded). We then construct two-component photoionization models including LyC leaking emission and compare them with 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.

2. Optical emission line properties of LyC emitters

In the present work we study eleven z ∼ 0.3 − 0.4 Lyman continuum emitters discovered with the COS spectrograph on-board HST and whose properties have been amply discussed in several papers (Izotov et al. 2016a,b, 2018a,b; Schaerer et al. 2016, 2018; Verhamme et al. 2017; Gazagnes et al. 2018; Chisholm 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.

These sources represent currently the best sample of local LyC emitters, for which a wide range of observations is available, including LyC, UV, and optical spectroscopy. Furthermore, Gazagnes et al. (2018, 2020) provided 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. Here we 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.

2.1. Observational data

2.1.1. 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,b, 2018a,b). We recently obtained new VLT spectra of 5 of the 11 leakers (J0925+1403, J0901+2119, J1011+1947, J1154+2443 and J1442−0209) with the X-shooter 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.

2.1.2. 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 (Guseva et al. 2019), and for which the [O III] λ4363 line is detected with an accuracy better than 4σ, allowing direct abundance determinations using the Te 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.

2.2. 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, that is, 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 for example by a high [O III]/Hβ ratio. However, at a finer level, the LyC emitters show some interesting distinctive properties, which we now discuss.

thumbnail Fig. 1.

O3Hβ–S2Hα diagram of the LyC leakers (open circles) and comparison sample with metallicities 12 + log(O/H) = 7.6 − 8.2 (gray dots). The plain blue line shows the ridge line of our SDSS sample, and the black dashed line is the ridge line from Wang et al. (2019) extrapolated to high [O III]/Hβ values. Many leakers, especially the strongest ones, show a deficiency in [S II] with respect to this latter ridge line.

thumbnail Fig. 2.

O32−O13 diagram.

2.2.1. A [S II] deficiency in LyC emitters

In Fig. 1 we show a classical diagnostic (O3Hβ–S2Hα) involving the [S II]/Hα ratio, which was used by Wang et al. (2019) to select LyC-emitter candidates among more metal-rich and therefore lower excitation (lower [O III]/Hβ) sources. Clearly several, if not a majority, of the eleven confirmed leakers show fairly weak [S II] λλ6716,6731 emission, or [S II]/Hα ratios close to or below the typical value for galaxies from our comparison sample at the same [O III] λ5007/Hβ ratio. 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. This 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 ten galaxies. Each bin was then fitted by a Gaussian distribution of a center x0 and a standard deviation σ. The ridge line defined by the x0 was then fitted with an eight-order polynomial,


where y = [O III] λ5007/Hβ, and the coefficients αn are reported in Table 1. Table 2 presents the coefficients of the ridge lines derived in other diagnostic plots x vs. y = [O III]/[O II].

Table 1.

Coefficients of the ridge lines in classical diagrams (Figs. 1, 4 and 8), determined for the entire galaxy sample and described by Eq. (1).

Table 2.

Coefficients of the ridge lines in oxygen and sulfur diagrams (Figs. 2, 5 and 9), described by Eq. (1).

For subsequent discussions, we also calculate Δ[S II], the lateral shift from the ridge line,


to quantify deviations of individual galaxies from the average observed line ratio of the entire comparison sample. With this definition Δ[S II] < 0 indicates a deficit in [S II]/Hα with respect to the average at the same [O III] λ5007/Hβ ratio. We note that the sulfur abundance ratio (S/O) is constant and shows a fairly small scatter over a wide range of metallicities, even beyond the range of O/H considered here (see, e.g., 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, five of our LyC emitters are within 1σ of the ridge line, whereas the rest are systematically offset to lower [S II]/Hα ratios, that is, they show an [S II] deficit, an empirical 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 are discussed below, providing a more quantitative description of the relationship (Sect. 4.4).

2.2.2. 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 [O III]/[O II] and [O I]/[O III], as shown in Fig. 2 (denoted the O32−O13 diagram hereafter). Indeed, naively one would expect to find a deficit of [O I] emission in LyC emitters as they should be highly ionized and possibly lacking neutral regions where [O I] is emitted. This is also shown by simple density-bounded photoionization models (see Figs. 5 and 6, and Stasińska et al. 2015). However, as also noted by Plat et al. (2019), [O I] is detected in 9 out of the 11 leakers of this study, and the observed [O I]/[O III] ratios of the LyC emitters are surprisingly high, shifted to even higher than average values when compared to the ridge line defined by the comparison sample in 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.

2.2.3. 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 Cols. 2 and 3 in Table 1). In this plane, we see a clear shift from the SDSS comparison sample of presumably nonleaking galaxies. Interestingly, the three most leaking galaxies (the three reddest circles, fesc ≳ 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 fesc ≲ 10%) exhibit a moderate deficiency in [S II] (−0.2 < Δ[S II] < 0) and [O I] excesses (Δ[O I] ≳ 0). However, we note that three galaxies in our sample with a large deficit in [S II] (Δ[S II] ≲ −0.3) have intermediate (J1011+1947, fesc = 11,4%) to low escape fractions (J1248+4259, J1333+6246, labeled on the plot). One of them (J1333+6246) strikingly stands apart with a nondetection of the [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.

thumbnail Fig. 3.

Δ[O I] vs. Δ[S II] derived from classical diagrams O3Hβ–O1Hα and O3Hβ–S2Hα. The three strongest leakers (J1256+4509, J1154+2443 and J1243+4646, by increasing order of fesc) are reprented by orange to brown circles in the lower part of the plot. The other galaxies discussed in this paper are labeled. The LyC emitters show a clear displacement in this plane compared to the comparison sample.

2.2.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. Likely related to this, Guseva et al. (2020) 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.

thumbnail Fig. 4.

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).

In another diagram discussed further below, relating [O III]/[O II] to [S II] λλ6716,6731/[S III] λ9068 (Fig. 5), we note that four out of five LyC emitters, for which we have near-IR measurements of [S III] λ9068, appear also shifted from the average [S II] λλ6716,6731/[S III] λ9068 ratio in a similar fashion as for [O I]/[O III] λλ4959,5007. More measurements will show if this behavior is systematic for LyC leakers or not.

thumbnail Fig. 5.

O32−O13 and O32−S23 diagnostics with BOND models. Photoionization models are represented by the colored symbols: the colorbar represents the variation of the escape fraction and the size of symbols is proportional to log(U). The shapes of model dependencies reflect the variations either with metallicity or stellar age and are shown and labeled in the top and bottom panels, respectively.

Finally, Izotov et al. (2017) proposed that LyC-emitting galaxies show peculiar line ratios of weak He I lines. As shown by Guseva et al. (2020), our X-shooter observations confirm the trend suggested by this study.

To gain further insight into the observed behavior of the emission line ratios we now proceed to comparisons with photoionization models, first with simple and then more complex geometries.

3. 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.

3.1. The models: BOND models from the 3MdB database

To examine the power of simple 1D photoionization models, we 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 specificities of the BOND code is that it is independent of any assumption on relations between O/H and N/O. Therefore, 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 (SED) 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 (IMF) 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 also studied the effect of LyC leakage. In the current study, we use models from the 3MdB_17 database2 (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 a fixed density of 100 cm−3, and 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., [O I], [S II]) with respect to higher ionization potential species (e.g., [O III], [S III]). Specifically, this would induce a lateral shift to the left-hand side in the O32−O13 and O32−S23 diagnostics plots (Figs. 5 and 9), which we examine in Sect. 3.2.

To describe density-bounded models, we used a parameter Hβfrac which cuts the nebula at different radii corresponding to a range of 10−100% of the total Hβ luminosity3. 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, density-bounded model to the Strömgren radius of its ionization-bounded analog. We directly translate Hβfrac into the corresponding LyC escape fraction by defining fesc(Hβ) = 1 − Hβfrac which 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.

3.2. 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 above-mentioned O32−O13, and other diagrams.

3.2.1. Classical diagrams

Figure 4 shows the LyC emitters, the comparison sample, and the one-zone models in the classical diagrams from Veilleux & Osterbrock (1987), including the [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). 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 (fesc = 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, that is, models with fesc = 0 or clearly only for low values of fesc <  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 strongly decreases 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].

3.2.2. 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 therefore the different zones of the ionized ISM in the LyC emitters. By focusing on different ionization stages of the same element, we also avoid dependencies on elemental abundance ratios and metallicity, at least to first order. These O32−O13 and O32−S23 diagnostics highlight two main points that motivate the introduction of more complex models in the Sect. 4:

(1) Although the oxygen line ratios (O32−O13) predicted by simple ionization-bounded models reproduce the SDSS sequence with ionizing spectra of young populations (∼1 − 2 Myr) quite well, the models are not able to account for galaxies with a certain [O I] λ6300 excess, that is, a higher-than-average [O I]/[O III] ratio at a given [O III]/[O II] 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 from the fact that the models do not take into account a DIG component, which can significantly contribute to the emission of [O I] and [S II] (Sanders et al. 2017, 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) because of problems with atomic data for sulfur (Izotov et al. 2006; Kewley & Dopita 2002; Kewley et al. 2019), or others. The sulfur line ratio [S II]/[S III] is also sensitive to more subtle effects including stellar winds, outflows, and photoevaporation of photo-dissociation regions (PDRs) 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 represents the [O III]/[O II] and [O III]/[O I] ratios predicted by the models as a function of their escape fraction. This diagram emphasizes the discrepancy between [O III]/[O I] and fesc, showing that [O III]/[O I] is overestimated by several orders of magnitude. This clearly shows that moderate to high escape fractions cannot be produced by simple density-bounded one-zone models without significantly decreasing the weak line emission of [O I] and [S II].

thumbnail Fig. 6.

Observed and predicted [O III]/[O II] (left) and [O III]/[O I] (right) line ratios as a function of the LyC escape fraction, comparing the observed LyC emitters and the one-zone photoionization models with log(U)≳ − 2.5. We note the very high range of the predicted [O III]/[O I] ratio.

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 observed in LyC emitters. On the other hand, while ionization-bounded models can describe the overall properties of most galaxies (the comparison sample) fairly well – 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 for 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 therefore examined composite photoionization models based on the combination of two distinct zones with different properties.

4. Composite “two-zone” photoionization models

4.1. Motivation and method

As demonstrated above (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) from Sect. 3.2.2 (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) from Sect. 3.2.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 model, was also suggested based on other observations. Gazagnes et al. (2018, 2020) and Chisholm et al. (2018) 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 that are 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 purposes, their model is basically equivalent to our two-zone model. As we 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.

thumbnail Fig. 7.

Schematic illustration of our two-zone photoionization models for LyC emitters, combining an ionization region with density-bounded regions chosen to have different ionization parameters (“high” and “low” U). In the “density-bounded channels” scenario, ionizing photons escape from the low-U region, while in the “density-bounded galaxy” scenario, they escape from the high-U region with fesc = 0 for the nonleaking region and fesc between 0% and 80% for the leaking region.

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 Eq. (4) in Stasińska et al. 2015). For ionization-bounded 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, α(Te) 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 s1/3 for Te = 104 K), which means that both low-density and clumpy mediums (ϵ ≪ 1) will 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 ionization-bounded 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 Eq. (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 Sect. 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-ionization parameter Uh and one with a low-ionization parameter Ul. The high- and low-ionization parameters were chosen “by hand” so that the nonleaking combined models accurately reproduce 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 (Figs. 5 and 9). We used Uh = 10−1 and Ul = 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. As 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 Qh and Ql, respectively. This is simply achieved by scaling the predicted line luminosities of the Ul model with the factor s = Qh/Ql. The two components are then combined with relative weights, ω for the Uh component and (1-ω) for the Ul 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 fesc of the two-zone model, is computed consistently from the relative weight of zones ω, by defining


As only one component is leaking in each scenario, Eq. (6) simply becomes in the DBC scenario and in the DBG scenario. In the absence of leakage (fesc = 0), the two scenarios are trivially identical.

The resulting emission line ratios predicted from our two-zone models 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). Nonleaking models are represented by the yellow symbols with 0% escaping photons. In this case the two scenarios (DBC and DBG) are identical by construction, and their symbols (yellow triangles on the first row, and yellow diamonds on the second row) are thus superposed. As mentioned before, we chose the ionization parameter of each component such that the combined nonleaking models adequately cover the upper part of the SDSS sequence.

thumbnail Fig. 8.

Classical diagrams (same as Fig. 4) showing the predictions from the two-zone photoionization models. Upper and lower panels: represent DBC and DBG models, respectively.

thumbnail Fig. 9.

O32−O13 and O32−S23 diagrams (same as Fig. 5) showing the predictions from the two-zone photoionization models. Upper and lower panels: DBC and DBG models, respectively.

4.2. [O I] and [S II] λλ6716,6731 emission in normal galaxies explained

As we can see from Fig. 8, the nonleaking 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 fesc) can be covered with the two-zone models combining two different ionization parameters. In particular, there seems to be no difficulty in explaining 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 Figs. 5 and 9). The inclusion of a second ISM component with a lower ionization parameter therefore seems both essential and sufficient to solve the above-mentioned discrepancies of classical models when compared to normal galaxies.

4.3. 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 Figs. 1 and 2), which characterizes the leakers, and which cannot be understood by one-zone models (see Sects. 2.2 and 3.2). Furthermore, the two-zone models account for the escape of LyC photons, and are thus consistent with this important observational property. Below we further explore the consistency of different line ratios and the extent to which 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 fesc 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 density-bounded 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 fesc 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, that is, the one leaking in the DBC scenario (upper panel). This explains why these line ratios are barely affected in the DBG scenario (lower panel), even for high escape fractions. On the contrary, leakage in the low-ionization component (DBC scenario) rapidly decreases [S II] and [O I]: only low escape fractions (fesc ≲ 10%) are compatible with the observed [O I] emission of the leakers in this scenario.

The same trend is visible on the O32−O13 plot (Fig. 9). The SDSS sequence seems bounded by fesc ∼ 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 strongly decreasing 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 leaky galaxy J1154+2443 (red circle, fesc = 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 high excitation, as indicated by their high O32 ratio: decreasing the contribution ω of the Uhigh region below a certain threshold produces two-zone models that are dominated by the low ionization (Ulow) component, and hence are unable to reproduce the observed line ratios of the leakers.

4.4. Links with fesc

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. Figure 10 also shows that LyC emitters with very high fesc ≳ 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.

thumbnail Fig. 10.

Same as Fig. 6 showing the predictions for the two-zone photoionization models.

Figure 11 (left) illustrates the evolution of fesc 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 nonleakers), as also pointed out by Wang et al. (2019). However, using Δ[S II] to determine the escape fraction is 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 a 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 fesc ≳ 30% all exhibit a 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 percent. This could mean that the observed fesc for these galaxies is underestimated, for example due to the fact 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 whether or not any of these sources indeed emit 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 fesc from the [S II] deficit remains challenging.

thumbnail Fig. 11.

LyC escape fraction as a function of the deficiency Δ[S II] in [S II] (left) and Δ[O I] (right). The gray-shaded areas represent dispersion of the SDSS comparison sample at 1, 2, and 3σ. The three strongest leakers (J1256+4509, J1154+2443 and J1243+4646, in increasing order of fesc) are represented by orange to brown circles in the upper part of the plots. The other galaxies that are discussed in Sect. 4.4 are labeled in blue.

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 most of the galaxies with small fesc reside, including the [S II]-deficient galaxies selected by Wang et al. (2019; green color along the contours of circles). However, we note that the two galaxies with the strongest [S II] emission (Δ[S II] close to zero) stand outside of the contours of the DBC scenario. This indicates that leakage powered exclusively by channels of low ionization is unable to simultaneously produce 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, Sect. 2.2). Positive Δ[O I] corresponds to galaxies showing an excess in [O I]/Hα compared to the SDSS median line, while a negative Δ[O I] marks an [O I] deficiency. Interestingly, we notice a visible trend of fesc relative to Δ[O I] in the observations: while the most [O I]-deficient galaxies and the two undetected galaxies (|Δ[O I]| < 1σ) exhibit the largest escape fractions, the observed fesc decreases for increasing Δ[O I]. This could indicate that [O I] emission, 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. However, we notice that two galaxies (J1248+4259 and J0901+2119) that are [O I]-deficient exhibit very small escape fractions (fesc ∼ 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 fesc 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. For the same reasons as for [S II] emission, leakage through low-ionization channels tends to be associated with low [O I] emission, that is, Δ[O I] ≲ 0. The galaxies that are not [O I] deficient are not well reproduced by the DBC scenario, especially the galaxy standing on the extreme right-hand side whose [O I] emission exceeds that of all the models considered here. This might suggest that the physics of [O I] emission is not well taken into account in our two-component models, and might require a boost from an additional component (see discussion in Sect. 5). However, we note that apart from the systematic offset toward lower [O I] emission, the DBC scenario qualitatively reproduces the observed trend of increasing fesc with decreasing Δ[O I]. The DBG scenario accurately reproduces all the Δ[O I] and fesc of the leakers (except for the strongest [O I] emitter) but Δ[O I] is independent from the escape fraction in this scenario; the observed variations in Δ[O I] only come from the decrease of [O III] which vertically moves the models closer to the ridge line but without affecting [O I] emission.

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 fesc, and other diagnostics sensitive to density boundaries are to be preferred (see, e.g., Figs. 5 and 9). It seems that [O III]/[O II], Δ[S II], and [O I] excess are good hints for LyC leakage, but there still seems to be other independent dimensions, which probably incorporate several physical parameters such as the nature 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).

5. Discussion

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

5.1. 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 escape 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 tracers, which 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 subgrid of our combined models, which includes the ones presented in Sect. 4. We defined a four-dimensional grid with the following free parameters: Uhigh and Ulow (varied from 10−4 to 10−1), the weight of the component of high ionization ω, and the escape fraction of the leaking component (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 prescription for 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 (Gazagnes et al. 2018; Chisholm et al. 2018; Schaerer et al. 2018). Hereafter, we summarize some interesting trends that merit further examination:

– The models minimizing χ2 always favor two distinct values for the ionization parameters Uhigh and Ulow. Inspection of their marginalized probability distribution functions shows that typically Uhigh ∼ 10−1 − 10−2, and Ulow ∼ 10−3 − 10−4. We interpret this as 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 that of 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 structure of the galaxy unchanged. However, we note that the two scenarios produce very different covering fractions. We note that for galaxies with a low escape fraction (fesc ≲ 10%), the ω obtained from the best-fitting model has a value comparable to the covering fractions, Cf, obtained from analysis of absorption lines (Gazagnes et al. 2018, 2020; Chisholm et al. 2018). 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 Cf ∼ (1 − ω), which is a better fit for some of these galaxies. On the other hand, for strong leakers (fesc ≳ 10%), we find Cf <  ω. 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 therefore overly 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. The best-fit models predict fesc ≈ 0 for all the weak leakers but one, and nonzero escape fractions (fesc ∼ 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 fesc is a factor of approximately two below the observed value.

Briefly, 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 robustly distinguishing between the two scenarios based only on emission lines is not easy. Furthermore, accurate predictions of fesc 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.

5.2. Low ionization line diagnostics: progress and caveats

In Sect. 3, we discussed the inadequacy of simple models to simultaneously reproduce 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 latter 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), Mingozzi et al. (2020) 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. However, the exact nature, spatial distribution, and other properties of such a component remain 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- and ionization-bounded H II regions can be found that reproduce the three ionization stages of oxygen. How these models would affect the sulfur lines, and whether or not 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).

Regardless of the exact origin of emission of the low-ionization lines (such as [O I] λ6300 and [S II] λλ6716,6731) in LyC emitting galaxies, the UV observations of these galaxies clearly reveal the presence of both low-H I-column-density (NHI) lines of sight towards the observer allowing the escape of LyC radiation and lines of sight with high NHI, where other neutral and low-ionization metal absorption lines are also formed (Gazagnes et al. 2018, 2020; Chisholm et al. 2018). Therefore, their ISM contains at least two “zones” or regions with different conditions, even if we do not know the exact spatial scales on which this happens. Furthermore, it is very plausible that properties that 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, Uhigh and Ulow, for the LyC leakers. Whether or not 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 1D photoionization models cannot reproduce all the observables.

5.3. 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 star-forming galaxies, they are generally 1D 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 and so on. Such models, which typically assume a single source of ionizing radiation and a very simple geometry, in contrast 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 highly complex models have so far been used for comparisons with observations. Stasińska et al. (2015) for example explored the combination of stellar photoionization and other emission sources (AGN, shocks, and others) with variable but arbitrary contributions, and confronted them with a large sample of star-forming galaxies from the SDSS. Cormier et al. (2012, 2019) used one- and two-component models to consistently analyze emission lines originating in the ionized and neutral ISM, as well as in adjacent PDRs. 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. These simulations 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 nonleakers, which show 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 from those observed in low-z LyC emitters and in z ∼ 2 − 3 galaxies. Understanding the origin and importance of such a discrepancy is nevertheless not straightforward. Different and complementary modeling approaches should help to develop more robust diagnostics and should therefore provide better insight into 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 60 galaxies at z ∼ 0.2 − 0.4 with direct LyC observations. This sample of confirmed LyC emitters and also nonLyC 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 will ideally allow 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.

6. Conclusion

We analyzed the main optical emission lines of low-z LyC-emitting galaxies (leakers) discovered recently with the HST (Izotov et al. 2016a,b, 2018a,b; Wang et al. 2019) in the classical excitation diagrams involving the [O III]/Hβ, [N II]/Hα, [S II]/Hα, and [O I]/Hα line ratios, as well as other diagnostic diagrams involving three ionization stages of oxygen (using [O III], [O II], and [O I]) and two of sulfur (using [S II] and [S III] observations). In comparison with the SDSS sample of star-forming galaxies with similar metallicities (12 + log(O/H)∼7.6 − 8.2), we found several peculiarities of the LyC emitters, including an excess of the [O I]/[O III] ratio for a fraction of the leakers (as already pointed out by Plat et al. 2019), and an [O I] deficiency for those with the highest LyC escape fractions (fesc ≳ 20%). We also confirm a deficit in [S II] found recently by Wang et al. (2019), and suggest that combined measurements of the [O I] and [S II] lines could provide a good indicator for strong LyC emission.

We then compared the observations to simple one-zone and two-zone photoionization models accounting for the escape of ionizing radiation in order 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; this underprediction is strong for density-bounded models with LyC escape, and also exists for ionization-bounded models (see Sect. 3), although it is less strong.

– Two-zone models, combining regions with a high- and low-ionization parameter (Uhigh, Ulow), where one of which is density-bounded, 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 fesc <  10% can be explained by the presence of low-column-density, 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 fesc, 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 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. (2018, 2020). Our finding of two different scenarios, favored by LyC emitters with a low or 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-sector modeling, including some density-bounded components to successfully capture the complexity of the ISM structure that leads to photon leakage. We also suggest that low-excitation 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.


In all the article we will use [O III] for [O III] λ5007, [O II] for [O II] λλ3726,3728, [O I] for [O I] λ6300, [N II] for [N II] λ6584, [S III] for [S III] λ9068 and [S II] for [S II] λλ6716,6731.


The 3MdB website describes the 3MdB and 3MdB_17 databases at


Since the study of Vale Asari et al. (2016), the grid in Hβfrac has been refined and is now available for steps of 10%.


We acknowledge the use of photoionization models from Stephane de Barros at an earlier stage of this work. YII and NGG acknowledge support from the National Academy of Sciences of Ukraine by its priority project No.0120U100935 “Fundamental properties of the matter in the relativistic collisions of nuclei and in the early Universe”. JMV acknowledges support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the IAA (SEV-2017-0709) and project AYA2016-79724-C4-4-P. RA acknowledges support from FONDECYT Regular Grant 1202007.


  1. Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42 [NASA ADS] [CrossRef] [Google Scholar]
  2. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Bassett, R., Ryan-Weber, E. V., Cooke, J., et al. 2019, MNRAS, 483, 5223 [CrossRef] [Google Scholar]
  5. Bian, F., Fan, X., McGreer, I., Cai, Z., & Jiang, L. 2017, ApJ, 837, L12 [Google Scholar]
  6. Chabrier, G. 2003, ApJ, 586, L133 [Google Scholar]
  7. Chisholm, J., Gazagnes, S., Schaerer, D., et al. 2018, A&A, 616, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Cormier, D., Lebouteiller, V., Madden, S. C., et al. 2012, A&A, 548, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. de Barros, S., Vanzella, E., Amorín, R., et al. 2016, A&A, 585, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385 [NASA ADS] [Google Scholar]
  12. Fletcher, T. J., Tang, M., Robertson, B. E., et al. 2019, ApJ, 878, 87 [Google Scholar]
  13. Garnett, D. R. 1989, ApJ, 345, 282 [NASA ADS] [CrossRef] [Google Scholar]
  14. Gazagnes, S., Chisholm, J., Schaerer, D., et al. 2018, A&A, 616, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Gazagnes, S., Chisholm, J., Schaerer, D., Verhamme, A., & Izotov, Y. 2020, A&A, 639, A85 [CrossRef] [EDP Sciences] [Google Scholar]
  16. Guseva, N. G., Izotov, Y. I., Fricke, K. J., & Henkel, C. 2019, A&A, 624, A21 [CrossRef] [EDP Sciences] [Google Scholar]
  17. Guseva, N. G., Izotov, Y. I., Schaerer, D., et al. 2020, MNRAS, 497, 4293 [CrossRef] [Google Scholar]
  18. Hirschmann, M., Charlot, S., Feltre, A., et al. 2017, MNRAS, 472, 2468 [NASA ADS] [CrossRef] [Google Scholar]
  19. Hogarth, L., Amorín, R., Vílchez, J. M., et al. 2020, MNRAS, 494, 3541 [CrossRef] [Google Scholar]
  20. Izotov, Y. I., Stasińska, G., Meynet, G., Guseva, N. G., & Thuan, T. X. 2006, A&A, 448, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016a, Nature, 529, 178 [NASA ADS] [CrossRef] [Google Scholar]
  22. Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683 [NASA ADS] [CrossRef] [Google Scholar]
  23. Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2017, MNRAS, 471, 548 [NASA ADS] [CrossRef] [Google Scholar]
  24. Izotov, Y. I., Schaerer, D., Worseck, G., et al. 2018a, MNRAS, 474, 4514 [Google Scholar]
  25. Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018b, MNRAS, 478, 4851 [Google Scholar]
  26. Izotov, Y. I., Thuan, T. X., & Guseva, N. G. 2019, MNRAS, 483, 5491 [NASA ADS] [CrossRef] [Google Scholar]
  27. Jaskot, A. E., & Oey, M. S. 2013, ApJ, 766, 91 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kakiichi, K., & Gronke, M. 2019, ApJ, submitted [arXiv:1905.02480] [Google Scholar]
  29. Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902 [Google Scholar]
  30. Katz, H., Ďurovčíková, D., Kimm, T., et al. 2020, MNRAS, 498, 164 [CrossRef] [Google Scholar]
  31. Kehrig, C., Vílchez, J. M., Telles, E., Cuisinier, F., & Pérez-Montero, E. 2006, A&A, 457, 477 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kewley, L. J., Heisler, C. A., Dopita, M. A., & Lumsden, S. 2001, ApJS, 132, 37 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511 [Google Scholar]
  35. Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2019, ApJ, 887, 80 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lebouteiller, V., Péquignot, D., Cormier, D., et al. 2017, A&A, 602, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Leitherer, C., Hernandez, S., Lee, J. C., & Oey, M. S. 2016, ApJ, 823, 64 [NASA ADS] [CrossRef] [Google Scholar]
  38. McLeod, A. F., Dale, J. E., Ginsburg, A., et al. 2015, MNRAS, 450, 1057 [NASA ADS] [CrossRef] [Google Scholar]
  39. Mingozzi, M., Belfiore, F., Cresci, G., et al. 2020, A&A, 636, A42 [CrossRef] [EDP Sciences] [Google Scholar]
  40. Mollá, M., García-Vargas, M. L., & Bressan, A. 2009, MNRAS, 398, 451 [NASA ADS] [CrossRef] [Google Scholar]
  41. Morisset, C., Delgado-Inglada, G., & Flores-Fajardo, N. 2015, Rev. Mex. Astron. Astrofis., 51, 103 [Google Scholar]
  42. Naidu, R. P., Forrest, B., Oesch, P. A., Tran, K.-V. H., & Holden, B. P. 2018, MNRAS, 478, 791 [CrossRef] [Google Scholar]
  43. Nakajima, K., & Ouchi, M. 2014, MNRAS, 442, 900 [NASA ADS] [CrossRef] [Google Scholar]
  44. Nakajima, K., Ellis, R. S., Robertson, B. E., Tang, M., & Stark, D. P. 2020, ApJ, 889, 161 [CrossRef] [Google Scholar]
  45. Oey, M. S., Meurer, G. R., Yelda, S., et al. 2007, ApJ, 661, 801 [NASA ADS] [CrossRef] [Google Scholar]
  46. Pellegrini, E. W., Rahner, D., Reissl, S., et al. 2020, MNRAS, 496, 339 [CrossRef] [Google Scholar]
  47. Pérez-Montero, E., & Vílchez, J. M. 2009, MNRAS, 400, 1721 [NASA ADS] [CrossRef] [Google Scholar]
  48. Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978 [NASA ADS] [CrossRef] [Google Scholar]
  49. Sanders, R. L., Shapley, A. E., Zhang, K., & Yan, R. 2017, ApJ, 850, 136 [NASA ADS] [CrossRef] [Google Scholar]
  50. Sanders, R. L., Jones, T., Shapley, A. E., et al. 2020, ApJ, 888, L11 [CrossRef] [Google Scholar]
  51. Schaerer, D., Izotov, Y. I., Verhamme, A., et al. 2016, A&A, 591, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Schaerer, D., Izotov, Y. I., Nakajima, K., et al. 2018, A&A, 616, L14 [CrossRef] [EDP Sciences] [Google Scholar]
  53. Schmitt, H. R. 1998, ApJ, 506, 647 [NASA ADS] [CrossRef] [Google Scholar]
  54. Shapley, A. E., Sanders, R. L., Shao, P., et al. 2019, ApJ, 881, L35 [CrossRef] [Google Scholar]
  55. Shimizu, I., Inoue, A. K., Okamoto, T., & Yoshida, N. 2016, MNRAS, 461, 3563 [CrossRef] [Google Scholar]
  56. Stasińska, G., & Leitherer, C. 1996, ApJS, 107, 661 [Google Scholar]
  57. Stasińska, G., & Schaerer, D. 1999, A&A, 351, 72 [Google Scholar]
  58. Stasińska, G., Cid Fernandes, R., Mateus, A., Sodré, L., & Asari, N. V. 2006, MNRAS, 371, 972 [NASA ADS] [CrossRef] [Google Scholar]
  59. Stasińska, G., Izotov, Y., Morisset, C., & Guseva, N. 2015, A&A, 576, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123 [Google Scholar]
  61. Vale Asari, N., Stasińska, G., Morisset, C., & Cid Fernandes, R. 2016, MNRAS, 460, 1739 [CrossRef] [Google Scholar]
  62. Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41 [NASA ADS] [CrossRef] [Google Scholar]
  63. Vanzella, E., Nonino, M., Cupani, G., et al. 2018, MNRAS, 476, L15 [Google Scholar]
  64. Vanzella, E., Caminha, G. B., Calura, F., et al. 2020, MNRAS, 491, 1093 [NASA ADS] [CrossRef] [Google Scholar]
  65. Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295 [NASA ADS] [CrossRef] [Google Scholar]
  66. Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Wang, B., Heckman, T. M., Leitherer, C., et al. 2019, ApJ, 885, 57 [CrossRef] [Google Scholar]
  69. Westmoquette, M. S., Smith, L. J., Gallagher, J. S., & Walter, F. 2013, MNRAS, 428, 1743 [NASA ADS] [CrossRef] [Google Scholar]
  70. Wilkins, S. M., Lovell, C. C., Fairhurst, C., et al. 2020, MNRAS, 493, 6079 [CrossRef] [Google Scholar]

All Tables

Table 1.

Coefficients of the ridge lines in classical diagrams (Figs. 1, 4 and 8), determined for the entire galaxy sample and described by Eq. (1).

Table 2.

Coefficients of the ridge lines in oxygen and sulfur diagrams (Figs. 2, 5 and 9), described by Eq. (1).

All Figures

thumbnail Fig. 1.

O3Hβ–S2Hα diagram of the LyC leakers (open circles) and comparison sample with metallicities 12 + log(O/H) = 7.6 − 8.2 (gray dots). The plain blue line shows the ridge line of our SDSS sample, and the black dashed line is the ridge line from Wang et al. (2019) extrapolated to high [O III]/Hβ values. Many leakers, especially the strongest ones, show a deficiency in [S II] with respect to this latter ridge line.

In the text
thumbnail Fig. 2.

O32−O13 diagram.

In the text
thumbnail Fig. 3.

Δ[O I] vs. Δ[S II] derived from classical diagrams O3Hβ–O1Hα and O3Hβ–S2Hα. The three strongest leakers (J1256+4509, J1154+2443 and J1243+4646, by increasing order of fesc) are reprented by orange to brown circles in the lower part of the plot. The other galaxies discussed in this paper are labeled. The LyC emitters show a clear displacement in this plane compared to the comparison sample.

In the text
thumbnail Fig. 4.

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).

In the text
thumbnail Fig. 5.

O32−O13 and O32−S23 diagnostics with BOND models. Photoionization models are represented by the colored symbols: the colorbar represents the variation of the escape fraction and the size of symbols is proportional to log(U). The shapes of model dependencies reflect the variations either with metallicity or stellar age and are shown and labeled in the top and bottom panels, respectively.

In the text
thumbnail Fig. 6.

Observed and predicted [O III]/[O II] (left) and [O III]/[O I] (right) line ratios as a function of the LyC escape fraction, comparing the observed LyC emitters and the one-zone photoionization models with log(U)≳ − 2.5. We note the very high range of the predicted [O III]/[O I] ratio.

In the text
thumbnail Fig. 7.

Schematic illustration of our two-zone photoionization models for LyC emitters, combining an ionization region with density-bounded regions chosen to have different ionization parameters (“high” and “low” U). In the “density-bounded channels” scenario, ionizing photons escape from the low-U region, while in the “density-bounded galaxy” scenario, they escape from the high-U region with fesc = 0 for the nonleaking region and fesc between 0% and 80% for the leaking region.

In the text
thumbnail Fig. 8.

Classical diagrams (same as Fig. 4) showing the predictions from the two-zone photoionization models. Upper and lower panels: represent DBC and DBG models, respectively.

In the text
thumbnail Fig. 9.

O32−O13 and O32−S23 diagrams (same as Fig. 5) showing the predictions from the two-zone photoionization models. Upper and lower panels: DBC and DBG models, respectively.

In the text
thumbnail Fig. 10.

Same as Fig. 6 showing the predictions for the two-zone photoionization models.

In the text
thumbnail Fig. 11.

LyC escape fraction as a function of the deficiency Δ[S II] in [S II] (left) and Δ[O I] (right). The gray-shaded areas represent dispersion of the SDSS comparison sample at 1, 2, and 3σ. The three strongest leakers (J1256+4509, J1154+2443 and J1243+4646, in increasing order of fesc) are represented by orange to brown circles in the upper part of the plots. The other galaxies that are discussed in Sect. 4.4 are labeled in blue.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.