Characterizing Extreme Emission Line Galaxies II: A Self-Consistent Model of Their Ionizing Spectrum

Observations of high-redshift galaxies ($z>$ 5) have shown that these galaxies have extreme emission lines with equivalent widths much larger than their local star-forming counterparts. Extreme emission line galaxies (EELGs) in the nearby universe are likely analogues to galaxies during the Epoch of Reionization and provide nearby laboratories to understand the physical processes important to the early universe. We use HST/COS and LBT/MODS spectra to study two nearby EELGs, J104457 and J141851. The FUV spectra indicate that these two galaxies contain stellar populations with ages $<\sim$ 10 Myr and metallicities $\leq$ 0.15 Z$_\odot$. We use photoionization modeling to compare emission lines from models of single-age bursts of star-formation to observed emission lines and find that the single-age bursts do not reproduce high-ionization lines including [O III] or very-high ionization lines like He II or [O IV]. Photoionization modeling using the stellar populations fit from the UV continuum similarly are not capable of reproducing the emission lines from the very-high ionization zone. We add a blackbody to the stellar populations fit from the UV continuum to model the necessary high-energy photons to reproduce the very-high ionization lines of He II and [O IV]. We find that we need a blackbody of 80,000 K and $\sim$60-70% of the luminosity from the young stellar population to reproduce the very-high ionization lines while simultaneously reproducing the low- intermediate-, and high-ionization emission lines. Our self-consistent model of the ionizing spectra of two nearby EELGs indicates the presence of a previously unaccounted-for source of hard ionizing photons in reionization analogues.


INTRODUCTION
The first galaxies are likely responsible for producing the photons that ionized the universe during the Epoch of Reionization (EoR, Wise et al. 2014;Robertson et al. 2015;Madau & Haardt 2015;Stanway et al. 2016). With telescopes like the Hubble Space Telescope (HST), Atacama Large Millimeolivier.15@osu.edu * Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NASA 5-26555. ter/submillimeter Array (ALMA), and other telescopes, astronomers have discovered thousands of galaxies in the highredshift universe (z > 6) in an attempt to study these first galaxies (e.g., Bouwens et al. 2015;Finkelstein et al. 2015;Livermore et al. 2017;Atek et al. 2018;Oesch et al. 2018;Bouwens et al. 2020;Harikane et al. 2020). More direct observations of these first galaxies and their environments during the EoR are just on the horizon with the coming era of the James Webb Space Telescope (JWST) and extremely large telescopes (ELTs), in combination with next-generation H I surveys such as the Square Kilometer Array (SKA). Soon we will be able to observe the first galaxies and their effects on the neutral gas surrounding them.
Over the past decade, the rest-frame UV emission-line spectra have been observed from several very high redshift galaxies (z > 6) revealing strong, high-ionization nebular lines (e.g., Stark et al. 2015;Stark 2016;Mainali et al. 2018) that are indicative of very hard radiation fields in lowmetallicity environments. With this in mind, much can be learned from detailed studies of nearby galaxies with properties analogous to what we expect to find in the distant universe. In particular, recent studies have focused on metalpoor, high-ionization galaxies with high specific star formation rates both locally (e.g., Berg et al. 2016;Izotov et al. 2016;Jaskot et al. 2017;Senchyna et al. 2017;Izotov et al. 2018;Senchyna et al. 2020;Berg et al. 2019b,a) and at intermediate redshifts using lensed galaxies (e.g., Hainline et al. 2009;Bayliss et al. 2014;Christensen et al. 2014;James et al. 2014;Rigby et al. 2018a;Berg et al. 2018).
These recently discovered samples of Reionization Analogues have extreme equivalent widths (EWs) in their emission lines from high-ionization species like O III], C III], C IV, and He II. As a result of their large equivalent widths and hard radiation fields, extreme emission line galaxies (EELGs) in the local universe provide a laboratory for understanding the physical processes active in the first galaxies. The large EWs from O III], C IV, and He II in galaxies out to z ∼ 2 have been associated with low metallicity environments and very young stellar populations from recent bursts of star-formation producing highly ionized gas (Atek et al. 2011;van der Wel et al. 2011;Maseda et al. 2013Maseda et al. , 2014Rigby et al. 2015;Berg et al. 2016;Senchyna et al. 2017;Chevallard et al. 2018;Berg et al. 2019b,a;Tang et al. 2019Tang et al. , 2021. Exploring similarly low metallicity environments in local EELGs can help us understand those young stellar populations and the highly ionized gas. These EELGs often demonstrate strong He II emission along with their other extreme emission lines, but the mechanisms producing the high-energy photons at >54eV are rarely understood. A number of studies have explored the production mechanism for the extreme He II emission at λ1640 in the restframe UV and λ4686 in the optical in a number of these reionization analogues (e.g., Schaerer 2002;Kehrig et al. 2015Kehrig et al. , 2018Kehrig et al. , 2021Schaerer et al. 2019;Senchyna et al. 2017Senchyna et al. , 2020. Some of the postulated sources of the He II emission include: Wolf-Rayet stars at low metallicity, high-mass X-ray binaries, metal-poor massive stars, and shocks, but, when examined in more detail, these suggestions tend to be inadequate in their photon production in the energy range needed for extreme He II emission (Kehrig et al. 2015(Kehrig et al. , 2018Senchyna et al. 2020;Wofford et al. 2021).
Here we examine two of these reionization era analogue galaxies. The first of these galaxies is J104457 which has the largest C IV λλ1548,1550 Å EW (-6.71Å ,-2.83Å ) measured in the local universe (Berg et al. 2016). The second galaxy is J141851 which has the largest He II λ1640 Å EW (-2.82Å ) measured in the local universe (Berg et al. 2019b). These two galaxies are prime targets for studying the physical processes required to create high ionization gas and to explore the conditions present in the EoR. A first look at the high-energy ionizing photons in these galaxies was presented in Berg et al. (2019a), but these fascinating objects demand a more complete analysis. Here we present the second paper in a series of work to characterize extreme emission line galaxies. In Paper I, Berg et al. (2021), we presented a detailed examination of the extreme emission lines, nebular physical conditions, and chemical abundances in J104457 and J141851. Owing to the detection of a number of very high-ionization emission lines (e.g., C IV, He II, O IV, Si IV, Ne V, and Ar IV), we adopted a novel 4-zone ionization model and found both EELGs to be characterized by much higher ionization parameters than reported in previous works. The main result of Paper I highlights a persistent high energy ionizing photon production problem (HEIP 3 ) that motivates the self-consistent modeling of their ionizing spectra presented here in Paper II. Paper II is structured as follows: We present the optical and UV spectra of J104457 and J141851 in Section 2. In Section 3, we investigate the UV stellar continuum. The procedure used for the stellar continuum fits is discussed in § 3.1 and the resulting stellar population properties are presented in § 3.2. Additionally, we generalize our findings in order to characterize the spectra of extremely metal poor stellar populations in § 3.3. Section 4 then combines constraints from the stellar continuum fits of Section 3 with the emission line ratios from the optical and UV spectra using CLOUDY (Ferland et al. 2017) photoionization modeling. We explore several input ionizing sources in Section 4.2, including single-aged burst binary models ( § 4.2.1, the best-fit stellar population models ( § 4.2.2), and combined stellar plus blackbody models in ( § 4.2.3). The resulting best fit to the suite of observed emission lines is produced by the latter models, requiring an additional source of very-high-energy ionizing photons. We discuss this result and the possible physical origin of the very-high-energy photons in Section 5. Finally, we present our conclusions in Section 6.

EXTREME UV EMISSION-LINE GALAXY OBSERVATIONS
The main objective of this work is to investigate whether the suite of nebular emission lines from two EELGs can be self-consistently reproduced by photoionization models using the ionizing photons from the massive star populations. To analyze the interplay of the ionizing stellar population and nebular gas of J104457 and J141851, we first measure the nebular emission lines from their UV and optical spectra.  Berg et al. (2019a) from the UV HST/COS G160M spectra (see Figure 2). Equivalent widths are listed for C IV λλ1548,1550, O III λλ1661,1666, and He II λ1640. We list velocity widths (FWHM), ∆V , of O III] λ1666 and He II λ1640, and both blue and red components of the C IV λ1548 profiles. Since the individual C IV lines present as doublets, we list the V blue peak , V red peak , and Vsep. measured for the C IV λ1548 line. The final line of the table lists the spectral resolution from COS used in our continuum fitting. a Data catalogues are available from http://www.sdss3.org/dr10/spectro/ galaxy_mpajhu.php. The Max Planck Institute for Astrophysics/John Hopkins University(MPA/JHU) SDSS data base was produced by a collaboration of researchers(currently or formerly) from the MPA and the JHU. The team is made up of Stephane Charlot (IAP), Guinevere Kauffmann and Simon White (MPA), Tim Heckman (JHU), Christy Tremonti (U. Wisconsin-Madison − formerly JHU) and Jarle Brinchmann (Leiden University − formerly MPA).

HST/COS UV Spectra
In this work we used FUV spectra taken with the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST) that were first reported in Berg et al. (2016), Berg et al. (2019b), and Berg et al. (2019a) from the Hubble Space Tele-scope Project ID 15465. We briefly summarize the observational details below.
All observations were centered on the NUV flux of J104457 and J141851 using the ACQ/Image acquisition mode. Science exposures were then taken in the TIME-TAG mode using the 2.5 PSA aperture. The FP-POS = ALL setting was used to take four images offset from one another in the dispersion direction, increasing the cumulative S/N and mitigating the effects of fixed pattern noise. Spectra were processed with CALCOS version 3.3.4, using the standard EXTRACTSEGMENTBOXCAR aperture method, and binned by the 6 native COS pixels per resolution element such that the nominal δv = 13.1 km s −1 .
The G160M spectra used a central wavelength of 1589 Å and had total exposures of 6439 and 12374s for J104457 and J141851, respectively. This setup provides wavelength coverage that is rich with stellar and nebular features. In particular, these spectra cover several important FUV high-ionization emission lines: Si IV λλ1394, 1403, O IV λλ1401,1405,1407, S IV λλ1405,1406,1417,1424, N IV λλ1483,1487, C IV λλ1548,1550 λλ1661,1666. For the purpose of testing the production of nebular emission lines, not only are the Si IV, C IV, and He II features difficult to interpret due to their potentially complex nebular and stellar contributions, but Si IV and C IV are also resonant lines. Further, the S IV and N IV lines are rather faint or not detected in the UV spectra of J104457 and J141851, and thus, only the O III] and O IV lines, along with He II, since its velocity profile shows the emission is mostly nebular in origin, remain to robustly constrain our models. Therefore, in order to expand our suite of strong UV emission lines constraining the shape of the ionizing continuum, we also consider the low-resolution HST/COS G140L spectra for these galaxies, which provide access to the C III] λλ1907,1909 lines.
The NUV acquisition images, whose flux is dominated by the light from the young ionizing stellar clusters, are shown in the left panels of Figure 1, and highlight the simple, very compact nature of these galaxies. The UV light from J104457 and J141851 is easily contained within the COS 2.5 aperture (gold circle), corresponding to a physical size of 0.5 kpc at their redshifts of z ∼ 0.01. In the right panels of Figure 1, we plot the spatial light distribution perpendicular to the dispersion axis (red dashed line in acquisition images). The spatial trace is smoothed by a Gaussian 1D kernel with σ = 1 pixel and compared to a Gaussian approximation of the full width half maximum (FWHM) of the non-Gaussian instrumental profile of FWHM inst = 0.6 (gold line). Similarly, the spectral light trace is plotted in green. Comparing the observed spatial trace to the Gaussian profile estimates the spatial extent of each galaxy. J141851 has less extended emission and has a fitted FWHM= 0.12 . However, the spatial profile of J104457 is somewhat larger (FWHM= 0.7 ), resulting in a degradation of the spectral resolution.
We measured the resolution of our spectral observations by fitting Gaussian profiles to three Milky Way absorption lines in the spectra of each galaxy (Si II λ1526, C IV λλ1548,1550, and Al II λ1671). For each galaxy we averaged the FWHM of the three lines together to find the spectral resolution. For J104457 and J141851, we found that the FWHM related to the spectral resolution was 188 km s −1 and 104 km s −1 , respectively. These values are larger than the native COS resolution: nominally, the spectral resolution for the G160M grating is R ∼ 20, 000 or ∆v = 15 km s −1 , with the observed spectral resolution broadened due to the considerable spatial extent of the galaxies (Figure 1). Similarly, Berg et al. (2019a) report GFWHMs of the O III] nebular emission lines of ∆v = 87.3, 64.7 km s −1 for J104457 and J141851, respectively (Berg et al. 2019a). Some of this velocity width could be due to the physical properties of the gas, but a degradation of the spectral resolution is also present. Nonetheless, at these resolutions, the individual emission line fluxes and velocities of interest can be resolved and precisely measured.

Optical LBT/MODS Spectra
While the FUV spectra enable modeling of the massive star populations, the FUV contains a limited number of nebular emission lines of only a handful of ions. Recently, Paper I presented a more detailed examination of the extreme emission lines in J104457 and J141851 and their chemical abundances, using the HST/COS G160M spectra together with new optical spectra from the Multi-Object Double Spectrographs (MODS, Pogge et al. 2010) on the Large Binocular Telescope (LBT, Hill et al. 2010). We refer the reader to Paper I for a full description of the observations and Berg et al. (2015) for a detailed discussion of the MODS data reduction pipeline. For convenience, we present a brief overview below.
MODS optical spectra of J104457 and J141851 were obtained on the UT dates of 2018 May 19 and 18, respectively. We obtained simultaneous, moderate-resolution blue (R ∼ 1850) and red (R ∼ 2300) spectra. J104457 and J141851 were each observed using the 1 longslit for 45 min. Spectra were reduced with the MODS reduction pipeline 1 . Onedimensional spectra were corrected for atmospheric extinction and flux calibrated based on observations of flux standard stars (Oke 1990). A portion of the resulting spectra are plotted in Figure 2, showing the broad spectra coverage provided by this LBT/MODS setup, extending from 3200 − 10, 000 Å and including a number of very highionization emission lines: [Ne III] λ3869, [Fe V] λ4227, He II 1 http://www.astronomy.ohio-state.edu/MODS/Software/modsIDL/ λ4686, and [Ar IV] λλ4711,4740. While only the UV spectra can be used to determine the massive star populations, we can use the high-ionization UV and optical emission lines together to better constrain the shape of the ionizing spectral energy distribution.
3. THE STELLAR CONTINUUM

Fitting the UV Continuum
We use the fitting procedure introduced in Chisholm et al. (2019) to determine the properties of the stellar populations in these galaxies. We first mask out any features in the spectrum that are produced by non-stellar sources. For example lines produced by Milky Way absorption and ISM emission and absorption lines associated with the galaxies are masked. We used Leitherer et al. (2011) andBruhweiler et al. (1981) as references to identify the lines to mask. In addition to masking the spectral lines not associated with the stars, we masked the regions of the spectrum where the signal-to-noise was less than 1. These masked regions are shown as the gray bands in Figure 3. We also normalize the spectrum by the flux density at 1495 Å in order to do our fitting.
To obtain the age and metallicity of the stellar populations we compare the spectra to theoretical models of single-age bursts of star formation. We test single star evolution using the high mass-loss rate stellar evolution models (Meynet et al. 1994) from Starburst99 (S99; Leitherer et al. 1999) and binary stars from BPASS v2.2.1 (Eldridge et al. 2017;Stanway & Eldridge 2018) separately. We explore Kroupa initial mass functions (Kroupa 2002) with a broken power-law with a high and low-mass slope of 2.3 and 1.3, respectively, and both 100 M and 300 M high mass cut offs. Our model grids for both single stars and binary stars cover a range of metallicities from 0.05 Z to 2 Z (0.05, 0.2, 0.4, 1.0, and 2.0 Z ) and ages from 1 Myr to 15 Myr in eight non-uniform steps (1, 2, 3, 4, 5, 8, 10, and 15 Myr) selected to capture variations in the FUV stellar continuum features with age. We used these ranges of metallicities and ages because these galaxies appear to be young, metal-poor star forming objects based on their extreme Hβ equivalent widths and high starformation rates (Berg et al. 2016(Berg et al. , 2019b. We do not expect, nor do we find, that high metallicity or old stellar populations create the observed UV spectrum as it is dominated by the light from the most massive stars from young stellar populations, but we included the models to ensure we did not introduce this bias to our fitting method. We convolve the 0.4 Å resolution of the S99 single star models to the spectral resolution of the individual galaxy spectra (Leitherer et al. 1999;Meynet et al. 1994). For the binary stellar models from BPASS, the resolution is 1 Å , so our observed galaxy spectra were convolved to that larger resolution for these fits (Eldridge et al. 2017). The fitting results are shown in Section 3.2. The 2.5 COS aperture used for the UV spectra is shown as a gold circle. The spectrograph position angle is noted in the corner and is used to determine the spatial and spectral axes. In comparison, the 1 LBT/MODS slit (black lines) captures most of the NUV light, but may miss a significant fraction of the extended nebular emission. In the right hand column we plot the spatial and spectral light profiles captured within the COS aperture. The 0.6 instrumental profile of COS is shown in gold.
We fit the stellar continuum assuming that the FUV light is a linear combination of multiple single age and single metallicity bursts of star formation scaled by a linear coefficient. The final intrinsic stellar continuum fit is the sum of all models weighted by their coefficients. We then fit for the best fit dust attenuation value (E(B-V)) by reddening the intrinsic spectrum using a foreground dust screen model and the Reddy et al. (2016) attenuation law. We tested using a Calzetti et al. (2000) and a SMC dust attenuation law (Gordon et al. 2003), but did not find that the continuum fits and stellar parameters, such as age and metallicity, statistically varied (although the E(B-V) parameter changed). The final fit has 41 free parameters (40 linear coefficients for the stellar continuum models and 1 dust attenuation model). We derive light-weighted stellar population parameters (i.e., age and metallicity) by calculating a weighted average using the model parameters and the best fit linear coefficients. In Table 2, we give the light-weighted parameters as well as the ionizing photon production efficiency (ξ ion ; or the number of ionizing photons divided by the luminosity at 1500Å).
To find the errors on these parameters, we run the fitting code on 1000 randomized spectra. During each randomized iteration, we created the randomized spectra by multiplying the observed flux at a given wavelength by a random number with the distribution width dictated by the uncertainty on the spectrum at that wavelength. For each iteration, we estimated the light-weighted properties and then took the standard deviation of the distribution as the uncertainty on the stellar population properties.  , highlighting the blue and Wolf-Rayet bump regions (blue and purple inset windows, respectively). In the blue inset window, there are strong emission lines, but they appear narrow and nebular in origin, whereas there is a lack of any features in the purple inset window for J141851 but a possible C IV bump in J104457 which we discuss in Section 5.3.

Stellar Continuum Parameters
The fits to the observed stellar continua provide estimates of the massive star population of these two EELGs. The first two rows in Table 2 list the estimated light-weighted ages and metallicities of these two galaxies. We include four different sets of models that we used to fit the stellar populations: single-star evolution with a 100 M upper IMF cutoff from S99 (Leitherer et al. 1999;Meynet et al. 1994), binary-star evolution with a 100 M upper IMF cutoff from BPASS (Eldridge et al. 2017), single-star evolution with a 300 M cutoff from S99 (Leitherer et al. 1999;Meynet et al. 1994), and binary-star evolution with a 300 M upper IMF cutoff from BPASS (Eldridge et al. 2017). The fits imply that both galaxies have young ( < ∼ 10 Myr) and low-metallicity (< 0.15 Z ) stellar populations. Due to the modest signal-to-noise ratios, the inferred stellar population properties have relatively high errors, but the estimates using the four different model assumptions are consistent with each other within the 1σ un-certainties. The BPASS models are not systematically older than the S99 models for these two galaxies, likely because the fitting algorithm identified a feature corresponding to an older stellar population in J141851, pushing the S99 models to an older age.
As a measure of the total production of ionizing photons, we use the stellar continuum fits to determine production efficiency of ionizing photons (ξ ion = Q/L 1500 ) by integrating the total number of ionizing photons from the fitted stellar population model (Q) from 1 − 912 Å and dividing by the FUV luminosity density at 1500 Å. We find log(ξ ion ) values of 25.82 and 25.68 dex [photons Hz erg −1 ] for J104457 and J141851, respectively, with relatively large errors (0.8 dex). These values are larger than (but statistically consistent with) the ξ ion values inferred from the optical Hβ flux ratios (Berg et al. 2019b), suggesting that the stellar population fits are producing an approximately sufficient number of hydrogen-ionizing  . The UV spectra of J104457 and J141851 and the best-fit continuum models. The data are plotted in black, and the spectra are normalized by their flux at 1495 Å. We show the continuum fits from the Starburst 99 models (Leitherer et al. 1999) in orange and BPASS models (Eldridge et al. 2017) in pink for J104457 in the top plot and J141851 in the bottom plot. We also plot the residuals, continuum fit subtracted from the data, in the bottom panel of each plot. The zoom-in regions highlight detections of absorption lines from highly ionized metal absorption features -Fe V and S V -that are only present in the photospheres of the most massive and hottest stars.
photons to match the observed Balmer emission lines. The total ionizing photon production is consistent regardless of the assumptions made for the stellar models (single-star versus binary-star evolution and the high-mass cutoff). Since the S99 and BPASS models produce different parameters for their age and metallicity to recreate the observed spectrum, it is surprising that ξ ion is consistent between the two models (See Table 2). However, the combination of older binary stars from BPASS produces a similarly shaped ionizing continuum as the single-stars combined for S99, thus the consistency of the ξ ion values is less surprising. Our high-mass cutoff models do produce ∼0.1 dex more ionizing photons in our log ξ ion calculations than our 100 M models. This difference is similar to the expected difference between the 100 M and 300 M models from BPASS found in Figure 35 of Eldridge et al. (2017). However, the errors on our calculated numbers are 0.8 dex and thus the difference in the 300 M and 100 M models is not significant. The inclusion of the 300 M only changes the fraction of single stars by a few percent in the BPASS models and does not imply a significant number of very massive stars (Stanway & Eldridge 2019).

Characteristic Stellar Spectral Features
In addition to determining the stellar populations of these galaxies, the FUV observations of extreme emission line galaxies present a rare opportunity to study the massive star population properties at very low metallicity (<10% Z ). Figure 3 reveals largely featureless continua, with a very weak C IV 1550 Å wind P-Cygni wind feature observed in J141851 and a nearly flat C IV region in J104457 as can be seen in the residuals in Figure 3. These weak winds are consistent with theoretical predictions of radiatively driven stellar winds, where mass-loss rates are significantly lower at extremely low metallicity (Vink et al. 2001). While the continua are nearly featureless power-law spectra, at our signalto-noise there are three identifiable weak photospheric features: Fe V 1409+1440 Å and S V 1502 Å (see Table 3; Bruhweiler et al. 1981;Leitherer et al. 2011). The combination of the weak C IV P-Cygni features, small highly ionized metal photospheric features, and an otherwise featureless stellar continuum are consistent with the young, lowmetallicity stellar populations inferred from the continuum fits.
Fe V 1409+1440 Å and S V 1502 Å absorption lines are indicative of very massive stars as only hottest stars keep significant amounts of Fe and S in this high of ionization states. If we observe absorption in Fe V and S V, this implies the most massive stars are contributing most of the light for these galaxies, and the massive stars only exist as the youngest stars. In Table 3, we compare the observed Fe V and S V equivalent widths of our two EELGS to two samples of star-forming galaxies at z ∼ 0 (Heckman et al. 2015;Chisholm et al. 2019) and z ∼ 2 (Rigby et al. 2018a,b), as well as the stacked composites for both samples (bottom of each section). The S V equivalent widths of both galaxies are significantly weaker than the comparison samples, likely because the metallicities of the EELGS are a factor of 5-10 lower than the comparison samples. Interestingly, the Fe V 1440 Å detection in J104457, is among the strongest detections in the sample. While these two EELGS have low Fe-abundances Berg et al. (2021), the Fe V feature is actually stronger than the alpha-element tracer S V. This suggests that the iron opacity within the atmospheres of these lowmetallicity stars is non-negligible. The only other spectra with Fe V 1440 Å strongly detected are J0150+1260, which is the second youngest galaxy in the Chisholm et al. (2019) low-redshift sample at 2.20 ± 2.17 Myr, and S000451.7-010321, which is dominated by a very young stellar population with a small contribution from an older population, suggesting that Fe V 1440 Å is a robust indicator of very young stellar populations even at low metallicity. This is consistent with the slightly older stellar population fit to J141851 and the non-detection of Fe V 1440 Å. These detections of Fe and S photospheric lines suggest that the UV dominant stars in these two EELGS are extremely young, massive, and low-metallicity.

Adopted Nebular Property Constraints
Owing to the detection of a number of very high-ionization emission lines (e.g., C IV, He II, O IV, Si IV, Ne V, and Ar IV), we adopted a novel 4-zone ionization model in Paper I and recalculated the nebular properties and chemical abundances of J104457 and J141851. Therefore, for consistency, we adopt the results of Paper I as a better estimate of the nebular properties of J104457 and J141851 in our photoionization models. Here we briefly discuss the nebular properties that we use to constrain the shape of the ionizing spectral energy distributions of EELGs.

Suite of UV+Optical Nebular Emission Lines
In order to compare our observations to photoionization modeling, we first need to precisely measure a suite of UV and optical emission lines that probe a large range in ionization potential energy. We began by adopting the nebular emission line intensities measured in Paper I. The detailed information on the line measurement and dereddening processes can be found in Section 2.3 of Paper I and the subsequent emission line intensities in their corresponding Table 2. To briefly summarize, the UV and optical emission lines were measured in a uniform, consistent manner. The optical emission lines were measured from continuum- Here we report the light-weighted ages and metallicities for both galaxies. We used the Reddy et al. (2016) reddening law to determine the E(B-V). We calculated the number of ionizing photons and report them here as log(ξion).
Object  In Paper I, we used the suite of nebular emission lines to measure the ionization parameter logU for these two EELGs and found that the ionization structure is not well described by a single logU measurement from [O III]/[O II]. Classically, studies use 3 ionization zones to describe the ionization structure of a galaxy or an H II region. Using the 4zone ionization model, we were able to constrain the logU for three ionization zones of the nebula to better understand the high-, intermediate-, and low-ionization energy photons separately. As a result, we adopted these three ionization parameters from Table 3 of Paper I for comparisons with our photoionization modeling.

Photoionization Modeling
The richness and quality of our combined UV and optical spectra of J104457 and J141851 allowed us to perform the first detailed models of these spectra with the goal of simultaneously constraining the direct-method gas-phase properties, inferred ionizing stellar populations, and the resulting nebular emission lines they produce. To do so, we used CLOUDY version C17.01 (Ferland et al. 2017) to run models with radiation field inputs from (1) single-aged bursts of binary star formation from the BPASS 2.2 models (Eldridge et al. 2017;Stanway & Eldridge 2018) in §4.2.1, (2) the stellar population continuum fits from Section 3 in §4.2.2, and (3) combination models of the fitted stellar populations plus an extra high energy source modeled as a blackbody in §4.2.3 in order to recreate our very high energy ionization lines.

Single-Aged Burst Binary Models
For our fiducial model grid, we used BPASS binary, singleaged stellar population models with ages ranging from 1-10 Myr and spanning a metallicity range of 0.005Z < Z < 0.3Z (Eldridge et al. 2017;Stanway & Eldridge 2018). We matched the gas-phase metallicity to the stellar metallicity in the models, as expected for alpha-elements, but investigated sub-solar abundance patterns of carbon and nitrogen relative to oxygen following the empirical trends seen for metal-poor dwarf galaxies (e.g., van Zee & Haynes 2006;Berg et al. 2012;Berg et al. 2019a). Paper I investigated the possibility of enhanced α/Fe ratios in our EELGs but found no evidence that the inclusion of increased α-elements were capable of producing the high-energy photons that are missing from these galaxies. Specifically, we varied the fraction of C/O and N/O relative to solar from 0 to 1 in steps of 0.25, encompassing the measured abundances of N/O = 0.44 N/O and 0.32 N/O and C/O = 0.42 C/O and 0.35 C/O for J104457 and J141851, respectively. Finally, we used a variety of constant densities ranging from 10 cm −3 to 10,000 cm −3 in logarithmic steps. We also used power law density gradients of the form n(r) = n 0 (r 0 )(r/r 0 ) α where n 0 is the density at the illuminated surface of the cloud which we set to 10 cm −3 , 100 cm −3 , 1000 cm −3 , and 10,000 cm −3 , and where r 0 is the radius from the source to the illuminated surface of the cloud, and slopes of α =[-1,-2,-3].
We compare our models to the observed spectra of J104457 and J141851 by comparing the emission line ratios of different species. For our study we use eight emission lines from the optical spectrum relative Similar to Berg et al. (2018) and others, we start by first finding the best fit model from our fiducial grid by using a reduced-χ 2 , χ 2 red fit. We compare the line ratios of low ionization energy species from J104457 and J141851 and the sample of small star forming galaxies presented in Berg et al. (2019a) Figure 4 we have plotted the BPASS single age burst models from 1 Myr (solid lines) to 3 Myr (dashed lines) for metallicities ranging from 0.005 -0.3 Z . The points for each galaxy are color coded by their metallicity. The results of the reduced χ 2 of each fit (χ 2 red = χ 2 /d.o.f.), the logU, the metallicity, the age, the C/O, and the N/O ratios are listed in the section labeled 4.2.1 in Table 5 and Table 6 for J104457 and J141851 respectively. Although these models fit the low-and intermediateionization lines relatively well, the lines with ionization energies through [Ar III], the χ 2 red is very large for the high-and very high-ionization lines, lines with ionization energies for [O III] and above (> 35 eV). The single-burst models do not match the high-ionization emission lines.
When we compare the high and very-high ionization energy lines ratios, O III], [Ar IV], [Ne III], He II λ4686, He II λ1640, and [O IV], predicted by the single-age burst models to the observed ratios from J104457, J141851, and the Berg et al. (2019a) sample, we find an offset in the line ratios. The observed high-ionization lines are systematically higher than the ratios predicted by the single-age bursts. We attempted changing a number of additional parameters in our models including ionization parameter, density, and power-law density distributions in order to recreate the observations, but these changes were still not enough to reproduce the observed line ratios. We show this offset in the lower row of Figure 4. We conclude that single-age bursts of star formation as predicted by BPASS do not produce enough high energy photons to recreate the observed highly ionized gas.

Stellar Population Fits Models
The single-age bursts of star formation underpredict the high energy ionization lines as shown in Section 4.2.1. We have more knowledge of the stellar populations in J104457 and J141851 than most EELGs because we have observations of the FUV stellar continuum which we fit in Section 3. As a result, we can test whether the observed stars are capable of producing the high-energy ionization lines by using the ionizing spectrum from the stellar populations discussed in Section 3.2 as the input for our CLOUDY models. We fix our stellar parameters to those produced in Section 3.2 for these models but we vary the gas phase parameters to match the nebular properties of 104557 and J141851 published in Paper I. Just as was done for the fiducial models, both constant density and power-law density distributions were considered. The affects of a power-law density distribution are negligible for reproducing high ionization emission and thus we only present the results from our constant density models.
We also specialized our fits for J104457 and J141851 by adjusting the C/O and N/O ratios to the values reported in Paper I to see if these models would reproduce the highionization lines we observed. We set the metallicity of the gas to the metallicity measured in Paper I and then vary the metallicity with linear steps of 0.02 (2% Z ). The results of these models can be seen in the first (J104457) and third (J141851) columns of Figure 5. We plot the metallicity steps around the observed metallicity as colorful lines and the observed metallicity as the black line for each galaxy.
As we did with the models in Section 4.2.1, we fit the models to the observed line ratios from the 8 optical lines shown in Figure 4 and 2 UV lines. These results are shown in the Section 4.2.2 region of Table 5 and Table 6. Similar to the models from the previous section, the models that are specialized using the observed stellar population and gas phase abundances do a reasonable job of reproducing the low-and intermediate-ionization lines, but fail to reproduce the veryhigh ionization lines. The stellar populations we observed are able to match the high-ionization lines O III] and [Ar IV] that the single-age bursts could not match well, as seen by the lower χ 2 red values for the high ionization zone for Section 4.2.2 in comparison to Section 4.2.1 in Table 5 and Table 6. Although the stellar populations from Section 3.2 can recreate the low-, intermediate-, and high-ionization lines, these models are unable to produce the observed values for the very high-ionization energy lines like He II and [O IV].

Combination Stellar + Blackbody Models
The stellar population fit models from the previous section failed to reproduce the very high-ionization emission lines as seen in the He II and [O IV] panels of "Stars Only" columns of Figure 5. This problem is often seen for He II observations (e.g. Kehrig et al. 2015Kehrig et al. , 2018Kehrig et al. , 2021Senchyna et al. 2017Senchyna et al. , 2020Wofford et al. 2021), but we emphasize that the discrepancy between photoionization models using the stellar continuum as the input ionizing spectra and observations of emission lines is a problem for all very-high ionization lines, especially He II and [O IV], as seen in §4.2.1 and §4.2.2 above. Here we explore what is required to reconcile the observed nebular emission lines with the models. In the absence of new, specialized models of low metallicity stellar populations, we can provide a first bench mark to improve these models using the addition of a blackbody (BB).
We combine the models from Section 4.2.2 with a hot BB in order to reproduce the unexplained excess flux we observe in the very high-ionization lines. We added BBs with temperatures of 60,000 K, 80,000 K, and 100,000 K because these temperatures span the range between the energies dominated by X-ray emission and the energies observed in the UV. Senchyna et al. (2020) demonstrated that similar objects do not have significant X-ray emission so we avoid adding high temperature BBs or a power-law contribution that would be detected in X-rays.
We added these three BBs (T=60kK, 80kK, 100kK) in varying strengths by including BBs with 10% steps of the total fractional luminosity produced by the stellar population fit in Section 3. We refer to this quantity as the luminosity of the young stars (L young ). We add this BB to the input spectrum of the young stars and then vary the logU of the photoionization models using the combined spectra as the shape of the SED for the ionizing source.
We ran photoionization models with the same gas metallicity steps as those in Section 4.2.2; this includes the observed gas-phase metallicity from Paper I and the 2%Z steps around that measured metallicity. We used the abundances from the 4-zone ionization structure as estimated by Paper I to fix our C/O and N/O ratios.
We plot the effects of the BB on each metallicity model in the columns labeled "Stars + BB" of Figure 5. The lower bound (dashed line) of each metallicity shows the 60kK 10% BB while the upper bound (dotted line) plots the 100kK 80% BB so each band covers the entire range produced by our different BB possibilities. The black band on this plot shows the models run with the metallicity fixed at the observed metallicity.
We see that the added BB increases the range of line ratios each metallicity model is capable of producing. For the observed metallicity (black band), we see that similar to the models from the previous two sections, we are still able to reproduce the emission from the low-, intermediate-, and highionization lines, which demonstrates that the added BB does not greatly affect the lower energy parts of the ionizing spectrum. Additionally, the very high-ionization lines now are included within the band the models produce. This shows us that the additional BB produces enough ionizing photons at energies >54 eV to reproduce the observed He II emission.
The shape of the BB is determined by the temperature and percent contribution to the luminosity. We show the effects of the BB temperature on the models in Figure 6, each temperature is a different color and the shaded region spans from a lower bound of 10% BB to an upper bound of 80% BB. We plot the observed line ratios for J194457 in the top panels and J141851 in the bottom panels. The low-and intermediateionization lines have the ability to give us more information about the temperature of the BB. High ionization parameters with large line ratios can rule out certain BB tempera-

Stars Only
Stars Only Stars + BB Stars + BB J141851 J104457 Figure 5. Line ratio diagrams for 5 emission lines spanning the different ionization zones with the observed line ratio plotted as the point for J104457 and J141851. The "Stars Only" column shows the CLOUDY photoionization models where we used the ionizing continuum from the stellar population we fit to these galaxies in Section 3.2 as the input ionizing spectrum. Each line is a different metallicity from low metallicity at 1% Z in yellow to 11% Z in purple for J104457 and for J141851 the range starts with 3% Z in yellow and goes to 13% Z in purple. The black dashed line shows the results of a model run with the observed gas-phase metallicity of each galaxy, 5.8% Z for J104457 and 8.7% Z for J141851. In the "Stars + BB" column the models are plotted in bands with the same metallicity color-coding system, but the bands now include the effects of the blackbodies we modeled. The band is bordered by a dotted line showing the 0% Lyoung (no BB) model and a dashed-dotted line showing the 80% Lyoung from a 100kK blackbody. The black hatched region shows the models run with the observed gas-phase metallicity of the galaxy.

Model
Ionization Zone J104457 = 0.058 Z Section and (lines used) χ 2 logU Z (Z ) Age (Myr) tures. The high-ionization lines provide more information about the percent contribution. The shapes of the BB models allow for any temperature to produce the very-high emission lines. At high line ratios, higher temperature BBs can produce the same emission with a smaller percent contribution. To demonstrate this for the He II/Hβ plot we show dashed lines of each additional 10% in the percentage contribution to luminosity for the 100kK BB. This shows that with larger percent contributions the BB can reach increasingly extreme line ratios. The very-high ionization zone includes the most information about the BB component of the fit as the He II and [O IV] lines are most sensitive to the high energy photons. J104457 requires an 80kK BB with 60% of the young stellar light while J141851 requires an 80kK BB producing 70% of the young light. We can see for J104457 in Figure 6 all three BB temperatures are capable of reproducing the observed line ratios for our low-and intermediate-ionization lines. For J141851 in Figure 6, we see that the [Ar III] line starts to fall outside the band of line ratios produced by the 100kK BB, indicating a 100kK BB may not produce the emission we see from J141851. In J104457, the low and intermediate zone fits do not contain much information that constrains the BB components, but in J141851 these lines may contain some temperature information.
Our solutions employ a number of different selections of the 8 optical lines and 2 UV lines we've been using in this paper. The fit results for J104457 and J141851 are shown in Table 5 and Table 6 respectively. We report the temperature of the BB and the percentage of total light that the BB makes up in this table in addition to the columns discussed in Section 4.2.1 and Section 4.2.2. The χ 2 red values for the combined fit from all 10 emission lines are high because the combined fits do not fit all 10 lines well simultaneously, particularly [O III] increases our χ 2 red values. The combined fit does favor the BB parameters fit using only the very-high ionization lines which gives us confidence the BB information from all 10 lines is accurate, since the BB information matches the parameters found by the lines with the most information about the BB. When we compare the fits that include this BB to the fits from Section 4.2.1 and Section 4.2.2, we see that the added BB allows for a better fit when using lines from the very-high ionization zone and in the total 10 line fit of all ionization zones. This demonstrates that including a high energy photon source strongly improves the predicted nebular emission structure of very-high ionization emission line galaxies.

The High-Energy Ionizing Photon Production Problem
J104457 and J141851 are two EELGs that demonstrate the inability of traditional photoionization models to accurately reproduce the very high-energy ionization lines. We show this discrepancy in Section 4.2 with the very high-ionization lines of He II and [O IV]. In addition to the excess He II and [O IV] emission we see compared to the photoionization models without an additional source of high energy photons, J104457 and J141851 show this discrepancy in their [Fe V] emission. Paper I detailed the ionization correction factors (ICFs) for these two galaxies and found that the observed [Fe V] emission was not matched by the ICF provided from photoionization modeling with only the observed stars producing the photons. This inability of the stellar populations to reproduce the high-energy photons through photoionization indicates there is a high-energy ionizing photon production problem (HEIP 3 ) as named by Paper I.
The HEIP 3 is observed as excess He II emission in these two EELGs but also in the extreme He II emission observed in other blue compact dwarf galaxies (e.g., Kehrig et al. 2015Kehrig et al. , 2018Kehrig et al. , 2021Berg et al. 2018;Senchyna et al. 2017Senchyna et al. , 2020Stanway et al. 2020). Paper I showed that the HEIP 3 results in underestimated ionization parameters which, when corrected, could increase the He II emission produced by photoionization models used to compare with the extreme He II emission in blue compact dwarf galaxies. We demonstrated that even this increased ionization parameter could not account for the He II emission we observe in our two EELGs when only the stars are ionizing the gas, so it is unlikely that correcting the ionization parameter would reconcile the excess He II emission observed in other samples of blue compact dwarf galaxies or EELGs. We have found that the addition of a BB to the observationally-motivated model of the stars can solve the HEIP 3 for J104457 and J141851. This additional source of very-high energy photons can reproduce the very-high ionization lines that are not reproduced as a result of the HEIP 3 as is described in Section 5.2.

Self-Consistent Modelling Success
The emission lines from EELGs trace the ionizing spectrum of these galaxies, and, by using photoionization modeling, we are able to compare different forms of stellar population synthesis to the observed emission lines in order to understand what is ionizing these galaxies. We tested a number of input spectra to use in our CLOUDY photoionization modeling to reproduce the observed emission lines from these two galaxies. Using only the canonical models of single-age bursts of star formation, we were able to recreate the emission lines from low-energy ionization species, but we were unable to produce the high-energy ionization lines with these input as seen in Figure 4. Figure 4 compares emission lines from not only J104457 and J141851, but also a number of compact dwarf galaxies in the nearby universe. Notably, the models of high-ionization lines fall short of the observations from all galaxies we plotted. This indicates that models with a single-age burst of star formation do not produce the correct ionizing continuum shape required to reproduce the observations of EELGs. We then tested the observed stellar populations, as we measured them from Section 3, as the input spectrum. This moved us from using sets of likely stellar populations to the stellar population that we observed in each galaxy. These stellar populations were determined using models with and without binaries, as well as with an extended high-mass cutoff. The shapes of the ionizing continuum from the singlestar and binary models were similar, so we used the binary model version as it has better wavelength resolution over the wavelength range we used to calculate the ionizing spectrum. With this observationally-motivated ionizing spectrum, we were able to reproduce the observed emission lines up to ionization energies of ∼35 eV, like OIII], which can be seen in the "Stars Only" columns of Figure 5. The stellar model fits to the observed stars in these galaxies produce enough photons to recreate the high-energy emission lines but not the very high-energy emission lines like He II and [O IV]. This inability to reproduce the highest energy emission lines demonstrates the HEIP 3 as introduced in Paper I.
Since the models fit to the observed stars did not reproduce the highest energy emission lines, we decided to add a BB to the input spectrum in order to get an idea of the temperature and percent contribution of the source needed to recreate these lines. The ad hoc BB has a peak with a temperature less than what is observable in the X-rays because Senchyna et al. (2020) and others show that EELGs like our sample do not tend to have significant X-ray detections. The BB temperature must also be high enough to not disturb the low-energy ionization lines that are already produced by the observed stellar population. With these restrictions, we fit using 8 optical lines from all 4 ionization zones and found that the best fit temperature of the BB was 80,000 K and the flux contribution for the BB was 60-70% of the total luminosity of the young stars, L young . Full details of the fits can be seen in Table 5 and Table 6. In order to understand more about our BB, we fit the models to multiple smaller sets of emission lines. The low-and intermediate-ionization lines carry little information about the BB but do favor ionization parameters indicative of their ionization zones. In Figure 6, we see that the high-and very high-ionization lines hold most of the information about the BB and thus determine the BB parameters for the full 8 line fit of all ionization zones. This added BB successfully reproduces the low-and intermediateionization lines already produced with observed stars as well as the very-high energy lines that were missed.
The addition of a BB to understand the the ionizing radiation needed to produce the very-high ionization lines we observe in our EELGs is not a new idea. Steidel et al. (2014) used a similar BB to model the ionizing radiation of starforming galaxies at z ∼ 2 − 3. We compare our method to previous models by plotting our models of the observed stellar population and ad hoc BB at fixed metallicity of 0.058 Z for J104457 as a grid on a UV diagnostic diagram as seen in Figure 7. The UV BPT diagram has been suggested as a way to determine if a galaxy is star-forming or AGN driven for galaxies at higher redshift where the UV diagnostic lines are the only observed diagnostics (Feltre et al. 2016). We use a diagnostic diagram consisting of C IV λλ1548,50/He II λ1640 versus O III] λ1666/He II λ1640 to distinguish the star-forming and AGN driven objects. The general range of star-formation models produced in Mainali et al. (2017) is shown as the teal shaded box on the upper right corner of the plot, while the AGN driven models are shown in the shaded purple box on the lower left (Feltre et al. 2016). Our models are plotted as a grid where the observed stellar population only is the rightmost line and each vertical line to the left is a 20% jump in relative flux contribution of the BB. The horizontal lines are the logU of the model which start at logU=−1.0 at the top and each horizontal line is a step of 0.5 in logU to the bottom-most line of logU=−3.5. These models move across the gap from star-forming into the AGN driven regions of parameter space which covers the observations of these galaxies well.
The observation for J104457 lies above the model grid in Figure 7. This discrepancy between the models and the observations stems from the method for calculating the abundances based on the ionization correction factor (ICF). We used abundances calculated by Paper I who used continuous star-formation models to determine the ICF for carbon in our galaxies. Continuous star-formation models underpredict the ICF when compared to models with bursts of starformation, which leads to lower carbon abundances. With increased C/O values, our models would increase along the C IV λλ1548,50/He II λ1640 axis and thus encompass the observation of J104457, as indicated by the arrow in Figure 7. In similar diagnostic diagrams using C IV λλ1548,50/He II λ1640 versus C III] λ1666/He II λ1640 we see the same effects; our models do not produce enough C emission and by increasing the C/O abundance our best-fit model would match the observations for these galaxies.
Our combined stellar population and ad hoc BB models span the gap between AGN and star-formation driven models. The space between the classical AGN and star-formation models is primarily where we find EELGs like J104457 and J141851. Since local EELGs lie between AGN and star-formation models, observations of high-redshift galaxies likely fall in this parameter space as well. We caution that reality for high-redshift galaxies will not conform to only starformation or AGN models and going forward a composite spectrum similar to the ones developed here should be used.
In summary, the models of the observed stellar continuum are combined with an additional hot BB component in order to simultaneously fit the emission lines from all 4 ionization zones. This combination model reproduces the emission of He II and [O IV] from the very high-ionization zone as well as [O II] and [N II] from the low-ionization zone and lines from the zones in between. This observationally-motivated stellar continuum model and ad hoc BB are capable of reproducing the emission lines from all 10 emission lines used in our fit. Figure 7. The UV BPT diagram plotted with C IV λλ1548,50/He II λ1640 versus O III] λ1666/He II λ1640. The teal and purple shaded boxes highlight the regions typical of photoionization models for star-formation only and AGN only, respectively. In comparison, we plot the models for J104457, where the right-most line maps the models produced by the stellar population fits for a range of ionization parameters, extending from logU = −1 at the top to logU = −3.5 at the bottom with horizontal lines marking steps of 0.5 in logU. Models adding contributions from an 80kK blackbody are plotted as narrow lines, where each subsequent model to the left increases the contribution from the blackbody by 20%. The arrow indicates the direction the models should move to account for the discrepancy between abundances calculated from ICFs using continuous and burst star-formation models.

What is the source of the high energy photons?
Our BB and observationally-motivated stellar continuum combined model successfully recreates the high energy line ratios observed in our two EELGs. The ionizing continuum of the stars, BB, and combined stars and BB are shown in Figure 8. The green line shows the ionizing continuum from the best fit stellar population using BPASS models as described in Section 3 scaled to account for BB contribution. The black spectrum shows the observed data, while the teal line shows the BB, and the blue line shows the stars and BB combined. The circles show the intersection of each model with the wavelength required to produce He II and C IV emission and are color-coded by model.
The BBs shown in Figure 8 contribute a significant portion of the flux required to reproduce the high energy ionization lines and create a distinctly different shape in the far UV for these galaxies. The sources(s) of the high energy photons are, to date, unidentified. We have explored the literature and found a number of candidates that could cause the very-high ionization emission lines we observe, and we discuss them and their feasibility below. Binaries. In Section 3.2, we tested binaries in our stellar population calculation by using the BPASS models and found that binary star models recreate the observed UV continuum with similar precision but at slightly older ages. Even with the binary population models we cannot reproduce the BB emission; however binary models at lower metallicity provide more uncertainty in their models as a result of the unknown binary fractions as a function of metallicity, initial periods, and mass ratios which are measured from the local universe but not in low metallicity environments Stanway et al. (2020).
Low-metallicity stellar atmospheres. These galaxies are very metal poor. As reported in Berg et al. (2016) for J104457 and Berg et al. (2019b) for J141851, the gas in these galaxies is only 5.8% and 7.1% Z respectively and the stars themselves only 9-14% and 6% Z as reported in Table  2. The stellar population synthesis models for low metallicity stars are theoretical or extrapolated to lower metallicities since there are not observations of individual stars with Z < 0.1Z (Lejeune et al. 1997;Gräfener et al. 2002;Smith et al. 2002;Lanz & Hubeny 2003;Topping & Shull 2015). However, stellar atmospheres are complicated structures that non-linearly depend on the temperatures, densities, mass-loss rates, and chemical abundances. Slight variations in these parameters could produce dramatically different amounts of ionizing photons. Theoretical work on Population III stars indicates that they could provide more high-energy photons if they experienced significant mass-loss (Schaerer 2002). Observations of massive stars with Z < 0.1Z will provide insight as to whether low metallicity stars are capable of producing the emission we model as a BB in this work. Telford et al. (2021, accepted) shows that very low-metallicity stars are not capable of producing enough high-energy photons to reproduce extreme emission lines like He II.
Very massive stars. Another option to produce our BB radiation is very massive stars, stars that are >100M , which could be more plausible in the low metallicity regime. Very massive stars have been detected in SBS0335-052E by Wofford et al. (2021) who claim the models of very massive stars reproduce the UV and optical emission lines they have observed. Additionally, Senchyna et al. (2020) used a sample of 10 star-forming galaxies from the Sloan Digital Sky Survey to conclude that the discrepancy between stellar population models and the observed UV lines could only be rectified with additional very massive stars. When we tested both S99 and BPASS models with stellar populations up to 300M our ionizing continua did not sufficiently change. The ionizing continuum was larger for the 300M models, but was not sufficient in reproducing the ionizing photons required in our BB.
Stripped stars. Another hot candidate for the source of the high energy photons is stripped stars in binaries as suggested by Götberg et al. (2019Götberg et al. ( , 2020. These stripped cores of stars are extremely hot and produce significant He II ionizing photons, but they require an older stellar population. As shown in Table 2 the stellar populations found using the UV spectra are only 1 -10 Myr old. This does not preclude an older stellar population existing in these galaxies, but the young stellar population is clearly dominant in the emission lines we observe. WISE photometry of J104457 indicates the infrared SED is dominated by dust heated by the starburst and does not clearly show an older stellar population (Brammer et al. 2012).
Wolf-Rayet Stars. Another metallicity dependent candidate are Wolf-Rayet (WR) stars. WR stars drive their envelopes off of the stellar core using line-driven winds, but at low metallicities there are few metals to absorb the momentum required to accelerate the wind. In addition, we do not observe any WR bumps in our spectrum for J141851, but there is a possible C IV bump that could signify WRs present in J104457, as evidenced in the inset panels of Figure 2 as the emission lines from most of those species appear narrow and nebular in origin. The existence of some WRs in J104457 could help partially solve our missing high-energy photon problem, but with only one WR emission line signature, it is unlikely we have a significant enough population to fully solve our problem. Finally, even when WR stars are observed in metal poor galaxies, like IZw18 and SBS0335-052E, they do not produce enough ionizing photons at energies high enough to produce He II emission as shown in Kehrig et al. (2015Kehrig et al. ( , 2018. With all of these things considered, WRs are not likely the source of our missing high-energy photons. Active Galactic Nuclei (AGN). Since these high energy emission lines mimic those produced by AGN, a reasonable candidate for our BB is a low luminosity AGN. However, we ran a photoionization model with an AGN spectrum and found that this spectrum still overproduces the very-high ionization lines when scaled so that the low-, intermediate-, and high-ionization lines match the observations, implying the spectrum is too hard to reproduce our suite of emission lines. The hardness of the AGN spectrum also implies the existence of hot gas that could be detected in the X-rays, but these dwarf galaxies are rarely observed with X-ray detections (Senchyna et al. 2017). The stellar masses of these galaxies are only 10 6.5−7 M which according to relations of the host mass to central black hole mass like the one produced in Reines & Volonteri (2015) would imply an intermediate mass black hole to power the AGN. Intermediate mass black holes are notoriously difficult to observe; however Xray surveys including the Chandra COSMOS-Legacy survey have allowed studies to find AGN in dwarf galaxies out to z ≤ 2.4 (Mezcua et al. 2018). Similar radio surveys have observed compact radio sources consistent with AGN activity in dwarf galaxies (Reines et al. 2020), and through resolved emission line studies from the optical (Moran et al. 2014;Mezcua & Domínguez Sánchez 2020). Mezcua et al. (2018) suggest that the AGN fraction decreases to ∼ 0.4 percent for dwarf galaxies with stellar masses in the range of 10 9 < M * ≤ 3 × 10 9 M with respect to both X-ray luminosity and stellar mass, which makes our ∼ 10 6.7 M galaxies unlikely to have AGN driving their emission lines.
Ultra-Luminous X-ray Sources (ULXs). Black holes could still cause the He II emission we see if they are in the form of high-mass X-ray binaries or ULXs where a black hole or neutron star is in a binary with an O or B type star and accretes material from the stellar wind of the companion. These accreting compact objects could produce the photons necessary to excite He II as suggested in Schaerer et al. (2019); however the spectral shape of high-mass X-ray binaries is debated. Senchyna et al. (2020) finds the ULX spectrum to be too hard for most small metal-poor galaxies as these galaxies do not tend to show X-ray detections in the high-energy bands which are usually produced by ULXs. Recently, however, Simmonds et al. (2021) performed a similar analysis using ULX spectra based on observed ULXs in nearby galaxies and found that the softer spectra are capable of reproducing He II emission, but they include the caveat that ULX spectra as a whole are not well understood. In particular, Schaerer et al. (2019) proposed the He II emission from IZw18 was produced by X-ray binaries, however, when Kehrig et al. (2021) explored the feasibility of this claim, they found a separation of ∼200 pc between the high-mass X-ray binary and the He II emission, again suggesting the X-ray binary is not responsible for the He II emission observed.
Supersoft X-ray Sources (SSSs). As ULXs are too hard, perhaps the younger sibling of compact object accretion, supersoft X-ray sources, where a white dwarf accretes material, could be the physical mechanism behind the high energy photons we require. SSSs have temperatures of ∼10-100 eV and bolometric luminosities of ∼ 10 36 − 10 38 erg s −1 and so, could possibly provide the very-high energy photons we need (Greiner et al. 1991;Kahabka et al. 1994). SSSs have not been observed in similar galaxies to-date; however the spectrum of these accreting WDs should fall into a similar parameter space as our BB. SSSs are rarely observed in nebulae similar to those found in J104457 and J141851 (Woods & Gilfanov 2016), and, as they are often viewed as the progenitors for Type Ia SNe, it is important to note that Type Ias are rarely found in low-metallicity environments (Kobayashi et al. 1998). Additionally, the presence of WDs implies an older stellar population than we observe in the UV, and as mentioned with the stripped stars, the IR infor-mation from WISE does not provide enough information to conclude whether there is an older stellar population residing in these galaxies. Even with an older stellar population we would need ∼50,000 and ∼10,000 SSSs for J104457 and J141851 respectively to produce the luminosity of our BB (J104457: 5.1 × 10 42 ergs −1 ; J141851: 1.1 × 10 42 ergs −1 ). The spatial distribution of these SSSs would be random in the location where the older stellar population formed, so a spatial analysis of the He II emission could help determine if SSSs provide the very-high energy photons we require.
Shocks. Shocks are another suggestion to produce the high-energy photons we need. Similar galaxies, IZw18 and SBS0335-052E, have soft X-ray emission that could correspond to shocks, but when the soft X-ray emission is integrated, it does not account for enough energy to produce the He II observed in those galaxies (Kehrig et al. 2015(Kehrig et al. , 2018(Kehrig et al. , 2021. Shocks should also produce emission lines from species like [Ne V] up to a line ratio of [Ne V] λ3426/Hβ < 0.3 for a shock propagating through a medium with density N e = 1 cm −3 (Izotov et al. 2012 Izotov et al. (2012). Additionally, Berg et al. (2018) tested models with pure shock emission to explain the extreme He II emission from the lensed star-forming galaxy SL2S J021737-051329 at z ∼ 2 and found that shocks generally produced too little C III] emission in these EELGs.
Catastrophic Cooling. Finally, suppressed super winds causing catastrophic cooling could cause very-high energy ionization line emission as suggested by Gray et al. (2019) and Danehkar et al. (2021). In this case, the rapid gas cooling through lines like O VI creates extreme emission beyond the capabilities of the stellar population. Gray et al. (2019) predict significant He II emission at the outer boundary of their nebula and Danehkar et al. (2021) predict O VI emission from the central parts of the nebula. Future work with spatially resolved He II emission could help determine if these confined winds are the source of the additional high-energy photons we observe. Additionally, future observations of these galaxies should look for emission lines indicative of this catastrophic cooling such as Si III λ1206 and O VI λ1037. We examined existing low resolution spectra for O VI for J104457 and J141851 and find that O VI is not significantly detected. This is an incomplete list of possible candidates that could generate the high energy photons we infer. Future work can help determine which mechanism provides the very-high energy photons. Spatial resolution of the He II emission and observations of the X-rays produced by J104457 and J141851 will help rule out some of the options provided in this list. We hope that by providing a BB with temperature and percentage of the young stellar population L young we are laying down a foundation for building a physical model of the source of the high-energy photons.

CONCLUSIONS
We used the HST/COS and LBT/MODS spectra presented in Paper I to study the stellar populations and extreme emission lines from two extreme emission line galaxies (EELGs), J104457 and J141851. These EELGs have the largest reported C IV and He II EWs in the local universe which suggests these are likely reionization analogue galaxies. We fit the stellar populations for these galaxies using the UV continuum from the HST/COS spectrum and find these are very young stars with ages ≤ 10 Myr and metal-poor with metallicities ≤0.15 Z . We used photoionization modeling to check if (a) single-age bursts of star formation, (b) the stellar populations from the UV continuum fits, and (c) an additional BB component combined with the observationallymotivated stellar populations were capable of reproducing the emission lines observed in these EELGs.
Our specific results from this work are: 1. The stellar populations that dominate the continuum in UV spectra are extremely young and metal-poor. We used the methods set up in Chisholm et al. (2019) to measure the age and metallicity of the stellar populations in J104457 and J141851 and found ages of 4.01 ± 1.13 Myr and 3.37 ± 1.24 Myr respectively. The metallicities were measured at 0.09±0.05 Z and 0.06±0.05 Z for J104457 and J141851 respectively. We used both single-star and binary-star models and models with high-mass cutoffs at 300M in addition to models with high-mass cutoffs at 100M . The use of very massive star models did increase the ionizing photon production rate but did not account for the increase needed to produce the extreme He II we observe.
2. Absorption features in the UV continuum also indicate young stellar populations. We found Fe V and S V absorption features in the UV continuum for both J104457 and J141851. These features are indicative of young stellar populations in higher metallicity galaxies, but their observation in these low metallicity stellar populations confirms the presence of very hot stellar atmospheres which require young massive stars.
3. We used photoionization modeling with CLOUDY to test whether single-age bursts of star-formation could reproduce the observed emission lines from our EELGs. We find that BPASS single-age burst models match the emission lines from species in the lowand intermediate-ionization zones, but do not produce enough high-ionization photons to recreate the high-and very-high ionization lines, including O III], [Ar IV], [Ne III], He II, and O IV. These models are frequently used to model the stellar populations found in EELGs, but they are insufficient to reproduce the highionization emission lines from EELGs. 4. We used the stellar populations fit in Section 3 as inputs for our photoionization modelling with CLOUDY in an attempt to reproduce all emission lines with the stars observed in the galaxy. Although our observationally-motivated stellar models were able to reproduce the low-, intermediate-, and high-ionization lines, the very-high ionization lines like [Ne III], He II, and O IV were still out of reach for these models. The inability of the observed stars to reproduce the observed emission indicates there is an additional source of very-high energy ionizing photons produced in these galaxies.
5. We added our observationally motivated stellar population from Section 3 to a BB of temperature 80kK that contributes 60% and 70% of the light from the young stars for J104457 and J141851 respectively. This ad hoc BB produces the very-high energy photons required to reproduce the very-high ionization lines like [Ne III], He II, and O IV. The combination of BB and observed stellar population models is capable of reproducing the entire suite of emission lines, from lowionization [O II] up to very-high ionization O IV.
6. The high energy ionizing photon production problem (HEIP 3 ) is often seen in EELGs where the extreme emission lines like He II and O IV cannot be reproduced by the stellar populations expected to reside there or by other combinations of physical phenomena like collections of WR stars or ULXs. The ad hoc BB we used to model the high-energy photons needed to produce the observed He II solves the HEIP 3 for the two EELGs studied in this paper by simultaneously reproducing the suite of emission lines across all ionization zones from low to very-high.
7. The BB we add to our observationally-motivated stellar models could come from a number of physical sources, but we are unsure which sources would produce this set of ionizing photons. Some of the previously suggested solutions from the literature include shocks, WR stars, very massive stars, and ultra luminous X-ray sources. We address these previous physical suggestions in Section 5.3, but none of the possible solutions stand out as a clear explanation of our highenergy photons. We require further observations at additional wavelengths to test more of these suggestions.
We have listed numerous possibilities of physical sources that could produce the high-energy photons needed to ex-plain the very-high ionization lines in J104457 and J141851. Further exploration of these galaxies in the X-ray could help limit the options by exploring the range of photons at highenergies. Additional information from spatially resolved He II emission would provide information that could help determine if the source of the ionizing photons is centrally located or if it exhibits a larger spatial distribution. Further UV spectroscopy could be useful in searching for O VI signatures related to some of the possible suggestions. We leave these explorations of J104457 and J141851 to future work.
GMO is thankful for the support for this program HST-GO-15465, that was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. DKE is supported by the US National Science Foundation (NSF) through Astronomy & Astrophysics grant AST-1909198. This paper made use of the modsIDL spectral data reduction reduction pipeline developed in part with funds provided by NSF Grant AST-1108693 and a generous gift from OSU Astronomy alumnus David G. Price through the Price Fellowship in Astronomical Instrumentation.