EMPRESS. VII. Ionizing Spectrum Shapes of Extremely Metal-poor Galaxies: Uncovering the Origins of Strong He ii and the Impact on Cosmic Reionization

Hiroya Umeda; Masami Ouchi; Kimihiko Nakajima; Yuki Isobe; Shohei Aoyama; Yuichi Harikane; Yoshiaki Ono; Akinori Matsumoto

doi:10.3847/1538-4357/ac602d

1. Introduction

Understanding the epoch of reionization (EoR) is one of the most important goals for astronomy today. Many studies have suggested that galaxies are the dominant sources of ionizing photons at the EoR (e.g., Wise et al. 2014; Madau & Haardt 2015; Stanway et al. 2016). From hydrodynamical simulations, high-redshift young galaxies at the early formation phase are expected to have low metallicities, low stellar masses, and high specific star formation rates (e.g., Wise et al. 2012). High-redshift galaxies are studied with various large telescopes including the Hubble Space Telescope (e.g., Bouwens et al. 2015; Harikane et al. 2016; Livermore et al. 2017; Atek et al. 2018; Oesch et al. 2018; Kikuchihara et al. 2020). Although there are spectroscopic studies for very bright or lensed galaxies at high redshifts (e.g., Stark et al. 2015; Mainali et al. 2018; Shibuya et al. 2018), dwarf galaxies at high redshift are too faint to be detected with current observational instruments. For example, even with the upcoming James Webb Space Telescope, we can detect Hα for galaxies with stellar masses below 10⁶ M_⊙ only at z ≲ 2 without gravitational lensing (Isobe et al. 2022).

Owing to observational difficulties, many studies use local dwarf galaxies as analogs of high-redshift young dwarf galaxies (e.g., Berg et al. 2019). Among the local dwarf galaxies, local extremely metal-poor galaxies (EMPGs) are recently gaining attention (e.g., Izotov & Thuan 1998; Thuan & Izotov 2005; Sánchez Almeida et al. 2016; Kojima et al. 2020). EMPGs are defined as galaxies with metallicities below 10% of the solar metallicity Z_⊙. So far, EMPGs with metallicities down to 1.6% Z_⊙ have been identified using data obtained with the Subaru/Hyper Suprime Cam (HSC; Kojima et al. 2020). Local EMPGs have low masses and high specific star formation rates (sSFRs), both of which are comparable to those of high-redshift young galaxies (e.g., Wise et al. 2012). Local EMPGs are therefore regarded as local analogs of high-redshift young galaxies.

The importance of relationships between local EMPGs and high-redshift galaxies is further emphasized because very strong high-ionization nebular emission lines (i.e., He ii, C iii], C iv) are detected from local EMPGs (Guseva et al. 2000; Shirazi & Brinchmann 2012; Senchyna et al. 2017; Berg et al. 2019; Kojima et al. 2021) and high-redshift galaxies (Stark et al. 2015; Laporte et al. 2017; Schmidt et al. 2017; Mainali et al. 2018). Most of the EMPGs have detections of optical recombination lines He ii λ4686 that need high-energy ionizing photons (>54.4 eV; López-Sánchez & Esteban 2010; Shirazi & Brinchmann 2012; Senchyna et al. 2017). A statistical study of local galaxies shows that fluxes of He ii λ4686 increase as the metallicity decreases (e.g., Schaerer et al. 2019). These results indicate the presence of hard ionizing sources in metal-poor environments.

The origin of strong He ii emission lines is still under debate. Xiao et al. (2018) examined contributions from stellar radiation using the photoionization code cloudy (Ferland et al. 2013) and the Binary Population and Spectral Synthesis (BPASS) code (Stanway et al. 2016; Eldridge et al. 2017). BPASS binary models incorporate contributions from hard radiation from Wolf–Rayet stars and stripped binary stars. However, a line ratio He ii λ4686/Hβ predicted in Xiao et al. (2018) falls short of the observed value by more than 0.5 dex, suggesting that stars are not major contributors to the He ii ionizing photons. Addition of extra stars (e.g., metal-poor, very massive, or Wolf–Rayet stars) could boost the He ii ionizing photons, but several studies suggest that these stars cannot fully account for the observed strength of He ii in EMPGs (e.g., Kehrig et al. 2018; Olivier et al. 2021). Pérez-Montero et al. (2020) suggest that ordinary stellar synthesis models with a nonzero ionizing photon escape fraction can explain He ii λ4686/Hβ, but their models require very high escape fraction values (∼70%) that are not compatible with the typical values reported for z ≲ 3 galaxies (e.g., Rutkowski et al. 2016; Steidel et al. 2018; Alavi et al. 2020; Izotov et al. 2021b; Saxena et al. 2022).

Nonthermal ionizing sources have been proposed to explain the strong nebular He ii from EMPGs. Active galactic nuclei (AGNs) are one of the candidates of nonthermal emission sources (e.g., Berg et al. 2018; Plat et al. 2019). Berg et al. (2018) found that AGN models can reproduce observed He ii λ1640, while other emission lines, such as C iv λλ1548,1551 and C iii]λ1909, cannot be reproduced simultaneously. However, their AGN models do not consider metal-poor (<0.25 Z_⊙) AGNs. Shocks are another candidate for the ionizing sources of He ii (e.g., Thuan et al. 2005; Berg et al. 2018; Plat et al. 2019; Isobe et al. 2021). Berg et al. (2018) also examined contributions from radiative shocks using radiative shock models of Allen et al. (2008). Calculations with radiative shock models also suggest that radiative shocks cannot explain the strong He ii λ1640 emission for a given intensity of other emission lines such as C iv λλ1548,1551 and C iii]λ1909 (see also Izotov et al. 2021a).

Alternatively, other hard radiation sources are also proposed. Schaerer et al. (2019) show that X-ray luminosities have a positive correlation with line ratios He ii λ4686/Hβ with X-ray binary population models, indicating the contributions of the He ii ionizing photons from high-mass X-ray binaries (HMXBs). Lebouteiller et al. (2017) demonstrated that HMXBs can be the sources of He ii ionizing photons in well-studied EMPG I Zw 18 through detailed modeling of radiation fields and interstellar mediums. Kehrig et al. (2021) argue that the ultraluminous X-ray source (ULX) in I Zw 18 cannot fully account for all of He ii ionizing photons. However, new integral field spectroscopic observations of I Zw 18 suggest that the ULX still remains the possibility (see Rickards Vaught et al. 2021). Moreover, while Senchyna et al. (2020) argue that HMXBs cannot be the main contributors of He ii in EMPGs through the photoionization models representing HXMB spectra by simple multicolor disk models, Simmonds et al. (2021) examined the photoionization models using several detailed spectrum models for HMXBs/ULXs and found that some of the models can reproduce He ii λ4686/Hβ and other line ratios, such as [O iii]λλ4959,5007/Hβ and [O iii]λλ4959,5007/[O ii]λ3727. The studies of Schaerer et al. (2019), Simmonds et al. (2021), and Lebouteiller et al. (2017) claim that HMXB/ULX can explain strong He ii emission lines (see Saxena et al. 2020 and Senchyna et al. 2020), but whether a simple physical model can explain all the emission lines from a galaxy in a self-consistent manner remains uncertain. There are also attempts to explain emission lines from EMPGs by the nonuniform interstellar medium (ISM; e.g., Ramambason et al. 2020; Berg et al. 2021; Izotov et al. 2021a). For example, Ramambason et al. (2020) showed that ordinary stellar population synthesis models can reproduce high-ionization emission lines [O iii]λλ4959,5007 simultaneously with low-ionization emission lines [O ii]λ3727 by introducing the so-called "two-zone" models.

Variations of models presented by these previous studies raise difficulties in reproducing strong He ii emission line fluxes for given low-ionization line fluxes. Instead, we aim to determine the self-consistent general spectrum shape by a parametric method. By obtaining the general spectrum shapes first, we can infer the possibility of countless ionization sources without running costly calculations using photoionization code each time. Moreover, we expect to uncover the spectral features that were not previously mentioned.

In this work, we search for the self-consistent spectrum shapes that can reproduce observed emission line fluxes of three EMPGs. Similar work such as that by Olivier et al. (2021; hereafter, O21) uses blackbody and BPASS models to search for the spectrum shape. Another similar study uses combinations of BPASS models and ULX models as ionizing spectra (Simmonds et al. 2021). With the minimum assumptions, we construct generalized spectral models that cover spectrum shapes examined in O21 and Simmonds et al. (2021). We use a generalized spectrum as an ionizing spectrum in a uniform ISM to compute emission line fluxes using the photoionization code cloudy. We then search for the self-consistent spectrum shapes using the Markov Chain Monte Carlo (MCMC) method.

Our sample and methodology are described in Section 2. We present our results for the general spectrum shapes and reproduction of emission lines in Section 3. In Section 4, we evaluate several candidates for the hard radiation sources using the general spectrum shape. We also examine the contributions of the hard radiation sources to cosmic reionization. In Section 5, we summarize our main results. Throughout this paper, magnitudes are in the AB system, and we assume a standard Λ CDM cosmology with parameters of (Ω_m, Ω_Λ, H₀) = (0.3, 0.7, 70 km s⁻¹ Mpc⁻¹) and the solar metallicity scale of Asplund et al. (2009), where 12 + $\mathrm{log}({\rm{O}}/{\rm{H}})$ = 8.69.

2. Sample and Methods

2.1. Sample

We investigate general ionizing spectrum shapes for a sample of three galaxies, J1631+4426, J104457, and I Zw 18 NW, whose properties are summarized in Table 1. All three galaxies are extremely metal-poor (Z < 10% Z_⊙), nearby (z ∼ 0), and low-mass (≲10⁶ M_⊙) galaxies. J1631+4426 is the most metal-poor galaxy detected, with Z = 1.6% Z_⊙ (Kojima et al. 2021; hereafter, K21). Thus, J1631+4426 is a good test-bed for investigating the nature of EMPGs. J104457 is a metal-poor galaxy, whose optical spectrum has been studied previously, and it has been shown that BPASS models alone cannot explain the high-energy ionization lines (Berg et al. 2021; hereafter, B21). Olivier et al. (2021) have also studied J104457 and have shown that the combination of a blackbody and BPASS model can explain the multilevel emission lines of J104457. Therefore, we can compare our results to those of Olivier et al. (2021). I Zw 18 NW is one of the most well-known EMPGs that has optical and X-ray data (Thuan & Izotov 2005; hereafter, TI05). This allows us to check the consistency between the ionizing spectrum shapes beyond 54.4 eV and the X-ray luminosity data for I Zw 18 NW. All three galaxies in our sample have optical spectroscopic data available. In Figure 1, we present the locations of all three galaxies on a diagram of [N ii]λ6583/Hα and [O iii]λ5007/Hβ emission line ratios, the so-called the Baldwin–Philips–Terlevich diagram (BPT diagram; Baldwin et al. 1981). The black dots represent the emission line ratios of z ∼ 0 star-forming galaxies (SFGs) and AGNs from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7; Abazajian et al. 2009). [N ii]λ6583/Hαand [O iii]λ5007/Hβ emission line ratios for all three galaxies can be explained by only stellar radiation (Kauffmann et al. 2003). However, because Kewley et al. (2013) suggest that metal-poor gas heated by AGN can produce emission line ratio values that fall on the SFG region, we cannot conclude the contributions from AGNs from this classification.

**Figure 1.** Our three galaxies on the BPT diagram. The red, blue, and orange dots represent the emission line ratios of J1631+4426, J104457, and I Zw 18 NW, respectively. The bars represent the errors of the emission line ratios. The red arrow indicates an upper limit with 1σ level for [N ii]λ6583/Hα of J1631+4426. The black dots represent z ∼ 0 SFGs and AGNs derived from the emission line catalog of SDSS DR7. The solid curve depicts the maximum photoionization models that can be reproduced solely by stellar radiation (Kauffmann et al. 2003). The region left (right) to the solid curve is the SFG (AGN and composite) region.
Download figure:
Standard image High-resolution image

Table 1. Sample Properties

Name of the Galaxy	z	i	log(M_⋆)	log(SFR)	Z	F(Hβ)	EW₀(Hβ)	L_{X,0.5−10 keV}	References
		(mag)	(M_⊙)	(M_⊙ yr⁻¹)	(Z_⊙)	(10⁻¹⁴ erg s⁻¹ cm⁻²)	(Å)	(10³⁹ erg s⁻¹)
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
J1631+4426	0.031	22.5	5.9	−1.28	0.016	0.13	123.5	⋯	(a)
J104457	0.013	18.7	6.8	−0.85	0.058	0.95	⋯	⋯	(b), (c)
I Zw 18 NW	0.002	20.1	5.5	−2.70	0.030	3.19	81.2	1.4–1.6	(d), (e), (f), (g)

Note. (1): Redshift. (2): i-band magnitude. The magnitudes are the cmodel magnitudes of Subaru/HSC (SDSS) for J1631+4426 (J104457 and I Zw 18 NW). (3): Stellar mass. (4): Star formation rate (SFR). The SFR of I Zw 18 NW is derived from the Hα flux in the same manner as Isobe et al. (2021). (5): Gas-phase metallicity derived by the direct T_e method. (6): Hβ flux. (7): Rest-frame equivalent width of an Hβ emission line. (8): X-ray luminosity at 0.5–10.0 keV obtained by Chandra X-ray observations. (9): References for values in columns (1)–(8): (a) Kojima et al. (2020), (b) B21, (c) SDSS DR2 (Abazajian et al. 2004), (d) SDSS DR6 (Adelman-McCarthy et al. 2008), (e) TI05, (f) Lebouteiller et al. (2017), (g) Thuan et al. (2004).

Download table as: ASCII Typeset image

In this work, we use observed optical emission lines from hydrogen, helium, oxygen, and sulfur ions to determine the best-fit model parameters. The optical emission lines and their flux ratios are listed in Table 2. We avoid using Hα for the fitting of J1631+4426 because the reported Hα flux is not consistent with the Balmer decrements. We also avoid using Hα for the fitting of I Zw 18 NW because the value of Hα flux has been obtained using a different spectrometer from the one used to obtain the fluxes of the other lines. We use line ratios between sulfur doublets [S ii]λ6316/λ6731 instead of line ratios between sulfur and hydrogen to minimize effects from the uncertainty of the sulfur to oxygen abundance ratio. Another reason to use [S ii]λ6316/λ6731 is that this line ratio is sensitive to the gas density in H ii regions (Osterbrock & Ferland 2006). We avoid including emission line fluxes from other elements to minimize systematics of abundance ratios. We note that He ii λ4686 lines of J1631+4426 and J104457 are narrow and have no broad-line features (Berg et al. 2021; Kojima et al. 2021). There is a report of a Wolf–Rayet bump centered at 4645 Å for I Zw 18 NW, but no underlying broad component is detected for He ii λ4686 (Legrand et al. 1997). Therefore, we treat He ii λ4686 lines of our three galaxies as the nebular origin. All the line ratios in Table 2 use the dust-corrected fluxes.

Table 2. Fluxes and Flux Ratios

Ion	J1631+4426	J104457	I Zw 18 NW
F_λ/F_Hβ

[O ii]λλ3727,3729	50.12 ± 2.66	26.37 ± 0.27	37.66 ± 0.61
He i λ4026	⋯	1.72 ± 0.03	1.03 ± 0.09
Hδ	27.53 ± 0.65	25.59 ± 0.42	26.52 ± 0.42
Hγ	46.88 ± 0.50	46.55 ± 0.67	47.53 ± 0.71
[O iii]λ4363	8.18 ± 0.48	13.51 ± 0.21	6.74 ± 0.14
He i λ4471	⋯	3.97 ± 0.09	2.76 ± 0.08
He ii λ4686	2.32 ± 0.38	1.80 ± 0.04	3.44 ± 0.09
Hβ	100.00 ± 0.37	100.0 ± 1.4	100.00 ± 1.45
[O iii]λ4959	55.76 ± 0.34	143.0 ± 1.5	67.15 ± 0.97
[O iii]λ5007	170.92 ± 0.38	427.5 ± 4.3	200.62 ± 2.89
Hα	⋯	296.7 ± 4.2	⋯

Line Ratios

[S ii]λ6316/λ6731	1.41 ± 0.40^a	1.23 ± 0.04	1.30 ± 0.61

Notes. Dust-corrected emission line fluxes and emission line ratios for all galaxies in our sample. The fluxes are all scaled with Hβ = 100, and taken from K21, B21, and TI05, except for [S ii] λ6316/λ6731 of J1631+4426.

^aThis line ratio of J1631+4426 is measured with the spectrum of Kojima et al. (2020) in the same manner as K21.

Download table as: ASCII Typeset image

2.2. Photoionization Models

We use the spectral synthetic code cloudy (version 17.02; Ferland et al. 2017) to calculate the emission line fluxes from a gas cloud ionized by a central ionizing source. With cloudy, we simulate physical conditions within a gas cloud to predict the emission line fluxes. To obtain the emission line fluxes, we use photoionization models with free parameters of ionizing spectra and nebulae. We describe these free parameters as well as the fixed parameters below.

2.2.1. Geometry and Density

In our photoionization model, we assume a spherical shell of the gas cloud surrounding the central ionizing source, in accordance with the default setting of cloudy.⁶ We fix the inner radius (outer radius) of the gas cloud R_in (R_out) at the default value. The default value of R_out is at R_out = 10³⁰ cm, and this large R_out value produces an effectively plane-parallel geometry. For simplicity, we assume the gas cloud with the constant hydrogen density n_H. To cover a range of typical n_H values for the hydrogen density of SFGs, we apply the hydrogen density ranging in 10^0.5 ≤ n_H ≤ 10⁵ cm⁻³.

2.2.2. Ionizing Spectra

Previous studies have examined ionizing spectra of blackbody, stellar, and hard radiation sources (e.g., ULXs), and combinations of these sources (e.g., Fernández et al. 2022; Olivier et al. 2021; Simmonds et al. 2021). However, for all the soft and hard radiation models, whether a single self-consistent model can explain the emission line fluxes from a single galaxy remains unclear. To solve this, we use a generalized ionizing spectrum composed of blackbody and power-law radiation. The flux of this spectrum at frequency ν is given as

$\begin{eqnarray}&&{F}_{\nu }=B(\nu ,\,{T}_{\mathrm{BB}})+{C}_{\mathrm{mix}}P(\nu ,\,{\alpha }_{{\rm{X}}}).\end{eqnarray} \tag{ 1 }$

B(ν, T_BB) and P(ν, α_X) are given as

$\begin{eqnarray}&&B(\nu ,\,{T}_{\mathrm{BB}})=\displaystyle \frac{h{\nu }^{3}}{\exp (h\nu /{{kT}}_{\mathrm{BB}})-1},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&P(\nu ,\,{\alpha }_{{\rm{X}}})={\nu }^{{\alpha }_{{\rm{X}}}}\exp (h\nu /{E}_{\mathrm{lc}})\exp (-h\nu /{E}_{\mathrm{hc}}),\end{eqnarray} \tag{ 3 }$

where T_BB, C_mix, α_X, h, and k are the blackbody temperature, the mixing parameter, the power-law index, the Planck constant, and the Boltzmann constant, respectively. For the power-law component, we define the lower (higher) cut energy E_lc (E_hc). We set E_lc = 0.1 Ryd (E_hc = 10,000 Ryd) to avoid strong free–free (pair-creation and Compton) heating. We limit the blackbody temperature and power-law index in the ranges of $4\leqslant \mathrm{log}{T}_{\mathrm{BB}}\,\leqslant 6$ and −3 ≤ α_X ≤ 1, respectively. For C_mix, we introduce the alternative parameter a_mix, defined as

$\begin{eqnarray}&&{a}_{\mathrm{mix}}\equiv {C}_{\mathrm{mix}}\displaystyle \frac{P(\nu =1\,\mathrm{Ryd}/h,{\alpha }_{{\rm{X}}})}{B(\nu =1\,\mathrm{Ryd}/h,{T}_{\mathrm{BB}})}.\end{eqnarray} \tag{ 4 }$

In our model calculations, we specify a value for a_mix instead of a value for C_mix to match the specification of cloudy. We allow a_mix to vary in the range of $-3\leqslant \mathrm{log}{a}_{\mathrm{mix}}\leqslant 3$ . Our generalized spectral models cover various spectrum shapes. For example, (T_BB, α_X, a_mix) = (10⁵ K, −1.8, −0.1) can mimic the FUV spectrum shape of an AGN model generated by cloudy "AGN" continuum command with the default parameters.

2.2.3. Chemical Abundance

We use the solar abundance ratios of Grevesse et al. (2010; i.e., GASS10) as an initial condition of our cloudy model calculations. We scale all the abundances linearly with the oxygen abundance (hereafter referred to as metallicity) Z, except for those of helium, carbon, and nitrogen. We assume that the helium and carbon abundances follow the nonlinear relationships that are given by Dopita et al. (2006). For the scaling of nitrogen with metallicity, we use the formula given by López-Sánchez et al. (2012). We adopt a metallicity range of $-3\leqslant \mathrm{log}\,Z/{Z}_{\odot }\leqslant 0$ . Note that we do not use emission lines from carbon and nitrogen, so the uncertainty for these ions' abundances will only have minor effects on our results.

2.2.4. Ionizing Spectrum Intensity

The intensity of the ionizing spectrum is specified with the ionization parameter U. The ionization parameter is defined as

$\begin{eqnarray}&&U\equiv \displaystyle \frac{Q({{\rm{H}}}^{0})}{4\pi {{R}_{{\rm{S}}}}^{2}{n}_{{\rm{H}}}c},\end{eqnarray} \tag{ 5 }$

where Q(H⁰), R_S, and c are the intensity of ionizing photons above the Lyman limit (≤912 Å), the Strömgren radius, and the speed of light, respectively. Because R_S is nearly equal to R_in in the geometry of our model, we approximate U by substituting R_S = R_in. We impose a flat prior for U in the range of $-5\leqslant \mathrm{log}\,U\leqslant -0.5$ .

2.2.5. Other Model Specifications

Our models are truncated at a neutral hydrogen column density N_{H I} at N_{H I} = 10²¹ cm⁻². We normalize the output emission line fluxes relative to the model Hβ flux for convenience of calculation.

2.3. Markov Chain Monte Carlo Methods

To estimate the best-fit parameters, we combine our cloudy photoionization models with MCMC methods. To run the MCMC methods, we use emcee (Ferland et al. 2013), a Python implementation of an affine invariant MCMC sampling algorithm. We use the log-likelihood Gaussian function in our MCMC methods. Our log-likelihood function $\mathrm{ln}{ \mathcal L }$ is

$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & -\displaystyle \frac{1}{2}\displaystyle \sum _{\lambda \in {\rm{\Lambda }}}\left[{\left(\displaystyle \frac{{F}_{\lambda ,\mathrm{obs}}-{F}_{\lambda ,\mathrm{mod}}}{{\sigma }_{\lambda ,\mathrm{obs}}}\right)}^{2}\right.\\ & & \left.+\mathrm{ln}(2\pi {\sigma }_{\lambda ,\mathrm{obs}}^{2})\right],\end{array}\end{eqnarray} \tag{ 6 }$

where Λ, F_λ,obs, ${F}_{\lambda ,\mathrm{mod}}$ , and σ_λ,obs are the set of emission lines used, the observed emission line flux at wavelength λ, the model emission line flux at wavelength λ, and the observed error of the emission line at wavelength λ, respectively. ${F}_{\lambda ,\mathrm{mod}}$ in Equation (6) is given by

$\begin{eqnarray}&&{F}_{\lambda ,\mathrm{mod}}={N}_{{\rm{H}}\beta }\displaystyle \frac{{F}_{\lambda ,\mathrm{cloudy}}}{{F}_{{\rm{H}}\beta ,\mathrm{cloudy}}},\end{eqnarray} \tag{ 7 }$

where F_λ,cloudy and N_Hβ are the output of cloudy for the emission line flux at the wavelength λ and the normalization factor for an Hβ emission line, respectively. F_λ,cloudy values are normalized at F_Hβ,cloudy = 1 by default in cloudy. We introduce N_Hβ as a new free parameter to scale the model fluxes to the observed fluxes. We allow N_Hβ (i.e., ${F}_{{\rm{H}}\beta ,\mathrm{mod}}$ by definition) to vary within a range of 3σ_Hβ,obs from F_Hβ,obs. Here, F_λ,obs values are normalized at F_Hβ,obs = 100 as shown in Table 2. We have a total of seven free parameters in our photoionization models. As a prior probability distributions for each model parameters, we adopt a simple uniform distributions in the range permitted for each variable. We summarize the prior distributions for all free parameters in Table 3.

Table 3. Prior Distributions for Free Parameters

	Parameter	Prior Range
(1)	$\mathrm{log}{a}_{\mathrm{mix}}$	[−3, 3]
(2)	$\mathrm{log}{T}_{\mathrm{BB}}/{\rm{K}}$	[4, 6]
(3)	α_X	[−3, 1]
(4)	$\mathrm{log}U$	[−5, −0.5]
(5)	$\mathrm{log}{n}_{{\rm{H}}}/{\mathrm{cm}}^{-3}$	[0.5, 5]
(6)	$\mathrm{log}Z/{Z}_{\odot }$	[−3, 0]
(7)	N_Hβ	within 3σ_Hβ,obs

Note. (1): Ratio of the blackbody flux to the power-law flux at 1 Ryd. (2): Blackbody temperature. (3): Power-law index. (4): Ionization parameter. (5): Hydrogen density. (6): Gas-phase metallicity. (7): Normalization factor for an Hβ emission line (i.e., model emission line flux for Hβ). Here, the observed Hβ emission line flux values are normalized at 100.

Download table as: ASCII Typeset image

We initialize the positions of the parameters by randomly selecting values from the prior distributions. We use 40 walkers and let each chain run for 1000 steps. We define the "best-fit" parameter set as a parameter set with the maximum likelihood value in the sampled parameter sets. The uncertainty is defined by the range between the minimum and maximum values of the parameter sets that satisfy a condition

$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & \geqslant & \mathrm{ln}{{ \mathcal L }}_{3\sigma }\\ & \equiv & -\displaystyle \frac{1}{2}\displaystyle \sum _{\lambda \in {\rm{\Lambda }}}\left[{\left(\displaystyle \frac{3{\sigma }_{\lambda ,\mathrm{obs}}}{{\sigma }_{\lambda ,\mathrm{obs}}}\right)}^{2}+\mathrm{ln}(2\pi {\sigma }_{\lambda ,\mathrm{obs}}^{2})\right].\end{array}\end{eqnarray} \tag{ 8 }$

Note that this uncertainty is only a rough standard, and the actual uncertainty could be more strictly given.

3. Results

3.1. Posterior Probability Distribution Functions

We obtain the best-fit parameters for all three galaxies in our sample using the method described in Section 2. We present the posterior probability distribution function (PDF) triangles for J1631+4426, J104457, and I Zw 18 NW in Figures 2–4, respectively.

**Figure 2.** Posterior PDF triangles for J1631+4426. One-dimensional PDFs are shown along the diagonals, while two-dimensional PDFs are shown on the off-diagonals. The darker regions on the two-dimensional PDFs indicate the higher density of the parameter sets. The red lines (black dashed lines) indicate the values of the best-fit parameters (uncertainty ranges). The yellow dots are the sampled parameter sets within the uncertainty ranges.
Download figure:
Standard image High-resolution image

**Figure 3.** Posterior PDF triangles for J104457. One-dimensional PDFs are shown along the diagonals, while two-dimensional PDFs are shown on the off-diagonals. The darker regions on the two-dimensional PDFs indicate the higher density of the parameter sets. The red lines (black dashed lines) indicate the values of the best-fit parameters (uncertainty ranges). The yellow dots are the sampled parameter sets within the uncertainty ranges.
Download figure:
Standard image High-resolution image

**Figure 4.** Posterior PDF triangles for I Zw 18 NW. One-dimensional PDFs are shown along the diagonals, while two-dimensional PDFs are shown on the off-diagonals. The darker regions on the two-dimensional PDFs indicate the higher density of the parameter sets. The red lines (black dashed lines) indicate the values of the best-fit parameters (uncertainty ranges). The yellow dots are the sampled parameter sets within the uncertainty ranges.
Download figure:
Standard image High-resolution image

One-dimensional and two-dimensional PDFs are plotted on the diagonal and off-diagonal, respectively. The gray scale in each of the two-dimensional PDFs represents the probability density of the sampled parameters. The darker region on the two-dimensional PDFs indicates a higher density of the sampled parameter sets. The red lines (the black dashed lines) indicate the values of the best-fit parameters (uncertainty ranges). The yellow dots are the sampled parameter sets within the uncertainty ranges. We summarize the best-fit parameters and their uncertainties in Table 4.

Table 4. Best-fit Parameters

	Parameters	J1631+4426	J104457	I Zw 18 NW
Best-fit Parameters

(1)	$\mathrm{log}{a}_{\mathrm{mix}}$	$-{2.05}_{-0.84}^{+0.74}$	$-{2.81}_{-0.08}^{+0.83}$	$-{1.18}_{-0.35}^{+0.14}$
(2)	$\mathrm{log}{T}_{\mathrm{BB}}/{\rm{K}}$	${4.54}_{-0.14}^{+0.02}$	${4.74}_{-0.11}^{+0.03}$	${4.45}_{-0.12}^{+0.05}$
(3)	α_X	$-{0.53}_{-0.71}^{+0.70}$	${0.26}_{-0.80}^{+0.09}$	$-{1.37}_{-0.20}^{+0.45}$
(4)	$\mathrm{log}U$	$-{1.55}_{-0.68}^{+0.05}$	$-{1.95}_{-0.07}^{+0.34}$	$-{1.67}_{-0.34}^{+0.22}$
(5)	$\mathrm{log}{n}_{{\rm{H}}}/{\mathrm{cm}}^{-3}$	${1.74}_{-0.39}^{+2.48}$	${2.82}_{-0.63}^{+0.19}$	${2.12}_{-0.82}^{+0.97}$
(6)	$\mathrm{log}Z/{Z}_{\odot }$	$-{1.61}_{-0.05}^{+0.19}$	$-{1.19}_{-0.03}^{+0.05}$	$-{1.48}_{-0.05}^{+0.07}$
(7)	N_Hβ	${99.9}_{-0.4}^{+0.9}$	${102.3}_{-3.5}^{+1.9}$	${101.9}_{-1.1}^{+1.4}$

Note. The best-fit parameters are the sampled parameters that maximize $\mathrm{log}{ \mathcal L }$ in Equation (6). The uncertainties for best-fit parameters are derived from the condition shown in Equation (8). (1): Ratio of blackbody flux to power-law flux at 1 Ryd. (2): Blackbody temperature. (3): Power-law index. (4): Ionization parameter. (5): Hydrogen density. (6): Gas-phase metallicity. (7): Normalization factor for an Hβ emission line. (i.e., model emission line flux for Hβ). Here, the observed Hβ emission line flux values are normalized at 100.

Download table as: ASCII Typeset image

In Figures 2–4, the prior ranges are wide enough to reveal the overall shapes of the PDFs that are not distorted by the prior boundaries except for that of $\mathrm{log}{a}_{\mathrm{mix}}$ of J104457. However, the uncertainty ranges for $\mathrm{log}{a}_{\mathrm{mix}}$ of J104457 are within its prior range. The PDF for some of the parameters has multimodal distribution. In contrast, parameter sets within the uncertainty ranges (yellow dots) are concentrated around the best-fit parameters. This suggests that peaks other than those of the best-fit parameters are just the local maxima of Equation (6). The best-fit values for the metallicity of our three galaxies are comparable (within 0.2 dex) with the literature values of metallicity.

3.2. Reproduction of Emission Lines

We compare the observed emission line fluxes with the model emission line fluxes reproduced with the best-fit parameters. In Figure 5, we present the normalized differences between the observed and the model fluxes for emission lines in the black dots. We also plot the 1σ_λ,obs and 3σ_λ,obs ranges for the emission lines in the red and yellow bars, respectively. From Figure 5, we find that the best-fit models simultaneously reproduce all the emission lines within around 3σ_λ,obs for all of our three galaxies. Interestingly, [O iii]λ4363 is the most deviant line for J1631+4426. This deviance is probably because our current model for the MCMC method could be weighting too much on reproducing [O iii]λλ4959,5007. The 1σ values for [O iii]λλ4959,5007are below 1% of the flux values, while the 1σ values for [O iii]λ4363 are about 5% of the flux value.

3.3. Ionizing Spectrum Shapes

In Figure 6, we depict the general ionizing spectrum shapes reproduced from the best-fit parameters for our three galaxies . Hereafter, we call these ionizing spectra as the "best-fit" ionizing spectra. For comparison, we overplot the ionizing spectrum of the BPASS binary model at the stellar age 10 Myr and the stellar metallicity Z_* at 0.05 Z_⊙. We use BPASS v2.2.1 (Eldridge et al. 2017; Stanway & Eldridge 2018) with a Kroupa initial mass function (IMF) and a high mass cutoff at 300 M_⊙. All ionizing spectra are normalized at 910 Å. We also plot the uncertainties for the spectral fluxes at the ionization energies of the ions used in our parameter search. The uncertainty bars are given by the maximum and the minimum spectral flux values reproduced with the sampled parameter sets satisfying Equation (8). We note again that these uncertainties are only rough standards, and the actual uncertainties could be be more strictly given.

Between the ionization energies of H⁰ (13.6 eV) and He⁺ (54.4 eV), we find that the best-fit ionizing spectra have two types of general spectrum shapes. The best-fit ionizing spectrum for J104457 has a convex upward spectrum shape like the one for BPASS binary models, except that most of the He ii ionizing photons are absorbed by the stellar atmosphere in the BPASS binary models. On the other hand, the best-fit ionizing spectra for J1631+4426 and I Zw 18 NW have convex downward spectrum shapes that BPASS binary models cannot reproduce. These convex downward spectrum shapes have dominant contributions from the power-law radiation around the ionization energy of He⁺ (54.4 eV). These two types of general spectrum shapes indicate a diversity of the ionizing spectrum shapes for EMPGs.

3.4. Comparison with Other Models

O21 claim that the combination of a BPASS binary model and blackbody radiation can reproduce emission lines from J104457 in a self-consistent manner. We compare our self-consistent model for J104457 with the model for J104457 presented in O21. We construct a mock model for the model presented O21 using the model parameters listed in the bottom column of Table 5 in the same paper. We use the same BPASS version and initial mass function as O21 for the mock model. In Figure 7(a), we present a spectrum of the mock model by the blue line, the spectrum from O21 by the green line, and our result for J104457 by the red solid line. Only the spectrum for below 54.4 eV is shown for the model in O21. We find that the spectrum in O21 has relatively higher fluxes around 54.4 eV than the spectrum of our result.

**Figure 7.** Comparison between our result for J104457 and other models. (a) The blue, green, red solid, and red dashed lines are the spectrum of the mock model for O21, the model in O21, our result for J104457, and our result when we include [Ne v]λ3426/[Ne iii]λ3869, respectively. Only the spectrum below 54.4 eV is shown for the model in O21. The other symbols are the same as those in Figure 6. (b) Emission line fluxes and emission line ratios from two models. The blue (back) dots are the emission line fluxes and emission line ratio from the mock model of O21 (our result for J104457). The other symbols are the same as those in Figure 5.
Download figure:
Standard image High-resolution image

Table 5. Best-fit Parameters with [Ne v] for J104457

	Parameters	Best-fit Values
(1)	$\mathrm{log}{a}_{\mathrm{mix}}$	−2.97
(2)	$\mathrm{log}{T}_{\mathrm{BB}}/{\rm{K}}$	4.80
(3)	α_X	0.37
(4)	$\mathrm{log}U$	−2.05
(5)	$\mathrm{log}{n}_{{\rm{H}}}/{\mathrm{cm}}^{-3}$	2.84
(6)	$\mathrm{log}Z/{Z}_{\odot }$	−1.16
(7)	N_Hβ	102.8

Note. (1): Ratio of the blackbody flux to the power-law flux at 1 Ryd. (2): Blackbody temperature. (3): Power-law index. (4): Ionization parameter. (5): Hydrogen density. (6): Gas-phase metallicity. (7): Normalization factor for an Hβ emission line. The best-fit values listed in this table are obtained with the analysis described in Section 2 but taking [Ne v]λ3462/[Ne iii]λ3869 into consideration.

Download table as: ASCII Typeset image

In Figure 7(b), we show the emission line fluxes and emission line ratios for the mock model (our model) as blue (black) dots. The mock model of O21 reproduces the observed He ii λ4686 fairly well. However, other emission lines, especially the low-ionization emission line fluxes such as [O ii]λλ3727,3729, He i λ4026, and He i λ4471, are systematically underproduced, unlike our result. One possible explanation for this underproduction by the mock model is that the general spectrum shape for the mock model has a relatively lower ratio of low-energy (13.6–24.6 eV) fluxes to high-energy (35.1–54.4 eV) fluxes than that of our result.

Next, we consider how taking [Ne v] emission lines into consideration affects our result. [Ne v] emission lines are detected in some of the He ii emitting EMPGs including J104457 (e.g., Izotov et al. 2004, 2021a; Senchyna et al. 2017; Isobe et al. 2022). Because [Ne v] is a very-high-ionization line with an ionization potential of 97.1 eV, its relative flux to Hβ can be used to constrain the contribution of very hard (i.e., X-ray) radiation to the ionizing spectrum shapes (e.g., Lebouteiller et al. 2017; Senchyna et al. 2020; Simmonds et al. 2021). Simmonds et al. (2021) show that some ULX models can reproduce He ii λ4686 and [Ne v]λ3426 flux values that are comparable to those from observations.

We conduct the same method as described in Section 2, but add [Ne v]λ3426 into consideration. We only examine for J104457 because [Ne v] emission lines are not detected in the other two galaxies, most likely due to the lack of sensitivity. To avoid uncertainty in the abundance ratio, we use [Ne v]λ3426/[Ne iii]λ3869 as we have done to [S ii]λλ6316,6731. The observed [Ne v]λ3426/[Ne iii]λ3869 value and its 1σ error for J104457 are 0.00347 ± 0.00315 (B21). We obtain similar best-fit values for each parameter from MCMC. In Figure 7, we have overplotted the ionizing spectrum produced from the parameter in Table 5 as the red dashed line. Although the blackbody component has a higher temperature, the newly obtained ionizing spectrum generally has a similar spectrum shape as the one we obtain without considering [Ne v]λ3426/[Ne iii]λ3869. The [Ne v]λ3426/[Ne iii]λ3869 value produced from new ionizing spectrum shapes is produced within around 3σ of the observed value. Moreover, this new ionizing spectrum shape produces other emission lines within around 3σ. Therefore, we confirm that the general spectrum shape for J104457 is consistent with [Ne v]λ3426, which requires very-high-energy (>97.1 eV) ionizing photons.

4. Discussion

4.1. Contributions from Stellar Radiation

Our best-fit ionizing spectra have blackbody temperature around (3–5) × 10⁴ K. B-type (O-type) stars have surface temperatures of (1–3) × 10⁴ K (≥3 × 10⁴ K). This suggests that blackbody radiation in our best-fit ionizing spectra, especially that of J104457, can be explained by B- and O-type stars. However, the convex downward spectra of J1631+4426 and I Zw 18 NW suggest that contributions from stars to hard photons around 35.1–54.4 eV might not be as significant as expected in BPASS binary models. Stellar models including fast-rotating stars (e.g., Kubátová et al. 2019) and Population III stars (e.g., Grisdale et al. 2021) also have convex upward spectra around 13.6–54.4 eV, unlike our results. These mismatches suggest that hard radiation from stellar sources cannot solely explain He ii ionizing photons for the two of our galaxies with convex downward spectrum shapes.

We also note that our general spectrum shapes do not take account stellar absorption. Stellar absorption weakens the intrinsic stellar radiation above ionization energies of H⁰ (13.6 eV), He⁰ (24.6 eV), and He⁺ (54.4 eV). This raises the question of whether stellar models with stellar absorption can really reproduce ionizing spectra of our three galaxies, including the blackbody-dominated spectra of J104457. More sophisticated spectral modeling is needed to further examine the stellar contributions, but this is beyond the scope of our paper.

4.2. Contributions from Black Holes

We examine the contributions from black holes (BHs) such as in AGNs and HMXBs/ULXs. BHs can produce very hard radiation through nonthermal radiation from Comptonization/synchrotron radiation. Because nonthermal radiation has a power-law-like shape, radiation from a BH could explain the power-law-like component seen in the ionizing spectrum of J1631+4426 and I Zw 18 NW. BH radiation also has hot (≳10⁵ K) thermal radiation from their accretion disks. Temperatures of accretion disks are roughly proportional to ${M}_{\mathrm{BH}}^{-1/4}$ , where M_BH is the mass of BH. Following the relation between the accretion disk temperature and M_BH, intermediate-mass BHs (10² ≤ M_BH ≤ 10⁵ M_⊙) and less massive ones would not produce sufficient ionizing photons in the soft energy range of <54.4 eV. For this reason, we will first examine if an AGN can solely account for the general spectrum shapes of our three galaxies by comparing them with spectra of AGNs. Then, we will investigate whether the combinations of a BH and stellar population can explain the general spectrum shapes of our three galaxies.

4.2.1. Only BHs

The inner-disk temperature T_in of the accretion disk is around 10⁴⁻⁵ K for AGNs. The blackbody radiation from the inner disk of an AGN gives a big blue bump (BBB) feature in the UV range of the ionizing spectrum. The BBB feature combined with nonthermal radiation can produce a spectrum shape like the one reproduced from the combination of blackbody and power-law radiation. We compare the shapes of four different AGN spectra with the general spectrum shapes of our three galaxies in Figure 8. These four AGN spectra are taken from Ferland et al. (2020), and they represent mean spectra of observed Type-1 AGNs grouped by Eddington ratios (Highest, High, Intermediate, and Low). As shown in Figure 8, all four AGN spectra have power-law like shapes between 13.6 and 54.4 eV, unlike the ionizing spectra of our three galaxies. The deformation of a disk blackbody by Comptonization in observed AGNs might be the reason why AGN spectra do not have blackbody-like shapes in 13.6–54.4 eV (see Ferland et al. 2020). From this comparison, we cannot conclude that AGNs can solely explain the general spectrum shapes for our three galaxies. However, we note that the ionizing spectrum shapes of AGNs remain unclear especially around the Lyman limit to soft X-rays because of intervening absorption.

**Figure 8.** Comparison of our result with AGN spectra. The dashed–dotted lines represent the mean AGN spectra at different Eddington ratios (Highest, High, Intermediate, and Low) taken from Ferland et al. (2020). The other symbols are the same as those in Figure 6.
Download figure:
Standard image High-resolution image

4.2.2. BHs and Stars

We will explore the combinations of BHs and stellar sources to see if they can account for the general spectrum shapes of our three galaxies. We choose ULXs as nonstellar sources because they are popular candidates for the origin of He ii. Moreover, there are detections of a ULX in I Zw 18 through X-ray observations (e.g., Kaaret & Feng 2013). As a ULX spectrum, we use a spectral model described in Gierliński et al. (2009). According to Simmonds et al. (2021), this ULX spectral model reproduces emission line ratios in metal-poor galaxies better than other ULX spectral models. We produce the ULX spectrum models from numerical calculations. To run the numerical calculations, we have referred to the code of the DISKIR model provided in xspec (Arnaud 1996). As presented in Simmonds et al. (2021), ULX spectral models cannot produce enough photons in UV region; thus we combine a ULX spectral model with a stellar spectral model (hereafter ULX + BPASS single model). We use spectra reproduced from BPASS single-star models as stellar spectra. We assume a single-aged stellar population with the same IMF as the BPASS binary model we presented in Figure 6. We do not use BPASS binary models because adding extra hard components from a ULX would not alter the convex upward spectrum shape of the spectra reproduced from BPASS binary models. We use fluxes at the ionization energies of H⁰, He⁰, O⁺, and He⁺ for the best-fit ionizing spectra to restrict the ULX + BPASS single models for our three galaxies. We allow a stellar age t and a ratio between 0.5 and 8 keV of the X-ray luminosity L_X to the SFR (L_X/SFR). An SFR for L_X/SFR is calculated from the UV luminosity at 1500 Å using the relationships presented in Kennicutt (1998). For each of our three galaxies, we fix the value of the stellar metallicity Z_* at the best-fit value of the gas-phase metallicity.

In Table 6, we summarize the parameters of the ULX + BPASS single models that reproduce the general spectrum shapes for our three galaxies. We compare the general spectrum shapes with the ULX + BPASS single-star models in Figure 9. The spectrum of ULX + BPASS single-star models for J1631+4426, J104457, and I Zw 18 NW are represented by the red, blue, and orange thin lines, respectively. We successfully reproduce the general spectrum shape of the best-fit ionizing spectra between 13.6 and 54.4 eV for our three galaxies, despite the diversity in the general spectrum shapes. For both types of the spectrum shapes, the L_X/SFR values are ≲10⁴¹ erg s⁻¹/M_⊙ yr⁻¹, which are comparable to the values obtained in Simmonds et al. (2021). The blackbody-like spectrum shape of J104457 can be explained by prominent stellar contributions relative to that of the ULX, while the convex downward spectrum shapes of J1631+4426 and I Zw 18 NW can be explained by the less prominent stellar contributions around 35.1–54.4 eV than that of J104457. Moreover, the stellar component for J104457 has a stellar age around 2 Myr, while the stellar components for J1631+4426 and I Zw 18 NW have a stellar age above 5 Myr. This suggests that the diversity of the spectrum shapes can arise from the stellar age difference of the stellar components. The stellar contributions become less prominent at ≳5 Myr because massive stars with masses ≳40 M_⊙ already reach their lifetimes (Xiao et al. 2018).

**Figure 9.** ULX+BPASS single models for J1631+4426 (red thin line), J104457 (blue thin line), and I Zw 18 NW (orange thin line). The other symbols are the same as those in Figure 6.
Download figure:
Standard image High-resolution image

Table 6. Parameters for ULX + BPASS Single Models

Name of the Galaxy	t	$\mathrm{log}{L}_{{\rm{X}}}/\mathrm{SFR}$
	(Myr)	(erg s⁻¹/M_⊙ yr⁻¹)
	(1)	(2)
J1631+4426	8.2	40.7
J104457	2.2	41.5
I Zw 18 NW	8.5	41.0

Note. (1): Stellar age. We assume a single-aged stellar populations. (2): 0.5–8 keV X-ray luminosity to SFR.

Download table as: ASCII Typeset image

4.3. Other Possibilities

We have only examined the ULX as a nonstellar component when considering the combination with the stellar component, but it does not limit the possibilities of other sources such as low-luminosity AGNs and supersoft X-ray sources. Photoionization by shocks is still not ruled out from the candidate by our analysis. The ionization of He ii by soft X-rays from cluster winds and superbubbles also remains a possibility (see Oskinova & Schaerer 2022). Further X-ray observations for J1631+4426 and J104457 may enable us to constrain contributions from X-ray emitting sources.

4.4. Multiwavelength Comparison and Relation with Cosmic Reionization

We investigate the consistency between our ionizing spectrum shapes and the observed values of I Zw 18, where multiwavelength data in the X-rays to optical photometry data are available. We assume that I Zw 18 has the same ionizing spectrum shape as our result for I Zw 18 NW to compare the model with the observed values of I Zw 18 as a whole. We use the X-ray luminosity obtained by the X-ray Multi-Mirror Mission (Kaaret & Feng 2013; Lebouteiller et al. 2017). For optical photometry, we use the SDSS i'- and z'-band magnitudes from DR16 (Ahumada et al. 2020) for I Zw 18 NW and I Zw 18 SE to compare the fluxes of the continuum. We do not use photometry data in other bands, where many strong emission lines are detected. To compare the optical photometry data, we replace the blackbody component of the ionizing spectrum with the BPASS single-star model reproduced with the parameters for I Zw 18 NW in Table 6. We normalize the luminosity by adopting a distance of 18.2 Mpc (Aloisi et al. 2007), an Hβ luminosity of 10^39.43 erg s⁻¹, and a relation between the number of hydrogen ionizing photons produced per second and the Hβ luminosity presented in Ono et al. (2010). We show the comparison between the multiwavelength data and the comparison model in Figure 10. The black diamonds represent unfolded X-ray luminosities and the leftward (rightward) black circle represents the i'- (z'-) band magnitude. The incident (net transmitted) spectrum of the comparison model is shown in the yellow (green) line. Here, the net transmitted spectrum is the sum of attenuated incident radiation and the outward portion of the emitted continuum and lines. The orange line is the same ionizing spectrum for I Zw 18 NW as the one presented in Figure 6. From Figure 10, we find that the optical photometry data is slightly below what our comparison model predicts. However, the net transmitted spectrum of our comparison model is consistent with the X-ray luminosities.

**Figure 10.** Comparison between the multiwavelength data and our models for I Zw 18. The black diamonds represent the unfolded X-ray luminosities. The leftward (rightward) black circle represents the i'- (z'-) band magnitudes. The yellow (green) line represents the ionizing spectrum (net transmitted spectrum) of the comparison model. We also overplot the ionizing spectrum of I Zw 18 NW in Figure 6 (BB+PL, orange line) under an assumption that the whole system has the same spectral shape. The other symbols are the same as those in Figure 6.
Download figure:
Standard image High-resolution image

As described in Section 1, a nearby EMPG can be a probe for young galaxies at the EoR. We evaluate the contribution of X-ray photons from the young galaxies to cosmic reionization, assuming that the young galaxies have the same ionizing spectrum shapes as the spectrum of I Zw 18. We adjust the ionizing photon escape fraction of 20% to match the standard picture of cosmic reionization (e.g., Ouchi et al. 2009; Robertson et al. 2013). To do so, we consider the model spectrum that linearly combines the incident ionizing spectrum and the net transmitted spectrum with their luminosities multiplied by 0.2 and 0.8, respectively. In Figure 11, we depict the cumulative fraction of the hydrogen ionizing photon number at each energy in the red line. As shown in Figure 11, about 90% (97%) of the ionizing photons have energy below 24.6 eV (54.4 eV). Moreover, the contributions of the high-energy (>54.4 eV) photons will be further reduced because the photoionization cross-section has a photon energy dependency of approximately E⁻³ (Osterbrock & Ferland 2006). Thus, we conclude that the hard ionizing radiation does not affect the standard picture of cosmic reionization (Robertson 2021).

**Figure 11.** Cumulative fraction of the hydrogen ionizing photon numbers at energies (red line) above the ionization energy of H⁰ (13.6 eV) for our spectrum model. The dotted lines mark the ionization energies of H⁰, O⁰. He⁰, O⁺, and He⁺. The cumulative fraction of the ionizing photon numbers at the ionization energy of He⁰ (24.6 eV) already reaches ∼90% (dashed line). The cumulative fraction at the ionization energy of He⁺ (54.4 eV) exceeds 97% (dotted–dashed line).
Download figure:
Standard image High-resolution image

5. Conclusions

We have searched for the general spectrum shapes of ionizing sources in EMPGs that can explain major emission lines including He ii. We have prepared parameterized ionizing spectra with generalized shapes and use these ionizing spectra in photoionization models with free parameters of ionizing spectra and nebulae. For three EMPGs with strong He ii λ4686, we have determined seven parameters of the ionizing spectra and nebulae using MCMC methods.

Our main findings can be summarized as follows:

1.
We have determined the best-fit parameters for our three galaxies. All the best-fit parameters are comparable to ordinary values. With the best-fit parameters, we successfully reproduced all major emission lines within ≲3σ errors.
2.
The ionizing spectrum produced from the best-fit parameters for J104457 has a blackbody-dominated shape around 13.6–54.4 eV with blackbody temperature 5 × 10⁴ K. In contrast, the ionizing spectra for J1631+4426 and I Zw 18 NW have significant contributions from power-law radiation around 54.4 eV, with relatively low-temperature (3 × 10⁴ K) blackbody radiation. As a result, the general spectrum shapes for J1631+4426 and I Zw 18 NW have convex downward shapes around 13.6–54.4 eV, which are fundamentally different from the shapes predicted by BPASS binary models. These two types of spectrum shapes indicate a diversity in ionizing spectrum shapes in EMPGs.
3.
We compare the obtained general shapes of ionizing spectra with those of the ionizing spectra from stars, AGNs, and ULXs. We find that ULX + BPASS single-star models can explain the general shapes of the obtained ionizing spectra for our three galaxies. Furthermore, we show that the diversity of the spectrum shapes can be attributed to the differences in the stellar age.
4.
We confirm that our result is consistent with the X-ray and optical photometry data for I Zw 18, for which the multiwavelength data are available. We also find that the hard ionizing radiation inferred from our general spectrum shape of I Zw 18 does not contradict the standard picture of cosmic reionization.

Spectroscopic observation of J1631+4426 and I Zw 18 NW for optical [Ne v] emission lines with high sensitivity would also be helpful in terms of constraining contributions of ionizing photons ≳100 eV in He ii emissions. Our new analysis method would provide new insights into the statistical nature of EMPGs if applied to a sample with a large number of EMPGs.

We thank Wako Aoki, Yoshihisa Asada, Kohei Inayoshi, Akio Inoue, Carolina Kehrig, Tomoya Kinugawa, Haruka Kusakabe, Vianney Lebouteiller, Daniel Schaerer, Yuta Tarumi, and Masayuki Umemura for sharing their data and figures and giving us helpful comments. This research is based in part on data collected at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical, and natural significance in Hawaii. This work is supported by the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, as well as KAKENHI Grant-in-Aid for Scientific Research (A) (20H00180, and 21H04467) through the Japan Society for the Promotion of Science (JSPS). This work was supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo. Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

EMPRESS. VII. Ionizing Spectrum Shapes of Extremely Metal-poor Galaxies: Uncovering the Origins of Strong He ii and the Impact on Cosmic Reionization

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction