EDP Sciences
Free Access
Issue
A&A
Volume 602, June 2017
Article Number A45
Number of page(s) 38
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201629675
Published online 07 June 2017

© ESO, 2017

1. Introduction

Star formation in primordial (or quasi-primordial) gas is a fundamental process taking place in the first galaxies that are not yet enriched with elements produced by stellar nucleosynthesis. Star formation proceeds when a cloud is gravitationally bound, dense, and cold enough to be subject to the Jeans instability (e.g., Krumholz 2012). Thermal pressure is removed by lowering the heating from UV photons through H2 self-shielding or absorption by dust particles (and conversion to infrared radiation; IR). Furthermore, the presence of metals, even in small amounts, significantly cools down the gas through radiative transitions such as [C ii] 157μm, [O i] 63μm, and [Si ii] 34μm in the neutral atomic medium, or CO in the molecular medium. In the diffuse interstellar medium (ISM), metal cooling is expected to become dominant over H2 cooling when the metallicity is 1/10Z (Glover & Clark 2014). The cooling rate from metals in the neutral phase and the abundance of H2 are therefore two critical parameters to understand the prerequisites for star formation in low-metallicity environments. At the same time, it is essential to identify the main heating mechanisms at work, especially in the neutral ISM, in order to establish the relationship between the thermal tracers and the star formation process. Infrared cooling lines are, for instance, widely used tracers to probe star formation at potentially all redshifts (e.g., De Looze et al. 2014), despite the lack of precise knowledge concerning the heating mechanisms.

The class of blue compact dwarf (BCD) galaxies contains some of the most metal-poor star-forming galaxies known. Apart from a subcomponent of SBS 0335-052 and a low star formation rate (SFR) BCD recently discovered through a blind H i survey (AGC 198691; Hirschauer et al. 2016), I Zw 18 is the nearby star-forming galaxy with the lowest metallicity known, i.e., 12 + log (O/H) = 7.22 or 1/30Z1, as measured by optical emission lines in the H ii regions (Searle & Sargent 1972; Skillman & Kennicutt 1993; Kunth et al. 1994; Garnett et al. 1997; Izotov & Thuan 1998). Observations of the neutral atomic medium probed by far-ultraviolet (FUV) absorption lines toward the massive stars suggest that the H i region might be even more metal poor (Kunth et al. 1994; Aloisi et al. 2003; Lecavelier des Etangs et al. 2004; Lebouteiller et al. 2013). According to Lebouteiller et al. (2013), it is possible that as much as 50% of the H0 mass in I Zw 18 is pristine. Blue compact dwarfs and in particular the well-studied galaxy I Zw 18, thus represent important probes of the thermal balance of the ISM in primitive environments.

The main heating mechanism in the ionized gas of H ii regions of star-forming galaxies is photoionization2 of H, He, and sometimes He+. In I Zw 18, Stasinska & Schaerer (1999) proposed that the H ii regions may be heated by other energy sources as well (shocks, conductive heating at the interface of an X-ray plasma), mainly owing to the supposedly too large electron temperature observed (see also Kehrig et al. 2016). However, Péquignot (2008, hereafter P08) later performed a detailed modeling of the I Zw 18-NW region using the code Nebu (Péquignot et al. 2001); these authors concluded that photoionization by hot stars could satisfactorily explain the entire optical line spectrum, provided that the H ii region topology is equivalent to an incomplete radiation-bounded shell embedded in a diffuse low-density matter-bounded medium of filling factor unity. In essence, the lower density of the diffuse ionized gas leads to a smaller fraction of H0, which is a dominant cooling agent in low-metallicity H ii regions, and therefore to a higher electron temperature.

The heating of H i regions is comparatively much less understood. While the main heating mechanism in the neutral ISM of our Galaxy is due to the photoelectric effect on polycyclic aromatic hydrocarbons (PAHs) and dust grains (e.g., Weingartner & Draine 2001b), both the low dust-to-gas mass ratio (D/G) and the low PAH abundance observed in BCDs (e.g., Wu et al. 2006; Rémy-Ruyer et al. 2014) may lead to important differences as compared to more metal-rich objects. P08 introduced in his model of I Zw 18-NW the heating of the H i region by the soft X-ray source observed in this galaxy and was able to account to order of magnitude for the low-ionization fine-structure lines then recently detected by Spitzer (in particular [Si ii] 34.8μm and [Fe ii] 26.0μm). P08 also made tentative predictions for the far-infrared (FIR) lines [C ii] 157μm and [O i] 63μm, which are the most important coolants in the neutral atomic medium. According to the models, these lines are mainly produced in an H i region of moderate ionization and temperature, that is, in an X-ray dominated region (XDR), using the terminology introduced in the framework of the physics of active galactic nuclei (e.g., Tine et al. 1997). The study of P08 implied that an X-ray source could provide an effective heating mechanism in the neutral ISM of low-metallicity BCD galaxies, and therefore a possible alternative to the traditional photoelectric effect heating. Thanks to the Herschel Space Observatory (Pilbratt et al. 2010) and, in particular, the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010), it is now possible to compare observations and models for [C ii] and [O i].

More recently, using Hubble/COS, Lebouteiller et al. (2013) observed I Zw 18-NW in the FUV absorption-lines C iiλ1334.5 and C ii* λ1335.7, arising from the ground level and fine-structure level of C+, respectively, and observed against the FUV continuum provided by the UV-bright stars in the stellar cluster. The authors roughly estimated an electron fraction ne/nH ~ 0.1% and attempted to substantiate the assumption of photoelectric effect on dust and PAHs, but with mitigate success, which may be viewed retroactively as evidence in favor of X-ray heating in the H i region, as proposed by P08.

Here, building on the model of P08, the heating by photoelectric effect, X-rays, and by other processes is examined. Consequences may pertain to other metal-poor star-forming galaxies, as there is growing evidence that ultraluminous X-ray sources (ULXs; LX ≳ 1039 erg s-1) are more numerous and more luminous in low-metallicity galaxies (e.g., Kaaret et al. 2011; Kaaret & Feng 2013; Brorby et al. 2014, 2015; Basu-Zych et al. 2016, and references therein). These ULXs are thought to be associated with high-mass X-ray binaries (HMXBs), involving either a stellar-mass or intermediate-mass black hole or, rather unexpectedly pulsating neutron stars (Bachetti et al. 2014; Fürst et al. 2016; Israel et al. 2016, 2017). The effects of X-rays are numerous, since they can photoevaporate small molecules and PAHs (while heating larger grains), and at the same time penetrate deep inside the H i region where they can ionize atomic and molecular hydrogen. The most obvious hallmark of X-ray photoionization of the ISM by luminous X-ray sources is the presence of highly ionized species such as He iiλ4686 recombination radiation. This effect has first been observed by Pakull & Angebault (1986) for the luminous black hole candidate LMC X-1. In the case of the ULX Holmberg II X-1 the detection of a X-ray ionized nebula has furthermore allowed an independent estimate of the total luminosity of the X-ray source from the He iiλ4686 emission (Pakull & Mirioni 2002). By modeling the observed ionization structure with Cloudy (Ferland et al. 1998) these observations imply a largely isotropic X-ray emission and largely exclude any significant beaming of the ULX into our line of sight.

Understanding the origin of [C ii] or [O i] in metal-poor galaxies and the impact of X-rays is an important challenge not only to constrain the neutral gas heating mechanism but also to evaluate the possible reservoir of molecular gas. The apparent lack of molecular gas in I Zw 18 (Vidal-Madjar et al. 2000; Wu et al. 2006; Leroy et al. 2007) is at variance with the present vigorous starburst episode. While it is possible that the earliest stages of star formation occur in the cold atomic gas, with molecular gas forming only at the onset of the star-forming cloud collapse (e.g., Glover & Clark 2012; Krumholz 2012), a significant reservoir of molecular gas that is not traced by CO may still exist, i.e., the so-called CO-dark gas (e.g., Tielens & Hollenbach 1985; Maloney & Black 1988; van Dishoeck & Black 1988; Grenier 2005; Wolfire et al. 2010). The low dust abundance in metal-poor galaxies results in a smaller photodissociated CO core while H2 is self-shielded, resulting in a CO-free molecular gas layer with abundant C+ and leading to an enhanced [C ii]/CO ratio for the global cloud emission (e.g., Poglitsch et al. 1996; Madden et al. 1997). I Zw 18 provides an opportunity to examine the origin of [C ii] and its hypothetical association with molecular gas.

A summary of relevant properties of I Zw 18 is provided in Sect. 2. Observations are described in Sect. 3. A topologically significant model of the NW region is then obtained using the photoionization and photodissociation code Cloudy (Sects. 4, 5). Various models, which are shown to be relevant to the full observed IR emitting region, are explored in Sect. 6. The presence of molecular gas and physical conditions in the diffuse gas are investigated in Sect. 7. Implications of X-ray heating of the H i gas are examined in Sect. 8. Conclusions are found in Sect. 9. Details about Herschel, Spitzer, and X-ray data treatments are provided in Appendices A, B, and C respectively.

2. Characteristics of I Zw 18

Some of the main characteristics of I Zw 18 are listed in Table 1. The most important properties are described in the following. Gas and dust masses are discussed separately in Sect. 5.

Table 1

Main I Zw 18 properties used in this study.

2.1. Distance

The distance to I Zw 18 has been the subject of much debate. Early determinations fell in the range 10−13 Mpc (Östlin 2000; Izotov & Thuan 2004). The distance was then revised to 18.2 Mpc when using the red giant branch tip (Aloisi et al. 2006) and 19.0 Mpc using Cepheids (Fiorentino et al. 2010; Marconi et al. 2010). The distance used by P08 for modeling the NW region was 13 Mpc although he briefly considered one model at 18.3 Mpc. A distance of 18.2 Mpc is adopted in our models. An update to the P08 model using this distance is presented in Sect. 4.2 for consistency and for comparison with the present study.

2.2. Constituents

Although I Zw 18 has been often considered to be a young galaxy, possibly showing its first episode of star formation, several studies have identified an old (>1 Gyr) stellar population (e.g., Aloisi et al. 2006; Annibali et al. 2013). The current onset of star formation could be due to the merging of dwarfs or sub-damped Lyman α systems, as suggested by the somewhat disrupted H i morphology observed, for instance, by Lelli et al. (2012).

I Zw 18 contains a main body and a secondary body. The main body contains two massive stellar clusters, NW and SE (Fig. 1), associated with giant H ii regions and surrounded by an irregular and filamentary halo of diffuse ionized gas (e.g., Izotov et al. 2001). Although the secondary body is gravitationally bound to the galaxy (Petrosian et al. 1997; van Zee et al. 1998), it is disconnected from the main body and contains stars that are older on average (Contreras Ramos et al. 2011). Both NW and SE contain a young stellar population but NW has been the more active recently (Contreras Ramos et al. 2011).

For comparison, the diameter of I Zw 18, as observed at 21 cm (e.g., Lelli et al. 2012), is close to that of the Large Magellanic Cloud, i.e., 7 kpc (e.g., Kim et al. 1998; Staveley-Smith et al. 2003), and the diameter of the NW H ii region, (400 pc), is about twice that of the LMC-30 Dor region.

thumbnail Fig. 1

H i column density contours from Lelli et al. (2012), with 2′′ resolution (beam size in the bottom left). Contours are drawn for 3 (dashed), 6, 9, 12, and 15 × 1021 cm-2. The largest red circle shows the PACS beam at 157μm ([C ii]) and the smallest red circle shows the beam at 63μm ([O i]). Both beams are centered at the emission centroid derived by the PACS optimal extraction method (Sect. 3.1.1). The green cross shows the location of the X-ray point source (Sect. 3.3). The background image is HST/ACS F555W. The NW region coincides with an H i hole. The H i column density peak lies between NW and SE.

Open with DEXTER

2.3. Star formation rate

Legrand et al. (2000) found that a low SFR (~10-4M yr-1) over the Hubble time could explain the metal enrichment of I Zw 18. The instantaneous SFR derived from Hα is 0.1−0.2M yr-1 (e.g., Dufour & Hester 1990; Petrosian et al. 1997; Cannon et al. 2002, 2005). De Looze et al. (2014) calculated a similar SFR using the combination of FUV and 24μm, providing 0.06M yr-1. De Looze et al. (2014) also investigated the applicability of several FIR lines for tracing SFR and obtained 0.02M yr-1 using preliminary measurements of [C ii], [O i], and [O iii] with Herschel/PACS. The radio continuum emission, which consists of the combination of thermal free-free emission and synchrotron radiation, provides yet another independent SFR estimate. In I Zw 18, the 1.4 GHz emission is dominated by synchrotron radiation (Cannon et al. 2005). Using L1.4 GHz = 8 × 1019 W Hz-1 (Hunt et al. 2005) and the SFR calibration from Bell (2003), we obtain 0.13M yr-1.

Annibali et al. (2013) examined color-magnitude diagrams (CMD) from deep Hubble/ACS images and found a larger value of 1 M yr-1 over the last 10 Myr in the most crowded regions (including NW). Part of the discrepancy could be explained by somewhat different timescales probed by each tracer (~100 Myr for the FUV vs. ~10 Myr for Hα or CMD). Some SFR determinations are sensitive to the escape of ionizing photons from the galaxy, but, considering the large amount of surrounding neutral gas, it is unlikely that this could explain the scatter of the different determinations. Various SFR values are considered to scale the cosmic ray (CR) ionization rate in Sect. 6.4.

2.4. Chemical abundances

The H ii region abundances are adopted (Table 1). A discontinuity seems to exist between these values measured from optical emission lines in the H ii regions and those determined from FUV absorption lines in the diffuse H i region (Aloisi et al. 2003; Lecavelier des Etangs et al. 2004; Lebouteiller et al. 2013). According to Lebouteiller et al. (2013), the oxygen abundance may be slightly lower by 0.18 ± 0.16 dex (2σ error bar) in the H i region. A greater discontinuity might exist for C and Si, but absorption-line saturation prohibits a reliable estimate. The origin of this discontinuity, which is in fact much larger in BCDs that are more metal rich than I Zw 18 (see summary in Lebouteiller et al. 2009), is subject to debate. Local self-enrichment of the H ii regions was initially proposed by Kunth & Sargent (1986), but contamination by metal-poor gas along the lines of sight was the favored explanation in Lebouteiller et al. (2013). We use hereafter the H ii region abundances while keeping in mind that all abundances in the H i region may be slightly lower.

The choice of using individual observed elemental abundances (as opposed to the solar abundance pattern scaled to the metallicity of I Zw 18) has some impact on the IR line ratio interpretation. The C/O abundance ratio in I Zw 18 is 2.5 times lower than the solar ratio. In general, C/O tends to decrease with decreasing metallicity in BCDs (e.g., Garnett et al. 1995), which is consistent with enrichment by massive stars at low metallicity. It can also be noted that the Si/O abundance ratio in I Zw 18 is about solar (in both the ionized gas and neutral gas, as discussed in Lebouteiller et al. 2013), indicating that both elements are produced in the same massive stars, that silicon is not significantly depleted on dust grains, and that [Si ii] may therefore be an important gas coolant (see Sect. 5). Iron is not depleted either (P08). Thus, there is no sign of depletion on dust grains in I Zw 18, which is consistent with the low dust-to-metal ratio in this galaxy (Rémy-Ruyer et al. 2015).

3. Observations

We present here observational data that were either not used or unavailable in P08, namely the dust mass and spectral energy distribution (SED), the [C ii] 157μm, [O i] 63μm, and [O iii] 88μm line fluxes from Herschel, the suite of Spitzer lines remeasured, recent X-ray observations, and the H0 mass. A summary of the observational constraints used in the models is provided in Sect. 5.

thumbnail Fig. 2

Herschel/PACS map of [C ii] (top), [O i] (middle), and [O iii] (bottom) emission in I Zw 18. The 25 spaxels of the PACS footprint are overplotted on the F555W HST/ACS image. The contours show the H i column density at 5″ resolution (Lelli et al. 2012). For display purposes, only the fits are shown for each spaxel with detection level >2σ, and the number indicates the detection level in σ. The red circle shows the beam size and the red cross shows the emission centroid as calculated by the optimal extraction (Sect. 3.1.2). Individual spaxel spectra are presented in Appendix A.1.

Open with DEXTER

3.1. Herschel/PACS

3.1.1. Datasets

I Zw 18 was observed by Herschel as part of the Dwarf Galaxy Survey Key Program (DGS; Madden et al. 2013). Observations are described in detail in Cormier et al. (2015), we describe in the following specific information relevant to the observation of I Zw 18. The PACS spectroscopy observations were performed in two steps. The [C ii] 157μm line was observed first in May 2011 (OBSID 1342220973) as part of the SHINING program (PI E. Sturm, KPGT_esturm_1) for 3.7 ks. The [O i] 63μm (OBSID 1342253757) and [O iii] 88μm (OBSID 1342253758) lines were then observed in October 2012 as part of the DGS (PI S. Madden, OT2_smadde01) for 13.8 ks and for 4.3 ks respectively. The input coordinates for the [C ii] observation were slightly different than the [O i] and [O iii] observations.

As explained in Cormier et al. (2015), the projection of the PACS array on the sky is a footprint of 5 × 5 spatial pixels (“spaxels”), corresponding to a 47″ × 47″ field of view. Each spaxel is 9.4″ in size. A single footprint observation was performed since I Zw 18 appears smaller than the footprint size. According to the PACS Observer’s Manual3, the point spread function (PSF) full width at half maximum (FWHM) ranges from 9.5″ (0.8 kpc at the adopted distance of 18.2 Mpc) between 55μm and 100μm to about 14″ at 200μm (1.2 kpc). The spectral resolution is about 90, 125, and 240 km s-1 for [O i], [O iii], and [C ii] respectively.

The data reduction was performed in HIPE 12.0 (Ott et al. 2010) using the default chop/nod pipeline script. The level 1 product (calibrated in flux and in wavelength, with bad pixels masks according to the HIPE reduction criteria) was then exported and processed by our in-house PACSman tool (Lebouteiller et al. 2012) for empirical error estimates and line flux extraction.

Figure 2 shows the footprint and line detections. Each line is well detected (>5σ) in at least one spaxel, and it is always unresolved in velocity. The observed spectra and the line fits are shown in Appendix A.1. Another, independent, observation of the [O iii] line was performed as a small map, but with a lower integration time (see Appendix A.2), so we decided to use only the pointed observation described here.

3.1.2. Line fluxes and spatial distribution

The line profile in each spaxel is adjusted with a Gaussian component and a flat baseline (Appendix A.1). The line width is fixed, constrained by the spectral resolution of the instrument, in particular for the spaxels for which the line is not detected.

Although each line is well detected in at least one spaxel, the low signal-to-noise ratio (S/N) in spaxels corresponding to the wings of the PSF together with the low spatial sampling of the PSF usually prevent an accurate determination of the emission spatial centroid. More specifically, for all lines, most of the spaxels around the brightest spaxel have a detection level 3σ (Fig. 2), so it is difficult to pinpoint the peak position or the source spatial shape with an accuracy smaller than the spaxel size (9.4′′). The optimal extraction algorithm of PACSman was used to obtain a more accurate estimate of the source centroid. The optimal extraction compares the spatial profile of the source with that of the instrument PSF, accounting for the uncertainties in the line flux measurements in all spaxels (see Appendix A.1). We find that the emission is point-like for [C ii], [O i], and [O iii] with an intrinsic extent 6″ (530 pc at the adopted distance of 18.2 Mpc). We also find a remarkable agreement between the centroids, despite the different map position angle and pointing coordinates (Fig. 2). The centroid location and the compact appearance both indicate that the line emission originates within the main body of I Zw 18. The [O i] and [O iii] observations, with a slightly higher spatial resolution than for [C ii], suggest that the centroid is closer to NW than SE (Fig. 2). Overall, the low spatial resolution of PACS observations unfortunately prevents us from disentangling the emission of the NW and SE regions. The fluxes we can derive therefore correspond to the global emission (implications for models are discussed in Sect. 5.6).

Cormier et al. (2015) provide flux determinations for all DGS objects, including I Zw 18, using various methods. We review these methods in the following and examine their applicability to the I Zw 18 observation in detail.

Method F1 scales the flux in the brightest spaxel by applying a point-source correction (from 25% to 67% between 50−220μm). This method is valid for a point-like source exactly centered in a spaxel. Any deviation from this hypothesis results in underestimating the flux determination. Method F1 provides the best S/N since it uses only the brightest spaxel, but an additional systematic uncertainty exists because of the pointing issues and possible deviation from a point source.

The second method uses the footprint subarray of 3 × 3 spaxels centered on the brightest spaxel and adds the line fluxes either from all spaxels (F3 × 3) or only from spaxels with >3σ detections (). A point-source correction factor is also required, although it is much smaller (from 4% to 17% for the range 50−220μm) than for F1. The 3 × 3 methods are more reliable than F1 when the source is not well centered in any spaxel, but it may increase the error bar on the flux determination by including spaxels with a low S/N. The method uses an incomplete sampling of the PSF and therefore results in a lower limit on the flux determination4.

Finally, method Fopt performs an optimal extraction by scaling the normalized instrument point-spread function. This is in principle the best method since it reaches a compromise between S/N and the ability to recover the total flux from a source that is not well centered in a spaxel. The flux calibration remains accurate as long as the source is point-like. Details on optimal extraction are given in Appendix A.1, where it can be seen that the emission in I Zw 18 appears point-like.

Table 2 lists the various flux determinations. We consider F1 to be a lower limit because the emission can never be perfectly centered in any spaxel and we use Fopt for our final fluxes. Cormier et al. (2015) used the F3 × 3 method for I Zw 18 as part of the global and systematic DGS analysis and our revised measurements agree within errors (Table 2). For all lines, we verified that the relatively large error bars of F3 × 3 encompass the Fopt determination. We discuss, in Sect. 5.1, how PACS fluxes are normalized for comparison with the other tracers used in this study.

Table 2

Herschel/PACS line flux determinations.

3.2. Spitzer/IRS

We used archival data from the Infrared Spectrograph (IRS; Houck et al. 2004) on board the Spitzer Space Telescope (Werner et al. 2004) data to measure in a consistent way the suite of lines originating mostly in the ionized gas (e.g., [Ne ii] 12.8μm and [Ne iii] 15.5μm), but also [Si ii] 34.8μm and [Fe ii], which partly originate from the neutral gas. For the low-resolution spectrum (SL and LL modules; R = λ/ Δλ ~ 57−126 over 5−14μm and 14−36.5μm, respectively), the deepest observation available was used (AORkey 16205568). For the high-resolution spectrum (SH and LH; R ~ 600 over 10−19.5μm and 20−36.5μm respectively), AORkey 16205568 was used together with 9008640 and 12622848, which are shallower but less affected by bad pixels in some spectral regions. These observations were all performed in staring mode, in which the source is observed in two nod positions. The PSF FWHM ranges between 2″ at 5μm to 11″ at 38μm.

The investigation of the spatial profiles of AORkey 16205568 in the cross-dispersion direction (Figs. 3 and 4) shows that the source appears somewhat extended in SL and SH (about 6″ FWHM). One can distinguish two components in the SL profile. One component is located at 09h34m02.29s/+55°1427.52, coinciding with NW, with a prominent [S iv] line and relatively shallow continuum, while the other component is located at 09h34m02.31s/55°1422.67, coinciding with SE, with a much weaker [S iv] line and a relatively steeper continuum. The NW component is responsible for 75% of the total [S iv] and 68% of the total H i recombination line Huα 12.37 μm (see fluxes in Appendix B). These values are in good agreement with the Hβ fraction originating from NW, 78% (Skillman & Kennicutt 1993). The NW component appears extended in the SL module while the SE component appears quasi-point-like (intrinsic broadening of 5″ and 1″, respectively). The NW and SE components are not distinguishable in the LL and LH modules because of the relatively lower spatial resolution.

thumbnail Fig. 3

Cross-dispersion profiles along the SL slit of the Spitzer/IRS observation 16205568. The histogram shows the emission as a function of the spatial position along the slit, in pixel units (1 px = 1.8′′). The spatial profile is modeled by two slightly extended sources, one of which corresponds to NW (blue, left) and the other to SE (red, right). The profile is shown for the entire integrated spectral order of module SL1 (~8−14μm; top) and for the [S iv] 10.5μm line only (bottom). Left corresponds to north in Fig. 4.

Open with DEXTER

thumbnail Fig. 4

Spatial positions of the two components seen in the IRS SL slit. The slit is shown with the white rectangle and the circles indicate the locations of the sources shown in Fig. 3. The size of the circles corresponds to the total FWHM of each source (i.e., including both the instrument PSF and the intrinsic broadening). The background image is the B band from the DSS2 survey.

Open with DEXTER

Line measurements are described in Appendix B, where a comparison is performed with Wu et al. (2007, hereafter W07). P08 already compared photoionization models to the W07 fluxes and found an overall good agreement. However, P08 noted that the [S iii] 33.5 μm flux measured by W07 is more than a factor of 2 larger than predictions; based on a similar discrepancy for the [Si ii] 34.8μm line arising in the same IRS module as [S iii] 33.5μm, P08 proposed that the fluxes of these two lines were overestimated. We find new [S iii] 33.5μm and [Si ii] fluxes that are lower by factors of 2.3 and 2.1, respectively, as compared to W07, thereby confirming the hypothesis of P08. Our revised measurement of both [S iii] lines at 18.7μm and 33.5μm confirms that the gas density in the ionized gas is well below 1000 cm-3. Furthermore, the [S iv] flux predicted by P08 is significantly larger than in W07. Our measured value confirms the low value of W07, hinting that the problem may be due to a doubtful S2+ di-electronic recombination coefficient (see P08 for more details).

The upper limit on the PAH emission was calculated using the deep low-resolution observation 16205568. We fitted the latter spectrum with the model of Galliano et al. (2011). The PAH optical properties are from Draine & Li (2007). Rather than fitting individual PAH emission bands, we have adjusted a template with two PAH components (neutral and ionized) with fixed properties. The upper limit on the PAH emission in the range 6−15μm is 1.9 × 10-17 W m-2 with little influence from the neutral/ionized PAH mixture. We review in Sect. 5.1 the extracted line fluxes from Spitzer and Herschel and how they are used as constraints for the models.

3.3. X-ray observations

Following P08, X-rays are to be considered a promising heating mechanism in the H i region of I Zw 18. We describe here the various X-ray measurements, and in particular the most recent observation by XMM-Newton.

I Zw 18 was observed with ROSAT with two instruments, first with PSPC (Position Sensitive Proportional Counters) in 1992 and then with the High Resolution Imager (HRI) in 1997. Fourniol et al. (1996) reported the PSPC detection of an unresolved source and calculated an unabsorbed X-ray luminosity of  erg s-1 between 0.1−2.4 keV. Martin (1996) independently examined the same observation, calculating a lower limit of LX ≳ 1 × 1039 erg s-1 with only Galactic absorption considered. The higher spatial resolution enabled by the subsequent HRI observations led Bomans & Weis (2002) to conclude to the existence of both a point source located in NW and a diffuse component.

I Zw 18 was later observed with Chandra in 2000 and X-ray Multi-Mirror Mission (XMM-Newton) in 2002. The Chandra observation, first reported in Bomans & Weis (2002), was analyzed in detail by Thuan et al. (2004). As noted by Thuan et al. (2004), the X-ray emission from I Zw 18 is dominated by a single point source associated with the NW region5 (see position in Fig. 1). The point-source luminosity is 3 × 1039 erg s-1 in the 0.5−10 keV range (value renormalized to a distance of 18.2 Mpc). Faint diffuse emission was also detected but contributed to at most 4% of the point-source flux.

Kaaret & Feng (2013) analyzed the XMM-Newton observation and found a larger flux, 1.4 × 1040 erg s-1, as compared to the Chandra observation. Based on this increased flux, Kaaret & Feng (2013) suggested that the X-ray point source emission is likely dominated by a single X-ray binary, as already proposed by Thuan et al. (2004). The X-ray source in I Zw 18 is located precisely within the NW cluster, according to the relative position with several quasars in the XMM-Newton images, whose positions are known with high accuracy (Pakull et al., in prep.; Fig. 1).

Diffuse X-ray emission might be present, possibly associated with a supernova (SN) cavity between the NW and SE regions (Thuan et al. 2004). The luminosity of the extended X-ray component measured by Thuan et al. (2004) is about 1038 erg s-1 (value renormalized to a distance of 18.2 Mpc). Although their studies are based on the same Chandra data, Ott et al. (2005) and Kaaret et al. (2011) both argue against the presence of detectable diffuse X-ray emission. In the following, only the X-ray point source is considered. The use of the X-ray observation as a constraint to our models is discussed in Sect. 5.3.

3.4. Photometry

Photometry data is used in our study to examine the predicted SED from the models, in particular in the IR range. We use the photometry data measurements from radio to IR (with 2MASS, Spitzer/IRAC, WISE, Herschel/PACS, and Herschel/SPIRE) compiled in Rémy-Ruyer et al. (2015). Figure 5 shows the Spitzer and Herschel dust emission maps. The measured photometry corresponds to the main body emission, and the NW and SE regions cannot be disentangled. Optical to FUV measurements (within a 20′′ aperture), taken from the NED archival data, are only used for illustrative purposes.

3.5. H  I observations

While P08 explored the H i region heating in I Zw 18, the H0 mass was not explicitly considered. The H0 mass, notably responsible for the [C ii] emission, provides an important constraint to our models (Sect. 5.5). Here the interferometric H i21 cm observations of Lelli et al. (2012) are adopted.

4. Model description

Our objective is to build a photoionization model for the incomplete H ii+H i region shell of I Zw 18-NW that takes into account new observational constraints presented in Sect. 3. In the following we describe the model components (Sect. 4.1) and present the Nebu and Cloudy models (Sect. 4.2).

4.1. Model topology

In the description of I Zw 18-NW by P08, the primary radiation sources are (1) the central young star cluster (responsible for the H ii region shell) and (2) the point-like X-ray source (responsible for a partially ionized warm H i region beyond the ionization front of the H ii region), resulting in a relatively simple geometry with a central UV+X source. The geometry is open, with 65% of the ionizing photons actually escaping the NW region through two different (FUV-)matter-bounded sectors (deprived of ionization fronts), while all photons (except for the hardest ones) are absorbed in the radiation-bounded sector (Fig. 6a). Both the Nebu and Cloudy calculations are performed in spherical symmetry (with no absorption by the far side) and the output parameters are post-processed using the covering factor of each sector. The model chemical composition is given in Table 1.

thumbnail Fig. 5

Photometric maps of I Zw 18 with Spitzer and Herschel (Rémy-Ruyer et al. 2015). The first panel shows the HST/ACS F555W image for reference, with the NW and SE regions circled in white. The dashed black circle indicates the beam size (from 2′′ for IRAC 8μm to 12′′ for PACS 160μm). The 8 μm band is not dominated by PAH emission in I Zw 18 but by stochastically heated small grains and warm grains in thermal equilibrium with the interstellar radiation field.

Open with DEXTER

thumbnail Fig. 6

Description of the sectors used in our models. The first panel shows the topology used in P08 and in most of our models (Table 7). The arrows reaching out of sectors #2 and #3 illustrate the fact that the fraction of escaping photon is significant. Optical depths at 1 Ryd are 1.18 and 0.05 for sectors #2 and #3 respectively. “CF” stands for covering factor. The sectors are drawn according to their respective covering factors but in reality the lines of sight are intermixed.

Open with DEXTER

Based on Hα morphology, and following P08, the inner and outer radii of the H+ shell are set to 1.5″ and 2.5″, respectively (130 pc and 215 pc at 18.2 Mpc). Although the H i 21 cm emission extends over several kiloparsecs from the stellar cluster (e.g., van Zee et al. 1998), the line and continuum IR emission seems to be concentrated within ~6″ (Sects. 3.1.2, 3.2), i.e., 530 pc at 18.2 Mpc. In the models, the H i region extends from 215 to 240 pc. The model diameter thus obtained is on the same order as the observed IR extent. In other words, the P08 model, focused on the NW H ii region, turns out to be suited for the IR emission as well since the latter is intimately linked to the H ii region (see also Sect. 5.6).

As in P08, the radial density profile is obtained from a thermal pressure law as a function of radial optical depth at 1 Ryd. We refer to P08 for details about model convergence. The same density profile is used for all sectors, but the calculation is stopped in each matter-bounded sector at a given optical depth. The density then levels off in the H i zone of the radiation-bounded sector, which ends at a given temperature.

4.2. New Nebu and Cloudy models

The H ii and H i region model selected in P08 (M2X) needs to be updated6 before it can be used to build an equivalent model using Cloudy. Some important input parameters and free parameters are summarized in Table 3. Several parameters are introduced here in order to manage the new heating processes and diagnostics, which are studied by means of Cloudy (Sects. 6, 7). We use Cloudy version c13.037 (Ferland et al. 2013).

The distance is now 18.2 Mpc instead of 13 Mpc with the primary luminosity and initial/final H ii region radii increased accordingly. The helium abundance by number is now He/H = 0.085 (instead of 0.080), as in one trial calculation of P08.

Table 3

Model summary and updates from P08.

Owing to the larger distance, the radial density profile, nH(r), must be updated. After convergence, the new parameters in expression (1) of P08 are as follows: Pin = 2.8, Pout = 24 (in /k/ 105 CGS), and τc = 4.2. The nH(r) obtained from Nebu is numerically introduced step by step in the Cloudy computation.

As in P08, the FUV radiation field is described as the sum of two blackbodies of similar power, with temperatures 4 × 104 K and 8 × 104 K, respectively. Unlike in P08, a realistic stellar continuum is implemented below 1 Ryd. Optical and UV photons contribute to the photoelectric heating on dust grains and to the photoionization of several neutral species (C0, Si0...), thereby influencing chemistry processes. The optical+UV continuum prescription used by Lebouteiller et al. (2013) is adopted. The radiation fields used in P08 and in the present study are compared in Fig. 7.

thumbnail Fig. 7

Blue curves show the input radiation field used in P08, comprised of three blackbodies with temperature 4 × 104, 8 × 104, and 2 × 106 K. The dashed curve is for LX = 4 × 1039 erg s-1 (P08 model M2X) and the solid curve is for LX = 8 × 1039 erg s-1 (P08 model M2X2). The red curves show the radiation field used in this study, with additional UV-optical contribution (from Lebouteiller et al. 2013) and with an improved X-ray spectrum prescription as compared to P08. The dashed curve is for LX = 1.4 × 1040 erg s-1 (luminosity inferred from observations) and the solid curve is for LX = 4 × 1040 erg s-1 (adopted standard). The black diamonds show the unfolded XMM-Newton spectrum (see Sect. 5.3 for more details).

Open with DEXTER

Although the usual signatures of Wolf-Rayet (WR) stars are relatively discrete in I Zw 18-NW (e.g., Brown et al. 2002; Kehrig et al. 2015), P08 considered that the observed nebular He ii recombination emission was essentially due to hot stars. It is now admitted that many low-metallicity massive stars could evolve into either very hot WR stars with weak winds (e.g., Crowther & Hadfield 2006) or else into chemically homogeneous transparent wind ultraviolet intense stars (“TWUIN stars”; Szécsi et al. 2015); these are both almost undetectable in the optical, but emit a plethora of radiation above 4 Ryd. Nonetheless, no star cluster synthetic model could pretend to predict the FUV continuum of I Zw 18-NW around 4 Ryd with any certainty, in particular the amplitude of the discontinuity expected at 4 Ryd. In the P08 and present photoionization models, the low-energy tail of the assumed X-ray source spectrum (Appendix C) contributes somewhat to He ii, and the discontinuity at 4 Ryd is empirically adjusted so that He ii is exactly fitted8. P08 provides an acceptable H ii region model that is useful to our H i region models, but the adopted FUV continuum may not be unique. The important point to emphasize here is that the assumed 1−6 Ryd FUV continuum has strictly no impact on the properties of the present H i region modeling.

Instead of the coarse representation of the intrinsic X-ray emission as a single blackbody at 2 × 106 K (P08), we take advantage of satisfactory fits to the observed XMM-Newton data, obtained by Kaaret & Feng (2013), who assumed either Kerr black hole, cutoff power-law, or diskbb (distribution of blackbodies from an accretion disk) models. Here we adopt the diskbb model spectrum (see Fig. 7 and Sect. 5.3).

In order to check the computations, we have compared the Nebu and Cloudy results in similar conditions, that is only with heating by UV and X-ray photoionization (no dust and no CR). Overall most significant line fluxes agree to better or much better than 15% (Table 4). This agreement between the Cloudy and Nebu model results is impressive, especially as the Cloudy computation was not performed in fully self-consistent conditions. Comparison between photoionization codes in standard conditions most often reveals larger discrepancies (e.g., Péquignot et al. 2001). This success may be partly explained by the fact that, due to the very low metallicity and high ionization in the I Zw 18 NW H ii region shell, the energy balance, dominated by H and He, is much simplified and the physical conditions approach the theoretical limit allowed for H ii regions.

Table 4

Comparison between observed extinction-corrected optical and UV line fluxes and models.

5. Summary of the observational constraints

The main observational constraints considered in this study are the optical and IR emission lines arising in the ionized gas, the IR lines [C ii], [O i], and [Si ii] arising in the H i region, and the H0 mass. The dust mass, dust SED, and X-ray luminosity are also used, but they are as well input parameters of the models and are left some freedom for different reasons explained in this section. The suitability of these constraints to the NW model is also discussed.

5.1. Homogenization of infrared and optical line fluxes

Optical and IR line fluxes are used to constrain the physical conditions of both the H ii and H i region in a consistent manner. The Herschel and Spitzer line fluxes (Sects. 3.1.2, 3.2) first need to be normalized9 to be compared to the optical tracers. All fluxes are scaled to Hβ = 1000.

For the normalization of Herschel/PACS measurements, we calculate the total deredenned Hβ flux in the PACS footprint from the Hα map (de Paz et al. 2003), assuming Hα/Hβ = 2.8 and E(BV) = 0.09 (Schlegel et al. 1998; Schlafly & Finkbeiner 2011). We find F(Hβ) ≈ 11.5 × 10-17 W m-2, implying a scaling factor of 1.15 × 10-19, which is the factor by which the PACS fluxes should be divided by to normalize to Hβ = 1000.

Since we consider the P08 predictions for the H ii region lines as robust, another, independent and informative, estimate of the PACS scaling factor is provided by the ratio between the observed and predicted ionized gas tracer [O iii] 88μm. We find 1.2 × 10-19, i.e., in good agreement with the geometrical factor derived directly from the Hα observation. The same normalization is used in our study for all PACS tracers [O iii], [C ii], and [O i], with a factor 1.2 × 10-19 (Table 5).

The Spitzer/IRS fluxes measured in each module do not strictly correspond to the emission within the corresponding apertures, as a fraction of the source emission is lost outside the aperture owing to the PSF size. For this reason, we cannot simply use the Hα fraction falling inside the IRS apertures. For point sources, the aperture correction is accounted for by the regular IRS flux calibration. In the case of I Zw 18, the emission is somewhat extended, requiring a specific extraction method to recover the total flux (Sect. 3.2). With this method, we expect that our Spitzer/IRS fluxes correspond to the total emission from the main body.

Similar to the scaling for PACS data, we can also test reference ionized gas tracers in the IRS range to estimate the required normalization by comparing the observed flux to the H ii region model predictions from P08. Results are shown in Table 5 where it can be seen that the scaling factors for all modules agree well with each other, implying that the total flux was recovered in the small apertures of the IRS10. The IRS scaling factors also agree well with those derived for PACS, implying that, as expected, the IRS line fluxes are well recovered by our extraction method.

The final normalized fluxes are provided in Table 6. The present study focuses on the H i region and does not aim to improve the H ii region model. Our adopted fluxes for the IR ionized gas tracers are close to those used in P08 and we use the same optical line fluxes as in P08. While our models are adapted to the NW region, we bear in mind that some of the tracers may be contaminated by the SE region (see discussion in Sect. 5.6).

Table 5

Factors to normalize fluxes in W m-2 to Hβ = 1000.

Table 6

Infrared line fluxes scaled to Hβ = 1000.

5.2. H  I region cooling lines

The maximum density in the H i region, the outer radius (or else the minimum temperature) at which the model calculation is stopped, and the X-ray source parameters can be constrained to some extent by the [C ii], [O i], [Si ii], and [Fe ii] fluxes measured in this study (Sect. 6.1). We expect these lines to be optically thin in I Zw 18 based on the low [O i] 145μm/63μm ratio found in other DGS sources (Cormier et al. 2015).

The [C ii] 157μm and [O i] 63μm lines are the dominant coolants in the neutral ISM as long as the metallicity is above some critical value (10-3.5Z; e.g., Santoro & Shull 2006). The [Si ii] 34.8μm line has attracted relatively less attention, partly because its wavelength falls at the edge of the Spitzer/IRS coverage and also because most of the silicon is usually depleted onto dust grains. In low-metallicity environments, however, depletion is weaker and [Si ii] (and to a lesser extent [Fe ii]) is an important coolant. Given the lack of depletion in I Zw 18 (Sect. 2.4), we expect [Si ii] to be an important constraint to the models; [Fe ii] is observationally and theoretically less reliable.

5.3. X-ray intrinsic spectrum

X-rays are a fundamental ingredient in the models, but the X-ray luminosity and spectrum shape as derived from observations may not correspond directly to the radiation absorbed in the H i region. The fact that the Spitzer observations (December 2005) and Herschel observations (May 2011 and October 2012) do not coincide in time with the X-ray observations (Sect. 3.3) is of little consequence since the light travel timescale across the shell is of order 103 yr. More importantly, however, the models indicate that typical cooling timescales in the temperature range 100−300 K are ~104 yr and ~105 yr for densities 500 cm-3 and 100 cm-3, respectively. Then, the X-ray spectrum we are considering in the models is nothing but an average over at least several 104 yr, which may readily differ from the present observations.

Moreover, the observed X-ray spectrum is itself subject to caution. There is observational evidence that the X-ray flux may vary considerably over a few years in I Zw 18 (Sect. 3.3; Kaaret & Feng 2013), which is compatible with state transitions in HMXBs between a low-luminosity hard spectrum state and a higher luminosity state with varying hardness. Such state transitions seem to have been observed in I Zw 18 (Kaaret & Feng 2013) and another BCD (VII Zw 403; Brorby et al. 2015).

The X-ray spectrum shape seen by the gas also bears some uncertainty, in particular for the soft X-rays that are absorbed in the H i region. Because of the degeneracy induced by this absorption, whose value is not accurately known along the X-ray source line of sight, the observed soft X-ray spectrum is poorly constrained (see Appendix C). Nonetheless, high-ionization lines such as [Ne v], with the help of detailed photoionization models, can considerably reduce the uncertainty on the intrinsic soft X-ray flux (Appendix C).

Despite these uncertainties, and ignoring possible strongly anisotropic X-ray emission (Pakull & Mirioni 2002; Kaaret et al. 2004; see however Bachetti 2016; King et al. 2001; Körding et al. 2002), it is relatively safe to assume that, within a factor of a few, the presently observed X-ray luminosity, 1.4 × 1040 erg s-1, should represent the (average) luminosity seen by the H i region in I Zw 18. In order to appraise the generality of the conclusions, combinations of blackbodies that are compatible with the X-ray observations are considered together with the diskbb spectrum (see Appendix C).

5.4. Dust mass

There is evidence that, globally, the D/G in I Zw 18 is much lower than the value assuming a simple scaling with metallicity (Galliano et al. 2008; Herrera-Camus et al. 2012; Rémy-Ruyer et al. 2014; Fisher et al. 2014; Rémy-Ruyer et al. 2015). The dust mass derived by Rémy-Ruyer et al. (2015) for the entire galaxy is robust and used here for reference. The authors derived two values assuming the carbon-rich dust component is described either by graphite or amorphous carbon. To remain consistent with the dust prescription in Cloudy, we use the dust mass derived with the “standard” dust composition of Galliano et al. (2011), where the carbon component is made of graphite, i.e., M. To put this value in perspective, the giant H ii region LMC-N 11, despite a smaller physical size (150−200 pc) as compared to the I Zw 18 NW region, harbors ~3 × 104M of dust (Gordon et al. 2014). The total gas mass in the main body is 108M (Lelli et al. 2012), leading to D/G ≈ 1/1000 D/GMW, where D/GMW is the Milky Way value11.

The global D/G value must be regarded with caution. On the one hand, the dust emission may be associated with only a fraction of the main body gas mass. On the other hand, there is evidence that the extended H i in I Zw 18 may contain as much mass as the main body (Lelli et al. 2012), in which case the global D/G may be driven to even lower values if the extended H i is dust-free. Based on the non-detections of CO (e.g., Leroy et al. 2007), the mass of H2 is ignored, and we keep in mind that D/G could be even lower if a significant fraction of gas exists as a CO-dark molecular gas. In our models we test several values for D/G. Apart from the dust mass, the dust SED shape is a qualitative constraint to our models.

5.5. H0 mass

The H0 mass in the main body of I Zw 18 is 108M, which includes the H i located in NW and in the high H i column density cloud between NW and SE (Lelli et al. 2012). Our observations show that the [C ii] and [O i] emission is compatible with a compact source within the main body (Sect. 3.1.2), so that the corresponding H0 mass should be less than 108M.

We use the 2″ resolution H i map from Lelli et al. (2012) to estimate the H0 mass associated with either the NW region itself (defined as a 430 pc diameter region centered on the NW stellar cluster) or with the slightly different region observed with IR tracers (Sects. 3 and 5.6). For the NW region, the H i column density lies between 3−6 × 1021 cm-2, which translates into a mass 5 × 106M. If we now consider instead the H0 mass associated with the H i column density peak between NW and SE (Lelli et al. 2012), we find a larger mass, 2 × 107M. In the following we consider that H0 masses around 0.5−2 × 107M are acceptable.

5.6. Applicability of the tracers to the NW region

The models are built for the NW H ii region of I Zw 18 and its surroundings. Owing to limited spatial resolution of IR observations, some degree of contamination by SE (gas and dust emission) and more diffuse regions is unavoidable. Nevertheless, our approach can be justified on both observational and theoretical grounds:

  • The emission in all IR gas tracers is compact with a size similar tothe NW region (6″; Sect. 3.1.2). The centroid for[O i] and [O iii] seems to be locatedcloser to NW, possibly coinciding with a dust-rich ionized gasshell located between NW and SE (labeled NW-D3 in Cannonet al. 2002), near theH i peak column density.

  • The convergence to a solution produces a kind of an average model because the (normalized) IR observations are rather global. Nonetheless, since the NW nebular emission is three times stronger than the corresponding SE emission at most wavelengths (e.g., Hβ, Huα, and [S iv]; Sect. 3.2) and given the location of the X-ray source within NW, it can be safely surmised that the model more closely reflects (average) properties of the NW region.

  • The ionization structure and metallicity do not vary much between NW and SE (Legrand et al. 2000; Skillman & Kennicutt 1993; Kehrig et al. 2016). Since our line measurements are normalized to Hβ, this means that line constraints should not be significantly affected by SE.

  • All H i regions will be subject to theX-ray radiation, whether they are associated with the NW or SEgiant H ii regions. Since theionization state is stable in theH i gas, the fine-structure lineemissivity depends essentially on local temperature, which iscontrolled by the local flux of X-rays. Geometrical details havelittle consequence.

  • The H i gas can only exist where the ionizing photons from the star cluster are exhausted, that is beyond an H ii region. Diffuse Hβ emission exists out of the main H ii regions of I Zw 18 and may be associated with diffuse [C ii] and [O i] emission. Trial calculations show that the ratio of strong lines from the H ii region, such as Hβ and [O iii] 5007 Å, to [C ii] and [O i] lines is a slow function of assumed parameters (density and distance to the source) over a large range of conditions. This further suggests that the NW shell model is representative of the more global emission (see also end of Sect. 7.2).

We conclude that, in practice, IR gas tracers can be safely assumed to arise in NW. A possible contamination by regions outside of NW (SE or diffuse) would not significantly alter our results. The H0 mass, on the other hand was calculated specifically for NW, or more precisely for the IR-emitting region that appears to coincide with NW (Sect. 5.5). Finally, the dust mass inferred from observations is the total dust mass of the galaxy (Sect. 5.4), so we keep in mind that the dust mass computed in the NW model can be smaller than the observed value. Several D/G values are explored in the models.

In summary, we consider that the observational constraints presented here can be applied to the NW model. The specific contamination of SE on the dust SED and dust mass is discussed further in Sect. 6.3.2.

Table 7

Cloudy model predictions for NW.

6. H  II + H  I region modeling

6.1. Modeling strategy

Armed with a satisfactory Cloudy model equivalent to the Nebu models (Sect. 4) and with new or updated observational constraints (Sect. 5), our strategy consists in exploring several parameters to evaluate their importance in the H i gas heating. Free parameters for the gas are the maximum density and the minimum temperature assumed in the radiation-bounded sector (Table 3). The other free parameters examined are the X-ray luminosity (Sect. 6.2), the photoelectric effect through D/G (Sect. 6.3), and the CR ionization rate (Sect. 6.4). Other processes such as mechanical heating are not considered.

For each set of free parameters, we monitor the predictions for the main observational constraints, namely the [C ii], [O i], and [Si ii] line fluxes but also the dust mass and H0 mass (Sect. 5.5). In most cases, the minimum temperature in the radiation-bounded sector is constrained by the observed [C ii] flux, while the maximum density is constrained by the [O i]/[C ii] ratio. In short, given the X-ray source, a larger maximum density implies a smaller H0 mass and smaller line fluxes, where [C ii] and [O i] are more sensitive to change than [Si ii] and [Fe ii]. A smaller minimum temperature implies a larger H0 mass and stronger [C ii] and moderately stronger [O i] lines. Increasing the X-ray luminosity and thus the H i zone temperature, all IR line intensities are increased; [Si ii] and [Fe ii] are selectively enhanced, while smaller H0 masses (larger minimum temperature) are required.

Relevant properties of our models are shown in Table 7. Since the combination of maximum density and depth (hence temperature) is constrained by the observed [C ii] and [O i] fluxes, the models are expected, by design, to reproduce both lines. This is not the case, however, for models with no X-rays (ℳ0a and ℳ0b, without and with dust, respectively), which severely underestimate the line fluxes and the H0 mass. The [Si ii] flux is comparatively better reproduced in these models since the line is partly emitted in the H ii region. The reasons why the H i region line fluxes are underestimated are twofold, (1) the UV luminosity provided by the central cluster is not able to heat a large enough H0 mass, even with a large D/G; and (2) we have only considered the OB stellar cluster as the heating source.

All other models include the X-ray source. In order to isolate the effect of each process, we explore models with no dust and no CR (ℳ1), with dust but no CR (ℳ2), and with CR but no dust (ℳ3). The objective is not to find the best combination but to test the impact of each parameter. The full model (ℳ4) is the only one that combines all physical processes. For each model we test different values of D/G and the CR rate, and for each set of values we allow the X-ray luminosity to vary. We then examine the output quantities to compare to observations, namely [C ii], [O i], [Si ii], the dust mass and SED shape, and the H0 mass. We monitor the relative importance of the different heating processes for each satisfactory model.

6.2. Exploratory model 1: X-ray heating

As originally proposed by P08, the heating provided by XR photoionization is a natural explanation to the H i region line emission. The intrinsic X-ray spectrum used in our models is described in Sect. 4.1 (see also Appendix C).

Figure 8 shows the suite of models with only the X-ray heating considered. The predicted [O i]/[C ii] line ratio increases with increasing X-ray luminosity and with increasing temperature. Both effects are due to the fact that [O i] dominates the cooling at the inner edge of the H i region while [C ii] dominates inside. The predicted [Si ii] emission tends to be somewhat weaker than the observed one. Improving [Si ii] using a larger LX is at the expense of [O i], which becomes too large. In addition to observational uncertainties, this moderate discrepancy may suggest that either the adopted Si abundance was underestimated or [O i] is affected by self absorption.

thumbnail Fig. 8

Exploration of models with X-ray heating only (no dust and no CR). The predicted line fluxes (model/observation) are plotted against each model parameter with from top to bottom, the dust mass and H0 mass (log M), minimum temperature (K) and maximum density (log  cm-3) in the radiation-bounded sector, and the X-ray luminosity (log L). Each point represents a model, and the symbol size is inversely proportional to the residuals between model predictions and observations for [C ii], [O i], and [Si ii]. The light gray zone corresponds to the observed values and their associated error bar. The dark gray zone corresponds to the best models.

Open with DEXTER

The model that reproduces the observations (line emission and H0 mass) best is labeled ℳ1. As expected, the [C ii], [O i] 63μm, and [Si ii] transitions are the strongest coolants, followed by [O i] 145μm and [Fe ii], in the H i region. The required X-ray luminosity is 4 × 1040 erg s-1, which is three times larger than the observed luminosity. According to the analysis of Sect. 5.3, such a difference is acceptable.

In summary, model ℳ1 shows that X-ray heating alone is already sufficient to explain the FIR line emission, through a reasonable alteration of the observed X-ray spectrum. Inasmuch as X-ray heating dominates in the H i gas, the H i region cooling lines provide a reliable average of the average soft X-ray luminosity over the last few 104 yr in I Zw 18. We now explore models involving alternative heating mechanisms.

6.3. Exploratory models 2: photoelectric effect heating

In this section we examine the photoelectric effect on dust grains and PAH molecules in addition to the X-ray heating. The objectives are to reproduce qualitatively the observed dust SED shape, obtain a reasonable dust mass, and quantify the photoelectric heating rate. Since the observed dust properties correspond to the entire galaxy (Sect. 5.4), models that overestimate the dust mass or the dust SED are not deemed acceptable.

thumbnail Fig. 9

Models with X-ray and photoelectric effect heating for D/G observed (left) or scaled with metallicity (right). See Fig. 8 for the plot description.

Open with DEXTER

For lack of an accurate description of the dust properties in I Zw 18 (or in any other extremely metal-poor environment), we use the Small Magellanic Cloud opacity files in Cloudy that are based on the size distribution given in Weingartner & Draine (2001a). Our results concerning the gas line emission are barely changed if we use grain properties reflecting the size distribution in the diffuse ISM of the Milky Way. Considering the lack of evidence for elemental depletion on dust grains and the low dust-to-metal ratio in I Zw 18 (Sect. 2.4), the gas-phase abundances are the same for all models. We test several choices for D/G and we calculate a grid with varying X-ray luminosity for each choice.

6.3.1. Uniform D/G (2a, 2b, 2c)

First, a uniform D/G value is assumed across the galaxy. Model ℳ2a uses a global D/G of 1/1000 D/GMW, corresponding to the observed value as constrained by Herschel (Sect. 5.4). Model ℳ2b uses 1/50 D/GMW, corresponding to the metallicity scaling. Figure 9 shows the results for both models.

As shown in Table 7, the solution for model ℳ2a is almost the same as that of model with only X-ray heating (ℳ1), which is due to the fact that the photoelectric heating rate is much smaller than the soft X-ray ionization heating rate (Fig. 10). The photoelectric effect represent only 4% of the total heating in the H i region in model ℳ2a (Table 7), as compared to 58% and 30% for ionization by soft X-rays of He0 and H0, respectively. The dust SED is, however, very underestimated (Fig. 10) and the dust mass somewhat underestimated too. This can be at least partly ascribed to the fact that we model the NW region while the observed dust SED and dust mass correspond to the entire galaxy (Sect. 5.4).

On the other hand, model ℳ2b provides relatively more photoelectric effect heating (57% of the total heating across the H i region; Table 7), which results in a lower required X-ray luminosity as compared to model ℳ2a. Nevertheless, the dust mass and dust SED are now significantly overestimated (Fig. 10) and [Si ii] is not as well reproduced as in model ℳ2a.

The photometry data at 100μm and the upper limit at 160μm provide an upper limit on the uniform D/G. We find that, with the current topology, a radiation-bounded cloud necessarily produces too much dust emission if the D/G is larger than 1/300 D/GMW (model ℳ2c; Table 7). The corresponding dust mass is in relatively good agreement with what observations suggest, but the SED is not well reproduced shortward of 100μm (Fig. 11), which could be due to a missing dust component (not contributing significantly to the dust mass). Model ℳ2c has similar physical conditions in the H i region (density, temperature) as model ℳ2a and the heating fraction contributed for by the photoelectric effect is only 10%.

thumbnail Fig. 10

Model results for ℳ2a (D/G observed, left) and ℳ2b (D/G scaled with metallicity, right). The plots on top show, from top to bottom, the total heating and cooling rate, the contributions to the heating rate, the contributions to the cooling rate, and the cumulative intensities of cooling lines for the radiation-bounded sector as a function of the depth within the cloud (the H ii-H i transition lies at 213 pc). For clarity, the heating and cooling mechanisms in the H ii region are not shown. The predicted dust SED is shown in the bottom panels for both models. The 20 cm measurement is thought to be dominated by synchrotron emission (see text). The vertical lines indicate the range used to calculate the total infrared luminosity (3−1100μm) and the green shaded area shows the modeled emission in this range. The blue shaded area shows the input radiation field. The model ℳ2b shows weak silicate emission bands at 10μm and 20μm that are not seen in the IRS spectrum.

Open with DEXTER

In summary, the H i region cooling lines can be reproduced in models with X-rays together with a wide range of uniform D/G value, but the photoelectric effect is never a significant heating mechanism in the H i region. Heating by soft X-rays dominate and the required X-ray luminosity is only twice as low when D/G increases from 1/1000 D/GMW to 1/50 D/GMW. Still, while a large D/G ratio does increase the photoelectric heating to significant amounts, the dust SED becomes overestimated if the D/G is larger than 1/300 D/GMW. In any case, there is no solution that reproduces well the entire dust SED under a uniform D/G assumption. This leaves us with three possibilities: modify the grain size distribution, modify the topology of the model (e.g., inclusions of dust-rich ionized gas with a relatively large ionization parameter), or keep the topology unchanged but modify the D/G spatial distribution. We investigate in Sect. 6.3.2 only the latter two possibilities because of the lack of constraints on the grain size distribution in such an extreme environment as I Zw 18.

6.3.2. Non-uniform D/G (2d, 2e)

In Sect. 6.3.1 we assumed that the D/G was uniform across the galaxy. The measurement of Rémy-Ruyer et al. (2014) only provides an average D/G. One way to reproduce the observed dust SED could be to keep the topology unchanged whilst modifying the D/G along some lines of sight. For instance, the dust SED assuming a radially decreasing D/G (ℳ2d) agrees better with the observations than models with uniform D/G (Fig. 11). The photoelectric effect in model ℳ2d represents only 12% of the H i heating (Table 7). Even then, model ℳ2d still somewhat overestimates the cold dust emission and therefore the total dust mass.

thumbnail Fig. 11

Attempts to reproduce the dust SED, ℳ2c (D/G maximum assuming standard topology), ℳ2d (radial D/G), and ℳ2e (extra sector with a relatively larger D/G). See Fig. 10 for the plot description.

Open with DEXTER

Another way to reproduce the observed SED is to alter the model topology. While in principle one could simply change the topology assuming a uniform D/G, in practice no satisfactory solution is found. We modify model ℳ2a (1/1000 D/GMW) to include an extra sector with 1/50 D/GMW and with a relatively small covering factor (ℳ2e). The extra sector can represent hypothetical dust-rich inclusions within the H ii region. Cannon et al. (2002) found that several dust-rich regions were located between the SE and NW regions and also close to ionizing sources in the SE region (finding based on variations of the observed Hα/Hβ ratio, even though part of these variations may be ascribed to collisional excitation of H0 in the H ii region; see P08 for more details). A close look at the Spitzer and Herschel photometry maps reveals that dust emission peaks in between the SE and NW regions (Fig. 5). In the nebular and stellar emission maps shown in Cannon et al. (2002), one can notice that compact nebular emission is associated with several stellar sources in SE12 and the corresponding dust emission could peak at relatively short wavelengths because of the compactness of the region. Whether the extra sector may correspond to the SE region or to dust-rich inclusions in the NW region makes little difference in our approach, and the only important constraints are the matter-bounded versus radiation-bounded nature of the new sector, the total Hβ luminosity, and to a lesser extent, the radio free-free emission. If the X-ray source illuminates this new sector (i.e., sector within NW), then it must be matter bounded, otherwise too much cold dust emission would be produced. On the other hand, if the impinging X-ray flux is fainter (i.e., clouds in SE), the extra sector does not necessarily need to be matter bounded.

In model ℳ2e, the new sector has a uniform gas density of 20 cm-3, extends from 50 pc to 100 pc from the stellar cluster, and has a D/G of 1/50 D/GMW and covering factor of 6%. The dust mass in this new sector is small, i.e., only 17 M. The total dust mass is 200M, i.e., somewhat lower than the determination by Rémy-Ruyer et al. (2015). The global D/G (all sectors combined) is 1/1000 D/GMW. The dust-rich sector parameters were tuned to reproduce the observed dust SED and the solution is not unique. Still, while the dust SED is relatively well reproduced (Fig. 11), the presence of this new sector or its origin in NW or SE has little importance for [C ii] and [O i] since these lines mostly come from the radiation-bounded region dominated by X-ray heating.

Overall, considering the lack of observational constraints on the grain size distribution and on the small-scale distribution of D/G in I Zw 18, we cannot prefer any of the models that are able to reproduce the observed dust SED (in particular the mid-IR data points). At any rate, the evidence points toward a non-uniform D/G, with a D/G globally lower than 1/50 D/GMW. The various attempts at reproducing the observed dust SED shape actually have little impact on the predicted [C ii] and [O i] line emission and on the required input parameters (in particular the X-ray luminosity). For all models that satisfactorily reproduce the SED shape and dust mass, the fraction of gas heating provided by the photoelectric effect is always 12%.

6.4. Exploratory models 3: cosmic ray heating

We now consider CR ionization as an additional parameter to the initial model ℳ1. For this test, dust is ignored. While the dependence of the CR heating rate on metallicity is uncertain, the CR ionization rate is expected to be proportional to the SN rate because CR are produced in SN remnants and increase their energy in SN shock fronts. The CR ionization rate is thus expected to scale roughly with the SFR (e.g., Wolfe et al. 2003; Abdo et al. 2010). There is a generally good correlation between the IR luminosity and the radio emission in star-forming galaxies (e.g., Yun et al. 2001), which seems to hold even in metal-poor BCDs (Wu et al. 2008b). Although I Zw 18 is an outlier in the IR-radio relation, possibly hinting at an unusually large radio emission, this is also easily explained by the low D/G and/or by a low fraction of UV photons absorbed by dust. In fact, the SFR derived from 1.4 GHz and from Hα agree within a factor of two (Sect. 2.3) and since the integrated 1.4 GHz emission is dominated by synchrotron radiation (Cannon et al. 2005), this suggests that the CR ionization rate may indeed scale with the SFR in I Zw 18. Nevertheless, the choice of SFR for I Zw 18 is not straightforward. Most SFR tracers (Hα, FUV+24μm, or 1.4 GHz) agree with 0.05−0.1M yr-1, but a CMD analysis indicates a larger value of ~1M yr-1 over the last 10 Myr (Annibali et al. 2013) (Sect. 2.3).

Considering the lack of constraints on the actual CR ionization rate, several values are tested in the following. As for the test of the photoelectric effect (Sect. 6.3), we calculate a grid with varying X-ray luminosity and find the best model that reproduces the observations. The description of the CR heating properties in Cloudy is assumed to hold in I Zw 18. We first consider a Galactic background ionization rate (2 × 10-16 s-1; e.g., Indriolo et al. 2015) and let Cloudy calculate the heating rate self-consistently (model ℳ3a in Table 7). The CR heating becomes significant only at large depths into the cloud. The required X-ray luminosity is marginally lower than the model with no CR, which is due to the additional heating from CR ionization. The heating rate contribution due to CR ionization in the H i region is 13%.

In model ℳ3b we test a CR rate that is five times the Galactic background rate. The heating rate in the H i region due to CR ionization becomes 65% and the required X-ray luminosity decreases to 1.25 × 1040 erg s-1. Although the H0 mass is fairly well reproduced, we notice that as the CR ionization rate increases it becomes more and more difficult to find a satisfactory solution for [C ii], [O i], and [Si ii], with either [O i]/[C ii] or [Si ii]/[C ii] increasing with the CR ionization rate.

In order to isolate and test the heating efficiency of CR ionization alone, we briefly considered a model with no X-ray but with a CR ionization rate scaled to reproduce the observed [C ii] flux. For a Galactic background rate (ℳ3c), the H0 mass required to reproduce [C ii] and the other H i region lines is significantly larger than observations (Table 7). For a rate five times larger (ℳ3d), observations are relatively well reproduced, although the prediction for [Si ii] is somewhat low. This shows that the heating provided by the X-ray source in I Zw 18 is in effect similar to what CR ionization could provide with a flat rate of 10-15 s-1 across the H i region. Since the characteristics of the X-ray source in I Zw 18 are well established, we argue that the large ionization fraction observed in the H i region is the consequence of the X-ray source.

6.5. Full model 4a

The full model ℳ4a combines all heating mechanisms (photoelectric effect, soft X-ray and CR ionization). We tentatively scale down the CR ionization rate by a factor of 10 as compared to the Galactic background (i.e., assuming a scaling with SFR) and keep in mind that the global heating in the H i region due to CR ionization remains under 13% even if a Galactic background rate is chosen (Sect. 6.4). We include an additional dust-rich sector to reproduce the dust SED (Sect. 6.3.2). This model therefore reproduces well the observed [C ii], [O i], and [Si ii] fluxes, the suite of Spitzer and optical H ii region lines, and the H0 mass. The dust mass is somewhat underestimated and the dust SED is qualitatively reproduced.

thumbnail Fig. 12

Radial density and temperature in the radiation-bounded sector of the full model ℳ4a. The plots are shown for the radiation-bounded sector.

Open with DEXTER

thumbnail Fig. 13

Heating and cooling rates for the full model (ℳ4a). The plots are shown for the radiation-bounded sector. See Fig. 10 for the plot description.

Open with DEXTER

The optical and infrared ionized gas line intensities predicted by model ℳ4a are in good agreement with observations (Tables 4, 8) and with Nebu (Sect. 4.2). The presence of dust has little impact on the predicted optical line intensities. The observed He i intensities are affected by stellar line contamination (see P08 for more details). The [S iv] and [Ar iv] lines are also comparatively not well reproduced. The ionization balances of both S and Ar are thought to be affected by uncertain dielectronic recombination rates and we refrain at this stage from elaborating on these discrepancies.

Interestingly, the flux of [Ne v] 14.3μm predicted by model ℳ4a is only a factor of 5 below the current observed upper limit. Since the [Ne v] emission traces the soft X-ray radiation, deeper observations with the James Webb Space Telescope (JWST) may be useful to constrain the X-ray source emission properties in a hardly accessible spectral range (see also Appendix C). The JWST will also be able to test model predictions for other faint neutral gas tracers, such as [Ar ii] and [Fe ii]. Nonetheless present collisional rates with H0 are inaccurate for these lines.

Figure 12 shows the density and temperature profiles for the radiation-bounded sector of model ℳ4a. The electron fraction ranges from ~0.5% down to ~0.05% and the electron temperature remains above 50 K throughout the H i region. Figure 13 shows the cooling and heating contributions. One can see that, as expected, [C ii], [O i], and [Si ii] are the main coolants in the H i region, with [O i] and [Si ii] dominating at the surface of the cloud and [C ii] dominating deeper. The optical depth values of [C ii], [O i], [Si ii], and [Fe ii] 26.0μm are 0.5, 2.3, 0.3, and 0.3, respectively; in a static model, i.e., no internal motions and assuming only thermal broadening for the lines.

The observed low D/G (Sect. 5.4) and low PAH abundance in I Zw 18 (e.g., Wu et al. 2006; Rémy-Ruyer et al. 2014) already suggest that dust grains and PAHs do not contribute significantly to the gas heating. This is confirmed by our model in which the main heating source in the H i region is the photoionization of He0 and H0 by the X-ray source. The photoelectric effect heating increases deeper into the cloud but never contributes to more than ~10% of the local heating rate.

Our model is helpful to understand the origin of the total infrared (TIR) emission. About 85% of the TIR in the model comes from the ionized phase in the H ii region (i.e., the two matter-bounded sectors, the additional dust-rich sector, and the radiation-bounded sector before the ionization front). With regard to the physical processes accounting for TIR, about 60% of the TIR in the model arises from free-free emission in the ionized phase, with the remaining 40% arising from dust (mostly in the ionized phase). The important contribution from the ionized phase in TIR implies that the ratio  = ([C ii]+[O i])/TIR cannot be used to probe the photoelectric effect heating efficiency in the H i region (see also discussion in Sect. 8.2).

Table 8

Comparison between the observed IR line fluxes and models.

7. Model discussion

7.1. Molecular gas content

Our models, describing the H i region physical properties of I Zw 18, provide useful constraints to the molecular gas content. We consider two cases: molecular gas in the uniformly distributed gas of the H i region model (Sect. 7.1.1) and in putative small dense clumps (Sect.7.1.2).

7.1.1. Molecular gas in the absence of dense clumps (variant 4b)

Molecular gas has not been detected yet in I Zw 18, neither from diffuse H2 absorption lines (Vidal-Madjar et al. 2000), H2 near-IR emission lines (Izotov & Thuan 2016), H2 mid-IR emission lines (Wu et al. 2006; this study), nor CO(1-0) (Leroy et al. 2007). Also, PAH emission has not been detected (Wu et al. 2006, 2008a; this study). This apparent lack of molecular gas is remarkable considering the recent formation of massive stars. While low observed CO luminosity leads to the hypotheses that either the molecular gas is more efficiently converted into stars or that CO does not trace the molecular gas, the overall scarcity of molecules could also imply that molecular gas may in fact not be a prerequisite for star formation to proceed (e.g., Glover & Clark 2012). Following Krumholz (2012), star formation could proceed in atomic gas if the temperature is cold enough, with molecular gas forming but only at latest stages, shortly before it is eventually destroyed by the FUV radiation from the newly formed massive stars. Alternatively, molecular gas could have been formed shortly before the star formation episode and then destroyed by UV radiation and shock from the newly born stars.

In model ℳ4b (Table 7), molecular gas is uniformly distributed in the relatively diffuse gas (200 cm-3) of the entire radiation-bounded sector (i.e., using the same sector configuration as model ℳ4a; Fig. 6). The molecular gas that could be hidden in small clumps is discussed separately in Sect. 7.1.2.

The “large” H2 molecule model is used in Cloudy (described in Shaw et al. 2005), with fully consistent formation and destruction rates. Several iterations are performed in order for the line optical depth to be firmly established. The main processes for H2 formation are the H route (H + H → H2 + e) and the dust route (e.g., Le Bourlot et al. 2012). The efficiency of the H route is a strong function of the X-ray luminosity. If the X-ray source were absent, free electrons would be provided mostly by neutral carbon and iron ionization. The H route prevails over the dust route in the following conditions:

  • The electron fraction in theH i region is significant, providingfree electrons for producing H through the reaction H + e → H + γ (e.g., Jenkins& Peimbert 1997).

  • The medium density is 100  cm-3, below which the H photodetachment (H + γ → H + e) starts to dominate (e.g. Kamaya & Hirashita 2001; Chuzhoy et al. 2007). Furthermore, regulation by electron recombination prevents a large molecular gas fraction through the H route unless the density becomes much larger than ~104  cm-3 (for a metallicity adapted to I Zw 18; Omukai et al. 2005).

  • The D/G falls below some critical value (see Glover 2003),

The formation of H2 in dense clouds can therefore be X-ray-induced (e.g., Haiman et al. 1996), with a significant number of free electrons provided by primary and secondary ionization13. The competing photodissociation of H2 by Lyman-Werner photons may, however, still partly control the molecular gas fraction if the X-ray and UV sources are correlated through similar origin and timescale (e.g., Machacek et al. 2003). The positive feedback of X-rays could manifest itself mostly if the UV field is weak but X-rays are still able to ionize interiors of dense clouds. The timescales are critical, and if the X-ray source survives long enough after the UV field weakens, H2 formation in dense clouds may be enhanced, thereby providing an important cooling.

Model ℳ4b predicts an H2 column density of 4 × 1013 cm-2 that is lower than the observed upper limit 1015 cm-2 measured by Vidal-Madjar et al. (2000) with FUV absorption lines. The authors also calculated a theoretical column density of 2 × 1012 cm-2 but underestimated the importance of the H2 formation process via the H route by assuming that free electrons were only provided by ionization of carbon. In model ℳ4b, the H formation route contributes to ~60% of the H2 mass. H2 is mostly formed through the H route at the inner edge of the H i shell, while the dust route eventually dominates deeper in the neutral cloud.

We therefore propose that the diffuse (as opposed to the clumps discussed in Sect. 7.1.2) H2 column density may be close to a few 1013 cm-2 based on our models, bearing in mind that the latter may not account for an accurate description of the photodissociation in a clumpy medium. Even gas with extremely low molecular gas fraction can produce a significant amount of self-shielding with a large enough total column density. Figure 1 of Draine & Bertoldi (1996) shows the self-shielding function without any dust. One can see that the self-shielding factor becomes important (0.5; i.e., half of the flux in the Lyman and Werner bands is absorbed) for N(H2) ≈ 2.5 × 1014 cm-2. In I Zw 18, the column density is thus barely sufficient, if at all, to enable H2 self-shielding unless it is distributed into dense clumps (Sect. 7.1.2).

The H2 mass distributed within the radiation-bounded sector of model ℳ4b is only ~0.1 M (corresponding to a molecular gas mass fraction of ~5 × 10-9). In the sector configuration of ℳ4b, the [C ii] and [O i] emission lines thus trace an almost purely atomic medium.

The CO(1-0) intensity predicted by model ℳ4b is 2 × 10-30 W m-2, to be compared to the upper limit by Leroy et al. (2007) of 8 × 10-22 W m-2 (with a beam of 380 pc). Both the sensitivity and the beam filling factor could explain in part why the current observed upper limit is so large as compared to the model.

7.1.2. Dense molecular clumps (variant 4c)

As mentioned in Sect. 7.1.1, absorption-line studies did not succeed in detecting H2 in the diffuse ISM along the lines of sight toward the massive stars. This is in agreement with model ℳ4b. Nonetheless, H2 may still exist in significant amounts out of the available lines of sight toward the massive stars in NW, which are located within an H i hole (Fig. 1), far from the dust-rich regions (Sect. 6.3.2). For example, H2 may belong to undetected, relatively small, and dense clumps (Vidal-Madjar et al. 2000). Thuan et al. (2005) proposed the existence of such clumps in the BCD SBS 0335052 based on the lack of diffuse H2 together with the detection of warm H2 emission in the near IR. Rubio et al. (2015) observed small CO clumps with ALMA in the low-metallicity galaxy WLM (12 + log (O/H) = 7.8) of 1.5 to 6 pc in size, which begs the question of whether such clumps can exist in I Zw 18.

Similar to Cormier et al. (2012), who modeled photodissociation regions (PDRs) in the dwarf galaxy Haro 11, we consider dense clumps in I Zw 18 with a small covering factor as compared to the H i shell (Fig. 6). In other words, only a small fraction of the radiation-bounded sector reaches the molecular phase. Since there are no constraints on the mere presence of clumps, we bear in mind that the following models are only exploratory. Apart from the upper limits on CO, H2, and PAH emission, another observational constraint is the dust SED and in particular the long wavelength emission corresponding to cold dust.

thumbnail Fig. 14

Spectral energy distribution for model ℳ4c (model with molecular clump). See Fig. 10 for the plot description. The dust-rich sector is ignored for this model.

Open with DEXTER

We assume D/G = 1/1000 D/GMW across the molecular sector. The dust-rich sector (added to reproduce the observed SED shortward of 100μm; Sect. 6.3.2) is ignored for this test since we wish to study the specific contribution of the clump sector. The density profile across the clump sector is calculated self-consistently by constraining a constant pressure (P/k = 8 × 106 CGS units). The density profile is thus similar to the radiation-bounded sector (#1) used in our previous modeling up to the photoionization front, but the density then eventually reaches ~106 cm-3, for which the medium becomes fully molecular (even with no dust grains).

Various combinations of covering factor and stopping column density can lead to significant H2 mass in the clump sector with more mass produced in variants with low covering factor and large stopping column density. As a test, we consider a sector reaching a total column density of 1025 cm-2 (AV ~ 3) with the covering factor scaled to accommodate the observed constraints (the CO upper limit, the observed dust SED, and the corresponding dust mass).

In practice, the CO upper limit is the most stringent observational constraint. It is reached for a clump covering factor of 0.1% (Table 8). For such a covering factor, the dust SED remains compatible with observations (Fig. 14) and the dust mass is also well reproduced (Table 7). The resulting upper limit on the molecular volume is about 143 pc3 (~25 × 25 pc2 surface area and ~5 pc deep). Considering the depth scale and assuming spherical clumps, this should correspond to a couple dozen clumps at most.

thumbnail Fig. 15

Top: abundance fraction of C0, C+, CO, and H2 plotted against the optical depth at 1 Ryd for the clump sector of model ℳ4c. Bottom: cumulative fraction of [C i] 609μm, [C ii], CO(1-0), H22.12μm, and H2 mass normalized to the maximum value.

Open with DEXTER

The main heating mechanisms in the clump sector are the photoionization of H and He at small depths, and the collisions with H2, photoelectric effect, and CR ionization at larger depths. The clump sector is almost fully (98%) molecular in mass. As shown in Fig. 15, most of the [C ii] emission in the clump sector occurs at depths for which the gas is molecular (75% of [C ii] is associated with a molecular gas fraction >50%). However, [C ii] does not trace well the H2mass, since most (80%) of the H2 mass is produced at depths where C is into C0 or CO. Comparatively, [C i], CO(1-0), or H22.12 emissions follow the H2 mass profile more closely, although they emit mostly at depths where the H2 mass is still small.

While [C ii] emission is associated with H2 in the clump sector, the fraction of [C ii] arising from this sector is only about 15% of the total [C ii] in the model (the same fraction is found for [O i]). Most of the emission arises from the relatively more diffuse H i region in sector #1 where [C ii] traces H0.

The H2 mass in model ℳ4c is ~4 × 107M, i.e., on the same order as the H0 mass. From CO observations and using a point-source luminosity, Leroy et al. (2007) found M(H2) <7.4 × 105M14. Assuming the star formation efficiency is the same as in the Milky Way, we can argue that XCO is significantly larger. If XCO is 100 times larger than the Galactic value as proposed by Leroy et al. (2007) similar to the factor extrapolated from the trend by Schruba et al. (2012) we obtain an upper limit on the H2 mass on the same order as our estimate in ℳ4c. A large XCO might be partly related to low cloud extinction due to low metallicity (i.e., enhanced photodissociation) rather than to a direct abundance effect (Glover & Low 2011).

The molecular volume and clump size we estimate are upper limits. First, these values are calculated based on the CO(1-0) observed upper limit. Second, our model does not account for internal motions. The CO emission could be different if we considered turbulence or a velocity gradient, in which case both the molecular volume and clump size could in fact be much smaller.

We conclude that it is possible that significant amounts of molecular gas exist in I Zw 18 if it is distributed within a few small dense clumps of size 5 pc (0.05′′ at 18.2 Mpc). Parsec-size clouds are difficult to observe in absorption because it requires a line of sight intersecting the clump toward a FUV background source. Furthermore, absorption-line measurements of a group of FUV background sources privilege low-extinction lines of sight. While such clumps could be directly detected with high spatial resolution CO observations, it is also possible to observe relatively strong signatures from [C i] and from the warm molecular layer in the mid-IR H2 lines in particular S(0) (Table 8) as long as the molecular clouds remain excited by the X-ray source and by the stellar clusters (hence providing a significant gas temperature). The near-IR H2 line 1-0 S(1) at 2.12μm is a few orders of magnitude below the current observed upper limit (Izotov & Thuan 2016), suggesting that 8 m class telescopes could possibly reach such levels predicted by our models (e.g., Vanzi et al. 2011). Clumps may also exist at significant distances from the clusters, where they can be near isothermal equilibrium with the CMB and CR, in which case only cold molecular gas can be detected. Another possible observational constraint to such dense clumps is provided by the dust emission longward of 500μm (Fig. 14).

The existence of molecular clumps in general in I Zw 18 would shed a new light on the importance of the molecular gas into the baryonic budget of low-metallicity dwarf galaxies; such galaxies have been generally thought to be dark-matter dominated although this is a disputed issue (e.g., Swaters et al. 2011; Ferrero et al. 2012; Kuhlen et al. 2013; Sawala et al. 2015). Such clumps are reminiscent of recurring and self-gravitating cold H2 clumps that may exist in a hierarchical structure of various sizes (Pfenniger & Combes 1994a,b; Combes & Pfenniger 1997; Combes 2005).

7.2. C+ cooling rate in the dense versus diffuse medium

While our models focus on the relatively dense gas associated with the star-forming region of I Zw 18, a significant fraction of the H0 mass of the galaxy is distributed within the large H i envelope (Lelli et al. 2012). In this section we compare the cooling rates observed in both phases.

The C+ cooling rate can be measured either using the [C ii] 157 μm emission line or the C ii* λ1335.7 FUV absorption line doublet (Fig. 16). In the Milky Way, the 2P3/2 fine-structure line population was first observed in absorption toward Galactic stars (Pottasch et al. 1979). The 157 μm emission line was first observed shortly after in NGC 2024 and the Orion nebula (Russell et al. 1980). A comparison between the cooling rate inferred from the absorption and the emission was first performed for our Galaxy by Stacey et al. (1985), indicating a rough agreement. The same comparison in extragalactic targets is challenging as in most sources the material observed in emission cannot always be easily associated with that observed in absorption (e.g., quasar lines of sight). Furthermore, most IR-bright galaxies are difficult to observe in the FUV because of absorption by dust. Inversely, most FUV-bright galaxies are difficult to observe in the IR because of the lack of dust. The advent of Herschel made it possible to detect relatively dust-poor galaxies at FIR wavelengths.

In that vein, the comparison between our results concerning [C ii] and those of Lebouteiller et al. (2013) concerning C ii* is informative and is, to our knowledge, the first comparison made in an extragalactic target. Lebouteiller et al. (2013) examined with the FUV absorption lines C iiλ1334.5 and C ii* λ1335.7, arising from the ground level and fine-structure level of C+, respectively. The observations of C ii*, performed with Hubble/COS, probe lines of sight toward bright young stars withing the 2.5″ aperture. The authors inferred an electron fraction upper limit of 2% in the H i region with a best guess of ~0.1%, which is within the range we have determined with Herschel/PACS and Spitzer/IRS line ratios (Sect. 6.5) and with our models (Sect. 6.5). Furthermore, the cooling rate as determined by the C ii* column density was found to be incompatible with the photoelectric effect on dust, which is again consistent with the results in the present study.

thumbnail Fig. 16

Main levels of the C+ ion. The dot-dashed lines correspond to weak transitions.

Open with DEXTER

We do not expect a priori an agreement between the C+ cooling rate calculated in the present study and the results of Lebouteiller et al. (2013). Firstly, only the most diffuse medium should be seen in absorption because the most dust-rich lines of sight suffer from UV absorption. However, we bear in mind that such a bias is likely negligible for I Zw 18 because of its low dust content. Secondly, and most importantly, the H i region seen in emission corresponds to the entire main body while the H i region seen in absorption by COS probes a region of 2.5″ across centered on NW and coinciding in fact with an H i hole (Fig. 1).

To quantify further the comparison between the cooling rate inferred in emission and absorption, we calculate the total cooling rate in erg s-1. From the observed [C ii] 157 μm line intensity, 10-17 W m-2, we find a total rate of 3.9 × 1031 W. For the measure in absorption, we first use the cooling rate per H atom obtained by Lebouteiller et al. (2013), 8.3 × 10-28 erg s-1 (H-1), which we then multiply by the total mass of H0 in the main body, 108M, assuming that the cooling rate per H atom is uniform across the H i body. We find a rate of 1038 erg s-1 = 1031 W, i.e., a factor of 4 difference with the rate measured in emission. Another way to show this discrepancy is to compare the COS measurement 8.3 × 10-28 erg s-1 (H-1) to the cooling rate predicted by the model (10-24−10-23 erg s-1; Fig. 13). Considering the volume density in the C+-emitting region in the model, the cooling rate becomes 10-26−10-25 erg s-1 (H-1), implying a factor of 10 larger than the COS rate.

thumbnail Fig. 17

Observed COS cooling line rate per H atom is shown for different temperature, density, and various electron fraction (x) values.

Open with DEXTER

The discrepancy we observe cannot be solely attributed to the relatively uncertain H0 mass of the main body used for the calculation and we argue that the cooling rate is genuinely lower in the H i region probed toward the NW cluster. Taking at face value the cooling rates measured in absorption and in emission and the H0 mass of the main body, it seems that most of the C+ cooling in I Zw 18 occurs in the star-forming region with a relatively smaller contribution from the diffuse H i envelope seen in absorption. The lower cooling rate in the diffuse H i region could be at least partly due to the integrated lower density of the ISM probed in absorption as compared to the H i shells probed with Herschel of about 300  cm-3. Lebouteiller et al. (2013) derived an upper limit 30  cm-3 on the density of the medium seen in absorption. We review this measurement by comparing the observed cooling rate to the theoretical value for various temperatures, densities, and electron fractions (Fig. 17). For a wide range of temperatures and electron fractions, we find densities on the order of 15−50 cm-3.

It is thus plausible that the C+ cooling rate measured in emission and absorption do not correspond to the same physical regions. On the one hand, the [C ii] 157 μm emission line probes the neutral shell of the radiation-bounded H ii region with an H i column density of 1022.2 cm-2. On the other hand, the C ii* absorption line appears to probe the more diffuse H i region possibly extending to (or simply located) several kiloparsecs away from the stellar clusters, and providing only little emission in the [C ii] line. The temperature is unfortunately difficult to measure from absorption lines in extragalactic objects because fine-structure absorption lines other than C ii* are usually not detected and because absorption components corresponding to individual clouds cannot be disentangled. The total H i column density seen in absorption is 1021.34 cm-2 (Lebouteiller et al. 2013). The lower chemical abundances measured in absorption by Lebouteiller et al. (2013) in the H i envelope as compared to the H ii region abundances only reinforce the hypothesis of two distinct physical regions while at the same time providing another explanation for the relatively lower cooling rate measured in absorption. Even though the carbon abundance measured with COS is uncertain, Lebouteiller et al. (2013) concluded there is a general metal deficiency in the H i envelope with abundances that are lower by a factor of 2 as compared to the H ii region. Lebouteiller et al. (2013) also hypothesized the presence of pristine gas along the lines of sight that would dilute abundances. If such metal-free gas does exist in the outskirts of the galaxy, it would also explain the measurement of the lower cooling rate per H atom.

The diffuse H i gas could, in principle, be considered an extra sector in our models, but we are limited by the unknown topology of this phase and by the lack of knowledge on the associated heating mechanisms. While clouds lying further out from neutral shells in the H ii region are supposedly shielded from the UV arising from the NW cluster, the clumpiness of the ISM and the scattering complicate this simple picture. Furthermore, the H i structures seen several kiloparsecs away from the stellar clusters may be heated by an old stellar population and by CR, both of which are not well constrained. Other fine-structure absorption lines in the FUV, such as O i* or Si ii* would be required to measure the physical conditions in greater detail and provide useful model constraints.

Nevertheless, we have briefly explored a model in which a low-density (10−50 cm-3) H i gas is located ~1 kpc away from the NW cluster. Such a model is in fact also motivated by Hα emission detected at very large distances, which is likely due to photoionization by radiation escaping from the main star-forming regions (NW+SE). Then H i zones beyond the ionization front of these far away (low-excitation) de facto “H ii regions” should exist, again probably heated by X-ray radiation. Zones deprived of significant H ii gas may as well be heated by X-rays. Interestingly, we find that the [C ii] and [O i] line fluxes relative to Hβ are not strongly affected. Thus, the inference from this tentative example (many other runs are needed to be more precise) is that the [C ii] and [O i] emission arising from the XDR is more or less proportional to Hα at large distances. [C ii] and [O i] may be relatively enhanced due to shielding, which would primarily reduce the relative size of the H+ region.

In summary, the chemical and physical conditions of the medium inferred from absorption lines point at a different origin from the neutral filaments observed in emission around the NW H ii region. Detailed properties of this more diffuse phase would be of interest to investigate, as it must play an important role in the evolution of the galaxy, notably through gas infall, future widespread star-formation, or star-formation quenching from chemical heating.

8. Implications of dominant X-ray heating

In this section we examine the consequences of our main result, which is chiefly that the main heating mechanism in the H i region of I Zw 18 is probably photoionization by soft X-rays.

8.1. Relationship between [C  II] and SFR

The [C ii] luminosity is often used as a SFR tracer (e.g., De Looze et al. 2014, and references therein), which usually has the implicit assumption that the gas cooling down through [C ii] is heated by the UV field arising from young stars through the photoelectric effect. The use of [C ii] as a SFR tracer may become questionable when the heating mechanism of the H i region is different.

The SFR in I Zw 18 calculated from [C ii] is only a few factors below other measurements (Sect. 2.3), while a larger discrepancy may have been expected considering that the heating is not dominated by the photoelectric effect. Nevertheless, we show below that there is a correlation between the X-ray luminosity (responsible for [C ii]) and the SFR for star-forming galaxies.

The X-ray emission in I Zw 18 corresponds to an ultraluminous X-ray source (ULX, defined as LX> 1039 erg s-1). Ultraluminous X-ray sources are not associated with the nuclei of galaxies and are thought to be intermediate-mass black holes or to involve a HMXB system with either a neutron star and/or a stellar mass black hole (see, e.g., Gilfanov et al. 2004a; Berghea et al. 2013; King & Lasota 2016). As they are late stages of dying massive stars, the HMXBs are expected to be related to star formation, with progenitors 8M, and with evolution timescales on the order of 106−7 yr. A correlation between SFR and LX has been shown observationally by Grimm et al. (2003) based on what seems to be a universal HMXB luminosity function in nearby star-forming galaxies (see also Gilfanov et al. 2004b). The SFR-LX relation is linear for high-SFR regime (>4.5M yr-1) and non-linear below. Contaminations include active galactic nuclei and low-mass X-ray binaries (LMXBs), the latter being a more important contaminant for galaxies with low SFR. The non-linear behavior is thought to be a metallicity effect, as proposed by Douna et al. (2015), reflecting the enhancement of the HMXB population at low-metallicity. This enhancement has been expected from theory (e.g., Linden et al. 2010; Fragos et al. 2013; Prestwich et al. 2013; Mapelli et al. 2010; Mapelli 2013) and is supported observationally by the rising evidence of ULXs in nearby metal-poor BCDs (e.g., Kaaret et al. 2011; Kaaret & Feng 2013; Brorby et al. 2014, 2015 and references therein). The measured values of SFR and LX in I Zw 18 are consistent with such a metallicity dependence15. Besides, the dependence with metallicity is inferred indirectly from the deficit of ULX in metal-rich luminous IR galaxies (Luangtip et al. 2014)

In summary, the SFR-LX relationship explains why [C ii] traces relatively well the SFR in I Zw 18. It also explains why the value of the often-used diagnostic ratio ([C ii]+[O i])/TIR is not unusual (Sect. 8.2).

8.2. Relationship between [C  II] and TIR

The ([C ii]+[O i])/TIR (ϵTIR; with TIR calculated between 3−1100μm) ratio is often used to estimate the photoelectric effect heating efficiency16, in which case it may also be used to calculate the physical scale of star formation (Stacey et al. 2010). It is also used to probe the so-called “[C ii] deficit” in ultraluminous IR galaxies (e.g., Díaz-Santos et al. 2013, and references therein). Inversely, metal-poor galaxies show an elevated ϵTIR (Cormier et al. 2015). A larger ϵTIR ratio is traditionally attributed to a larger photoelectric effect heating efficiency, which may be considered a consequence of a lower grain charging due to the relatively lower absorption of UV photons in a dust-poor environment (Israel et al. 1996; Israel & Maloney 2011) or to a smaller dust grain size distribution (Galliano et al. 2005). Our result that the photoelectric effect heating is negligible in I Zw 18 allows us to test diagnostics held by ϵTIR.

For I Zw 18 we find ϵTIR = 0.6%. If we tentatively assume that half of the TIR emission arises in the H ii region (approximate fraction derived from the SED fitting; Galliano et al. 2008), we derive the following values for the PDR/H i region component alone ϵTIR = 1.2%. For comparison, Lebouteiller et al. (2012, 2013) found ϵTIR = 0.4−0.7% in the PDR-dominated regions of the LMC-N11B nebula. In Haro 11, Cormier et al. (2012) found 0.18% for the entire galaxy and 1.1% for the PDR component alone. In Galactic PDRs, Okada et al. (2013) found 0.1−0.9%. We can also compare to other measurements of integrated galaxies (not corrected for the H ii region contribution from TIR). Cormier et al. (2015) find 0.2−1.0% in dwarf galaxies of the DGS while Croxall et al. (2012) find 0.2−0.7% in two super-solar metallicity galaxies.

In summary, ϵTIR in I Zw 18 is similar to the values in other galaxies, although somewhat larger. As shown in Sect. 6.3, this cannot be explained in terms of photoelectric efficiency for this galaxy as photoionization of the H i region by X-rays can solely explain the emission of [C ii] and [O i]. In order to understand why ϵTIR is not extreme in I Zw 18, we must take a closer look at the tracers involved in the quantity ϵTIR.

First, the free-free contribution in the TIR range, TIRff, could become important when D/G decreases to extremely low values. In I Zw 18, using the number of ionizing photons in P08 (scaled at 18.2 Mpc), we find TIRff ≈ 6.9 × 106L, i.e., only about four times lower than the observed TIR luminosity 2.9 × 107L (Rémy-Ruyer et al. 2015). The observed luminosity of [C ii]+[O i] calculated at the same distance is 1.8 × 105L, which results in ϵTIR,ff ~ 2.6%. Hence even if I Zw 18 were virtually dust free, ϵTIR would remain on the order of a few percent, i.e., on the same order as the observed ratio in more metal-rich objects.

Apart from the free-free emission, the observed TIR also includes an important contribution from the dust in the ionized phase. Such dust population does not participate significantly to the heating of the H ii region, and evidently does not contribute to the H i region heating where [C ii] and [O i] originate. In our models, the fraction of the observed TIR that arises by dust in the H i region alone is only 15%. This shows that to calculate the actual photoelectric effect heating efficiency in the H i region, we would need to apply a significant correction to TIR.

The fact that ϵTIR does not reach extreme values in a dust-free environment can be surprising at first, as it suggests a physical connection between cooling lines and TIR, even though the heating mechanisms are different as compared to more metal-rich sources. While TIR scales with the UV field intensity (through dust and free-free emission), [C ii] and [O i] scale with the X-ray emission through photoionization heating. It turns out that both energy sources (UV and X-rays) are correlated (Sect. 8.1), implying that the cooling line emission and TIR scale together on first approximation.

Nevertheless, in a completely metal-free environment the free-free emission would dominate TIR, and [C ii]+[O i] would be null, leading in effect to ϵTIR = 0. Therefore, the fact that ϵTIR does not reach extreme values remains valid as long as the main H i region coolants are metallic species.

8.3. Star-formation quenching from X-ray vs. photoelectric effect heating

Forbes et al. (2016) proposed that the gas depletion timescale in dwarf galaxies may be dominated by the photoelectric effect heating and not by SN feedback. The photoelectric heating could then be an important quenching mechanism for star formation for isolated systems, implying a significant metallicity-dependence for the star formation history. Hu et al. (2017) find conflicting results, with SNe dominating the feedback, but the photoelectric effect still has some impact on the SFR.

The H i region heating in I Zw 18 is in fact dominated by X-ray ionization. Our result, together with the enhanced abundance of HMXBs in low-metallicity galaxies (Sect. 8.1), suggests that star formation may have been mostly suppressed, or at least influenced, by X-ray heating. Our hypothesis stands for a closed-box model of star formation, but interestingly, X-ray and EUV ionization generated by star formation have also been proposed to regulate star formation in galaxies through the reduced cooling rate of the halo gas to be accreted (Cantalupo 2010).

Overall, there seems to exist a compensation between the heating provided by HXMBs and the photoelectric heating from dust, where both heating mechanisms vary with metallicity (or D/G) in opposite ways. Although it is plausible that the photoelectric effect dominates the gas heating in most sources in the DGS sample the enhancement of the HMXB population at low metallicity may play a role in extreme environments such as I Zw 18. The key parameter is not metallicity but, to first order, the competition between D/G and LX (modulated by the topology and the photon absorption). For instance, our models for I Zw 18 predict that the photoelectric effect would be a significant heating mechanism (accounting for ~31% of the total heating in the H i region) if D/G were to be scaled with the metallicity (Sect. 6.3.1).

9. Conclusions

In this study the extremely metal-poor (1/30Z) galaxy I Zw 18 is examined to identify and quantify the main gas heating mechanisms in the H i region.

To this end, a detailed modeling of the H ii+H i region with the photoionization code Cloudy is undertaken, building on a previous model by P08. Our new Herschel/PACS observations ([C ii] 157μm, [O i] 63μm, and [O iii] 88μm) and updated Spitzer/IRS measurements (particularly [Si ii] 34.8μm) are used as constraints. Other new constraints include observations of the dust SED, dust mass, H0 mass, and X-ray spectrum.

The photoelectric effect, CR, and X-ray heating are all considered in different consistent models. The selected models reproduce well the optical to FIR emission lines, dust SED, and other data. The main conclusions are as follows:

  • The photoelectric effect is not a significant heating mechanism inthe H i region, whether we use theD/G constrained by Herschel or a larger D/G scaled withmetallicity.

  • To order of magnitude, the weak photoionization of H0 and He0 by X-rays can provide enough heating to explain the observed low-ionization fine-structure line intensities. The dominant coolants are [O i] and [Si ii] at the inner edge of the H i region and [C ii] deeper inside.

  • More precisely, the H i zone fluxes imply a luminosity LX ≈ 4 × 1040 erg s-1 that is about three times larger than the value inferred from a diskbb model fitting of the XMM-Newton observations. This difference is not critical, given that (1) the source may not emit isotropically; (2) the X-ray source is variable; and (3) the model X-ray flux corresponds to a time average over the H i gas cooling timescale of several 104 yr. This large average X-ray luminosity could imply the existence of a massive accreting black hole in I Zw 18-NW.

  • The lack of constraints in an environment such as I Zw 18 prevents from any definite conclusion concerning CR heating. Nonetheless, CR heating is negligible if the “standard” description in Cloudy is adopted.

  • The models predict only about 0.1M of H2 in the region traced by [C ii] and [O i] if the radiation-bounded sector has uniform properties (no clumps). Under such conditions, [C ii] and [O i] are therefore not (CO-dark) molecular gas tracers.

  • We also explore the possible existence of undetected molecular clumps. These H2 clouds should reach a density of 106 cm-3 and have a small covering factor (0.1%). The main constraints come from the CO emission and the dust SED. This dense sector can shelter 4 × 107M of H2, in which case H2 must belong to a few clouds of a few parsecs in size (0.05″). Most of the [C ii] and [O i] emission in the clump sector is associated with a large molecular gas fraction, but theses lines do not trace well the H2 mass. Even in the presence of clumps, most (85%) of the [C ii] and [O i] emission in the model is associated with atomic gas in the more diffuse sector.

  • The cooling rate derived from C ii* in absorption against UV-bright stars is somewhat less than the rate measured from [C ii] in the H i shell adjoined to the H ii region. In absorption, the lines of sight may intersect a rather diffuse neutral (~10 cm-3) up to several kiloparsecs away from the stellar cluster.

  • The dust SED is not well reproduced with uniform D/G values. A warm dusty component is required to fit the mid-IR part of the SED. Adding this component has no impact on the FIR line emission and on any other result.

  • More than half of the observed TIR flux is free-free emission from the ionized phase. The rest is dust emission (mostly in the ionized phase). Only 15% of the TIR comes from dust in the H i region. Diagnostics involving TIR, such as the photoelectric heating efficiency proxy [C ii]+[O i]/TIR, should be used with caution.

That the H i region heating is dominated by X-ray ionization has important implications for extremely metal-poor galaxies in general. There seems to be a competition between the photoelectric effect (scaling with metallicity) and the X-ray ionization, the latter becoming more and more important in low-metallicity star-forming galaxies through the larger occurrence and luminosity of HMXBs (expected theoretically and checked observationally). Some uncertainty pertains to the difficulty in recovering the X-ray energy distribution actually seen by the gas. We show that some of this uncertainty can be lifted through the use of high-ionization lines such as [Ne v] (Appendix C).

In environments with both an extremely low D/G and bright X-ray binaries, the physical conditions can no longer be easily derived from FIR diagnostics. In I Zw 18, we notice that the [C ii]+[O i]/TIR ratio is surprisingly normal, with a value similar to what is observed in more metal-rich sources. We argue that this ratio is kept stable because of a correlation between X-ray luminosity, which is responsible for the heating of the gas in which [C ii] and [O i] originate, and UV radiation field through the SFR, which is responsible for most of the TIR emission. This implies that the value of [C ii]+[O i]/TIR cannot be interpreted as a photoelectric effect heating efficiency proxy (or as an indirect metallicity tracer) in I Zw 18. For the same reason (i.e., the relationship between SFR and X-ray luminosity in star-forming galaxies), H i region cooling lines, and in particular [C ii], scale with the SFR even though the gas heating mechanism may be different.

Finally, we propose that X-ray heating may be an important quenching mechanism for star formation. This quenching could be efficient at large scales, as shown by the comparison of the C+ cooling rate measured in emission to that measured in absorption. Our models suggest that X-rays can heat the gas at large distances.

The study of I Zw 18 is interesting to put in perspective with the expected occurrence of binary systems made of compact objects in the first galaxies. The latter are star-forming H i rich objects that share similarities with nearby BCDs, and the heating of the H i region may be dominated by the same mechanisms. The number and luminosity of HMXBs is expected to be enhanced in low-metallicity galaxies, and it is therefore likely that such HMXBs contributed significantly to heating the early Universe (e.g., Mirabel et al. 2011). I Zw 18 provides an interesting testbed, but the fraction of escaping X-ray photons needs to be estimated. The ISM is expected to be more porous in dust- and metal-poor environments (e.g., Cormier et al. 2015), which would favor both the escape of FUV and X-ray photons.


1

We use the oxygen abundance obtained by Péquignot (2008) divided by the solar value from Asplund et al. (2009).

2

Outside the H ii regions, in the so-called warm diffuse ionized medium of disk galaxies, some extra heating exists in extremely low-density regions (e.g., Reynolds et al. 1999) and may be due, e.g., to photoelectric effect on dust or dissipation of interstellar turbulence.

4

Alternative methods consist of combining the spaxel spectra (either 3 × 3 or 5 × 5) before performing the line fitting. However, possible variations of the baseline between spaxels and co-addition of spurious features may introduce systematic errors. In any case, the fluxes measured this way, despite large systematic errors, are compatible with our final values.

5

Ott et al. (2005) reported a different position, located at the edge of the Hα nebula, somewhat offset from the NW stellar cluster. While the coordinates of the X-ray source reported in their Table 6 does in fact coincide with the NW cluster, their interpretation seems to rely on incorrect astrometry of HST images.

6

The secondary ionization by X-ray photoelectrons was inadvertently skipped in the P08 computation. Secondary ionizations result in more energy going into ionization, at the expense of the thermal energy available for line excitation in the H i region, so the predicted line fluxes become smaller than reported in P08 for a given primary X-ray luminosity.

7

The default atomic data of Cloudy version c13.03 are used except in that Abrahamsson et al. (2007), also used in P08, is preferred to Launay & Roueff (1977) for the critical collisional excitation O0(3P) + H0(2S). The more recent rate is many times larger than the old rate.

8

Fitting He ii and a number of other nebular lines is a prerequisite to any model of the H ii region (see P08).

9

We do not refer to the normalization of fluxes due to a possible contamination by regions outside NW (Sect. 5.6) but to the normalization between observations with different apertures.

10

The present absolute flux values are globally larger than those in W07, and the scaling factor reduces the fluxes to levels corresponding approximately to the P08 normalization. P08 did not want to consider all of the Hβ emission, but that arising from the main NW shell. However, P08 calibrated the different instruments, using Hα fluxes from comparable regions (not necessarily the region he wanted to model). In the same way, P08 derived 65% of ionizing photons escaping from the NW shell.

11

D/GMW is 1/148, close to the values in Zubko et al. (2004), 1/158, and Jones et al. (2013), 1/156.

12

This is in contrast with the NW region where most of the nebular emission seems to be offset from the stars (except in the NW D1 H ii region), in filaments forming an incomplete shell.

13

Number of secondary ionizations is 25 (/ 1 keV) (Shull & van Steenberg 1985).

14

Value rescaled at 18.2 Mpc. with a Galactic CO/H2 conversion factor.

15

I Zw 18 would fall closer to the linear relationship with SFR of ~1 M yr-1 as suggested by the CMD analysis from Annibali et al. (2013) (Sect. 2.3). Photon escape probability and the use of indirect tracers to infer SFR may therefore play some role and explain in part the different behavior of SFR-LX observed at low-metallicity (as the metallicity decreases, photon escape is presumably enhanced due to a lower D/G and to a more porous ISM).

16

ϵTIR represents the gas cooling over the power absorbed by dust. Considering the low efficiency values (typically less than 1%), the dust IR emission is assumed to approximate the UV photon energy absorbed by dust.

19

A small map was performed in SH (AORkey 21825792) with dedicated offset observations, but we do not consider it because of its low S/N.

20

Since the additional constraint provided by [Ne v] is related to the ionized gas irradiated by the X-ray source, we need to consider the X-ray luminosity around 4 × 1040 erg s-1 used in our modeling (Sect. 5.3).

Acknowledgments

We thank the anonymous referee for constructive feedback. V.L. was supported by a CEA/Marie Curie Eurotalents fellowship when this study was initiated. We acknowledge financial support from “Programme National de Cosmologie et Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. M.C. gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group, grant number KR4801/1-1. This research was partly supported by the Agence Nationale pour la Recherche (ANR) through the program SYMPATICO (Programme Blanc Projet ANR-11-BS56-0023). We wish to thank Jérôme Rodriguez and Giulia Migliori for their help with the XMM-Newton spectra. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Herschel/PACS spectral extraction

Appendix A.1: PACS optimal extraction of pointed observations

thumbnail Fig. A.1

Individual Herschel/PACS spaxel background-subtracted spectra for [C ii] 157 μm (top) and [O i] 63 μm (bottom). The observed spectrum is shown in black while the line fit is shown in orange. The number between parentheses indicates the spaxel coordinate within the footprint. The second column is shifted to reflect the PACS footprint projection on the sky.

Open with DEXTER

thumbnail Fig. A.2

Individual Herschel/PACS spaxel background-subtracted spectra for [O iii] 88 μm. See Fig. A.1 for the description.

Open with DEXTER

The Herschel/PACS spectra are shown in Figs. A.1 and A.2. We derive the following radial velocities, 795 ± 100 km s-1 for [C ii], 760 ± 40 km s-1 for [O i], and 795 ± 75 km s-1 for [O iii]. While the spectral resolution of Spitzer is too low to provide a useful comparison (>500 km s-1; Sect. 3.2), Lelli et al. (2012) modeled the H i 21 cm velocity distribution of I Zw 18 main body as a disk with a centroid close to the SE component and a central velocity of 765 km s-1. The velocity ranges from 735 km s-1 to 790 km s-1 from north to south. Interferometric Hα measurements are in agreement with these values (Petrosian et al. 1997).

The optimal extraction uses the “as-built” telescope monochromatic PSF17 projected on the PACS footprint. The flux is calibrated by ensuring that the encircled energy fraction is similar to that in the PACS photometer documentation. The algorithm does not take into account optical effects within the spectrometer. Nevertheless, the spectrometer PSF is to first order similar to the telescope PSF and it is considered to be a satisfactory model when used with a proper flux calibration. The flux is determined by scaling the normalized PSF to the data. The scaling is performed by using a multi-linear regression algorithm, following the method used for the Spitzer/IRS spectra in Lebouteiller et al. (2010). The spaxels located far from the source centroid do not significantly affect the optimal extraction since the latter follows the PSF profile. We have validated the algorithm on point sources from the HERUS sample (Farrah et al. 2013).

The optimal extraction for [C ii], [O i], and [O iii] in I Zw 18 is shown in Fig. A.3. From the spatial profiles, we can estimate whether a source is point-like or not. Practically, only two or three detected spaxels are sufficient to show evidence of significantly broadened sources (for a single source in the footprint). The emission appears point-like for the three PACS lines. The intrinsic source extent is estimated to be 6′′.

We can also quantify the contribution from a large-scale low surface brightness component by adding a flat emission component in the optimal extraction algorithm. The large-scale component is fitted simultaneously with the point source. The flux of this large-scale component represents 10%, 18%, and 7% of the total line flux for [C ii], [O i], and [O iii], respectively. We conclude that our measurements of the point source are not significantly affected by a potential low surface brightness component. This finding is also consistent the point-like appearance of I Zw 18 in the PACS photometry 70μm and 100μm bands (with a similar PSF size as the PACS spectroscopy observations) (Rémy-Ruyer et al. 2015).

thumbnail Fig. A.3

Optimal extraction of [C ii] (top), [O i] (middle), and [O iii] (bottom). For each line the top panel shows the observed footprint (with crossed spaxels indicating non-detection) and the template footprint (chosen from a library of precomputed projected PSFs). The footprint is shown as a regular grid for display purposes only; in reality, the spaxels are not uniform in size and one row is shifted with respect to the others. The bottom panel shows the spaxel fluxes ordered by flux. Observed values are indicated in orange. The theoretical distribution from the template footprint is shown with the black curve. The 3 brightest spaxels suffice to locate the peak.

Open with DEXTER

Appendix A.2: Unchopped [O  III] 88 μm raster map

I Zw 18 was observed independently with Herschel/PACS on April 6th, 2013 as part of program OT2_kcroxall_1 (PI K. Croxall), for 2.7 ks. The observation (OBSID 1342269442) is a spectral map of the [O iii] 88μm line. Contrary to our DGS observations in which observations were performed using the chop/nod mode (with a fast removal of the telescope baseline, Sect. 3.1.1), OBSID 1342269442 uses the unchopped grating scan mode in which an offset position was observed before and after the on-source exposure. The map consists of four raster positions shifted by half a spaxel, allowing a better spatial sampling of the line peak emission.

We reduced the data and analyzed it with PACSman. The PACS spectral map projection algorithm is explained in Lebouteiller et al. (2012). The final map is shown in Fig. A.4. The emission corresponds to a compact source located at 09h34m02.12s/55°1427.4, i.e., coinciding with the NW cluster (Fig. A.4). The source appears mostly point-like, with a slight elongation along the NW-SE axis. This result is consistent with the Spitzer/IRS results showing that most of the IR ionized gas line emission arises in the NW region and that a relatively lower fraction arises in the SE region (Sect. 3.2).

The integrated source flux in the reconstructed map is 2.05 ± 0.25 × 10-17 W m-2, close to the value derived from our observation (OBSID 1342253758; Table 2). For testing purposes we have also extracted the flux in each individual raster position of the map using the methods described in Sect. 3.1.2. Only Fopt (optimal extraction) could not be used because only one spaxel is significantly detected in each raster position. Despite the lower S/N as compared to OBSID 1342253758, we find that the various determinations are overall in good agreement (Table A.1). The error bars on the [O iii] flux from the spectral map observation may be underestimated because of the presence of systematic uncertainties in the baseline used for the fit. Such uncertainties are better canceled in chop/nod observations than in unchopped scan observations because of the fast removal of the telescope background (see PACS Observer’s manual18). In our analysis we only use the chop/nod [O iii] observations described in Sect. 3.1.1.

Table A.1

[O iii] 88μm line flux determinations from the spectral map.

thumbnail Fig. A.4

The [O iii] map from the program OT2_kcroxall_1 (contours) obtained with Herschel/PACS. The background is the Hubble/ACS F555W image. The beam is shown in the bottom-left.

Open with DEXTER

Appendix B: Spitzer/IRS line flux measurements

The apertures of the various Spitzer/IRS observations19 are overlaid in Fig. B.1. It can be seen that the LL (10.7″ × 168″) and LH (11.1″ × 22.3″) observations include most of the emission from the two main H ii regions. Concerning the short wavelength observations, AORkey 16205568 is better oriented in both SL (3.7″ × 57″) and SH (4.7″ × 11.3″) as compared to the other observations. We therefore consider AORkey 16205568 as our reference observation for SL and SH, while all observations can be used in principle for LL and LH.

thumbnail Fig. B.1

Overlay of apertures on an Hα image (from de Paz et al. 2003). Dashed colored rectangles correspond to the low-resolution apertures IRS SL (3.6″ × 57″) and LL (10.5″ × 168″). Solid colored rectangles correspond to the high-resolution apertures SH (4.7″ × 11.3″) and LH (11.1″ × 22.3″). The large black rectangle corresponds to the PACS footprint (47″ × 47″).

Open with DEXTER

We use the extracted spectra from the Combined Atlas of Spitzer/IRS Sources (CASSIS), which is described in Lebouteiller et al. (2011) for low-resolution data and Lebouteiller et al. (2015) for high-resolution data. For the low-resolution observations, the background sky emission was calculated via the image corresponding to the other nod position. For the high-resolution observations, no offset image was subtracted, but the background emission was removed during spectral extraction, as explained in Lebouteiller et al. (2015).

In the following, we compare regular extraction (i.e., full aperture for high-resolution data and tapered column for low-resolution data) and optimal extraction methods. The regular method integrates the flux within the extraction window and uses either a point-source calibration (accounting for light losses outside of the aperture due to the PSF size) or extended source calibration (assuming uniformly extended emission within and outside the aperture). In practice, we find that the source is quasi-point-like in the long wavelength observations (LL, LH) and that it is somewhat extended in the short wavelength observations (SL, SH). Despite the noticeable extent, and for lack of a partially extended source calibration, the point-source calibration is still much more adequate than the extended source calibration. We consider hereafter that the regular extractions of LL and LH observations (with a point-source calibration) reflect the total emission from both NW and SE. On the other hand, we consider that the regular extraction of SL and SH might miss a fraction of the total flux.

Optimal extraction is adapted to point sources only and, as such, it is a reliable method for the LL and LH observations for which only a minor correction due to the source extent may be necessary. The flux calibration already accounts for the aperture loss outside the slit. Therefore, the LL and LH optimal extraction are assumed to reflect the total emission from the NW and SE regions. For SL and SH, optimal extraction is only an approximation because the source is noticeably extended in these modules.

We have also used a custom optimal extraction to extract the NW and SE regions simultaneously in the low-resolution observation 16205568 (6th column in Table B.1) using the SMART/AdOpt manual extraction tool described in Lebouteiller et al. (2010). With this special algorithm, the PSF is computed at the source position and broadened to reflect the intrinsic source broadening, so that the extracted spectra are calibrated in a reliable way and reflect the total source flux. The corresponding measurement therefore provides to first order the total NW+SE H ii region emission.

thumbnail Fig. B.2

Spectral fits for the most important lines used in this study. Several spectra are available for each Spitzer/IRS line in both low- and high-resolution modes with line widths of ~3000 km s-1 and ~400 km s-1 respectively. We selected one spectrum per line for illustrative purposes. The diamonds show the data and the connected diamonds represent the range used for fitting the line. The many light gray curves show Monte-Carlo iterations for determining the line flux uncertainty. The red curve shows the final fit.

Open with DEXTER

Figure B.2 shows some line fits for illustration. We list in Tables B.1 and B.2 the line measurements for each observation. The fact that we obtain similar results in both resolution modes, despite the different slit/aperture sizes, implies that we are successfully recovering the total flux of the NW+SE H ii regions. The final adopted fluxes are given in Table B.3, along with a comparison to W07.

Table B.1

Spitzer/IRS low-resolution measurements.

Table B.2

Spitzer/IRS high-resolution measurements.

Table B.3

Spitzer/IRS line fluxes.

Appendix C: X-ray source intrinsic spectrum

In this section, we explore the possible intrinsic X-ray spectra that are compatible with the XMM-Newton observation, the observed optical and IR line fluxes (e.g., [Ne v], [C ii], [O i]), and the H0 column density. All tests are performed with xspec (Arnaud 1996). We are mostly interested in photons with energies between ~0.1−2 keV because they are absorbed in the H i region. Lower energy photons are absorbed in the H ii region and higher energy photons freely cross the H ii + H i regions. For this reason, one must bear in mind that the fits to the XMM-Newton observation are actually mostly constrained by the well detected hard X-ray component that is not absorbed. Recovering the soft-X-ray component heavily depends on the model chosen that ties the hard- and soft-X-ray components in a physically motivated way.

We first consider a diskbb model. For the attenuation in the Milky Way, we used the foreground Galactic extinction 0.091 (Schlegel et al. 1998; Schlafly & Finkbeiner 2011) and the AV-to-column density conversion of Guver & Ozel (2009). We find NH ≈ 2 × 1020 cm-2, in good agreement with the value measured at 21 cm (e.g., Stark et al. 1992; Kalberla et al. 2005). The total absorbing hydrogen column density in the Milky Way (including molecular hydrogen) is 2.7 × 1020 cm-2 (Willingale et al. 2013). Figure C.1 shows the diskbb model and the absorbed spectrum. Most of the absorption occurs in I Zw 18, with a column density of 1021 cm-2. The luminosity between 0.3−1 keV is 1040 erg s-1 and the temperature at the inner edge of the accretion disk is 1.04 keV. These results are similar to those found by Kaaret & Feng (2013). These values are unchanged when using the abundances measured in absorption by HST/COS rather than abundances measured in emission in the H ii regions (Sect. 2.4). The column density for the diskbb model is somewhat lower than expected based on H i absorption measurements. The X-ray source is located within the stellar cluster (Sect. 3.3), toward which the H i column density has been calculated to be 2.2 × 1021 cm-2 (Aloisi et al. 2003; Lecavelier des Etangs et al. 2004; Lebouteiller et al. 2013). The column density measured in absorption, however, represents an average value (toward the most FUV-bright stars within the 2″ COS aperture), which may differ from that derived along the line of sight of the X-ray point source.

thumbnail Fig. C.1

X-ray spectrum of I Zw 18. The X-ray emission is calculated with a diskbb model (distribution of blackbodies from an accretion disk). The diamonds show the unfolded XMM-Newton spectrum. The dashed red line shows the spectrum absorbed by I Zw 18 alone while the dotted red line shows the spectrum absorbed by the Milky Way alone. The top red solid line shows the intrinsic (unabsorbed) X-ray spectrum. The thick black curve shows the sum of the radiation field components (two blackbodies in the optical and the X-ray source).

Open with DEXTER

While the diskbb model was preferred in the present study, combinations of blackbodies were also considered. Such combinations allow us to fully explore the range of acceptable spectra without strong assumptions about the nature of the source. On the other hand, using blackbodies ensures that physically unrealistic sharp features and/or discontinuities are discarded from the analysis. The single 2 × 106 K blackbody used in P08 neglected the hard X-ray tail (2 keV) which barely impacts on the H i region properties, yet it was a coarse first approximation, given the X-ray data now available.

Using a sum of four blackbodies with temperatures between 0.05−1.20 keV, and column densities between 0−1 × 1022 cm-2 proved to provide sufficient flexibility (nine free parameters). Figure C.2a shows a broad selection of combinations providing the best χ2. Unlike for the diskbb model, owing to the lack of physical ties between the hard and soft X-ray components, the soft X-rays are little constrained and a degeneracy develops between the coldest blackbody fluxes and the (here poorly defined) column density. We constrain further the column density in the X-ray model assuming that the X-ray source lies within the NW stellar cluster, hence assuming a column density of 1−2.5 × 1021 cm-2 (Fig. C.2b). Two different spectral shapes emerge, one with a soft X-ray excess and one without. In Fig. C.2b, the X-ray luminosity range of the models is 1.1−2.5 × 1040 erg s-1. Below ~0.3 keV, the freedom left by the xspec approach with a blackbody combination is considerable.

Some of the degeneracy can be lifted using signatures of the X-ray emission on the surrounding ISM. Our standard photoionization model, in which the diskbb X-ray spectrum is adopted, predicts fluxes 3−5 times weaker than the observational upper limits for the high-ionization lines [Ne v] 14.3μm and [Ne v] 24.3μm (Table 8). Also, no detection has been reported for [Ne v] 3426 Å in I Zw 18. Producing [Ne v] in the H ii region requires photon energies 100 eV. In Fig. C.2b, the soft X-ray excess allowed by xspec around 100 eV can exceed the relatively flat diskbb model spectrum by more than one order of magnitude. Obviously, at least some of these blackbody combinations produce much stronger [Ne v] lines than diskbb. Figure C.2c shows the combination spectra allowed when the additional constraint brought by the upper limit on [Ne v] is applied to photoionization model results20.

This illustrates how detailed photoionization modeling of the ISM around X-ray sources may provide powerful constraints on the low-energy tail of X-ray source spectra and therefore on the nature and properties of such sources. A preliminary inference from the present analysis is that the soft X-ray spectrum of the I Zw 18-NW source cannot much exceed that corresponding to the diskbb model. Not only is the solution above ~0.4 keV well defined, but it turns out that, within a factor of ~2.5, the set of blackbody combinations follows the diskbb spectrum down to less than 0.1 keV. The fall of the blackbody combinations toward low energies, essentially an artefact of the fitting assumptions, is of no consequence since the stellar continuum then dominates.

thumbnail Fig. C.2

Test of the 4-blackbody combinations with xspec. The open diamonds show the diskbb model while the filled diamonds show the XMM-Newton unfolded spectrum. The set of best blackbody combinations with column density <1022 cm-2 are shown in panel a). A subset with column densities between 1−2.5 × 1021 cm-2 is shown in panel b). Only in the combinations left in panel c) are the associated photoionization models fulfilling the observed upper limit on the [Ne v] flux. For simplicity, we plot the X-ray spectra of panel c) with the same luminosity as for the other panels, but in reality an increased luminosity is used, according to Sect. 5.3. In each panel, the colors indicate the relative variation of the χ2 among the plotted combinations (lower χ2 in red, larger χ2 in blue).

Open with DEXTER

All Tables

Table 1

Main I Zw 18 properties used in this study.

Table 2

Herschel/PACS line flux determinations.

Table 3

Model summary and updates from P08.

Table 4

Comparison between observed extinction-corrected optical and UV line fluxes and models.

Table 5

Factors to normalize fluxes in W m-2 to Hβ = 1000.

Table 6

Infrared line fluxes scaled to Hβ = 1000.

Table 7

Cloudy model predictions for NW.

Table 8

Comparison between the observed IR line fluxes and models.

Table A.1

[O iii] 88μm line flux determinations from the spectral map.

Table B.1

Spitzer/IRS low-resolution measurements.

Table B.2

Spitzer/IRS high-resolution measurements.

Table B.3

Spitzer/IRS line fluxes.

All Figures

thumbnail Fig. 1

H i column density contours from Lelli et al. (2012), with 2′′ resolution (beam size in the bottom left). Contours are drawn for 3 (dashed), 6, 9, 12, and 15 × 1021 cm-2. The largest red circle shows the PACS beam at 157μm ([C ii]) and the smallest red circle shows the beam at 63μm ([O i]). Both beams are centered at the emission centroid derived by the PACS optimal extraction method (Sect. 3.1.1). The green cross shows the location of the X-ray point source (Sect. 3.3). The background image is HST/ACS F555W. The NW region coincides with an H i hole. The H i column density peak lies between NW and SE.

Open with DEXTER
In the text
thumbnail Fig. 2

Herschel/PACS map of [C ii] (top), [O i] (middle), and [O iii] (bottom) emission in I Zw 18. The 25 spaxels of the PACS footprint are overplotted on the F555W HST/ACS image. The contours show the H i column density at 5″ resolution (Lelli et al. 2012). For display purposes, only the fits are shown for each spaxel with detection level >2σ, and the number indicates the detection level in σ. The red circle shows the beam size and the red cross shows the emission centroid as calculated by the optimal extraction (Sect. 3.1.2). Individual spaxel spectra are presented in Appendix A.1.

Open with DEXTER
In the text
thumbnail Fig. 3

Cross-dispersion profiles along the SL slit of the Spitzer/IRS observation 16205568. The histogram shows the emission as a function of the spatial position along the slit, in pixel units (1 px = 1.8′′). The spatial profile is modeled by two slightly extended sources, one of which corresponds to NW (blue, left) and the other to SE (red, right). The profile is shown for the entire integrated spectral order of module SL1 (~8−14μm; top) and for the [S iv] 10.5μm line only (bottom). Left corresponds to north in Fig. 4.

Open with DEXTER
In the text
thumbnail Fig. 4

Spatial positions of the two components seen in the IRS SL slit. The slit is shown with the white rectangle and the circles indicate the locations of the sources shown in Fig. 3. The size of the circles corresponds to the total FWHM of each source (i.e., including both the instrument PSF and the intrinsic broadening). The background image is the B band from the DSS2 survey.

Open with DEXTER
In the text
thumbnail Fig. 5

Photometric maps of I Zw 18 with Spitzer and Herschel (Rémy-Ruyer et al. 2015). The first panel shows the HST/ACS F555W image for reference, with the NW and SE regions circled in white. The dashed black circle indicates the beam size (from 2′′ for IRAC 8μm to 12′′ for PACS 160μm). The 8 μm band is not dominated by PAH emission in I Zw 18 but by stochastically heated small grains and warm grains in thermal equilibrium with the interstellar radiation field.

Open with DEXTER
In the text
thumbnail Fig. 6

Description of the sectors used in our models. The first panel shows the topology used in P08 and in most of our models (Table 7). The arrows reaching out of sectors #2 and #3 illustrate the fact that the fraction of escaping photon is significant. Optical depths at 1 Ryd are 1.18 and 0.05 for sectors #2 and #3 respectively. “CF” stands for covering factor. The sectors are drawn according to their respective covering factors but in reality the lines of sight are intermixed.

Open with DEXTER
In the text
thumbnail Fig. 7

Blue curves show the input radiation field used in P08, comprised of three blackbodies with temperature 4 × 104, 8 × 104, and 2 × 106 K. The dashed curve is for LX = 4 × 1039 erg s-1 (P08 model M2X) and the solid curve is for LX = 8 × 1039 erg s-1 (P08 model M2X2). The red curves show the radiation field used in this study, with additional UV-optical contribution (from Lebouteiller et al. 2013) and with an improved X-ray spectrum prescription as compared to P08. The dashed curve is for LX = 1.4 × 1040 erg s-1 (luminosity inferred from observations) and the solid curve is for LX = 4 × 1040 erg s-1 (adopted standard). The black diamonds show the unfolded XMM-Newton spectrum (see Sect. 5.3 for more details).

Open with DEXTER
In the text
thumbnail Fig. 8

Exploration of models with X-ray heating only (no dust and no CR). The predicted line fluxes (model/observation) are plotted against each model parameter with from top to bottom, the dust mass and H0 mass (log M), minimum temperature (K) and maximum density (log  cm-3) in the radiation-bounded sector, and the X-ray luminosity (log L). Each point represents a model, and the symbol size is inversely proportional to the residuals between model predictions and observations for [C ii], [O i], and [Si ii]. The light gray zone corresponds to the observed values and their associated error bar. The dark gray zone corresponds to the best models.

Open with DEXTER
In the text
thumbnail Fig. 9

Models with X-ray and photoelectric effect heating for D/G observed (left) or scaled with metallicity (right). See Fig. 8 for the plot description.

Open with DEXTER
In the text
thumbnail Fig. 10

Model results for ℳ2a (D/G observed, left) and ℳ2b (D/G scaled with metallicity, right). The plots on top show, from top to bottom, the total heating and cooling rate, the contributions to the heating rate, the contributions to the cooling rate, and the cumulative intensities of cooling lines for the radiation-bounded sector as a function of the depth within the cloud (the H ii-H i transition lies at 213 pc). For clarity, the heating and cooling mechanisms in the H ii region are not shown. The predicted dust SED is shown in the bottom panels for both models. The 20 cm measurement is thought to be dominated by synchrotron emission (see text). The vertical lines indicate the range used to calculate the total infrared luminosity (3−1100μm) and the green shaded area shows the modeled emission in this range. The blue shaded area shows the input radiation field. The model ℳ2b shows weak silicate emission bands at 10μm and 20μm that are not seen in the IRS spectrum.

Open with DEXTER
In the text
thumbnail Fig. 11

Attempts to reproduce the dust SED, ℳ2c (D/G maximum assuming standard topology), ℳ2d (radial D/G), and ℳ2e (extra sector with a relatively larger D/G). See Fig. 10 for the plot description.

Open with DEXTER
In the text
thumbnail Fig. 12

Radial density and temperature in the radiation-bounded sector of the full model ℳ4a. The plots are shown for the radiation-bounded sector.

Open with DEXTER
In the text
thumbnail Fig. 13

Heating and cooling rates for the full model (ℳ4a). The plots are shown for the radiation-bounded sector. See Fig. 10 for the plot description.

Open with DEXTER
In the text
thumbnail Fig. 14

Spectral energy distribution for model ℳ4c (model with molecular clump). See Fig. 10 for the plot description. The dust-rich sector is ignored for this model.

Open with DEXTER
In the text
thumbnail Fig. 15

Top: abundance fraction of C0, C+, CO, and H2 plotted against the optical depth at 1 Ryd for the clump sector of model ℳ4c. Bottom: cumulative fraction of [C i] 609μm, [C ii], CO(1-0), H22.12μm, and H2 mass normalized to the maximum value.

Open with DEXTER
In the text
thumbnail Fig. 16

Main levels of the C+ ion. The dot-dashed lines correspond to weak transitions.

Open with DEXTER
In the text
thumbnail Fig. 17

Observed COS cooling line rate per H atom is shown for different temperature, density, and various electron fraction (x) values.

Open with DEXTER
In the text
thumbnail Fig. A.1

Individual Herschel/PACS spaxel background-subtracted spectra for [C ii] 157 μm (top) and [O i] 63 μm (bottom). The observed spectrum is shown in black while the line fit is shown in orange. The number between parentheses indicates the spaxel coordinate within the footprint. The second column is shifted to reflect the PACS footprint projection on the sky.

Open with DEXTER
In the text
thumbnail Fig. A.2

Individual Herschel/PACS spaxel background-subtracted spectra for [O iii] 88 μm. See Fig. A.1 for the description.

Open with DEXTER
In the text
thumbnail Fig. A.3

Optimal extraction of [C ii] (top), [O i] (middle), and [O iii] (bottom). For each line the top panel shows the observed footprint (with crossed spaxels indicating non-detection) and the template footprint (chosen from a library of precomputed projected PSFs). The footprint is shown as a regular grid for display purposes only; in reality, the spaxels are not uniform in size and one row is shifted with respect to the others. The bottom panel shows the spaxel fluxes ordered by flux. Observed values are indicated in orange. The theoretical distribution from the template footprint is shown with the black curve. The 3 brightest spaxels suffice to locate the peak.

Open with DEXTER
In the text
thumbnail Fig. A.4

The [O iii] map from the program OT2_kcroxall_1 (contours) obtained with Herschel/PACS. The background is the Hubble/ACS F555W image. The beam is shown in the bottom-left.

Open with DEXTER
In the text
thumbnail Fig. B.1

Overlay of apertures on an Hα image (from de Paz et al. 2003). Dashed colored rectangles correspond to the low-resolution apertures IRS SL (3.6″ × 57″) and LL (10.5″ × 168″). Solid colored rectangles correspond to the high-resolution apertures SH (4.7″ × 11.3″) and LH (11.1″ × 22.3″). The large black rectangle corresponds to the PACS footprint (47″ × 47″).

Open with DEXTER
In the text
thumbnail Fig. B.2

Spectral fits for the most important lines used in this study. Several spectra are available for each Spitzer/IRS line in both low- and high-resolution modes with line widths of ~3000 km s-1 and ~400 km s-1 respectively. We selected one spectrum per line for illustrative purposes. The diamonds show the data and the connected diamonds represent the range used for fitting the line. The many light gray curves show Monte-Carlo iterations for determining the line flux uncertainty. The red curve shows the final fit.

Open with DEXTER
In the text
thumbnail Fig. C.1

X-ray spectrum of I Zw 18. The X-ray emission is calculated with a diskbb model (distribution of blackbodies from an accretion disk). The diamonds show the unfolded XMM-Newton spectrum. The dashed red line shows the spectrum absorbed by I Zw 18 alone while the dotted red line shows the spectrum absorbed by the Milky Way alone. The top red solid line shows the intrinsic (unabsorbed) X-ray spectrum. The thick black curve shows the sum of the radiation field components (two blackbodies in the optical and the X-ray source).

Open with DEXTER
In the text
thumbnail Fig. C.2

Test of the 4-blackbody combinations with xspec. The open diamonds show the diskbb model while the filled diamonds show the XMM-Newton unfolded spectrum. The set of best blackbody combinations with column density <1022 cm-2 are shown in panel a). A subset with column densities between 1−2.5 × 1021 cm-2 is shown in panel b). Only in the combinations left in panel c) are the associated photoionization models fulfilling the observed upper limit on the [Ne v] flux. For simplicity, we plot the X-ray spectra of panel c) with the same luminosity as for the other panels, but in reality an increased luminosity is used, according to Sect. 5.3. In each panel, the colors indicate the relative variation of the χ2 among the plotted combinations (lower χ2 in red, larger χ2 in blue).

Open with DEXTER
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.