Reconciling escape fractions and observed line emission in Lyman-continuum-leaking galaxies
L. Ramambason, D. Schaerer, G. Stasińska, Y. I. Izotov, N. G. Guseva, J.M Vílchez, R. Amorín, C. Morisset
AAstronomy & Astrophysics manuscript no. 38634corr © ESO 2020October 9, 2020
Reconciling escape fractions and observed line emission inLyman-continuum-leaking galaxies (cid:63)
L. Ramambason , , D. Schaerer , , G. Stasi´nska , Y. I. Izotov , N. G. Guseva , J.M Vílchez , R. Amorín , , C.Morisset Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland AIM, CEA, Université Paris-Saclay, Université Paris Diderot, Université de Paris, F-91191, Gif-sur-Yvette, France. CNRS, IRAP, 14 Avenue E. Belin, 31400 Toulouse, France LUTH, Observatoire de Meudon, 92195 Meudon Cedex, France Bogolyubov Institute for Theoretical Physics, National Academy of Sciences of Ukraine, 14-b Metrolohichna str., Kyiv, 03143,Ukraine Instituto de Astrofísica de Andalucía, CSIC, Apartado de correos 3004, E-18080 Granada, Spain Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile Instituto de Astronomía, Universidad Nacional Autónoma de México, AP 106, 22800 Ensenada, B. C., MexicoReceived 11 June 2020 / Accepted 17 September 2020
ABSTRACT
Context.
Finding and elucidating the properties of Lyman-continuum(LyC)-emitting galaxies is an important step in improving ourunderstanding of cosmic reionization.
Aims.
Although the z ∼ . − . Methods.
To address these questions we construct one- and two-zone photoionization models accounting for the observed LyC escape,which we compare to the observed emission line measurements. The main diagnostics used include lines of [O iii ], [O ii ], and [O i ]plus sulfur lines ([S ii ], [S iii ]) and a nitrogen line ([N ii ]), which probe regions of di ff erent ionization in the ISM. Results.
We find that single (one-zone) density-bounded photoionization models cannot reproduce the emission lines of the LyCleakers, as pointed out by earlier studies, because they systematically underpredict the lines of species of low ionization potential, suchas [O i ] and [S ii ]. Introducing a two-zone model, with di ff ering ionization parameter and a variable covering fraction and where oneof the zones is density-bounded, we show that the observed emission line ratios of the LyC emitters are well reproduced. Furthermore,our model yields LyC escape fractions, which are in fair agreement with the observations and independent measurements. The [O i ] λ f esc > ∼ i ] λ ii ] λλ ii ] λλ f esc < ∼ Conclusions.
We conclude that two-zone photoionization models are su ffi cient and required to explain the observed emission lineproperties of z ∼ . − . Key words.
Galaxies: starburst – Galaxies: high-redshift – Cosmology: dark ages, reionization, first stars – Ultraviolet: galaxies
1. Introduction
Significant progress has been made in recent years in searchesfor and studies of star-forming galaxies, which emit copiousamounts of Lyman continuum (LyC) radiation. This ionizingradiation can escape the interstellar medium (ISM) and ionizethe surroundings, thus contributing to intergalactic ionizing ra-diation. If a su ffi cient population of such galaxies exists in theearly Universe, their presence could explain cosmic reionization, (cid:63) Based on data obtained with the European Southern ObservatoryVery Large Telescope, Paranal, Chile, under 0102.B-0942(A) which is known to have started shortly after the Big Bang and tohave ended ∼ z ∼ ff orts to search for LyC emitters or LyC-leaking galaxies (abbreviated to “leakers”), several groups haveidentified such galaxies from low redshift up to z ∼ f esc , ranging from a few percent up to ∼
72 % in some cases (Izotov et al. 2018a; Vanzella et al. 2020).Although di ff erent methods and criteria have been used to search Article number, page 1 of 16 a r X i v : . [ a s t r o - ph . GA ] O c t & A proofs: manuscript no. 38634corr for these galaxies, some of the distinct observational featuresshared by many LyC emitters are strong Ly α emission with adouble-peaked Ly α line profile (see Verhamme et al. 2015, 2017;Vanzella et al. 2020) and intense optical emission lines with ahigh [O iii ] / [O ii ] 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 valuesof [O iii ] / [O ii ] could pinpoint density-bounded galaxies, thatis, LyC-emitting galaxies, this simple prediction was put to testby Izotov and collaborators, who selected compact star-forminggalaxies with intense emission lines and [O iii ] / [O ii ] >
5, andwho were thus able to find 11 LyC emitters at z ∼ . − . z > f esc ∼ −
72 %), and are therefore good analogs ofthe sources of cosmic reionization (cf. Schaerer et al. 2016).Despite the fact that high [O iii ] / [O ii ] was used as a basicselection criterion for these searches, it nevertheless turns outthat f esc does not well correlate with the [O iii ] / [O ii ] ratio (seeIzotov et al. 2018b; Naidu et al. 2018; Bassett et al. 2019; Naka-jima 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 todetect LyC emission from two low- z galaxies with fairly low ex-citation ([O iii ] / [O ii ] ∼ . − . ii ] deficiency and theLyC escape fraction of the known leakers. These findings are notentirely surprising, as it is well known that the [O iii ] / [O ii ] ra-tio depends also on several other factors, including in particularthe ionization parameter and metallicity. Whether or not density-bounded models are applicable to known LyC emitters thereforeremains to be examined quantitatively.On the other hand, the detailed analysis of UV absorptionlines (including the Lyman series lines of hydrogen and otherlow-ionization lines) of the z ∼ . − . f esc < ∼ iii ] / [O ii ], be-cause the observed emission line ratios are not compatible withdensity-bounded H ii regions and therefore with a significant In all the article we will use [O iii ] for [O iii ] λ ii ] for [O ii ] λλ i ] for [O i ] λ ii ] for [N ii ] λ iii ] for[S iii ] λ ii ] for [S ii ] λλ ionizing photon escape. The same authors also stressed the im-portance of considering high- and low-excitation lines. Further-more, they conjectured that some leakage might occur becauseof a nebular covering factor smaller than unity, and that [O i ]emission could come from optically thick clumps embedded in adensity-bounded medium (Stasi´nska et al. 2015). Related to this,Plat et al. (2019) recently noted that low- z LyC emitters show ahigher [O i ] / [O iii ] ratio than their ionization-bounded photoion-ization models, in contrast to expectations from density-boundedmodels where [O i ] / [O iii ] should be reduced. To solve this dis-crepancy these latter authors invoke a contribution of hard ra-diation from active galactic nuclei (AGNs) or radiative shocks.However, no consistent model allowing for the observed leak-age and escape fraction, including the geometrical constraintsderived from the UV lines (as mentioned above), and explainingthe intensities of the major emission lines has yet been made.To progress further on these issues, here we compute bothone- and two-zone photoionization models with metallicities tai-lored to those of the observed LyC emitters to quantitatively ex-amine how such models compare to the various detailed obser-vational constraints, including measured LyC escape fractions.First, we re-examine the most important observed emission lineratios of the z ∼ . − . ff erent optical emission lines,their dependency on geometrical factors, leaking, and other fac-tors. Hopefully, these models will lead to the development of anew diagnostic of LyC leakage using rest-frame optical emissionlines that is applicable to other observations, especially of high- z galaxies, for which many spectra will be obtained in the nearfuture with the James Webb Space Telescope (JWST).The need for multi-zone (or multi-component) models is notonly of importance for LyC emitters. Indeed, numerous papershave pointed out the necessity to account for variations in den-sity, ionization parameter, or geometry in order to simultane-ously explain the observed strong emission line ratios of galaxiesincluding di ff erent elements and ionization stages (see Stasi´nska& Leitherer 1996; Stasi´nska et al. 2006). Similarly, detailed stud-ies of nearby galaxies with a large set of observations haveclearly shown the limitations of one-zone models and the needfor more complex and realistic representations (e.g., Stasi´nska& Schaerer 1999; Cormier et al. 2012, 2019; Lebouteiller et al.2017). Here we explore some of these aspects with a focus onLyC-emitting galaxies.This paper is organized as follows. In Sect. 2 we presentour observational samples. Section 3 shows the main opti-cal emission line diagrams used to compare to "classical"single-component photoionization models (both density- andionization-bounded). We then construct two-component pho-toionization models including LyC leaking emission and com-pare them with observations (Sect. 4). Our results are dis-cussed in Sect. 5. Section 6 summarizes and concludes the pa-per. Throughout this paper, we assume a solar chemical com-position following Asplund et al. (2009) where log( Z / Z (cid:12) ) = + log(O / H) − . Article number, page 2 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters [SII] 6716,6731A / H [ O III ] A / H Wang et al. 2019Ridge line - SDSS x , per bin of 0.1SDSS sampleLyC leakers Fig. 1.
O3H β -S2H α diagram of the LyC leakers (open circles) and com-parison sample with metallicities 12 + log(O / H) = . − . iii ] / H β values. Many leakers, especially the strongest ones,show a deficiency in [S ii ] with respect to this latter ridge line.
2. Optical emission line properties of LyC emitters
In the present work we study eleven z ∼ . − . + log(O / H) = . − .
16, with a median oxygenabundance of 12 + log(O / H) of 7.91.These sources represent currently the best sample of localLyC emitters, for which a wide range of observations is avail-able, including LyC, UV, and optical spectroscopy. Furthermore,Gazagnes et al. (2018) and Gazagnes et al. (2020) provided adetailed analysis of their UV absorption lines (including Lymanseries lines and other low-ionization absorption lines), whichprovides important geometrical constraints on the ISM of thesegalaxies. Here we analyze the optical emission line properties ofthese LyC emitters.Before discussing the observed emission line properties ofthe LyC emitters, we briefly summarize the observational dataused, including the main galaxies of interest (11 LyC emitters)and an appropriate comparison sample.
The emission line measurements from the SDSS spectra of theLyC emitters are reported in the discovery papers (Izotov et al.2016a,b, 2018a,b). We recently obtained new VLT spectra of5 of the 11 leakers (J0925 + + + + [OI] 6300A / [OIII] 5007A110 [ O III ] A / [ O II ] , A Ridge line - SDSS x , per bin of 0.1SDSS sampleLyC leakers Fig. 2.
O32-O13 diagram.
We use the improved emission line measurements from thesespectra, and in particular also the [S iii ] λ For comparison we use a sample of compact star-forming galax-ies from the Sloan Digital Sky Survey (SDSS) Data Release 14(Abolfathi et al. 2018) analyzed in earlier publications (Gusevaet al. 2019), and for which the [O iii ] λ σ , allowing direct abundance determi-nations using the T e method. In practice, from the list of 5607galaxies with direct metallicity measurements, we retain 4571galaxies in the metallicity range 12 + log(O / H) = . − .
2, forcomparison with the LyC emitters discussed here.
The emission line properties of the LyC emitters and the com-parison sample are shown in Figs. 1, 2, and 4. By selection, allthe galaxies studied here, including the LyC emitters, have emis-sion lines ratios compatible with those of star-forming galaxies,that is, they fall in the range of stellar photoionization sourcesin the classical diagrams by Veilleux & Osterbrock (1987) usedfor spectral classification (including the [N ii ] BPT diagram in-troduced by Baldwin et al. 1981). Also, by selection, the sourcesstudied here have a relatively low metallicity, which implies afairly high excitation, reflected for example by a high [O iii ] / H β ratio. However, at a finer level, the LyC emitters show some in-teresting distinctive properties, which we now discuss. Article number, page 3 of 16 & A proofs: manuscript no. 38634corr
Table 1.
Coe ffi cients of the ridge lines in classical diagrams (Figs. 1, 4and 8), determined for the entire galaxy sample and described by Eq. 1.Here y = [O iii ] λ / H β . x [N ii ] / H α [S ii ] / H α [O i ] / H αα -2.40 -2.00 -3.27 α α -5.74e2 -3.98e2 -4.33e2 α α -9.77e4 -6.40e3 -6.68e3 α α -1.94e4 -1.23e4 -1.28e4 α α -2.79e3 -1.74e3 -1.83e3 Table 2.
Coe ffi cients of the ridge lines in oxygen and sulfur diagrams(Figs. 2,5 and 9), described by Eq. 1. Here y = [O iii ] / [O ii ]. x [O i ] / [O iii ] [S ii ] / [S iii ] α -1.67 4.13e-1 α -1.07 -5.15e-1 α α α -9.67e-1 1.30 α -5.33 -9.20e1 α α α -1.84 3.31e1 ii ] deficiency in LyC emitters In Fig. 1 we show a classical diagnostic (O3H β -S2H α ) involv-ing the [S ii ] / H α ratio, which was used by Wang et al. (2019) toselect LyC-emitter candidates among more metal-rich and there-fore lower excitation (lower [O iii ] / H β ) sources. Clearly several,if not a majority, of the eleven confirmed leakers show fairlyweak [S ii ] λλ ii ] / H α ratios close to orbelow the typical value for galaxies from our comparison sampleat the same [O iii ] λ / H β ratio. To quantify this comparisonwe derive an SDSS “ridge line” which marks the median loca-tion of normal star-forming galaxies in the same range of O / H(metallicity) as our LyC emitters. This was obtained by binningthe data in logarithmic bins of 0.1 dex along the [O iii ] / H β axiswith each bin having at least ten galaxies. Each bin was thenfitted by a Gaussian distribution of a center x and a standarddeviation σ . The ridge line defined by the x was then fitted withan eight-order polynomial,log x = (cid:88) n = α n (log y ) n , (1)where y = [O iii ] λ / H β , and the coe ffi cients α n are reportedin Table 1. Table 2 presents the coe ffi cients of the ridge linesderived in other diagnostic plots x vs y = [O iii ] / [O ii ].For subsequent discussions, we also calculate ∆ [S ii ], the lat-eral shift from the ridge line, ∆ [SII] = log ([SII] / H α ) − f (cid:0) log ([OIII] / H β ) (cid:1) , (2)to quantify deviations of individual galaxies from the averageobserved line ratio of the entire comparison sample. With thisdefinition ∆ [S ii ] < ii ] / H α with respectto the average at the same [O iii ] λ / H β ratio. We note thatthe sulfur abundance ratio (S / O) is constant and shows a fairlysmall scatter over a wide range of metallicities, even beyond the [ O I ] J1333+6246J1248+4259 J1011+1947
SDSS sampleLyC leakers f e s c Fig. 3. ∆ [O i ] vs. ∆ [S ii ] derived from classical diagrams O3H β -O1H α and O3H β -S2H α . The three strongest leakers (J1256 + + + f esc ) are reprentedby orange to brown circles in the lower part of the plot. The other galax-ies discussed in this paper are labeled. The LyC emitters show a cleardisplacement in this plane compared to the comparison sample. range of O / H considered here (see, e.g., Guseva et al. 2020).Therefore, the mean trends described by this ridge line shouldbe quite independent of the exact S / O abundance.As shown in Fig. 1, five of our LyC emitters are within 1 σ ofthe ridge line, whereas the rest are systematically o ff set to lower[S ii ] / H α ratios, that is, they show an [S ii ] deficit, an empiri-cal finding already pointed out by Wang et al. (2019). Severalreasons could be invoked to explain such a behavior: S / H abun-dance variations; a decrease of [S ii ] due to the loss of the outer,neutral regions related to the escape of LyC radiation; a lowercontribution of di ff use ionized gas (DIG), which is in generalknown to “boost” [S ii ] emission (cf. Sanders et al. 2020); orpossibly other explanations. The [S ii ] deficit and its link withthe observed LyC escape fraction are discussed below, providinga more quantitative description of the relationship (Sect. 4.4). Another peculiarity of the LyC emitters is found in the diagramcomparing the three ionization stages of oxygen traced by theoptical line ratios [O iii ] / [O ii ] and [O i ] / [O iii ], as shown in Fig.2 (denoted the O32-O13 diagram hereafter). Indeed, naively onewould expect to find a deficit of [O i ] emission in LyC emit-ters as they should be highly ionized and possibly lacking neu-tral regions where [O i ] is emitted. This is also shown by sim-ple density-bounded photoionization models (see Figs. 5 and 6,and Stasi´nska 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 sur-prisingly high, shifted to even higher than average values whencompared to the ridge line defined by the comparison samplein Fig. 2. This empirical finding clearly shows the presence ofneutral oxygen in the LyC leakers, which already indicates that Article number, page 4 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters the ISM in these galaxies can probably not be described by sim-ple density-bounded models, and whose presence and intensityneeds to be explained quantitatively by models.
In Fig. 3, we show the displacement relative to the ridge linesin O3H β -O1H α and O3H β -S2H α diagrams (see columns 2 and3 in Table 1). In this plane, we see a clear shift from the SDSScomparison sample of presumably nonleaking galaxies. Interest-ingly, the three most leaking galaxies (the three reddest circles, f esc > ∼ ∆ [S ii ] and ∆ [O i ] ( < -0.2 dex). On the other hand, most ofthe small leakers in our sample (yellow circles with f esc < ∼ ii ] (-0.2 < ∆ [S ii ] <
0) and[O i ] excesses ( ∆ [O i ] > ∼ ii ] ( ∆ [S ii ] < ∼ -0.3) haveintermediate (J1011 + f esc = + + + ii ] line but the largest [O i ] excess in the sample. A moredetailed analysis of these shifts in [O i ] and [S ii ], including acomparison to models, is given in Sect. 4.4. In the [N ii ] BPT diagram (O3H β -N2H α , Figs. 4 and 8) it is pos-sible to note a slight o ff set of most of the leakers with respectto the average of the comparison sample. Likely related to this,Guseva et al. (2020) noted that the N / O abundance of LyC emit-ters is again higher than the average of star-forming galaxies atthe same O / H abundance. The origin of this behavior is currentlyunknown.In another diagram discussed further below, relating[O iii ] / [O ii ] to [S ii ] λλ / [S iii ] λ iii ] λ ii ] λλ / [S iii ] λ i ] / [O iii ] λλ i lines. As shownby Guseva et al. (2020), our Xshooter observations confirm thetrend suggested by this study.To gain further insight into the observed behavior of theemission line ratios we now proceed to comparisons with pho-toionization models, first with simple and then more complexgeometries.
3. Comparison to classical photoionization models
The first natural step in a comparison between the observedemission line properties of LyC emitters is to examine how sim-ple 1D photoionization models allowing for varying degrees ofionizing photon escape behave and compare with the data.
To examine the power of simple 1D photoionization models, weused Cloudy photoionization models assembled in the BOND(Bayesian Oxygen and Nitrogen abundance Determinations)grid of models (Vale Asari et al. 2016), which is available onthe 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 ni-trogen abundances in giant H ii regions. One of the specificitiesof the BOND code is that it is independent of any assumptionon relations between O / H and N / O. Therefore, the grid spans awide range in O / H, N / O, and ionization parameter U, and coversdi ff erent starburst ages and nebular geometries (spherical bub-ble or thin shell). The spectral energy distribution (SED) of theionizing radiation is obtained from the population synthesis codePopStar (Mollá et al. 2009) for a Chabrier (Chabrier 2003) stellarinitial mass function (IMF) between 0.15 and 100 M (cid:12) . Varyingstarburst ages describe variations in the ionizing radiation fieldhardness, which arise due to the evolution of the stellar popula-tions ionizing the H ii regions. A subset of the BOND model gridhas been used by Stasi´nska et al. (2015), who also studied thee ff ect of LyC leakage. In the current study, we use models fromthe 3MdB_17 database (ref = ’BOND_2’), which is an updatedversion of 3MdB using Cloudy v17.02 (Ferland et al. 2017) andcovering a somewhat refined grid.For our study, we used a sub-grid of BOND models, with pa-rameters appropriate for the LyC-emitting galaxies and the com-parison sample. Specifically, the metallicity range is limited to12 + log(O / H) = U varies between -4 and -1 (with a step of 0.25), the age between1 and 3 Myr (in steps of 1 Myr). The models are spherical filledbubbles with a fixed density of 100 cm − , and abundance ratiosof N / O = − and S / O = − . . The choice of a spherical filledbubble geometry influences the ionization structure of the mod-eled H ii regions: based on Stasi´nska et al. (2015), we expect thinshell models to produce fainter emission from low-ionization-potential species (e.g., [O i ], [S ii ]) with respect to higher ion-ization potential species (e.g., [O iii ], [S iii ]). Specifically, thiswould induce a lateral shift to the left-hand side in the O32-O13and O32-S23 diagnostics plots (Fig. 5 and 9), which we examinein Section 3.2.To describe density-bounded models, we used a parameter H β f rac which cuts the nebula at di ff erent radii corresponding toa range of 10–100% of the total H β luminosity. This parame-ter corresponds to the fraction of H β luminosity emitted by theentire modeled galaxy relative to the H β luminosity that wouldbe 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 β f rac into the corre-sponding LyC escape fraction by defining f esc (H β ) = H β f rac which we use as a proxy to estimate how much leakage is al-lowed by a given density-bounded model. We refer to the modelsused here as simple, one-zone models. We now compare the one-zone photoionization models to the ob-servations using the classical diagrams introduced by Veilleux& Osterbrock (1987) for spectral classification, the above-mentioned O32-O13, and other diagrams. The 3MdB website describes the 3MdB and 3MdB_17 databases at https://sites.google.com/site/mexicanmillionmodels/ Since the study of Vale Asari et al. (2016), the grid in H β f rac hasbeen refined and is now available for steps of 10%.Article number, page 5 of 16 & A proofs: manuscript no. 38634corr [NII] 6584A / H [ O III ] A / H log U=-1-1.5 -2 -2.5 -37.8 8.012+log(O/H)=8.2 [SII] 6716,6731A / H -1 -1.5 -2 -2.5 -37.67.8 8.08.2 [OI] 6300A / H -1 -1.5 -2 -2.5 -37.8 8.08.2 f e s c [NII] 6584A / H [ O III ] A / H log U=-1-1.5 -2 -2.5 -3Age=123 Kewley et al. 2001Ridge-line - SDSSLyC leakersSDSS sample [SII] 6716,6731A / H -1 -1.5 -2 -2.5 -3123 [OI] 6300A / H -1 -1.5 -2 -2.5 -3123 f e s c Fig. 4.
Classical diagrams with BOND models and observations. Photoionization models are represented by the colored symbols: the colorbarrepresents the variation of the escape fraction and the size of symbol is proportional to log( U ). The shape either represents the variation ofmetallicity (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 theSDSS galaxy sample in the restricted metallicity range, 12 + log(O / H) = . − .
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).
Figure 4 shows the LyC emitters, the comparison sample, andthe 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 av-erage SDSS galaxies in the corresponding plots (cf. Sect. 2.2).For illustration, dashed lines indicate the theoretical delimita-tions between star-forming galaxies (below) and AGN (above),according to the models of Kewley et al. (2001).Overall, the ionization-bounded ( f esc =
0) photoionizationmodels cover reasonably well the observed line ratios of the en-tire 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 examinedhere; discussing the origin of these o ff sets is beyond the presentscope. The position of the LyC leakers indicates high ionizationparameters and young stellar ages of the ionizing populations;the latter is in good agreement with the young ages found fromthe observed UV spectral features and from spectral modeling ofthese galaxies (see, e.g., Izotov et al. 2019).The second most important conclusion from Fig. 4 is that thedomain of the emission line ratios observed in the LyC emittersis only covered by ionization-bounded photoionization models,that is, models with f esc = f esc < .
1, which is in contradiction with the observed escape
Article number, page 6 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters fractions reaching up to ∼
70% for some of the leakers. As al-ready shown by Stasi´nska et al. (2015), increasing the escapefraction in these simple, one-zone 1D photoionization modelsstrongly decreases the [S ii ] and [O i ] emission, which resultsfrom trimming the outer parts of the H ii region. This showsthat overall the observed LyC emitters (or at least those with asignificant LyC escape fraction) cannot be reproduced by sim-ple density-bounded models, as these systematically underpre-dict lines of species with low ionization potentials, such as [O i ]and [S ii ]. In Fig. 5 we examine complementary and probably more pow-erful emission line diagnostics to probe the ionization structureand therefore the di ff erent zones of the ionized ISM in the LyCemitters. By focusing on di ff erent ionization stages of the sameelement, we also avoid dependencies on elemental abundanceratios and metallicity, at least to first order. These O32-O13 andO32-S23 diagnostics highlight two main points that motivate theintroduction of more complex models in the Sect. 4:(1) Although the oxygen line ratios (O32-O13) predictedby simple ionization-bounded models reproduce the SDSS se-quence with ionizing spectra of young populations ( ∼ − i ] λ i ] / [O iii ] ratio at a given [O iii ] / [O ii ] ratio, whichis observed for the LyC emitters and other sources (see Fig. 5left). For [S ii ] / [S iii ], shown on the right panel of this figure,the deficiency of the models is even more striking, as they failto reproduce the bulk of the observations. The two discrepan-cies might arise from the fact that the models do not take intoaccount a DIG component, which can significantly contributeto the emission of [O i ] and [S ii ] (Sanders et al. 2017, 2020).Alternatively, these discrepancies could also be attributed to thepresence of shocks, which are not included in our models (see,e.g., Plat et al. 2019) because of problems with atomic data forsulfur (Izotov et al. 2006; Kewley & Dopita 2002; Kewley et al.2019), or others. The sulfur line ratio [S ii ] / [S iii ] is also sen-sitive to more subtle e ff ects including stellar winds, outflows,and photoevaporation of photo-dissociation regions (PDRs) inGalactic dense gaseous environments (e.g., Westmoquette et al.2013; McLeod et al. 2015).(2) The density-bounded models allowing for LyC photonescape do not at all reproduce the observed oxygen and sul-fur line ratios of the leakers, as already noted above. On thecontrary, while leakers populate the upper part of the SDSS se-quence, the density-bounded models produce very low emissionin [O i ] and [S ii ], which drives them in the opposite directionto the bottom-left corner of our diagnostic plots. Figure 6 rep-resents the [O iii ] / [O ii ] and [O iii ] / [O i ] ratios predicted by themodels as a function of their escape fraction. This diagram em-phasizes the discrepancy between [O iii ] / [O i ] and f esc , showingthat [O iii ] / [O i ] is overestimated by several orders of magnitude.This clearly shows that moderate to high escape fractions cannotbe produced by simple density-bounded one-zone models with-out significantly decreasing the weak line emission of [O i ] and[S ii ].In short, although density-bounded models would explainthe Lyman-continuum-escaping photons, they clearly do notmatch the spectral signatures arising from the outer part of theH ii region and observed in LyC emitters. On the other hand,while ionization-bounded models can describe the overall prop-erties of most galaxies (the comparison sample) fairly well— except for sulfur line ratios—, they underpredict the emissionfrom low-ionization species (e.g., the [O i ] line) in the LyC emit-ters, and most importantly, they do not account for the escapeof LyC photons. These failures of simple, one-zone photoion-ization models clearly indicate that more complex models arerequired to describe the observed emission line properties of theLyC emitters. We therefore examined composite photoionizationmodels based on the combination of two distinct zones with dif-ferent properties.
4. Composite "two-zone" photoionization models
As demonstrated above (Sect. 3.2), single-component (one-zone) photoionization models fail to simultaneously reproducethe observed line ratios and escape fractions of LyC leakers. Toanswer our problem (1) from Sect. 3.2.2 (insu ffi cient emissionin low-ionization lines), we introduce a component of low ion-ization (low U ) that would mimic a DIG contribution in star-forming galaxies. The relative contribution of this componentis described by a factor ω that varies from one galaxy to an-other. In the case of LyC emitters, although they are expected tobe relatively poor in DIG (see Oey et al. 2007), this componentcould e ff ectively represent the e ff ect of low-ionization channelscarved into the ISM through which some LyC photons escape,and which may contain gas with a low-volume filling factor, suchas clumps or filaments. As an answer to problem (2) from Sect.3.2.2 (inconsistency between the escape fraction and the pres-ence of low-ionization lines), we include in our two-zone modelthe possibility for LyC photon escape from one region, while thesecond remains ionization-bounded.A two-zone, “picket-fence”-type model, was also suggestedbased on other observations. Gazagnes et al. (2018), Chisholmet al. (2018), and Gazagnes et al. (2020) analyzed the UV ab-sorption lines of hydrogen and low-ionization metal lines of theLyC emitters studied here. Their work clearly shows the pres-ence of absorption lines formed in the neutral gas, with columndensities that are too high for the escape of LyC photons. Fur-thermore, the depths of the saturated absorption lines show thatchannels of low column density must be present in these galax-ies, allowing LyC escape. For our purposes, their model is basi-cally equivalent to our two-zone model. As we show below, suchan idealized two-zone model can also capture and describe theobserved optical emission lines of the LyC emitters.In practice, we have considered two di ff erent scenarios thatare represented in Fig. 7. Each scenario is based on a linear com-bination of two models from the BOND database (see Sect. 3.1),one with a high-ionization parameter weighted by a factor ω , andone 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 thegalaxy which is highly ionized and homogeneous (filling factor (cid:15) ∼ U ( r = R S ) = Q π R S nc , and R S = (cid:32) Q πα ( T e ) n (cid:15) (cid:33) / , (3) Article number, page 7 of 16 & A proofs: manuscript no. 38634corr [OI] 6300A/[OIII] 5007A110 [ O III ] A / [ O II ] , A log U=-1 -1.5 -2 -2.5 -312+log(O/H)=7.88.0 8.2 A /[SIII] 9068 A -1-1.5 -2 -2.5 -37.67.88.08.2 f e s c [OI] 6300A/[OIII] 5007A110 [ O III ] A / [ O II ] , A log U=-1 -1.5 -2 -2.5 -3Age=123 Ridge-line - SDSSSDSS sampleLyC leakers A /[SIII] 9068 A -1-1.5 -2 -2.5 -3123 f e s c Fig. 5.
O32-O13 and O32-S23 diagnostics with BOND models. Photoionization models are represented by the colored symbols: the colorbarrepresents the variation of the escape fraction and the size of symbols is proportional to log( U ). The shapes of model dependencies reflect thevariations either with metallicity or stellar age and are shown and labeled in the top and bottom panels, respectively. where Q is the number of ionizing photons, n the density of hy-drogen, α ( T e ) the recombination coe ffi cient of hydrogen, and (cid:15) the filling-factor for an inhomogeneous medium (0 < ∼ (cid:15) < ∼ U ( r = R S ) = A ( Qn (cid:15) ) / , (4)where A is almost constant ( A ≈ . × − cm s / for T e = K), which means that both low-density and clumpy mediums( (cid:15) (cid:28)
1) will correspond to a model with a low input ioniza-tion parameter. We note that the density-bounded models cre-ated by imposing a cut in radius to the full ionization-boundedmodel do not necessarily have a mean ionization parameter thatmatches their input ionization parameter. The two componentsdo not have to share exactly the same Strömgren radius but onecan ensure that the physical sizes are compatible by adjusting thefilling factor of one component to match the size of the other. Wenow further describe the scenarios in Fig. 7.
Density-bounded channels (DBC) : Ionization-boundedgalaxy with density-bounded channels. The galaxy is ionization-bounded but LyC photons can escape through channels, whichare matter-bounded. Such holes could be carved in the ISM bysupernovae explosions or could result from hydrodynamical fragmentation e ff ects induced by turbulence or stellar feedback(e.g., Kakiichi & Gronke 2019; Hogarth et al. 2020). Weexpect such mechanisms to e ffi ciently clear out channels, eitherproducing low-density channels or almost emptied holes withsmall clumps of matter remaining within. In both cases, wemodel such channels by a density-bounded component withlow input ionization parameter (see Equation 4). Since the LyCemitters studied here are clearly dominated by very young stellarpopulations with ages < ∼ − ff ects to be the maindriver for the creation of these channels or holes, although theexplanation may be more complex (cf. Hogarth et al. 2020). Density-bounded galaxy (DBG) : Density-bounded galaxywith ionization-bounded clumps. The core of the galaxy is com-pletely ionized by massive young stars, which produce a largenumber of ionizing photons exceeding what can be absorbedby the available nebular gas. This is represented by a high- U ,density-bounded component, from which photons can directlyescape at the edge. We add a low-ionization component to rep-resent denser, ionization-bounded filaments or small clumps,which can naturally emit the [O i ] and [S ii ] emission lacking Article number, page 8 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters f esc [ O III ] A / [ O II ] , A log U=-2.5-2-1.5-1 12+log(O/H)=7.67.88.08.2 f esc [ O III ] A / [ O I ] A -2.5-2-1.5-1 7.67.88.08.2 LyC leakers f e s c 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 theobserved LyC emitters and the one-zone photoionization models with log( U ) > ∼ − .
5. We note the very high range of the predicted [O iii ] / [O i ]ratio. Fig. 7.
Schematic illustration of our two-zone photoionization modelsfor LyC emitters, combining an ionization region with density-boundedregions chosen to have di ff erent ionization parameters (“high” and“low” U ). In the "density-bounded channels" scenario, ionizing photonsescape from the low- U region, while in the "density-bounded galaxy"scenario, they escape from the high- U region with f esc = f esc between 0% and 80% for the leaking region. in one-zone models (see section 3) and where the observed UVlow-ionization lines are formed. Both components have the samechemical composition and the same ionizing spectrum, emittedby a central source.For a first approach, we combined two arbitrary photoion-ization models with di ff erent ionization parameters: one with ahigh-ionization parameter U h and one with a low-ionization pa-rameter U l . The high- and low-ionization parameters were cho-sen "by hand" so that the nonleaking combined models accu-rately reproduce the upper part of the SDSS sequence on our di-agnostic plots (see Sect. 4.2) and the tight SDSS sequence in theO32-O13 and O32-S23 plots (Fig. 5 and 9). We used U h = − and U l = − . . We briefly comment on this choice and an opti-mizing procedure below (Sect. 5.1).We computed both scenarios for an identical stellar popula-tion of 1 Myr, a metallicity 12 + log(O / H) =
8, and N / O = . 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 . Thetwo components are then combined with relative weights, ω forthe U h component and (1- ω ) for the U l component. Any ratio ofemission lines A and B of the combined model then becomes L A c L B c = ω L A h + (1 − ω ) L A l s ω L B h + (1 − ω ) L B l s , where s = Q h / Q l . (5)For each scenario, we allowed the LyC escape fraction of theleaking component to vary from 0% to 80% in steps of 10%,provided by the 3MdB_17 database. The total f esc of the two-zone model, is computed consistently from the relative weightof zones ω , by defining f cesc = ω f hesc + (1 − ω ) f lesc . (6)As only one component is leaking in each scenario, Eq. 6 simplybecomes f cesc = (1 − ω ) f lesc in the DBC scenario and f cesc = ω f hesc in the DBG scenario. In the absence of leakage ( f esc = i ] and [S ii ] λλ As we can see from Fig. 8, the nonleaking two-zone models fol-low a bended or sometimes straight “mixing line” from very highto moderate excitation (in [O iii ] λ / H β ) and covering a widerange in the [N ii ] / H α , [S ii ] / H α , and [O i ] / H α ratios. Despiteour somewhat arbitrary choice of models describing the individ-ual zones, we already see that observed line ratios of the bulkof the SDSS galaxy sample (which primarily hosts galaxies with Article number, page 9 of 16 & A proofs: manuscript no. 38634corr [NII] 6584A / H [ O III ] A / H [SII] 6716,6731A / H [OI] 6300A / H f e s c [NII] 6584A / H [ O III ] A / H [SII] 6716,6731A / H [OI] 6300A / H f e s c Fig. 8.
Classical diagrams (same as Fig. 4) showing the predictions from the two-zone photoionization models. The upper and lower panelsrepresent DBC and DBG models, respectively. very low or zero f esc ) can be covered with the two-zone mod-els combining two di ff erent ionization parameters. In particular,there seems to be no di ffi culty in explaining the relatively highobserved intensities of [O i ] λ ii ] λλ / [S iii ] λ ffi cient to solve the above-mentioned discrepan-cies of classical models when compared to normal galaxies. In contrast to the one-zone models, the combined two-zone mod-els are able to cover the entire range of emission line ratios of theLyC emitters, as shown by Figs. 8 and 9. Most importantly, the
Article number, page 10 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters [OI] 6300A/[OIII] 5007A110100 [ O III ] A / [ O II ] , A Ridge-line - SDSSSDSS sampleLyC leakers A /[SIII] 9068 A f e s c [OI] 6300A/[OIII] 5007A110100 [ O III ] A / [ O II ] , A Ridge-line - SDSSSDSS sampleLyC leakers A /[SIII] 9068 A f e s c Fig. 9.
O32-O13 and O32-S23 diagrams (same as Fig. 5) showing the predictions from the two-zone photoionization models. The upper and lowerpanels represent DBC and DBG models, respectively. models can—at least qualitatively—explain the observed [O i ]excess and the [S ii ] deficiency (see Fig. 1 and Fig. 2), whichcharacterizes the leakers, and which cannot be understood byone-zone models (see Sects. 2.2 and 3.2). Furthermore, the two-zone models account for the escape of LyC photons, and are thusconsistent with this important observational property. Below wefurther explore the consistency of di ff erent line ratios and the ex-tent to which the models also correctly reproduce the observedLyC escape fractions.Let us now examine how the predictions from the two-zonemodels in classical diagrams vary with LyC escape fraction. Inthe O3H β -N2H α diagram (Fig. 8 left), most of the models aresuperposed quite independently of their f esc value. This is due tothe fact that [N ii ] originates from the highly excited region inthe inner part of the H ii region, which remains una ff ected whentrimming the outskirts of the galaxy. Only the most density-bounded models from the DBG scenario, where the highly ion-ized component dominates almost all the galaxy ( ω = ff ect the [N ii ] / H α ratio when f esc exceeds 10%. Incontrast, the [S ii ] and [O i ] line ratios (shown in the middle and right panels) are very sensitive to LyC leakage. We note thatthese emissions primarily arise from the low-ionization compo-nent, that is, the one leaking in the DBC scenario (upper panel).This explains why these line ratios are barely a ff ected in theDBG scenario (lower panel), even for high escape fractions. Onthe contrary, leakage in the low-ionization component (DBC sce-nario) rapidly decreases [S ii ] and [O i ]: only low escape frac-tions ( f esc < ∼ i ] emis-sion of the leakers in this scenario.The same trend is visible on the O32-O13 plot (Fig. 9). TheSDSS sequence seems bounded by f esc ∼ −
10% for the DBCscenario (first row). Such values are su ffi cient to explain the ob-served escape fraction of some of the leakers, represented byyellow circles, but it fails to reproduce higher leakage of ioniz-ing photons, represented by the orange to brown circles. Mean-while, the DBG scenario can produce escape fractions up to 80%without strongly decreasing the emission of [O i ] and [S ii ]. Onthe other hand, the O32-S23 plan seems to be better covered bymodels from the DBC scenario, while the DBG scenario tends tooverestimate the [S ii ] / [S iii ] ratio. For the five leaking galaxies Article number, page 11 of 16 & A proofs: manuscript no. 38634corr with [S iii ] measurements, the observed escape fractions seemconsistent with predictions from the DBC scenario, even for themost leaky galaxy J1154 + f esc = ω < ∼ . ω of the U high region below a certain threshold producestwo-zone models that are dominated by the low ionization ( U low )component, and hence are unable to reproduce the observed lineratios of the leakers. f esc In Fig. 10 we now revisit the predicted O32 and O31 line ra-tios in models with LyC escape presented earlier for one-zonemodels (see Fig. 6). As already mentioned and clearly shownin this figure, the two-zone models solve the [O i ] problem ofsimple models, and they cover the entire range of observed lineratios of the LyC emitters. A comparison of the two panels alsoshows that, at least approximately and without adjustment of thecomponents of the two-zone models, the models are capable ofconsistently reproducing the two oxygen line ratios and the ob-served escape fraction. Figure 10 also shows that LyC emitterswith very high f esc > ∼
20% can only be reproduced by our DBGscenario (where the high U component is leaking), whereas bothscenarios are possible for the galaxies with lower escape frac-tions.Figure 11 (left) illustrates the evolution of f esc as a functionof ∆ [S ii ]. This plot suggests that a su ffi ciently strong [S ii ] de-ficiency (e.g., ∆ [S ii ] < ∼ − .
2) e ff ectively picks out LyC leak-ers, although not all of them. Indeed, several LyC emitters arealso located within the normal range of the comparison sample(presumably dominated by nonleakers), as also pointed out byWang et al. (2019). However, using ∆ [S ii ] to determine the es-cape fraction is not obvious. Compared to the dispersion of ourcomparison sample represented by the gray-shaded areas, 6 outof 11 leakers show a deficit in [S ii ] larger than 1 σ (5 out of 11with deficit larger than 2 σ ). If we include the 2 leaking galaxiesfrom Wang et al. (2019), this amounts to 8 out of 13 galaxieswith a clear [S ii ] deficit at 1 σ (6 out of 13 at 2 σ ). The threestrongest leakers (J1256 + + + f esc > ∼
30% all exhibit a large deficiency in [S ii ] ofat least 0.2 dex from the ridge line ( > σ ). J1011 + ii ]deficit. However, the most [S ii ]-deficient galaxy in our sample(J1248 + + ii ],have escape fractions of only a few percent. This could mean thatthe observed f esc for these galaxies is underestimated, for exam-ple due to the fact that LyC radiation escapes only along other,unobserved lines of sight. Interestingly, in the analysis of the UVabsorption lines of more than 30 low- z galaxies from Gazagneset al. (2020), J1248 + α emission (high equivalent width) despiteshowing a large covering fraction in the H i absorption lines.Demonstrating whether or not any of these sources indeed emitmore ionizing photons in other directions seems di ffi cult. In anycase, our analysis confirms that the method proposed by Wanget al. (2019) could be a promising way to identify new leaking galaxies, although the determination of f esc from the [S ii ] deficitremains challenging.From Fig. 11 (left) we notice once more that only the DBGscenario is able to reproduce the more extreme leakage above20%. The DBC scenario is able to cover the space where most ofthe galaxies with small f esc reside, including the [S ii ]-deficientgalaxies selected by Wang et al. (2019) (green color along thecontours of circles). However, we note that the two galaxies withthe strongest [S ii ] emission ( ∆ [S ii ] close to zero) stand outsideof the contours of the DBC scenario. This indicates that leakagepowered exclusively by channels of low ionization is unable tosimultaneously produce strong [S ii ] emission and escaping LyCphotons. This comes from the fact that in this scenario leakageis allowed exclusively by trimming the low ionization compo-nent, which always results in decreasing [S ii ] emission. Thiscould either suggest that density-bounded LyC could still be thepreferred mechanism even for rather small ( < ∼ 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 trendof f esc relative to ∆ [O i ] in the observations: while the most [O i ]-deficient galaxies and the two undetected galaxies ( | ∆ [O i ] |> σ )exhibit the largest escape fractions, the observed f esc decreasesfor increasing ∆ [O i ]. This could indicate that [O i ] emission,which traces the photo-dissociated interface between ionized andatomic gas or cold dense clumps in the H ii region might belinked to the disappearance of the remaining Lyman continuum,leading to lower escape fractions. However, we notice that twogalaxies (J1248 + + i ]-deficientexhibit very small escape fractions ( f esc ∼ ff ects are likely to cause the drop of Lyman contin-uum even in [O i ]- deficient environment. However, the overalltrend shown by the observations in Fig. 11 suggests that selecting[O i ]-deficient galaxies could also be an e ff ective way to targetgalaxies with high leakage and might correlate better with f esc than [S ii ] deficiency. On the other hand, the [S ii ] lines are gen-erally stronger than [O i ], making the former easier to measure.In Fig. 11 (right) the contours from the DBC scenario in-clude only the two galaxies with the lowest escape fractions.For the same reasons as for [S ii ] emission, leakage throughlow-ionization channels tends to be associated with low [O i ]emission, that is, ∆ [O i ] < ∼
0. The galaxies that are not [O i ] de-ficient are not well reproduced by the DBC scenario, especiallythe galaxy standing on the extreme right-hand side whose [O i ]emission exceeds that of all the models considered here. Thismight suggest that the physics of [O i ] emission is not well takeninto account in our two-component models, and might require aboost from an additional component (see discussion in Section.5). However, we note that apart from the systematic o ff set towardlower [O i ] emission, the DBC scenario qualitatively reproducesthe observed trend of increasing f esc with decreasing ∆ [O i ]. TheDBG scenario accurately reproduces all the ∆ [O i ] and f esc of theleakers (except for the strongest [O i ] emitter) but ∆ [O i ] is inde-pendent from the escape fraction in this scenario; the observedvariations in ∆ [O i ] only come from the decrease of [O iii ] whichvertically moves the models closer to the ridge line but withouta ff ecting [O i ] emission.We conclude that ∆ [S ii ] and ∆ [O i ] could be interestingproxies to estimate the escape fraction but only in the caseof ionization-bounded leakage (DBC scenario). In the case of Article number, page 12 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters f esc [ O III ] A / [ O II ] , A DBC scenarioDBG scenarioLyC leakers f esc [ O III ] A / [ O I ] A f e s c Fig. 10.
Same as Fig. 6 showing the predictions for the two-zone photoionization models. f e s c J1333J1248 J0901J1011
DBC scenarioDBG scenarioLyC leakersWang et al. 2019 f e s c J1333J1248 J0901J1011 f e s c 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 dispersionof the SDSS comparison sample at 1, 2, and 3 σ . The three strongest leakers ((J1256 + + + f esc ) are represented by orange to brown circles in the upper part of the plots. The other galaxies that are discussed in Section 4.4 are labeled inblue. density-bounded leakage (DBG scenario), likely producing thehighest escape fractions, [O i ] and [S ii ] seem uncorrelated with f esc , and other diagnostics sensitive to density boundaries are tobe preferred (see, e.g., Fig. 5 and 9). It seems that [O iii ] / [O ii ], ∆ [S ii ], and [O i ] excess are good hints for LyC leakage, but therestill seems to be other independent dimensions, which proba-bly incorporate several physical parameters such as the natureof the ionizing cluster, its evolutionary stage, the gas content ofeach galaxy, its star formation rate, and metallicity, all of whichcould 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 presentwork.
In this study, we used the observed LyC escape fractions and op-tical emission lines to build simple but consistent two-zone pho-toionization models of the LyC emitters, and to discriminate be-tween di ff erent leakage scenarios. An obvious further step wouldbe to use these two-component models to infer the expected es-cape fraction from the emission line observations only. Such in-direct methods are badly needed in order to optimally exploit up-coming JWST observations of rest-frame optical emission linesfor galaxies during the epoch of reionization.A first automated routine was implemented in an e ff ort toextract best-fit values for the parameters of the two-componentmodels. Although very promising, this method requires a robustimplementation and a careful choice of tracers, which we post-pone to future work. We present here the first trends obtainedfrom this study and discuss possible caveats of this method.We ran a χ -minimization routine on a subgrid of our com-bined models, which includes the ones presented in Sect. 4. We Article number, page 13 of 16 & A proofs: manuscript no. 38634corr defined a four-dimensional grid with the following free param-eters: U high and U low (varied from 10 − to 10 − ), the weight ofthe component of high ionization ω , and the escape fraction ofthe leaking component f holesesc (both varied between 0 and 1). Asobservational constraints we used six line ratios: the four ratiosused in the classical diagrams (Fig. 8), and two ratios from theO32-O13 diagrams (Fig. 9). We did not include upper limits.As our prescription for f combinedesc is fully determined by the pa-rameters of each model (see Eq. 6), we could also determinethe corresponding value of the total LyC escape fraction for the χ 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 fractionsdeduced from the Lyman series absorption lines (Gazagnes et al.2018; Chisholm et al. 2018; Schaerer et al. 2018). Hereafter, wesummarize some interesting trends that merit further examina-tion: – The models minimizing χ always favor two distinct valuesfor the ionization parameters U high and U low . Inspection oftheir marginalized probability distribution functions showsthat typically U high ∼ − − − , and U low ∼ − − − .We interpret this as confirmation of our conclusion in Sect.4 that two-component models do a better job at reproducingthe observed line ratios than single-component models. – Regardless of the considered scenario, the weight of thehigh-ionization component is always much higher than thatof the low-ionization one ( ω ∼ . ff erent cov-ering fractions. We note that for galaxies with a low escapefraction ( f esc < ∼ ω obtained from the best-fittingmodel has a value comparable to the covering fractions, C f ,obtained from analysis of absorption lines (Gazagnes et al.2018; Chisholm et al. 2018; Gazagnes et al. 2020). Thisis consistent with what we expect from the DBC scenario,where the covering fraction should be close to ω but cannotbe reconciled with the DBG scenario, where C f ∼ (1 − ω ) , which is a better fit for some of these galaxies. On the otherhand, for strong leakers ( f esc > ∼ C f < ω . Thisstrengthens the conclusion in Sect. 4.4 that the DBC scenariois unable to reproduce strong leakage, as it predicts too high acovering fraction, and therefore overly low escape fractions.However, we note that the DBC scenario sometimes betterreproduces the six line ratios we considered, even for strongleakers. – Finally, the escape fraction obtained for the best-fit modelsallows us to distinguish the four galaxies with the strongestleakage in our sample. The best-fit models predict f esc ≈ f esc ∼ . − + f esc is a factor of approx-imately two below the observed value.Briefly, this exercise supports our conclusion in Sect. 4.4 thatthe geometry of weak leakers resembles that of the DBC sce-nario, and strong leakers more the DBG scenario. However, we find that robustly distinguishing between the two scenarios basedonly on emission lines is not easy. Furthermore, accurate predic-tions of f esc from the optical emission line ratios do not appearto be straightforward. This will probably require a careful choiceof the relevant line fluxes or ratios to be used in such an analy-sis, a proper treatment of upper limits and uncertainties, a robuststatistical method to overcome the caveats of a simple χ min-imization, the exploration of SEDs allowing also for di ff erentages, and other improvements, which are beyond the scope ofthe present publication. In Sect. 3, we discussed the inadequacy of simple models tosimultaneously reproduce the high and low ionization lines, asrevealed for example in the O32-O13 and O32-S23 diagnostics(Fig. 5). Specifically, we emphasized the o ff set of BOND modelsin O32-S23, which systematically underestimate the [S ii ] / [S iii ]ratio. This o ff set of standard photoionization models has beenshown in previous work (e.g., Pérez-Montero & Vílchez 2009;Kehrig et al. 2006), in which such models fail to explain thehardening in ionizing radiation shown by most low-mass, low-metallicity H ii galaxies in the O32-S23 plane. Problems with thesulfur 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 andMappings models to reproduce the [S ii ] / [S iii ] ratio have for ex-ample recently been shown by Mingozzi et al. (2020). Consid-ering that this only seems to a ff ect sulfur lines, these latter au-thors suggest that this is due to limitations in stellar atmospheremodeling and / or the atomic data for sulfur. Based on the obser-vation that the discrepancy holds for single H ii regions (fromKreckel et al. 2019), Mingozzi et al. (2020) exclude the role ofdi ff use ionized gas as a cause of this issue. However, in our study,we see in Fig. 9 that adding a low-ionization component, whichprovides additional emission of [S ii ] without producing [S iii ],helps to solve this discrepancy. We conclude that for both ourreference sample of star-forming galaxies and the LyC-emittinggalaxies, the addition of such a low-ionization component is suf-ficient to reproduce the observed [S ii ] / [S iii ] ratios. Furthermore,this component also contributes to low ionization emission ofoxygen, simultaneously reproducing the observed line strengthof [O i ] λ v ], [Fe v ], [O iv ], He ii , andothers) and also of low-ionization lines (e.g., Schmitt 1998;Stasi´nska et al. 2015; Plat et al. 2019). The e ff ect of a contribu-tion from shock models with di ff erent magnetic field strengths(tuned to reproduce a typical He ii / H β ratio) on the O32-O13diagram is shown in Stasi´nska et al. (2015). Clearly, somemodel combinations of shock with both density- and ionization-bounded H ii regions can be found that reproduce the three ion-ization stages of oxygen. How these models would a ff ect the sul-fur lines, and whether or not they could also be consistently re-produced remains to be seen. Given di ff erent degeneracies anduncertainties it seems currently di ffi cult to exclude a contribu-tion from shocks (see Stasi´nska et al. 2015). Similarly, a weak Article number, page 14 of 16. Ramambason et al.: Consistent emission line modeling of Lyman continuum emitters
AGN component (contribution ∼
3% of the hydrogen ionizingphotons) cannot be excluded based on the observed emissionline diagrams for the galaxies of our reference sample (Stasi´nskaet al. 2015). Again, the sulfur line ratios used here may providean additional handle on these questions.Finally, the low ionization lines are also a ff ected by di ff useionized gas (DIG) in galaxies, as recently discussed by Mingozziet al. (2020), Shapley et al. (2019), Sanders et al. (2020), andothers. In the LyC emitters studied here, the contribution of theDIG would presumably be low, if we adopt the observed scalingof the DIG contribution with surface density of star-formation(e.g., Oey et al. 2007; Shapley et al. 2019), since the LyC emit-ters are compact by selection, and have a very high SFR surfacedensity (Izotov et al. 2018b).Regardless of the exact origin of emission of the low-ionization lines (such as [O i ] λ ii ] λλ i -column-density( N HI ) lines of sight towards the observer allowing the escape ofLyC radiation and lines of sight with high N HI , where other neu-tral and low-ionization metal absorption lines are also formed(Gazagnes et al. 2018; Chisholm et al. 2018; Gazagnes et al.2020). Therefore, their ISM contains at least two “zones” or re-gions with di ff erent conditions, even if we do not know the exactspatial scales on which this happens. Furthermore, it is very plau-sible that properties that influence the ionization parameter (e.g.,density, clumpiness, filling factor) di ff er between these regions.This justifies our two-zone approach using two di ff erent ioniza-tion parameters, U high and U low , for the LyC leakers. Whether ornot and to what extent such an approach is also motivated for“normal” galaxies (most of which are presumably not LyC emit-ters) remains to be seen. In any case, given the ubiquitous com-plexity of the ISM in resolved star-forming galaxies it is proba-bly not surprising that simple 1D photoionization models cannotreproduce all the observables. Some of the methodology used and some of the questions ad-dressed here certainly call for other applications and for fu-ture improvements. Indeed, although photoionization models arefrequently used to predict and analyze emission lines in star-forming galaxies, they are generally 1D spherical or plane paral-lel models computed with Cloudy or similar codes, which thenallow for exploration in a wide range of parameter space in termsof ionization parameters, radiation field, abundances and so on.Such models, which typically assume a single source of ionizingradiation and a very simple geometry, in contrast to observationsof integrated galaxy spectra, are obvious oversimplifications ofthe true complexity and diversity of the radiation field, ISM, andgeometry.Relatively few more highly complex models have so far beenused for comparisons with observations. Stasi´nska et al. (2015)for example explored the combination of stellar photoionizationand other emission sources (AGN, shocks, and others) with vari-able but arbitrary contributions, and confronted them with a largesample of star-forming galaxies from the SDSS. Cormier et al.(2012, 2019) used one- and two-component models to consis-tently analyze emission lines originating in the ionized and neu-tral ISM, as well as in adjacent PDRs. These have been com-pared to detailed multi-wavelength observations of 39 nearbydwarf galaxies. A very complex tailored model including upto 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 othermodels can and need to be envisaged, although su ffi cient inde-pendent observational constraints are necessary to avoid degen-eracies in a multi-dimensional parameter space.Another approach is to use high-resolution hydrodynamicsimulations of galaxies, ideally including radiation transfer anda proper treatment of the ISM, to predict the resulting emis-sion lines and compare them to observations. Such models, withfairly di ff erent degrees of sophistication have for example beenpresented 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 captureseveral observational trends and allow useful predictions, in par-ticular for high-redshift galaxies. However, direct and detailedcomparisons to observations are not trivial to undertake or to in-terpret. For example, Katz et al. (2020) predicts optical and IRemission lines of LyC leakers and nonleakers, which show goodagreement for the [O iii ] λλ ii ] λλ β emission lines when compared to observations. In contrast,their predictions for [S ii ] λλ / H α and [O iii ] / H β arequite di ff erent from those observed in low- z LyC emitters and in z ∼ − ff erentand complementary modeling approaches should help to developmore robust diagnostics and should therefore provide better in-sight 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 HSTLow- z LyC survey (PI A. Jaskot), will provide for the first time astatistical sample of more than 60 galaxies at z ∼ . − . ff erentapproaches. The first comparison with simple two-zone modelsand a small sample presented here is very encouraging, but in-complete. Nevertheless, there is hope that optical emission linescan to some extent be used to identify LyC emitters, and willideally allow determination (even approximate) of the escapefraction of ionizing photons from these galaxies. Clearly, futureimprovements on these lines will be very welcome and of greatinterest, especially in view of the upcoming launch of the JWST,which will routinely observe rest-frame optical emission linesout 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 classicalexcitation diagrams involving the [O iii ] / H β , [N ii ] / H α , [S ii ] / H α ,and [O i ] / H α line ratios, as well as other diagnostic diagrams in-volving 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 galaxieswith similar metallicities (12 + log(O / H) ∼ . − . i ] / [O iii ] ratio for a fraction of the leakers (as alreadypointed out by Plat et al. 2019), and an [O i ] deficiency for thosewith the highest LyC escape fractions ( f esc > ∼ ii ] found recently by Wang et al. (2019), andsuggest that combined measurements of the [O i ] and [S ii ] linescould provide a good indicator for strong LyC emission. Article number, page 15 of 16 & A proofs: manuscript no. 38634corr
We then compared the observations to simple one-zone andtwo-zone photoionization models accounting for the escape ofionizing radiation in order to construct, for the first time, consis-tent models which reconcile their observed emission lines andthe measured LyC escape fractions. The main results from thisanalysis can be summarized as follows: – Classical one-zone photoionization models underpredict thelow-excitation line emission of [O i ] and [S ii ] arising fromthe outskirts of H ii regions; this underprediction is strong fordensity-bounded models with LyC escape, and also exists forionization-bounded models (see Sect. 3), although it is lessstrong. – Two-zone models, combining regions with a high- and low-ionization parameter ( U high , U low ), where one of which isdensity-bounded, allowed us to resolve these discrepanciesand to reproduce the observed line ratios of the LyC emitters.Furthermore, the models show LyC escape fractions compa-rable to those measured directly from the UV LyC observa-tions. – From the two-zone models we also found a possible di-chotomy between galaxies with di ff erent LyC escape frac-tions. Indeed, the models indicate that galaxies with f esc <
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 es-cape. The observed [O i ] emission can originate from thesefilaments / channels, in what we called the DBC scenario(Fig. 7). For much higher f esc , our photoionization modelsshow that LyC leakage must occur from the region with high-ionization parameter (DBG scenario) in order to explain theobserved emission line ratios.Our two-zone models are inspired by and compatiblewith the “picket-fence’-type model, which successfully repro-duces the observed UV absorption lines of hydrogen and low-ionization metal lines of the LyC emitters, as demonstrated bythe studies of Gazagnes et al. (2018, 2020). Our finding of twodi ff erent scenarios, favored by LyC emitters with a low or highescape fraction, is also compatible with the recent analysis ofGazagnes et al. (2020) combining information from the Ly α lineprofiles and absorption lines. Our results are also consistent withthe numerical simulations of Kakiichi & Gronke (2019), whichsuggest that the mechanism powering LyC escape may evolve ina continuous way from the DBC scenario, producing low escapefractions, to the DBG scenario, which appears as a necessaryconfiguration to reach high levels of leakage.Finally, we have discussed the possibility of improved anal-ysis tools to predict escape fractions based on optical emissionlines and using two-zone photoionization models (Sect. 5). Weconclude that robust predictions would require multi-phase andmulti-sector modeling, including some density-bounded compo-nents to successfully capture the complexity of the ISM structurethat leads to photon leakage. We also suggest that low-excitationtracers 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 studythe morphology and dynamical state of leaking and nonleakinggalaxies. Acknowledgements.
We acknowledge the use of photoionization models fromStephane de Barros at an earlier stage of this work. YII and NGG acknowl-edge support from the National Academy of Sciences of Ukraine by its priorityproject No.0120U100935 "Fundamental properties of the matter in the relativis-tic collisions of nuclei and in the early Universe". JMV acknowledges supportfrom the State Agency for Research of the Spanish MCIU through the ”Centerof Excellence Severo Ochoa” award to the IAA (SEV-2017-0709) and project AYA2016-79724-C4-4-P. RA acknowledges support from FONDECYT RegularGrant 1202007.
References