Neutral gas heating by X-rays in primitive galaxies: Infrared observations of IZw18 with Herschel

TL;DR

This study investigates how X-ray radiation heats the neutral gas in extremely metal-poor galaxy IZw18, revealing X-ray binaries as key heating sources and challenging traditional diagnostics of star formation and gas content.

Contribution

It provides new evidence that X-ray heating dominates in primitive galaxies and assesses the diagnostic limitations of fine-structure lines in such environments.

Findings

X-ray binaries are primary heating sources in IZw18.

[CII] and [OI] lines indicate an X-ray luminosity of 4×10^40 erg/s.

H2 is negligible but could be hidden in dense clouds.

Abstract

(abridged) The dominant thermal mechanisms in the neutral interstellar medium, which acts as a star-forming gas reservoir, are uncertain in extremely metal-poor galaxies. Our objective is to identify the heating mechanisms in one such galaxy, IZw18, and assess the diagnostic value of fine-structure cooling lines. We also seek to constrain the mass of H$_2$, which, despite being an important catalyst and tracer of star formation, remains elusive in this object. Building on a previous photoionization model within a multi-sector topology, we provide additional constraints from the [CII] and [OI] lines and the dust mass recently measured with Herschel. The heating of the HI region appears to be mainly due to photoionization by radiation from a bright X-ray binary source, while photoelectric effect (PE) is negligible. The [CII] and [OI] lines imply an average X-ray luminosity of $4\times10^{40}$ erg s$^{-1}$, while the [NeV] upper limits bring strong constraints to the soft X-ray flux arising from the binary. A negligible amount of H$_2$ is predicted, but $\lesssim10^7$ M$_\odot$ of H$_2$ may be hidden in sufficiently dense clouds of order $\lesssim10$ pc in size. Regardless of the presence of significant amounts of H$_2$, [CII] and [OI] do not trace the so-called CO-dark gas, but the almost purely atomic medium. Although the [CII]+[OI]/TIR ratio is close to values found in more metal-rich sources, it cannot be safely used as a PE heating efficiency proxy. This ratio seems to be kept stable due to a correlation between the X-ray luminosity and the star-formation rate. We propose that X-ray heating is an important process in extremely metal-poor sources. The weak PE heating due to the low dust-to-gas ratio could be compensated for by the larger occurrence and power of X-ray binaries in low-metallicity galaxies. We speculate that X-ray heating may quench star formation.

Figures

Tables

Table 1. Table 1: Main I Zw 18 properties used in this study.

Parameter	Value	Reference
Coordinates (J2000)	$09$h$34$m$02.2$s, $+55$d$14$m$28.0$s
Distance	$18.2$ Mpc	(1)
M(H0) (main body)	${10}^8$ M⊙	(2,3)
M(H0) (NW)	$\approx 2\times {10}^7$ M⊙	(this study)
M(H+) (NW)	$3.4\times {10}^6$ M⊙	(4)
M∗	$1-20\times {10}^6$ M⊙	(5,6)
M∗ (NW)	$\sim 3\times {10}^5$ M⊙	(7)
M${}_{}{}^a{}_{\mathrm{dust}}{}^{}$	${562}_{-287}^{+586}$ M⊙, ${724}_{-308}^{+535}$ M⊙	(6)
D/Gb	$5_{-1.7}^{+9.3}\times {10}^{-4}$	(6)
LFIR	$2.88\pm 0.13\times {10}^7$ L⊙	(6)
LUV	$1.06\times {10}^8$ L⊙	(8)
SFR	$\approx 0.1$ M⊙ yr-1	(9)
Abundances	$12+\log \left(C/H\right)=6.59$	(4)
	$12+\log \left(N/H\right)=5.57$	(4)
	$12+\log \left(O/H\right)=7.22$	(4)
	$12+\log \left(\mathrm{Ne}/H\right)=6.41$	(4)
	$12+\log \left(\mathrm{Si}/H\right)=6.00$	(4)
	$12+\log \left(S/H\right)=5.63$	(4)
	$12+\log \left(\mathrm{Fe}/H\right)=5.79$	(4)
	$12+\log \left(\mathrm{Ar}/H\right)=5.79$	(4)

Table 2. Table 2: Herschel /PACS line flux determinations.

Method	[C ii] 157 $\mu$m	[O i] 63 $\mu$m	[O iii] 88 $\mu$m
$F_{\mathrm{opt}}$	$97\pm 23$	$76\pm 20$	$248\pm 31$
$F_1$	$(>)71\pm 7$	$(>)69\pm 12$	$(>)207\pm 14$
$F_{3\times 3}$	$108\pm 45$	$143\pm 82$	$267\pm 102$
$F_{3\times 3}^{\prime }$	$(>)81\pm 16$	$(>)80\pm 21$	$(>)217\pm 38$
Cormier et al. (2015)a	$106\pm 8$	$115\pm 21$	$284\pm 34$
Final	$97\pm 25$	$76\pm 25$	$248\pm 35$

Table 3. Table 3: Model summary and updates from P08 .

Parameter	Value/comment
Topology
Number of sectors	Fixed (3 to 5)
Covering factors	Fixed, except for dense clump sector
Inner radius	Fixed ($130$ pc)
Ionization front	Fixed ($220$ pc)
Physical conditions
Density	Pressure law, maximum value free
Stopping temperature	Minimum value free
(outer radius of H i region)	(maximum value free)
Radiation field
Blackbody	Fixed ($\log L=41.58$, $T=4\times {10}^4$ K)
Blackbody${}^{\text{a}}$	Fixed ($\log L=41.58$, $T=8\times {10}^4$ K)
X-rays	diskbb spectrum, luminosity free
Optical-UV	Added in present study (Fig. 7)
CR	Several ionization rate values tested
Chemical composition
Metal abundances	Fixed (see Table 1)
	He/H=$0.085$ (updated from P08)
D/G	Several values tested

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

Line	Obs.	Nebu${}^{\text{a}}$	$ℳ0a$	$ℳ4a$
H i 4861Å	1000	1000	1000	1000
H i 6563Å	2860	2841	2858	2854
He i 4471Å	21.4	36.6	36.9	35.2
He i 6678Å	25.3	27.0	28.1	26.7
He i 7065Å	24.4	23.9	22.4	21.8
He ii 4686Å	36.8	36.8	38.8	35.6
$[$N ii$]$ 6548+84Å	9.2	9.2	7.8	7.3
$[$O i$]$ 6300+63Å	8.5	8.6	9.9	9.0
$[$O ii$]$ 3726+9Å	238	236.8	207	203
$[$O ii$]$ 7320+30Å	6.3:	6.6	5.7	5.4
$[$O iii$]$ 4363Å	65.9	66.0	65.6	66.3
$[$O iii$]$ 5007+4959Å	2683	2679.4	2691	2653
$[$Ne iii$]$ 3869Å+968Å	191	191.2	224	219.9
$[$S ii$]$ 6716Å	22.5	18.9	20.7	19.0
$[$S ii$]$ 6731Å	16.9	14.0	15.5	14.4
$[$S ii$]$ 6716/6731	1.33	1.35	1.33	1.32
$[$S iii$]$ 6312Å	6.7	6.4	7.6	7.2
$[$S iii$]$ 9069+532Å	114	113	129.4	122
$[$Ar iii$]$ 7136+751Å	23.5	23.4	20	18.6
$[$Ar iv$]$ 4740Å	4.5	7.2	9.5	9.5
$[$Fe iii$]$ 4658Å	4.5	4.5	4.8	4.5
$[$Fe iii$]$ 4986Å	7.4	6.7	…	…
C iii$]$ 1909+07Å	467	465	430	411
Si iii$]$ 1882Å	164:	207	225	209
C iv 1549Å	512:	393.	413	…
$\log$ $L$(H$\beta$) (erg s-1)	39.43	39.43	39.36	39.43

Table 5. Table 5: Factors to normalize fluxes in W m -2 to H β = 1000 𝛽 1000 \beta=1000 .

Observation	Reference tracer(s)	Factor
Herschel	[O iii] 88 $\mu$m, aperture corrected H$\beta$	$\approx 1.25\times {10}^{-19}$
Spitzer
…SL	Hu$\alpha$ 12.37 $\mu$m, [Ne ii] 12.8 $\mu$m	$\sim 1\times {10}^{-19}$
…SH	[Ne iii] 15.6 $\mu$m	$\approx 1.3\times {10}^{-19}$
…LL	[Ne iii] 15.6 $\mu$m, [S iii] 18.7 $\mu$m	$\approx 1.5\times {10}^{-19}$
…LH	[S iii] 18.7 $\mu$m	$\approx 1.2\times {10}^{-19}$

Table 6. Table 6: Infrared line fluxes scaled to H β = 1000 𝛽 1000 \beta=1000 .

	Normalized		Adopted
Herschel
$[$C ii$]$ 157 $\mu$m	$76\pm 20$		$76\pm 20$
$[$O i$]$ 63 $\mu$m	$60\pm 20$		$60\pm 20$
$[$O iii$]$ 88 $\mu$m	$198\pm 30$		$200\pm 30$
Spitzer	SL/LL	SH/LH
H i Hu$\alpha$ $12.37$ $\mu$m	$8.5\pm 3$		$8.5\pm 3$
$[$O iv$]$ 25.9 $\mu$m	$40\pm 10$	$37.5\pm 8$	$38\pm 8$
	(+[Fe ii] 26.0 $\mu$m)
$[$Ne ii$]$ 12.8 $\mu$m	$6\pm 3$	$\sim 6$	$6\pm 3$
$[$Ne iii$]$ 15.5 $\mu$m	$48\pm 10$	$46\pm 7.5$	$47\pm 10$
$[$Ne iii$]$ 36.0 $\mu$m
$[$Ne v$]$ 14.3 $\mu$m
$[$Ne v$]$ 24.3 $\mu$m
$[$Si ii$]$ 34.8 $\mu$m	$63\pm 13$	$67\pm 25$	$65\pm 12$
$[$S iii$]$ 18.7 $\mu$m	$27\pm 7$	$29\pm 8$	$28\pm 7$
$[$S iii$]$ 33.5 $\mu$m	$37\pm 7$	$42\pm 8$	$40\pm 7$
$[$S iv$]$ 10.5 $\mu$m	$43\pm 5$	$53.8\pm 15.4$	$50\pm 10$
$[$Ar ii$]$ 7.0 $\mu$m		…
$[$Ar iii$]$ 9.0 $\mu$m		…
$[$Ar iii$]$ 21.8 $\mu$m
$[$Fe ii$]$ 17.9 $\mu$m
$[$Fe ii$]$ 26.0 $\mu$m	$40\pm 10$	$12.5\pm 8:$	$13\pm 8:$
	(+[O iv] 25.9 $\mu$m)
$[$Fe iii$]$ 23.0 $\mu$m
H2 S(0) 28.2 $\mu$m
H2 S(1) 17.0 $\mu$m
H2 S(2) 12.29 $\mu$m	…

Table 7. Table 7: Cloudy model predictions for NW.

			Input parameters			Output parameters
Model	Description	Conf.a	$L_X$	CRb	D/G	$n_{\max }^c$	$T_{\mathrm{stop}}^d$	$M\left(H^0\right)$	$M_{\mathrm{dust}}$	[C ii]	[O i]	[Si ii]	${\Gamma }_{\mathrm{tot}}{\left(\mathrm{HI}\right)}^f$
			($\log$ erg s-1)	($\times$ Gal.)	($\times$ Gal.)	(cm-3)	(K)	($\log$ M⊙)	($\log$ M⊙)	(M/Oe)	(M/Oe)	(M/Oe)	PE	CR
Exploratory models
$ℳ0a$	($\approx$P08’s $M2$)	(a)	$0$	$0$	$0$	$150$	${50}^g$	$3.7$	$0$	$0.01$	$0.01$	$0.27$	$0$%	$0$%
$\dots \dots b$	Dust (Z-scaled D/G)	(a)	$0$	$0$	$1/50$	$50$	${50}^g$	$4.7$	$2.6$	$0.02$	$0.06$	$0.35$	$12$%	$0$%
$ℳ1$	XR	(a)	$40.6$	0	0	$150$	$100$	$7.2$	$0$	$1.01$	$1.10$	$0.81$	$0$%	$0$%
$ℳ2a$	XR + dust (observed D/G)	(a)	$40.6$	$0$	$1/1000$	$200$	$95$	$7.2$	$2.2$	$1.00$	$1.07$	$0.77$	$4$%	$0$%
$\dots \dots b$	XR + dust (Z-scaled D/G)	(a)	$40.3$	$0$	$1/50$	$250$	$110$	$7.2$	$3.4$	$1.01$	$1.05$	$0.62$	$57$%	$0$%
$\dots \dots c$	XR + dust (maximum D/G)	(a)	$40.6$	$0$	$1/300$	$200$	$100$	$7.2$	$2.6$	$1.00$	$1.08$	$0.77$	$10$%	$0$%
$\dots \dots d$	XR + dust (radial D/G)	(a)	$40.6$	$0$	$1/\left(150\searrow 750\right)$	$200$	$105$	$7.2$	$3.0$	$1.02$	$1.13$	$0.80$	$12$%	$0$%
$\dots \dots e$	XR + dust (extra sector)	(b)	$40.6$	$0$	$1/(1000,50)$	$150$	$90$	$7.4$	$2.3$	$1.02$	$1.02$	$0.75$	$4$%	$0$%
$ℳ3a$	XR + CR	(a)	$40.6$	$1$	$0$	$200$	$110$	$7.2$	$0$	$1.01$	$1.13$	$0.81$	$0$%	$13$%
$\dots \dots b$	XR + CR	(a)	$40.1$	$5$	$0$	$250$	$125$	$7.2$	$0$	$1.06$	$1.00$	$0.65$	$0$%	$65$%
$\dots \dots c$	CR only	(a)	$0$	$1$	$0$	$25$	$230$	$8.1$	$0$	$1.04$	$1.23$	$0.80$	$0$%	$96$%
$\dots \dots d$	CR only	(a)	$0$	$5$	$0$	$250$	$150$	$7.3$	$0$	$1.02$	$1.01$	$0.66$	$0$%	$96$%
Full models
$ℳ4a$	Best parameter set	(b)	$40.6$	$0.1$	$1/(1000,50)$	$250$	$100$	$7.4$	$2.3$	$1.05$	$1.05$	$0.77$	$3$%	$0.7$%
$\dots \dots b$	$ℳ4a$ + molecules	(b)	$40.6$	$0.1$	$1/(1000,50)$	$250$	$80$	$7.3$	$2.2$	$1.03$	$1.00$	$0.71$	$3$%	$1$%
$\dots \dots c$	$ℳ4b$ + clumps	(c)	$40.5$	$0.1$	$1/1000$	$150$	$90$	$7.3$	$2.6$	$1.00$	$1.12$	$0.84$	$3$%	$2$%
													($23$%)	($17$%)

Table 8. Table 8: Comparison between the observed IR line fluxes and models.

Line	Observed	Nebu${}^{\text{a}}$	$ℳ4a$	$ℳ4b$	$ℳ4c$
Ionized gas
H i Hu$\alpha$ $12.37$ $\mu$m	$8.5\pm 3$	9.5	8.9	8.9	9.1
$[$O iii$]$ 88 $\mu$m	$200\pm 30$	207.9	188.0	188.1	188.7
$[$O iv$]$ 25.9 $\mu$m	$38\pm 8$	49.6	45.3	45.3	43.6
$[$Ne ii$]$ 12.8 $\mu$m	$6\pm 3$	3.0	2.0	2.0	2.1
$[$Ne iii$]$ 15.6 $\mu$m	$47\pm 10$	47.0	45.7	45.7	45.9
$[$Ne iii$]$ 36.0 $\mu$m		4.2	4.1	4.1	4.2
$[$Ne v$]$ 14.3 $\mu$m		0.7	0.8	0.8	0.6
$[$Ne v$]$ 24.3 $\mu$m		0.8	1.0	1.0	0.7
$[$S iii$]$ 18.7 $\mu$m	$28\pm 7$	28.3	23.3	23.3	24.3
$[$S iii$]$ 33.5 $\mu$m	$40\pm 7$	51.6	47.7	47.7	50.0
$[$S iv$]$ 10.5 $\mu$m	$49\pm 10$	117.7	107.7	107.7	103.2
$[$Ar ii$]$ 7.0 $\mu$m		1.4	0.5	0.5	0.5
$[$Ar iii$]$ 9.0 $\mu$m		8.1	6.3	6.3	6.6
$[$Ar iii$]$ 21.8 $\mu$m		0.6	0.5	0.5	0.5
$[$Fe iii$]$ 22.9 $\mu$m		3.2	2.8	2.8	3.1
Neutral gas
$[$C i$]$ 609 $\mu$m	…	…	$3\times {10}^{-4}$	$3\times {10}^{-4}$	$0.2$
$[$C i$]$ 370 $\mu$m	…	…	$6\times {10}^{-4}$	$6\times {10}^{-4}$	$0.8$
$[$C ii$]$ 157 $\mu$m	$76\pm 20$	$77.1$	79.0	77.0	77.9
$[$O i$]$ 63 $\mu$m	$60\pm 20$	$60.5$	63.3	63.2	68.3
$[$Si ii$]$ 34.8 $\mu$m	$65\pm 12$	$52.6$	49.9	49.9	54.7
$[$Fe ii$]$ 17.9 $\mu$m		$1.6$	1.1	1.1	1.1
$[$Fe ii$]$ 26.0 $\mu$m	$13\pm 8$	$9.6$	4.4	4.4	4.7
Molecular gas
H2 0-0 S(0) 28.2 $\mu$m		…	…	$7\cdot {10}^{-6}$	$2.7$
H2 0-0 S(1) 17.0 $\mu$m		…	…	$2\times {10}^{-4}$	$0.1$
H2 0-0 S(2) 12.3 $\mu$m		…	…	$2\times {10}^{-6}$	$5\times {10}^{-5}$
H2 1-0 S(1) 2.12 $\mu$m		…	..	$7\times {10}^{-4}$	$3\times {10}^{-3}$
CO(1-0)		…	…	$2\times {10}^{-11}$	$8\times {10}^{-3}$
PAH (6-15 $\mu$m)		…	…	200.0	212.4

Table 9. Table 9: [O  iii ] 88 88 88   μ 𝜇 \mu m line flux determinations from the spectral map.

Obs. & Method	$F_1$	$F_{3\times 3}$	$F_{3\times 3}^{\prime }$
Raster 1	$(>)144\pm 33$	$133\pm 104$	$(>)367\pm 104$
Raster 2	$(>)150\pm 23$	$335\pm 192$	$(>)411\pm 133$
Raster 3	$(>)236\pm 32$	$239\pm 173$	$(>)505\pm 108$
Raster 4	$(>)191\pm 22$	$270\pm 203$	$(>)251\pm 59$
Map	Aperture extraction: $205\pm 25$

Table 10. Table 10: Spitzer /IRS low-resolution measurements.

Method	Regular	Optimal	Regular	Optimal	Optimal custom${}^{\text{a}}$	Regular	Optimal
SL	9008640	9008640	16205568	16205568	16205568	12622848	12622848
	(Ap.)	(Ap., meth.)	(Pref.)	(Meth.)	(Pref.)	(Ap.)	(Ap., meth.)
H i Hu$\alpha$ $12.37$ $\mu$m			$8.2\pm 2.5$	$4.1\pm 1.5$	$2.6\pm 2.0,5.5\pm 1.8$	$7.0\pm 3.5$	$6.5\pm 3.3$
$[$Ne ii$]$ 12.8 $\mu$m				$4.1\pm 1.7$	$3.0\pm 1.6,5.2\pm 2.6$	$8.5\pm 3.5$	$6.0\pm 3.1$
$[$S iv$]$ 10.5 $\mu$m			$40\pm 3$	$19\pm 2$	$14\pm 2,34\pm 2$	$39\pm 4$	$29\pm 3$
$[$Ar ii$]$ 7.0 $\mu$m	$\dots$	$\dots$			$\dots$
$[$Ar iii$]$ 9.0 $\mu$m
LL	9008640	9008640	16205568	16205568	16205568	12622848	12622848
	(Pref.)	(Pref., meth.)	(Pref.)	(Pref., meth.)	(Pref.)	(Pref.)	(Pref., meth.)
$[$O iv$]$+$[$Fe ii$]$	$66\pm 19$	$52\pm 12$	$53\pm 22$	$43\pm 11$	$53\pm 11$	$79\pm 20$	$60\pm 15$
$[$Ne iii$]$ 15.5 $\mu$m	$71\pm 13$	$66\pm 17$	$83\pm 8$	$65\pm 8$	$76\pm 10$	$104\pm 41\ast$	$59\pm 13$
$[$Ne iii$]$ 36.0 $\mu$m	$\dots$				$46\pm 24\ast$	$135\pm 70\ast$	$60\pm 24\ast$
$[$Ne v$]$ 14.3 $\mu$m
$[$Ne v$]$ 24.3 $\mu$m
$[$S iii$]$ 18.7 $\mu$m	$39\pm 17$	$46\pm 15$	$43\pm 6$	$40\pm 6$	$47\pm 8$	$42\pm 17$	$56\pm 16$
$[$S iii$]$ 33.5 $\mu$m		$65\pm 16$	$60\pm 10$	$45\pm 13$	$48\pm 12$	$91\pm 30$	$57\pm 11$
$[$Ar iii$]$ 21.8 $\mu$m
$[$Fe iii$]$ 23.0 $\mu$m
$[$Si ii$]$ 34.8 $\mu$m		$96\pm 30$	$96\pm 14$	$50\pm 11$	$65\pm 12$	$130\pm 63$	$115\pm 24$
$[$Fe ii$]$ 17.9 $\mu$m
H2 S(0) 28.2 $\mu$m
H2 S(1) 17.0 $\mu$m

Table 11. Table 11: Spitzer /IRS high-resolution measurements.

Method	Regular (full)	Optimal	Regular (full)	Optimal
SH	9008640	9008640	16205568	16205568
	(Ap.)	(Ap., meth.)	(Pref.)	(Meth.)
H i Hu$\alpha$ $12.37$ $\mu$m	$12\pm 6,18\pm 10$
$[$Ne ii$]$ 12.8 $\mu$m	$17\pm 3\ast ,13\pm 3\ast$		$\dots ,\dots$
$[$S iv$]$ 10.5 $\mu$m	$56\pm 5,37\pm 4$	$14\pm 3,22\pm 4$	$66\pm 22,83\pm 25$	$23\pm 2,26\pm 1$
$[$Ne v$]$ 14.3 $\mu$m
$[$Ne iii$]$ 15.5 $\mu$m	$47\pm 4,44\pm 4$	$28\pm 1,28\pm 1$	$51\pm 7,78\pm 7$	$21\pm 3,28\pm 7$
$[$S iii$]$ 18.7 $\mu$m	$25\pm 3,32\pm 4$	$13\pm 2\ast ,13\pm 2\ast$	$37\pm 4,31\pm 4$	$21\pm 3,25\pm 2$
$[$Fe ii$]$ 17.9 $\mu$m
H2 S(1) 17.0 $\mu$m			$\dots ,\dots$	$\dots ,\dots$
H2 S(2) 12.3 $\mu$m
LH	9008640	9008640	16205568	16205568
	(Pref.)	(Pref., meth.)	(Pref.)	(Pref., meth.)
$[$O iv$]$ 25.89 $\mu$m	$35\pm 8,49\pm 11\ast$	$29\pm 6,25\pm 4$
$[$Ne iii$]$ 36.0 $\mu$m		$\dots ,\dots$
$[$Ne v$]$ 24.3 $\mu$m
$[$S iii$]$ 33.5 $\mu$m	$59\pm 23,46\pm 22$		$38\pm 11,52\pm 16$	$58\pm 8,43\pm 8$
$[$Ar iii$]$ 21.8 $\mu$m
$[$Fe iii$]$ 23.0 $\mu$m
$[$Si ii$]$ 34.8 $\mu$m	$64\pm 14,88\pm 16$	$65\pm 6,52\pm 10$	$111\pm 50,117\pm 32$	$95\pm 2\ast ,61\pm 10\ast$
$[$Fe ii$]$ 26.0 $\mu$m	$34\pm 10,33\pm 11$		$16\pm 7,13\pm 5$
H2 S(0) 28.2 $\mu$m	$\dots ,\dots$

Table 12. Table 12: Spitzer /IRS line fluxes.

	W07	Low-resolution		High-resolution
H i Hu$\alpha$ $12.37$ $\mu$m	…	$8.5\pm 3^{\text{a}}$	SL		SH
$[$O iv$]$ 25.89 $\mu$m	$49\pm 3$	$60\pm 15$	LL	$45\pm 10$	LH
		(+[Fe ii])
$[$Ne ii$]$ 12.8 $\mu$m	$9\pm 1$	$6\pm 3$	SL	$\sim 8^{\text{b}}$	SH
$[$Ne iii$]$ 15.5 $\mu$m	$46\pm 2$	$70\pm 15$	LL	$60\pm 10$	SH
$[$Ne iii$]$ 36.0 $\mu$m	…		LL		LH
$[$Ne v$]$ 14.3 $\mu$m	…		LL		SH
$[$Ne v$]$ 24.3 $\mu$m	…		LL		LH
$[$Si ii$]$ 34.8 $\mu$m	$158\pm 8$	$95\pm 20$	LL	$80\pm 30$	LH
$[$S iii$]$ 18.7 $\mu$m	$23\pm 2$	$40\pm 10$	LL	$35\pm 10$	LH
$[$S iii$]$ 33.5 $\mu$m	$120\pm 12$	$55\pm 10$	LL	$50\pm 10$	LH
$[$S iv$]$ 10.5 $\mu$m	$48\pm 3$	$43\pm 5$	SL	$70\pm 20$	SH
$[$Ar ii$]$ 7.0 $\mu$m	…		SL	…	…
$[$Ar iii$]$ 9.0 $\mu$m	…		SL	…	…
$[$Ar iii$]$ 21.8 $\mu$m	…		LL		LH
$[$Fe ii$]$ 17.9 $\mu$m	…		LL		SH
$[$Fe ii$]$ 26.0 $\mu$m	$34\pm 3$	$60\pm 15$	LL	$15\pm 5$	LH
		(+[O iv])
$[$Fe iii$]$ 23.0 $\mu$m	…		LL		LH
H2 S(0) 28.2 $\mu$m	…		LL		LH
H2 S(1) 17.0 $\mu$m	…		LL		SH
H2 S(2) 12.29 $\mu$m	…	…	SL		SH

Full text

*Aims. *Our objective is to identify the dominant thermal mechanisms in one such galaxy, I Zw 18 ($1/30$ Z⊙), assess the diagnostic value of fine-structure cooling lines, and estimate the molecular gas content. Even though molecular gas is an important catalyst and tracer of star formation, constraints on the molecular gas mass remain elusive in the most metal-poor galaxies.

*Methods. *Building on a previous photoionization model describing the giant H ii region of I Zw 18-NW within a multi-sector topology, we provide additional constraints using, in particular, the [C ii] $157$ $\mu$m and [O i] $63$ $\mu$m lines and the dust mass recently measured with the Herschel Space Telescope.

Results. The heating of the H i region appears to be mainly due to photoionization by radiation from a bright X-ray binary source, while the photoelectric effect is negligible. Significant cosmic ray heating is not excluded. Inasmuch as X-ray heating dominates in the H i gas, the infrared fine-structure lines provide an average X-ray luminosity of order $4\times 10^{40}$ erg s-1 over the last few $10^{4}$ years in the galaxy. The upper limits to the [Ne v] lines provide strong constraints on the soft X-ray flux arising from the binary. A negligible mass of H2 is predicted. Nonetheless, up to $\sim 10^{7}$ M⊙ of H2 may be hidden in a few sufficiently dense clouds of order $\lesssim 10$ pc ($\lesssim 0.1$″) in size. Regardless of the presence of significant amounts of H2 gas, [C ii] and [O i] do not trace the so-called “CO-dark gas”, but they trace the almost purely atomic medium. Although the [C ii]+[O i] to total infrared ratio in I Zw 18 is similar to values in more metal-rich sources ($\sim 1$%), it cannot be safely used as a photoelectric heating efficiency proxy. This ratio seems to be kept stable owing to a correlation between the X-ray luminosity and the star formation rate.

*Conclusions. *X-ray heating could be an important process in extremely metal-poor sources. The lack of photoelectric heating due to the low dust-to-gas ratio tends to be compensated for by the larger occurrence and power of X-ray binaries in low-metallicity galaxies. We speculate that X-ray heating may quench star formation.

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$ $\mu$m, [O i] $63$ $\mu$m, and [Si ii] $34$ $\mu$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 $\gtrsim 1/10$ Z⊙ (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({\rm O/H})=7.22$ or $1/30$ Z⊙111We use the oxygen abundance obtained by Péquignot (2008) divided by the solar value from Asplund et al. (2009). , 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 photoionization222Outside 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. 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$ $\mu$m and [Fe ii] $26.0$ $\mu$m). P08 also made tentative predictions for the far-infrared (FIR) lines [C ii] $157$ $\mu$m and [O i] $63$ $\mu$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 $\lambda 1334.5$ and C ii* $\lambda 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 $n_{e}/n_{\rm H}\sim 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; $L_{\rm X}\gtrsim 10^{39}$ 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 $\lambda 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 $\lambda 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 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 Section 2. Observations are described in Section 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 Section 6. The presence of molecular gas and physical conditions in the diffuse gas are investigated in Section 7. Implications of X-ray heating of the H i gas are examined in Section 8. Conclusions are found in Section 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 Section 5.

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 Section 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 $\alpha$ 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., $\approx 7$ kpc (e.g., Kim et al. 1998; Staveley-Smith et al. 2003), and the diameter of the NW H ii region, ($\approx 400$ pc), is about twice that of the LMC-30 Dor region.

2.3 Star formation rate

Legrand et al. (2000) found that a low SFR ($\sim 10^{-4}$ M*⊙* yr*-1*) over the Hubble time could explain the metal enrichment of I Zw 18. The instantaneous SFR derived from H$\alpha$ is $0.1-0.2$ M*⊙* 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$ $\mu$m, providing $0.06$ M⊙ yr*-1*. De Looze et al. (2014) also investigated the applicability of several FIR lines for tracing SFR and obtained $\approx 0.02$ M⊙ 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 $L_{\rm 1.4GHz}=8\times 10^{19}$ W Hz*-1* (Hunt et al. 2005) and the SFR calibration from Bell (2003), we obtain $0.13$ M⊙ yr*-1*.

Annibali et al. (2013) examined color-magnitude diagrams (CMD) from deep Hubble/ACS images and found a larger value of $\approx 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 ($\sim 100$ Myr for the FUV vs. $\sim 10$ Myr for H$\alpha$ 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 Section 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\pm 0.16$ dex ($2\sigma$ 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 $\approx 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$ $\mu$m, [O i] $63$ $\mu$m, and [O iii] $88$ $\mu$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 Section 5.

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$ $\mu$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$ $\mu$m (OBSID 1342253757) and [O iii] $88$ $\mu$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\times 5$ spatial pixels (“spaxels”), corresponding to a $\approx 47\arcsec\times 47\arcsec$ field of view. Each spaxel is $\approx 9.4\arcsec$ in size. A single footprint observation was performed since I Zw 18 appears smaller than the footprint size. According to the PACS Observer’s Manual333http://herschel.esac.esa.int/Docs/PACS/html/pacs_om.html, the point spread function (PSF) full width at half maximum (FWHM) ranges from $\approx 9.5\arcsec$ ($\approx 0.8$ kpc at the adopted distance of $18.2$ Mpc) between $55$ $\mu$m and $100$ $\mu$m to about $\approx 14\arcsec$ at $200$ $\mu$m ($\approx 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\sigma$) 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 $\lesssim 3\sigma$ (Fig. 2), so it is difficult to pinpoint the peak position or the source spatial shape with an accuracy smaller than the spaxel size ($\approx 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 $\lesssim 6\arcsec$ ($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 $F_{1}$ scales the flux in the brightest spaxel by applying a point-source correction (from $25\%$ to $67$% between $50-220$ $\mu$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 $F_{1}$ 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\times 3$ spaxels centered on the brightest spaxel and adds the line fluxes either from all spaxels ($F_{3\times 3}$) or only from spaxels with $>3\sigma$ detections ($F^{\prime}_{3\times 3}$). A point-source correction factor is also required, although it is much smaller (from $4$% to $17$% for the range $50-220$ $\mu$m) than for $F_{1}$. The $3\times 3$ methods are more reliable than $F_{1}$ 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 $F^{\prime}_{3\times 3}$ method uses an incomplete sampling of the PSF and therefore results in a lower limit on the flux determination444Alternative methods consist of combining the spaxel spectra (either $3\times 3$ or $5\times 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. .

Finally, method $F_{\rm opt}$ 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 $F_{1}$ to be a lower limit because the emission can never be perfectly centered in any spaxel and we use $F_{\rm opt}$ for our final fluxes. Cormier et al. (2015) used the $F_{3\times 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 $F_{3\times 3}$ encompass the $F_{\rm opt}$ determination. We discuss, in Section 5.1, how PACS fluxes are normalized for comparison with the other tracers used in this study.

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$ $\mu$m and [Ne iii] $15.5$ $\mu$m), but also [Si ii] $34.8$ $\mu$m and [Fe ii], which partly originate from the neutral gas. For the low-resolution spectrum (SL and LL modules; $R=\lambda/{\Delta\lambda}\sim 57-126$ over $\approx 5-14$ $\mu$m and $\approx 14-36.5$ $\mu$m, respectively), the deepest observation available was used (AORkey $16205568$). For the high-resolution spectrum (SH and LH; $R\sim 600$ over $\approx 10-19.5$ $\mu$m and $\approx 20-36.5$ $\mu$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 $\approx 2\arcsec$ at $5$ $\mu$m to $\approx 11\arcsec$ at $38$ $\mu$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\arcsec$ FWHM). One can distinguish two components in the SL profile. One component is located at 09h34m02.29s/+55∘14$\arcmin$27.52$\arcsec$, coinciding with NW, with a prominent [S iv] line and relatively shallow continuum, while the other component is located at 09h34m02.31s/55∘14$\arcmin$22.67$\arcsec$, coinciding with SE, with a much weaker [S iv] line and a relatively steeper continuum. The NW component is responsible for $\approx 75$% of the total [S iv] and $\approx 68$% of the total H i recombination line Hu$\alpha$ 12.37 $\mu$m (see fluxes in Appendix B). These values are in good agreement with the H$\beta$ fraction originating from NW, $\approx 78\%$ (Skillman & Kennicutt 1993). The NW component appears extended in the SL module while the SE component appears quasi-point-like (intrinsic broadening of $\approx 5\arcsec$ and $\approx 1\arcsec$, respectively). The NW and SE components are not distinguishable in the LL and LH modules because of the relatively lower spatial resolution.

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 $\mu$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$ $\mu$m line arising in the same IRS module as [S iii] $33.5$ $\mu$m, P08 proposed that the fluxes of these two lines were overestimated. We find new [S iii] $33.5$ $\mu$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$ $\mu$m and $33.5$ $\mu$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$ $\mu$m is $\leq 1.9\times 10^{-17}$ W m*-2* with little influence from the neutral/ionized PAH mixture. We review in Section 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 $L_{\rm X}=6.2_{-4.9}^{+2.3}\times 10^{39}$ erg s*-1* between $0.1-2.4$ keV. Martin (1996) independently examined the same observation, calculating a lower limit of $L_{\rm X}\gtrsim 1\times 10^{39}$ 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 region666Ott et al. (2005) reported a different position, located at the edge of the H$\alpha$ 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. (see position in Fig. 1). The point-source luminosity is $3\times 10^{39}$ 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\times 10^{40}$ 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 $10^{38}$ 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 Section 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 i $21$ 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 Section 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 $\approx 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.

Based on H$\alpha$ morphology, and following P08, the inner and outer radii of the H+ shell are set to $1.5\arcsec$ and $2.5\arcsec$, respectively ($\approx 130$ pc and $\approx 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 $\sim 6\arcsec$ (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 ($M2_{\rm X}$) needs to be updated777The 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. 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.03888The 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 ${\rm O}^{0}(^{3}P)+{\rm H}^{0}(^{2}S)$. The more recent rate is many times larger than the old rate. (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.

Owing to the larger distance, the radial density profile, $n_{\rm H}(r)$, must be updated. After convergence, the new parameters in expression (1) of P08 are as follows: $P_{\rm in}=2.8$, $P_{\rm out}=24$ (in $/k/10^{5}$ CGS), and $\tau_{\rm c}=4.2$. The $n_{\rm H}(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\times 10^{4}$ K and $8\times 10^{4}$ 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 Figure 7.

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 fitted999Fitting He ii and a number of other nebular lines is a prerequisite to any model of the H ii region (see P08). . 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\times 10^{6}$ 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.

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 manner111111We 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. . The Herschel and Spitzer line fluxes (Sects. 3.1.2, 3.2) first need to be normalized to be compared to the optical tracers. All fluxes are scaled to H$\beta=1000$.

For the normalization of Herschel/PACS measurements, we calculate the total deredenned H$\beta$ flux in the PACS footprint from the H$\alpha$ map (de Paz et al. 2003), assuming H$\alpha$/H$\beta=2.8$ and $E$(B-V)$=0.09$ (Schlegel et al. 1998; Schlafly & Finkbeiner 2011). We find $F({\rm H}\beta)\approx 11.5\times 10^{-17}$ W m*-2*, implying a scaling factor of $\approx 1.15\times 10^{-19}$, which is the factor by which the PACS fluxes should be divided by to normalize to H$\beta=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$ $\mu$m. We find $\approx 1.2\times 10^{-19}$, i.e., in good agreement with the geometrical factor derived directly from the H$\alpha$ 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\times 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$\alpha$ 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 IRS121212The 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$\beta$ emission, but that arising from the main NW shell. However, P08 calibrated the different instruments, using H$\alpha$ 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.. 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).

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$ $\mu$m/$63$ $\mu$m ratio found in other DGS sources (Cormier et al. 2015).

The [C ii] $157$ $\mu$m and [O i] $63$ $\mu$m lines are the dominant coolants in the neutral ISM as long as the metallicity is above some critical value ($\gtrsim 10^{-3.5}$ Z⊙; e.g., Santoro & Shull 2006). The [Si ii] $34.8$ $\mu$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 $10^{3}$ years. More importantly, however, the models indicate that typical cooling timescales in the temperature range $100-300$ K are $\sim 10^{4}$ years and $\sim 10^{5}$ years 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 $10^{4}$ years, 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\times 10^{40}$ 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., $724^{+535}_{-308}$ M⊙. To put this value in perspective, the giant H ii region LMC-N 11, despite a smaller physical size ($\approx 150-200$ pc) as compared to the I Zw 18 NW region, harbors $\sim 3\times 10^{4}$ M⊙ of dust (Gordon et al. 2014). The total gas mass in the main body is $10^{8}$ M⊙ (Lelli et al. 2012), leading to D/G$\approx 1/1000$ D/GMW, where D/GMW is the Milky Way value141414D/GMW is $1/148$, close to the values in Zubko et al. (2004), $1/158$, and Jones et al. (2013), $1/156$. .

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 $\approx 10^{8}$ M⊙, 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 $10^{8}$ M⊙.

We use the $2\arcsec$ resolution H i map from Lelli et al. (2012) to estimate the H0 mass associated with either the NW region itself (defined as a $\approx 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\times 10^{21}$ cm*-2*, which translates into a mass $5\times 10^{6}$ M⊙. 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, $\approx 2\times 10^{7}$ M⊙. In the following we consider that H0 masses around $\approx 0.5-2\times 10^{7}$ M⊙ 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 to the NW region ($\lesssim 6\arcsec$; Sect. 3.1.2). The centroid for [O i] and [O iii] seems to be located closer to NW, possibly coinciding with a dust-rich ionized gas shell located between NW and SE (labeled NW-D3 in Cannon et al. 2002), near the H 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$\beta$, Hu$\alpha$, 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$\beta$, this means that line constraints should not be significantly affected by SE.

• •

All H i regions will be subject to the X-ray radiation, whether they are associated with the NW or SE giant H ii regions. Since the ionization state is stable in the H i gas, the fine-structure line emissivity depends essentially on local temperature, which is controlled by the local flux of X-rays. Geometrical details have little 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$\beta$ 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$\beta$ 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 Section 6.3.2.

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 model are expected, by design, to reproduce both lines. This is not the case, however, for models with no X-rays ($\mathcal{M}0a$ and $\mathcal{M}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 ($\mathcal{M}1$), with dust but no CR ($\mathcal{M}2$), and with CR but no dust ($\mathcal{M}3$). The objective is not to find the best combination but to test the impact of each parameter. The full model ($\mathcal{M}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 $\mathcal{M}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 Section 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 $L_{\rm X}$ 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.

The model that reproduces the observations (line emission and H0 mass) best is labeled $\mathcal{M}1$. As expected, the [C ii], [O i] $63$ $\mu$m, and [Si ii] transitions are the strongest coolants, followed by [O i] $145$ $\mu$m and [Fe ii], in the H i region. The required X-ray luminosity is $4\times 10^{40}$ erg s*-1*, which is three times larger than the observed luminosity. According to the analysis of Section 5.3, such a difference is acceptable.

In summary, model $\mathcal{M}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 $10^{4}$ years in I Zw 18. We now explore models involving alternative heating mechanisms.

6.3 Exploratory models $\mathcal{M}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.

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 ($\mathcal{M}2a$, $\mathcal{M}2b$, $\mathcal{M}2c$)

First, a uniform D/G value is assumed across the galaxy. Model $\mathcal{M}2a$ uses a global D/G of $1/1000$ D/GMW, corresponding to the observed value as constrained by Herschel (Sect. 5.4). Model $\mathcal{M}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 $\mathcal{M}2a$ is almost the same as that of model with only X-ray heating ($\mathcal{M}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 $\approx 4$% of the total heating in the H i region in model $\mathcal{M}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 $\mathcal{M}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 $\mathcal{M}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 $\mathcal{M}2a$.

The photometry data at $100$ $\mu$m and the upper limit at $160$ $\mu$m provide an upper limit on the 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 $\mathcal{M}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$ $\mu$m (Fig. 11), which could be due to a missing dust component (not contributing significantly to the dust mass). Model $\mathcal{M}2c$ has similar physical conditions in the H i region (density, temperature) as model $\mathcal{M}2a$ and the heating fraction contributed for by the photoelectric effect is only $10$%.

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 Section 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 ($\mathcal{M}2d$, $\mathcal{M}2e$)

In Section 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 ($\mathcal{M}2d$) agrees better with the observations than models with uniform D/G (Fig. 11). The photoelectric effect in model $\mathcal{M}2d$ represents only $\approx 12$% of the H i heating (Table 7). Even then, model $\mathcal{M}2d$ still somewhat overestimates the cold dust emission and therefore the total dust mass.

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 $\mathcal{M}2a$ ($1/1000$ D/GMW) to include an extra sector with $1/50$ D/GMW and with a relatively small covering factor ($\mathcal{M}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$\alpha$/H$\beta$ 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 SE161616This 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. 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$\beta$ 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 $\mathcal{M}2e$, the new sector has a uniform gas density of $20$ cm*-3*, extends from $50$ pc to $\approx 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 $\approx 17$ M⊙. The total dust mass is $\approx 200$ M⊙, i.e., somewhat lower than the determination by Rémy-Ruyer et al. (2015). The global D/G (all sectors combined) is $\approx 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 $\lesssim 12\%$.

6.4 Exploratory models $\mathcal{M}3$: Cosmic ray heating

We now consider CR ionization as an additional parameter to the initial model $\mathcal{M}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$\alpha$ 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$\alpha$, FUV+$24$ $\mu$m, or $1.4$ GHz) agree with $0.05-0.1$ M⊙ yr*-1*, but a CMD analysis indicates a larger value of $\sim 1$ M⊙ 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\times 10^{-16}$ s*-1*; e.g., Indriolo et al. 2015) and let Cloudy calculate the heating rate self-consistently (model $\mathcal{M}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 $\approx 13$%.

In model $\mathcal{M}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 $\approx 65$% and the required X-ray luminosity decreases to $\approx 1.25\times 10^{40}$ 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 ($\mathcal{M}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 ($\mathcal{M}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 $\approx 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 $\mathcal{M}4a$

The full model $\mathcal{M}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.

The optical and infrared ionized gas line intensities predicted by model $\mathcal{M}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$ $\mu$m predicted by model $\mathcal{M}4a$ is only a factor of $\approx 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 $\mathcal{M}4a$. The electron fraction ranges from $\sim 0.5$% down to $\sim 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$ $\mu$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 $\sim 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 $\epsilon^{\prime}_{\rm TIR}=$([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).

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 $\mathcal{M}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 $\mathcal{M}4b$ (Table 7), molecular gas is uniformly distributed in the relatively diffuse gas ($\lesssim 200$ cm*-3*) of the entire radiation-bounded sector (i.e., using the same sector configuration as model $\mathcal{M}4a$; Fig. 6). The molecular gas that could be hidden in small clumps is discussed separately in Section 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 (${\rm H}^{-}+{\rm H}\rightarrow{\rm H}_{2}+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:

• $\bullet$

The electron fraction in the H i region is significant, providing free electrons for producing H- through the reaction ${\rm H}+e^{-}\rightarrow{\rm H}^{-}+\gamma$ (e.g., Jenkins & Peimbert 1997).

• $\bullet$

The medium density is $\gtrsim 100$  cm*-3*, below which the H- photodetachment (${\rm H}^{-}+\gamma\rightarrow{\rm 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 $\sim 10^{4}$  cm*-3* (for a metallicity adapted to I Zw 18; Omukai et al. 2005).

• $\bullet$

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 ionization. The number of secondary ionizations is $\approx 25$ for a $1$ keV photon181818Number of secondary ionizations is $25\ (E\gamma/1\,{\rm keV})$ (Shull & van Steenberg 1985). 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 $\mathcal{M}4b$ predicts an H2 column density of $4\times 10^{13}$ cm*-2* that is lower than the observed upper limit $\lesssim 10^{15}$ cm*-2* measured by Vidal-Madjar et al. (2000) with FUV absorption lines. The authors also calculated a theoretical column density of $2\times 10^{12}$ 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 $\mathcal{M}4b$, the H- formation route contributes to $\sim 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 Section 7.1.2) H2 column density may be close to a few $10^{13}$ 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({\rm H}_{2})\approx 2.5\times 10^{14}$ 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 $\mathcal{M}4b$ is only $\sim 0.1$ M⊙ (corresponding to a molecular gas mass fraction of $\sim 5\times 10^{-9}$). In the sector configuration of $\mathcal{M}4b$, the [C ii] and [O i] emission lines thus trace an almost purely atomic medium.

The CO(1-0) intensity predicted by model $\mathcal{M}4b$ is $2\times 10^{-30}$ W m*-2*, to be compared to the upper limit by Leroy et al. (2007) of $8\times 10^{-22}$ W m*-2* (with a beam of $\approx 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 $\mathcal{M}4c$)

As mentioned in Section 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 $\mathcal{M}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 0335$-$052 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({\rm 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.

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$ $\mu$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\times 10^{6}$ 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 $\sim 10^{6}$ 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 $10^{25}$ cm*-2* ($A_{V}\sim 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 $\approx 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 $14^{3}$ pc3 ($\sim 25\times 25$ pc2 surface area and $\sim 5$ pc deep). Considering the depth scale and assuming spherical clumps, this should correspond to a couple dozen clumps at most.

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 Figure 15, most of the [C ii] emission in the clump sector occurs at depths for which the gas is molecular ($\approx 75$% of [C ii] is associated with a molecular gas fraction $>50$%). However, [C ii] does not trace well the H2 mass, 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 H2 $2.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 $\mathcal{M}4c$ is $\sim 4\times 10^{7}$ M⊙, 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\times 10^{5}$ M⊙191919Value rescaled at $18.2$ Mpc. with a Galactic CO/H2 conversion factor.. Assuming the star formation efficiency is the same as in the Milky Way, we can argue that $X_{\rm CO}$ is significantly larger. If $X_{\rm CO}$ 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 $\mathcal{M}4c$. A large $X_{\rm CO}$ 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 $\lesssim 5$ pc ($\lesssim 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$ $\mu$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$ $\mu$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 $\mu$m emission line or the C ii* $\lambda 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 $\mu$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 $\lambda 1334.5$ and C ii* $\lambda 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\arcsec$ aperture. The authors inferred an electron fraction upper limit of $2$% in the H i region with a best guess of $\sim 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.

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 $\approx 2.5\arcsec$ 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 $\mu$m line intensity, $10^{-17}$ W m*-2*, we find a total rate of $3.9\times 10^{31}$ W. For the measure in absorption, we first use the cooling rate per H atom obtained by Lebouteiller et al. (2013), $8.3\times 10^{-28}$ erg s*-1* (H*-1*), which we then multiply by the total mass of H0 in the main body, $\approx 10^{8}$ M⊙, assuming that the cooling rate per H atom is uniform across the H i body. We find a rate of $10^{38}$ erg s${}^{-1}=10^{31}$ W, i.e., a factor of $\approx 4$ difference with the rate measured in emission. Another way to show this discrepancy is to compare the COS measurement $8.3\times 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 $\approx 10^{-26}-10^{-25}$ erg s*-1* (H*-1*), implying a factor of $\gtrsim 10$ larger than the COS rate.

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 $\lesssim 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 $\mu$m emission line probes the neutral shell of the radiation-bounded H ii region with an H i column density of $10^{22.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 $10^{21.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 $\approx 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 $\sim 1$ kpc away from the NW cluster. Such a model is in fact also motivated by H$\alpha$ 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$\beta$ 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$\alpha$ 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 $L_{\rm X}>10^{39}$ 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 $\gtrsim 8$ M⊙, and with evolution timescales on the order of $10^{6-7}$ yr. A correlation between SFR and $L_{\rm X}$ 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-$L_{\rm X}$ relation is linear for high-SFR regime ($>4.5$ M⊙ 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 $L_{\rm X}$ in I Zw 18 are consistent with such a metallicity dependence202020I Zw 18 would fall closer to the linear relationship with SFR of $\sim 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-$L_{\rm X}$ 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). . 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-$L_{\rm X}$ 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 ($\epsilon_{\rm TIR}$; with TIR calculated between $3-1100$ $\mu$m) ratio is often used to estimate the photoelectric effect heating efficiency212121$\epsilon_{\rm 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., 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 $\epsilon_{\rm TIR}$ (Cormier et al. 2015). A larger $\epsilon_{\rm 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 $\epsilon_{\rm TIR}$.

For I Zw 18 we find $\epsilon_{\rm 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 $\epsilon_{\rm TIR}=1.2$%. For comparison, Lebouteiller et al. (2012, 2013) found $\epsilon_{\rm 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, $\epsilon_{\rm TIR}$ in I Zw 18 is similar to the values in other galaxies, although somewhat larger. As shown in Section 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 $\epsilon_{\rm TIR}$ is not extreme in I Zw 18, we must take a closer look at the tracers involved in the quantity $\epsilon_{\rm 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 TIR${}_{\rm ff}\approx 6.9\times 10^{6}$ L⊙, i.e., only about four times lower than the observed TIR luminosity $2.9\times 10^{7}$ L⊙ (Rémy-Ruyer et al. 2015). The observed luminosity of [C ii]+[O i] calculated at the same distance is $1.8\times 10^{5}$ L⊙, which results in $\epsilon_{\rm TIR,ff}\sim 2.6$%. Hence even if I Zw 18 were virtually dust free, $\epsilon_{\rm 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 $\epsilon_{\rm 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 $\epsilon_{\rm TIR}=0$. Therefore, the fact that $\epsilon_{\rm 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 $L_{\rm X}$ (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 $\sim 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/30$ Z⊙) 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$ $\mu$m, [O i] $63$ $\mu$m, and [O iii] $88$ $\mu$m) and updated Spitzer/IRS measurements (particularly [Si ii] $34.8$ $\mu$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:

• $\bullet$

The photoelectric effect is not a significant heating mechanism in the H i region, whether we use the D/G constrained by Herschel or a larger D/G scaled with metallicity.

• $\bullet$

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.

• $\bullet$

More precisely, the H i zone fluxes imply a luminosity $L_{\rm X}\approx 4\times 10^{40}$ 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 $10^{4}$ years. This large average X-ray luminosity could imply the existence of a massive accreting black hole in I Zw 18-NW.

• $\bullet$

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.

• $\bullet$

The models predict only about $0.1$ M⊙ 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.

• $\bullet$

We also explore the possible existence of undetected molecular clumps. These H2 clouds should reach a density of $10^{6}$ cm*-3* and have a small covering factor ($\lesssim 0.1$%). The main constraints come from the CO emission and the dust SED. This dense sector can shelter $\lesssim 4\times 10^{7}$ M⊙ of H2, in which case H2 must belong to a few clouds of a few parsecs in size ($\lesssim 0.1^{\prime\prime}$). 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 ($\gtrsim 85$%) of the [C ii] and [O i] emission in the model is associated with atomic gas in the more diffuse sector.

• $\bullet$

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 ($\sim 10$ cm*-3*) up to several kiloparsecs away from the stellar cluster.

• $\bullet$

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.

• $\bullet$

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

Acknowledgements.

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

Appendix A Herschel/PACS spectral extraction

A.1 PACS optimal extraction of pointed observations

The Herschel/PACS spectra are shown in Figures 18 and 19. We derive the following radial velocities, $795\pm 100$ km s*-1* for [C ii], $760\pm 40$ km s*-1* for [O i], and $795\pm 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 $\approx 765$ km s*-1*. The velocity ranges from $\approx 735$ km s*-1* to $\approx 790$ km s*-1* from north to south. Interferometric H$\alpha$ measurements are in agreement with these values (Petrosian et al. 1997).

The optimal extraction uses the “as-built” telescope monochromatic PSF222222http://pacs.ster.kuleuven.ac.be/pubtool/PSF/ 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 Figure 20. 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 $\lesssim 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 $\approx 10$%, $\approx 18$%, and $\approx 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$ $\mu$m and $100$ $\mu$m bands (with a similar PSF size as the PACS spectroscopy observations) (Rémy-Ruyer et al. 2015).

A.2 Unchopped [O iii] 88 $\mu$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$ $\mu$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 Figure 21. The emission corresponds to a compact source located at 09h34m02.12s/55∘14$\arcmin$27.4$\arcsec$, i.e., coinciding with the NW cluster (Fig. 21). 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\pm 0.25\times 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 Section 3.1.2. Only $F_{\rm opt}$ (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 9). 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 manual232323http://herschel.esac.esa.int/Docs/PACS/html/pacs_om.html). In our analysis we only use the chop/nod [O iii] observations described in Section 3.1.1.

Appendix B Spitzer/IRS line flux measurements

The apertures of the various Spitzer/IRS observations252525A small map was performed in SH (AORkey 21825792) with dedicated offset observations, but we do not consider it because of its low S/N. are overlaid in Figure 22. It can be seen that the LL ($10.7\arcsec\times 168\arcsec$) and LH ($11.1\arcsec\times 22.3\arcsec$) 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\arcsec\times 57\arcsec$) and SH ($4.7\arcsec\times 11.3\arcsec$) 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.

We use the extracted spectra from Cornell 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 10) 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.

Figure 23 shows some line fits for illustration. We list in Tables 10.1 and 11.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 12.3. along with a comparison to W07.

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 $\sim 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 $A_{V}$-to-column density conversion of Guver & Ozel (2009). We find $N_{\rm H}\approx 2\times 10^{20}$ 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\times 10^{20}$ cm*-2* (Willingale et al. 2013). Figure 24 shows the diskbb model and the absorbed spectrum. Most of the absorption occurs in I Zw 18, with a column density of $10^{21}$ cm*-2*. The luminosity between $0.3-1$ keV is $\approx 10^{40}$ 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\times 10^{21}$ 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\arcsec$ COS aperture), which may differ from that derived along the line of sight of the X-ray point source.

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\times 10^{6}$ K blackbody used in P08 neglected the hard X-ray tail ($\gtrsim 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\times 10^{22}$ cm*-2* proved to provide sufficient flexibility (nine free parameters). Figure 25a shows a broad selection of combinations providing the best $\chi^{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\times 10^{21}$ cm*-2* (Fig. 25b). Two different spectral shapes emerge, one with a soft X-ray excess and one without. In Figure 25b, the X-ray luminosity range of the models is $1.1-2.5\times 10^{40}$ erg s*-1*. Below $\sim 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$ $\mu$m and [Ne v] $24.3$ $\mu$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 $\gtrsim 100$ eV. In Figure 25b, 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 25c shows the combination spectra allowed when the additional constraint brought by the upper limit on [Ne v] is applied to photoionization model results292929Since 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\times 10^{40}$ erg s*-1* used in our modeling (Sect. 5.3)..

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 $\sim 0.4$ keV well defined, but it turns out that, within a factor of $\sim 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.