# Resolving the Mechanical and Radiative Feedback in J1044+0353 with KCWI
Spectral Mapping

Crystal L. Martin Department of Physics
University of California Santa Barbara
Santa Barbara, CA 93106, USA  Zixuan Peng  Department of Physics
University of California Santa Barbara
Santa Barbara, CA 93106, USA  Yuan Li  Department of Physics
University of California Santa Barbara
Santa Barbara, CA 93106, USA

(Received November 22, 2023; Revised March 12, 2024; Accepted March 14, 2024)

###### Abstract

We present integral field spectroscopy toward and around J1044+0353, a rapidly
growing, low-metallicity galaxy which produces extreme [O III] line emission.
A new map of the O32 flux ratio reveals a density-bounded ionization cone
emerging from the starburst. The interaction of the hydrogen ionizing
radiation, produced by the very young starburst, with a cavity previously
carved out by a galactic outflow, whose apex lies well outside the starburst
region, determines the pathway for global Lyman continuum (LyC) escape. In the
region within a few hundred parsecs of the young starburst, we demonstrate
that superbubble breakthrough and blowout contribute distinct components to
the [O III] line profile, broad and very-broad emission line wings,
respectively. We draw attention to the large [O III] luminosity of the broad
component and argue that this emission comes from photoionized, superbubble
shells rather than a galactic wind as is often assumed. The spatially resolved
He II $\lambda 4686$ nebula appears to be photoionized by young star clusters.
Stellar wind emission from these stars is likely the source of line wings
detected on the He II line profile. This broader He II component indicates
slow stellar winds, consistent with an increase in stellar rotation (and a
decrease in effective escape speed) at the metallicity of J1044+0353. At least
in J1044+0353, the recent star formation history plays a critical role in
generating a global pathway for LyC escape, and the anisotropic escape would
likely be missed by direct observations of the LyC.

††journal: ApJ††facilities: WMKO(KCWI), HST(ACS,WFC3)††software:  Astropy
(Astropy Collaboration et al., 2013), BPASS (Eldridge et al., 2017), Cloudy
(Ferland et al., 2013), CWITools (O’Sullivan & Chen, 2020), the Infrared
Science Archive (https://irsa.ipac.caltech.edu/frontpage/), KDERP
(https://github.com/Keck-DataReductionPipelines/KcwiDRP), LMFIT (Newville et
al., 2021), MVT-binning (https://github.com/pierrethx/MVT-binning ) PyNeb
(Luridiana et al., 2015), and SExtractor (Bertin & Arnouts, 1996)

##  1 Introduction

Most of the hydrogen in the universe underwent a phase change, from neutral to
ionized, during the first billion years of galaxy assembly. Stars more massive
than roughly 20${\rm~{}M_{\odot}}$ produce Lyman continuum (LyC) photons at a
significant rate, so LyC leakage from star-forming galaxies is widely believed
to be the dominant source of the radiation that ionized the intergalactic gas
(Bouwens et al., 2015a; Finkelstein et al., 2019). How this LyC radiation
leaked out of galaxies, however, is still not well understood. Integral field
spectroscopy of the nearest analogs of the galaxies that likely drove cosmic
reionization can provide critical insight into how reionization proceeded.

At redshifts where the LyC can be directly observed, neither galaxy
metallicity nor stellar mass provide a reliable prediction of the LyC escape
fraction, $f_{esc}$. Two of the best indicators of high $f_{esc}$ are a high
star formation rate (SFR) surface density or a high nebular [O III] to [O II]
flux ratio (Flury et al., 2022). High [O III] / [O II] ratios, for example,
are consistent with density-bounded channels, where the column density of
neutral hydrogen is insufficient to absorb all the LyC radiation (Jaskot &
Oey, 2013; Izotov et al., 2018; Plat et al., 2019; Ramambason et al., 2020).
Similarly, supernova explosions may be more likely to open pathways for LyC
escape from starbursts with high SFR surface density since their measured
outflow speeds increase with SFR surface density (Kornei et al., 2012; Heckman
& Borthakur, 2016). Yet efforts to demonstrate correlations between these
quantities and $f_{esc}$ yield confusing results: (1) outflow velocities,
measured from the Doppler shifts of interstellar absorption lines, have not
shown a correlation with $f_{esc}$ (Henry et al., 2015; Jaskot et al., 2017),
and (2) density-bounded channels are not a unique explanation for high [O III]
/ [O II] ratios (Plat et al., 2019; Flury et al., 2022). The resulting tension
has been cited as evidence that winds might become ineffective at the highest
SFR surface densities due to catastrophic radiative cooling (Gray et al.,
2019; Danehkar et al., 2021). It is therefore important to understand whether
line-of-sight effects caused by anisotropic feedback might offer an
alternative explanation for these apparent contradictions.

We observed J1044+0353 with multiple pointings of the Keck Cosmic Web Imager
(KCWI), an integral field spectrograph on the Keck II telescope. The
properties of our target have many similarities to reionization-era galaxies
and are listed in Table 1. The [O III] equivalent width is 1400 Å in the
$\lambda 5007$ line alone, an extreme value at any redshift and an indication
of efficient ionizing photon production (Chevallard et al., 2018). The high
specific star formation rate, $4.1\times 10^{-8}$ yr-1 (Berg et al., 2022),
underlines this small galaxy’s rapid, recent growth. The far-ultraviolet (FUV)
spectrum of J1044+0353 shows strong, high-ionization emission lines, including
He II $\lambda 1640$ (Berg et al., 2016, 2022). Recombination of He+2 produces
He II $\lambda 1640$ and optical He II $\lambda 4686$ emission but requires
photons with energies $\geq\ 4$ Rydberg. Normal stellar populations, as
described by Binary Population and Spectral Synthesis v2.1 (BPASS) models
(Eldridge et al., 2017) for example, do not produce enough 54.4 eV photons to
explain the He II $\lambda 1640$ emission from J1044+0353 (Berg et al., 2019,
2021; Olivier et al., 2022). This high-energy ionizing photon production
problem has been observed in other local galaxies (Shirazi & Brinchmann, 2012;
Jaskot & Oey, 2013; Stasińska et al., 2015; Senchyna et al., 2017; Schaerer et
al., 2019; Mingozzi et al., 2022, 2023) and is a common feature reionization-
era galaxies (Stark et al., 2015; Mainali et al., 2017; Bunker et al., 2023;
Senchyna et al., 2023; Topping et al., 2024). The damped Ly$\alpha$ profile in
the same FUV spectrum of J1044+0353 (Hu et al., 2023), however, indicates a
low $f_{esc}$ value along our sightline. In known LyC leakers, the intensity
of He II relative to H$\beta$ appears to be driven by metallicity rather than
$f_{esc}$ (Marques-Chaves et al., 2022). J1044+0353 has a gas-phase oxygen
abundance of $12+\log(O/H)=7.45$ (Peng et al., 2023) but lacks the strong
nitrogen enhancement seen in GN-z11 (Senchyna et al., 2023). Massive stars of
similar metallicity are expected to have a harder ionizing spectrum on
empirical grounds; see Figures 2-4 of Sander (2022) for example. Green Pea
galaxies are a well-studied sample selected based on their high [O III]
equivalent widths (Cardamone et al., 2009), but their redshifts,
$0.112<z<0.360$, place them over ten times further away than J1044+0353.

Figure 1 illustrates the information gained from spatially resolving the
starburst environment. The compact starburst is off center in J1044+0353. In a
preliminary analysis of our KCWI data cube (Peng et al., 2023), we mapped the
star formation history across the surrounding galaxy and discovered a post-
starburst population roughly 1.3 kpc east of the starburst. Whereas the
interstellar absorption lines in the FUV spectrum had indicated a very weak
outflow (Xu et al., 2022), the global gas kinematics revealed a large, bipolar
outflow (Peng et al., 2023). Initially, it came as a surprise that the apex of
this galactic wind was outside the FUV spectral aperture, well separated from
the compact starburst circled in Figure 1. However, the post-starburst
population, where the stars are at least 10-20 Myr old, has produced more
supernovae than the compact starburst, which is less than 4 Myr old (Olivier
et al., 2022). This spatial mapping of the star formation history and gas
kinematics can help us understand LyC escape.

Figure 1:  Composite Hubble image of J1044+0353 constructed from 8 filters:
F125LP, F150LP, F165LP, F336W, F438W, F606W, F665N, and F814W. Using the pivot
wavelengths, filters were mapped into the optical range. The stretching method
is adapted from Lupton et al. (2004) but generalized to more than three
filters. The red circle of diameter 2$\farcs$5 ($\approx 350$ pc) identifies
the Hubble COS aperture used to study the starburst previously (Berg et al.,
2022; Xu et al., 2022; Hu et al., 2023). The three slicer pointings (white)
outline the footprint of the KCWI mosaic cube. We follow the naming convention
from Figure 1 of Peng et al. (2023): SB for the compact starburst at
10:44:57.797, +3:53:13.16, FEN for the group of star clusters 1.24 kpc to the
east 10:44:58.092, +03:53:12.58, and MEN and NEN for the complexes between
them at 10:44:57.964, +03:53:12.40 and 10:44:57.986, +03:53:13.81. The MEN and
NEN regions are resolved into reddish-colored clusters strung out along a
north-south direction between SB and FEN. The Balmer absorption is strong
across MEN, NEN, and FEN, so we refer to them as the post-starburst region.
The red spiral galaxy visible south of J1044+0353 at 10:44:57.874,
+03:53:09.39 has an emission-line redshift of $z=0.27522$.

A single, short burst of star formation fails to produce much LyC escape. The
main sequence lifetime of the lowest mass star undergoing core-collapse,
$\approx 8-10$ ${\rm~{}M_{\odot}}$, is roughly an order of magnitude longer
than the lifetime of the lowest mass star producing significant LyC radiation
(Leitherer et al., 1999). Considering the steep slope of the stellar initial
mass function (IMF), a short burst stops producing LyC radiation before most
of the supernova feedback has taken place, a situation we call the supernova
timing problem. Cosmological, zoom-in simulations produce channels for Lyman
escape several ways, including runaway OB stars (Kimm & Cen, 2014), radiative
feedback (Kakiichi & Gronke, 2021; Chevance et al., 2022), and frequent galaxy
mergers in the reionization era (Martin-Alvarez et al., 2023; Witten et al.,
2024). After reionization, the environment of dwarf galaxies evolves, yet
feedback continues to play a major role in their evolution (Agertz et al.,
2020; Rey et al., 2022). Most local, dwarf galaxies are very diffuse, inactive
systems (Weisz et al., 2012; Guo et al., 2016), so the high specific SFR in
J1044+0353 indicates that some mechanism recently drove gas inwards. The low
gas-phase metallicity across J1044+0353, relative to the mass - metallicity
relation, provides indirect evidence for this recent gas inflow (Peng et al.,
2023).

In this paper we map the ionization structure across the field shown in Figure
1 in order to better understand whether the supernova feedback produced by
J1044+0353’s recent star formation history creates low-density pathways
through the interstellar medium. To gain additional insight into the origin of
the hard ionizing radiation, we compare the structure of the He II nebula in
J1044+0353 to the distribution of young star clusters visible in Hubble images
and the gas kinematics in the starburst region. Our presentation is organized
as follows. In Section 2, we fully describe the observations and data
reduction, including the generation of the KCWI mosaic in Figure 2.1 and the
archival Hubble imaging in 2.2. Section 3 presents measurements of the
hydrogen and helium recombination lines; we resolve the He+2 nebula and
measure total recombination-line luminosities. Section 4 maps out spatial
variations in the gas density and ionization structure, and we illustrate the
LyC escape pathways using the [O III] / [O II] flux ratio. Section 5 describes
the mechanical feedback in the starburst region and identifies two sources of
high-velocity, emission-line luminosity. Section 6 discusses whether stellar
populations modeled with a standard IMF can match the measured ionizing photon
efficiency and inferred supernova power. Section 7 summarizes the results and
their implications for stellar feedback and $f_{esc}$ in the reionization era.

Throughout this paper, we correct for internal reddening using the ${\rm
A_{V}}=0.20\pm 0.10$, measured (from the Balmer decrement) in the starburst
region (Peng et al., 2023), and an SMC extinction curve (Gordon et al., 2003).
We correct fluxes for Galactic extinction (Fitzpatrick, 1999) using the the
dust map of Schlafly & Finkbeiner (2011a), which we accessed via the IRSA
interface111 https://irsa.ipac.caltech.edu/cgi-bin/bgTools/nph-bgExec. All
magnitudes are in the AB system, and all emission-line equivalent widths are
quoted in the rest frame. We adopt the total recombination rate coefficients
from Storey & Hummer (1995) and the effective recombination rate coefficients
from Pequignot et al. (1991), both at a fiducial electron temperature of
$2.0\times 10^{4}$ K (Peng et al., 2023). The optical emission lines in the
integrated KCWI spectrum indicate a spectroscopic redshift $z_{spec}=0.012873$
, consistent with SDSS and COS spectra (Berg et al., 2022). However, a model
of the local velocity field places J1044+0353 further away, at an effective
redshift ${\sc ZDIST}=0.01363\pm 0.00068$ (Willick et al., 1997). We adopt a
flat $\Lambda{\rm CDM}$ cosmology with $\Omega_{\Lambda}=0.7$,
$\Omega_{m}=0.3$, and ${\rm H_{0}}$= 70 km s-1 Mpc-1. The luminosity distance
is then 59.0 Mpc, and the angular scale is 278 pc/″.

##  2 Data

A decade before the launch of JWST, hundreds of reionization-era galaxies were
identified by their high [O III] equivalent width, which produces a
recognizable broad-band infrared color (Schaerer & de Barros, 2010; Labbé et
al., 2013; Roberts-Borsani et al., 2016; Smit et al., 2015). JWST spectroscopy
has now confirmed this rapid increase in [O III] emission line strength with
redshift (Matthee et al., 2023; Oesch et al., 2023).

We selected our target, J1044+0353, from the NASA-Sloan Atlas (NSA). The NSA
v0_1.2 catalog includes nearly 27,000 dwarf galaxies, defined here as galaxies
fainter than $\mbox{${\rm~{}M_{r}}$}=-18$ and with NSA stellar masses less
than $10^{9}$${\rm~{}M_{\odot}}$. We searched this sample for rest-frame [O
III] $\lambda 5007$ equivalent width exceeding 1000 Å because strong [O III]
emission is a common spectral feature of high-redshift galaxies (Labbé et al.,
2013; Smit et al., 2014, 2015; Stark et al., 2015; Roberts-Borsani et al.,
2016; Stark et al., 2017; Laporte et al., 2017; Hashimoto et al., 2018; Berg
et al., 2018; Tamura et al., 2019). We found only 35 dwarf galaxies with ${\rm
W}(\mbox{[O III]}\ \lambda 5007)>1000$ Å, hereafter the O3EWGT1000 sample.
Dwarf galaxies with extreme [O III] equivalent width are therefore rare at low
redshift, but the O3EWGT1000 sample includes a high fraction of narrow-line He
II-emitters, 17 of the 35 galaxies including J1044+0353. Considering that only
71 of the 27,000 NSA dwarf galaxy spectra are listed as He II emitters
(Shirazi & Brinchmann, 2012), galaxies in the O3EWGT1000 sample are 200 times
more likely than average to emit nebular He II $\lambda 4686$.

The presence of He II recombination lines requires an ionizing spectrum with a
significant photon luminosity at an energy of 54.4 eV. We are observing this
galaxy population with integral field spectroscopy in order to gain new
insight about the source of the hard radiation, the escape channels for
ionzing radiation, and the mechanical feedback. We present a case study of
J1044+0353 because our observations spatially resolve the very-high ionization
zone, and Hubble imaging provides new insight the number of ionzing sources.

###  2.1 Integral Field Spectroscopy with KCWI

The Keck Cosmic Web Imager (KCWI) is the blue channel of an integral field
spectrograph on the Keck II telescope (Morrissey et al., 2018). We observed
J1044+0353 (RA=10:44:57.80, DEC=+03:53:13.2) with KCWI on 2018 January 14
under clear skies and average seeing of 0$\farcs$9\. These observations were
made before commissioning was fully completed, and keywords in the FITS
headers needed to be manually synchronized with the observer’s log.

We configured KCWI with the small slicer, the BL grating, and a camera angle
that produced a central wavelength of 4600 Å. In this configuration, the slits
are 0$\farcs$35 wide, and the 24 slices subtend an 8$\farcs$4 $\times$
20$\farcs$4 field of view. The dichroic limited spectral coverage to
wavelengths less than 5600 Å. Falling sensitivity shortward of 3700 Å produced
a steep decline in sensitivity toward the atmospheric cutoff. The combination
of narrow slits plus the BL grating provided an average resolution $R\approx
3600$ over a bandpass which includes: (1) recombination lines from He II, He
I, and the hydrogen Balmer series, and (2) collisionally excited metal lines
from ions spanning a wide range of ionization states.

We obtained exposures of 1200 s and 300 s at three KCWI pointings. As
illustrated in Figure 1, the slit position angle was 90 degrees east of north,
and offsets were made in the direction perpendicular to the slits. The first
exposure was centered on the compact starburst, and we then offset roughly
half the slicer width north and south, thereby covering a 20$\farcs$4 by 16″
field of view. We made the offsets in half-slice multiples to improve the
point spread function (PSF) sampling. This strategy partially offsets the
coarse sampling of the slits relative to the finer pixel sampling along the
slit, 0$\farcs$146 pixel-1, providing nearly uniform sampling in both spatial
dimensions. To measure the sky spectrum, the slicer was chopped to a blank
field, and 600 s sky exposures were made between pairs of 1200 s and 300 s on
target exposures. This sequence provided 4500 s of integration time over the
central half of the final mosaic. These exposures detected the bright emission
from [O III] $\lambda\lambda 4959,5007$ over the entire field of view.

####  2.1.1 Cube Construction

We used the KCWI Data Extraction and Reduction Pipeline v1.1.0 (Morrissey et
al., 2018, KDERP) to reduce the raw frames and build calibrated data cubes. We
first produced bias subtracted and gain-corrected intensity frames with cosmic
rays and scattered light subtracted. We then mapped each detector pixel to a
slice, position, and air wavelength. We extracted the sky spectrum from the
paired sky observations, modeled the 2D sky spectrum for each integration, and
subtracted the sky counts prior to cube generation. Each cube was corrected
for differential atmospheric refraction based on the airmass of the
observation, the position angle of the slicer, and the parallactic angle. Each
stage of the pipeline produced a matching variance image or cube that
describes the detector noise and the shot noise from photon counting
statistics.

We observed the spectrophometric standard star Feige 34 with the same
instrument configuration during the observing run. We extracted the Feige 34
stellar spectrum, subtracted the sky spectrum recorded on the star-free region
of the slicer, and corrected the spectrum for atmospheric attenuation (Buton
et al., 2013) according to the airmass of the observation. We computed the
inverse sensitivity in narrow bands that excluded strong stellar absorption
lines. We fit the inverse sensitivity with a high-order polynomial across the
$\lambda 3650$ to $\lambda 5625$ bandpass. The extracted, flux-calibrated
spectrum of Feige 34 is consistent with its reference spectrum to within $\pm
3\%$ across this bandpass. We applied this inverse sensitivity function to
each cube. The measurements reported in this paper use this absolute flux
calibration.222 The steepening of the inverse sensitivity function shortward
of $\lambda 3650$ was not easily fit, as it tended to introduce high-order
wiggles at longer wavelengths. Inspection of the spectrum around [Ne V]
$\lambda 3425$ did not reveal an emission line, so we did not flux calibrate
the shortest wavelengths.

####  2.1.2 Construction of the Mosaic Cube

After cropping the KCWI cubes to the valid data region, we registered the
world coordinates of the compact starburst with its position in the SDSS
g-band image. For the offset frames that only covered half of the compact
starburst, we used the telescope offsets to register them relative to the base
pointing. Following O’Sullivan & Chen (2020), we coadded the individual cubes
using the CWITools Python package333https://github.com/dbosul/cwitools. This
procedure produces a mosaic cube with 0$\farcs$146 $\times$ 0$\farcs$146
$\times$ 0.5 Å voxels.

A similar procedure was applied to the variance frames using the -vardata True
flag. We compared the global noise properties of the intensity cube to the
variance cube and found the former was typically higher. We rescaled our
variance cubes by their ratio, a factor $\approx 1.5$. The resampling to a
regular grid introduces covariance among adjacent voxels of the data cube. We
ignore this covariance in our analysis because its effect on spatial scales
larger than the point-spread-function is negligible.

A few spaxels (out of over 15,000 in the entire mosaic) have saturated [O III]
lines. We identified saturated pixels with more than 60,000 counts in the raw
2D frames. In the 1200 s exposures, saturation clips the cores of the [O III]
$\lambda\lambda 4959,5007$ and H$\beta$ emission lines in slits that intersect
the spatially unresolved core of the starburst. No slits show saturated [O
III] $\lambda 4959$ or H$\beta$ emission in the short 300 s exposures.
Whenever a measurement requires integration over a masked spaxel, we simply
measure the H$\beta$ and [O III] $\lambda 4959$ line flux from the matched
mosaic of 300 s cubes. Measurements of the [O III] $\lambda 5007$ line in the
300 s mosaic are assigned a flux equal to three times the [O III] $\lambda
4959$ flux, consistent with the respective radiative decay rates of the two
transitions.

####  2.1.3 Spectral Extraction

The starburst spectrum shown in Figure 2 was extracted by summing 1-D spectra
within 10 spaxels of the brightest spaxel. The radius of this roughly circular
aperture is 1$\farcs$53, similar to the size of an SDSS fiber. We confirm the
strong emission lines from low, medium, high, and very-high ionization line
described previously by (Berg et al., 2021). As an external check on the flux
calibration, we compared the fluxed KCWI and SDSS spectra. We extracted
spectra at all spaxels within 1$\farcs$53 of the SDSS_PLUG position and added
them to provide a KCWI spectrum of the starburst. We corrected these spectra
for foreground Galactic reddening using the visual extinction from Schlafly &
Finkbeiner (2011b) and the Fitzpatrick (1999) average extinction curve. To
register the wavelength scales, we corrected the KCWI spectra to the
heliocentric reference frame, and we converted the vacuum SDSS wavelengths to
air. We found that the SDSS spectrum downloaded from the SkyServer had a
higher overall flux density than our aperture-matched KCWI spectrum. Further
investigation revealed that the absolute fluxes of SDSS spectra are rather
uncertain, and the spectra should be scaled to the fiber magnitudes following
Brinchmann et al. (2004). We multiplied the SDSS spectrum by the g, r, and i
throughput curves, and averaged the three scale factors to obtain a
multiplicative correction of 0.74. This produced a good overall agreement in
flux density between the SDSS and KCWI spectra at the red end of our bandpass.
Tremonti et al. (2004) reported that SDSS spectra are systematically 12% bluer
than SDSS photometry, and our comparison confirms a discrepancy of similar
magnitude in the part of the u-band covered by the SDSS spectrum. Therefore,
the external comparison to SDSS did not add information beyond the internal
calibration.

Figure 2: Example rest-frame KCWI spectrum extracted from a 3$\farcs$0
diameter aperture centered on the compact starburst. Note the prominent
emission lines from the very-high ionization zone: He II $\lambda 4686$, [Fe
V] $\lambda 4227.4$, [Fe IV] $\lambda 4906.6$, and [Ar IV] $\lambda\lambda
4711.35,4740.20$. The spectral resolution separates the [O II] $\lambda\lambda
3726.03,28.82$ doublet lines. Broad, low intensity wings are visible on the [O
III] and H$\beta$ lines.

We measured W([O III]) across a set of circular apertures centered on the
starburst. In a circular aperture of radius 1$\farcs$5, we measure a rest-
frame W($\lambda 5007$) of nearly 1400 Å. The measured EW grows rapdily with
increasing aperture radius out to a radius of 2$\farcs$2, where it reachs 1780
Å, and then slowly rises to 1930 Å within an aperture of radius of
7$\farcs$2\. Observing J1044+0353 with a large aperture that subtends several
physical kiloparsecs produces a larger EW than does an SDSS fiber spectrum,
which subtends a diameter of 830 pc. We conclude that observations of a
J1044+0353 clone at high-redshift, where the angular diameter distance
produces scales of 5 to 8 kpc per arcsecond, through a seeing-limited
1$\farcs$2 slitlet (6-10 kpc in width) would capture most of the nebular
emission. However, spectra obtained through narrow 0$\farcs$2 NIRSpec shutters
(1 to 1.6 kpc) would yield lower equivalent widths and miss much of the
nebular emission.

####  2.1.4 Extraction of Emission-Line Images

We extracted pure emission-line images from the data cube. For each emission-
line of interest, we defined a narrowband filter with a top-hat transmission
function and added the continuum-subtracted line flux over the included
wavelength slices. Each synthesized filter was centered on the observed-frame
wavelength of the line, and the filter widths ranged from 6.5 to 14.5 Å. The
continuum under the line was modeled using the median level measured within
two line-free windows, one on each side of the emission line. The continuum
uncertainty was set equal to the difference between the two medians and
propagated through to the narrowband variance image.

We adaptively binned the line images to improve sensitivity to low-surface
brightness features. The narrow-band images presented a special challenge
because some of the pixels had negative values, a result of background and
continuum subtraction. Our approach combined the robust treatment of negative
values developed for binning background-subtracted X-ray data (Sanders &
Fabian, 2001) and the flexible bin sizes of the Cappellari & Copin (2003)
algorithm. Using a weighted Voronoi tessellation and following (Diehl &
Statler, 2006), we adaptively binned the flux to a local signal-to-noise
requirement. This technique preserves spatial resolution in high S/N ratio
regions while improving sensitivity in low surface brightness regions at the
cost of reduced spatial resolution. Structures visible in the binned images
motivated the definitions of the apertures used for photometry, but the
photometry was performed on unbinned images, avoiding the need to model the
covariance introduced by the adaptive binning.

Figure 3 illustrates the ionization structure of the galaxy. The very-high
ionization zone emitting He II $\lambda 4686$, shown in violet, is spatially
resolved but localized around the young starburst. The high-ionization gas
traced by the [O III] emission is detected into the circumgalactic medium,
i.e., well beyond the faintest countours of continuum emission in the Hubble,
SDSS, or DESI images. Throughout the region east of the starburst, the Balmer
emission lines sit in broad Balmer absorption troughs. We fit the absorption
with a second component (having a Lorentzian profile) and extracted both a
pure emission-line image, i.e. the flux relative to the bottom of the trough,
and the absorption equivalent width.

Figure 3:  KCWI emission-line image of J1044+0353: red (H$\beta$), green ([O
III] $\lambda 5007$), violet (He II $\lambda 4686$). For scale, the red circle
has a diameter of 2$\farcs$5 (695 pc). The H$\beta$ to [O III] ratio is
highest along the northeast and southeast lobes of the global, bipolar outflow
(see Peng et al. (2023) Figure 4). The [O III] nebula is spectroscopically
detected across the entire field; the [O III] $\lambda 5007$ contour levels
are spaced by actors of two from $1.5\times 10^{-17}$ to $4.80\times
10^{-16}$ergs s-1 cm-2

###  2.2 High-resolution Morphology from Hubble Imaging

We retrieved Hubble images of J1044+0353 from the Mikulski Archive for Space
Telescopes (MAST, PID 16209 & PID 16763). The composite image in Figure 1
includes filters from the UV (F125LP, F150LP, F165LP), UV/optical (F336W,
F438W, FR462N), and red/infrared (F606W, F665N, and F814W). We corrected the
WFC3/UV frames for charge transfer effects. Cosmic rays were removed from all
frames. Frames in each filter were drizzled onto matched images with a pixel
scale of 0$\farcs$02 per pixel. We aligned each image with a g-band image from
DESI Legacy Imaging Survey (Dey et al., 2019). The COS and SDSS spectra of
J1044+0353 cover only the starburst region.

To investigate the star formation history within the starburst region, we
detected and deblended star clusters using photutils.segmentation. Different
parameter combinations were explored across various filters, but they all
identified at least four subregions. Figure 4 shows the regions obtained using
the F336W image and the following parameters: nlevels = $1024$, contrast =
0.005, and n_sigma = 3. Table 2 lists 8-band photometry for these regions.
Region 303, which is the brightest in both F814W and F665N, consists of two
clusters separated by just 24 pc. Region 304 is slightly brighter than 303 in
the UV filters. We fit these spectral energy distributions (SEDs) with stellar
population synthesis models using BEAGLE (Chevallard & Charlot, 2016). Our
fiducial model describes each region with a single-age stellar population and
adopts a SMC extinction curve (Gordon et al., 2003). The stellar metallicity,
gas-phase metallicity, and ionization parameter were fixed at the following
values: $Z_{*}=0.05\mbox{${\rm~{}Z_{\odot}}$}$,
$Z=0.05\mbox{${\rm~{}Z_{\odot}}$}$ , and $\log(U)=-1.75$, respectively. The
fitted ages are 3.0 Myr, consistent with the Olivier et al. (2022) spectral
fit. The estimated stellar masses of Regions 303 and 304 are $3.5\times
10^{5}$${\rm~{}M_{\odot}}$ and $3.2\times 10^{5}$${\rm~{}M_{\odot}}$,
respectively. In order to quantify the systematic errors on these masses and
ages, we also fit the same photometry using an alternative extinction curve
(Charlot & Fall (2000), CF00) and a constant star formation rate, in total
four settings. Table 2 shows that the results are not sensitive to either the
reddening or the star formation history. We argue that the dominant systematic
uncertainty is actually the initial mass function (IMF). The BEAGLE models
assume an IMF based on the Galactic disk (Chabrier, 2003), but the
interstellar metallicity is much lower in J1044+0353. A higher upper mass
limit would not change the estimated cluster masses significantly.

Figure 4:  Top: Hubble F150W and F814W images zoomed in on the starburst.
Small circles (0$\farcs$1 radius) mark the centroids of individual star
clusters within the four region indentified by the segmentation map. The box
marks the centroid of the He II emission. Middle: Hubble F665N image shows the
H$\alpha$ emission at high resolution. Star clusters are marked with black
circles. The northwest shell extends 395 pc (1$\farcs$42) from Region 303. The
brighter, southwest extended 341 pc (1$\farcs$23) from Region 304. Contours
compare the extended UV emission; levels are $6.6\times 10^{-17}$, $1.3\times
10^{-16}$, and $2.7\times 10^{-16}$ ergs s-1 cm-2 Å-1 per square arcsecond
(uncorrected for extinction). For scale, an 0$\farcs$2 radius circle is shown
at the center of the 2$\farcs$5 red circle. Bottom: KCWI He II $\lambda 4686$
image. The very high-ionization nebula surrounds the compact starburst. The
outer red circle is 2$\farcs$5 in diameter. Contour levels are $2.35\times
10^{-18}$, $2.35\times 10^{-17}$, and $2.35\times 10^{-16}$ ergs s-1 cm-2 per
square arcsecond. The centroid is at 10:44:57.797, +03:53:13.17 with a
positional uncertainty of one spaxel, $\pm 0\farcs 0729$ or 40 pc. This places
its origin near the two star clusters in Region 303. The nebula is visibly
elongated, and the position angle of $-12\deg$, following the line of young
star clusters. This morphology is consistent with multiple star clusters
contributing $h\nu>54.4$ eV photons.

Figure 4 resolves the H$\alpha$ morphology within the core of the seeing-
limited PSF. Two H$\alpha$ shells protrude northwest of Region 303 and
southwest of Region 304. The curved filaments defining the northwest shell
converge on the Region 303 star clusters. The southwest shell emerges from the
Region 304 cluster. The UV surface brightness distribution is not symmetric
about the starburst; a bright region fills the southwest shell.

We measured the far-UV flux of J1044+0353 directly from the ACS/SBC F150LP
image. The filter pivot wavelength is 1606 Å. We used two apertures: a
2$\farcs$2 by 2$\farcs$56 box around the starburst region, and a 8$\farcs$42
by 4$\farcs$74 box that included most of the UV light from the galaxy.
Uncorrected for extinction, the starburst and total magnitudes were m150 =
18.70 and 18.34, respectively. After correcting for Galactic and internal
extinction, see Table 1, the absolute magnitudes are M${}_{FUV}^{SB}=-16.32\pm
0.17$ and M${}_{FUV}^{Tot}=-16.68\pm 0.17$. These values, listed in Table 3,
supercede the GALEX values which are affected by object blending.

Finally, we measured the axis ratio of J1044+0353 from the F814W isophotes.
Despite the irregular morphology of J1044+0353, the interstellar medium likely
has a disk component; the outflow axis is perpendicular to the major axis
(Peng et al., 2023). The rotation of the gas disk would be just $\pm 15$ ${\rm
km~{}s}^{-1}$, the circular velocity estimated by Xu et al. (2022). The
violent kinematics of the bipolar outflow, and a possible infalling gas cloud
(Peng et al., 2023), mask the rotational velocity of a gas disk in the KCWI
cube. At the F814W contour where the length of the major-axis is 2.3 kpc
(8$\farcs$27), we find the minor-to-major axis ratio is 0.64. If the
underlying distribution of stars is disk-like, then we estimate a disk
inclination of $50^{\circ}$ in the thin disk limit, higher for a thick disk.
For example, modeling the starlight as an oblate spheroid with intrinsic
height-to-radius ratio of 0.5, increases the inclination to $70^{\circ}$. This
edge-on orientation suggests that the radial outflow speed is significantly
higher than the line-of-sight outflow speeds, which are 40 - 50 ${\rm
km~{}s}^{-1}$.

##  3 Ionizing Photon Luminosities, Strömgren Radii, and RMS Densities

In this section, we describe the hydrogen and helium emission from the ionized
gas in and around J1044+0353. We describe the nebular morphology, measure the
Strömgren radii of each ionization zone, derive the ionizing photon
luminosities absorbed by the ISM at 1 and 4 Rydberg. The luminosities and
sizes of the H+ and He+2 nebulae then determine their root-mean-square (RMS)
electron density.

###  3.1 Recombination Line Images

Prominent emission lines from [Ar IV] $\lambda\lambda 4711.35,4740.20$ and He
I emission are visible in Figure 2 near the He II line. Imaging through a
narrowband filter would not separate these individual lines. In the post-
starburst regions, a narrow He I emission component sits in a broad,
photospheric He I absorption trough, which must be fitted and subtracted to
form a He I narrowband image. Narrow-band imaging would not measure this He I
line flux accurately because it would not recognize the depressed, local
continuum level. Integral field spectroscopy offers the only accurate method
for mapping recombination-line nebulae.

####  3.1.1 He+2 Recombination

Nebular He II $\lambda 1640$ and $\lambda 4686$ emission is produced by the
recombination of He+2, and thus demonstrates the presence of hard ionizing
photons with energies $h\nu>54.4$ eV. The velocity width of these nebular
lines is narrow, and the majority of local He II\- emitters, once AGN are
excluded, do not present a broad emission-line component (Shirazi &
Brinchmann, 2012). However, spectra of some star-forming galaxies show a broad
He II spectral feature attributed to the optically thick stellar winds of
classical Wolf-Rayet stars (Crowther, 2007; Vacca & Conti, 1992). These narrow
and broad components are not exclusive. Spatially resolved studies confirm
that classical Wolf-Rayet stars power many of the He II nebulae in M33 (Kehrig
et al., 2011) for example. However, the association between narrow He II
emission and classical Wolf-Rayet stars does not extend to low metallicity
galaxies (Kehrig et al., 2008), where very massive stars (VMSs), defined as
those having masses above 100${\rm~{}M_{\odot}}$(Vink et al., 2011), produce
broad, stellar He II emission lines, sometimes underneath a narrow nebular
emission line (Smith et al., 2016; Leitherer et al., 2018; Wofford et al.,
2023). After describing the overall morphology of the He II emission in
relation to the young, massive clusters in J1044+0353, we use our data cube to
carefully examine the He II line profile.

The bottom panel of Figure 4 shows the He II line image. The centroid
coincides with Region 303. The elongation of the He II contours follows the
chain of young clusters extending from Region 301 to 304. The He II emission
is detected beyond the COS aperture, indicated by the red circle in both
Figures 4 and 3. A He I nebula is detected in the $\lambda 4471$ line image
(not shown) at all position angles around SB. In contrast to the He II nebula,
however, we also detect He I further to the east, around FEN, where no He II
emission is detected.

If recombination accounts for all the He II luminosity, then the shape of the
He II profile will be related to other nebular emission lines emitted by the
very-high ionization zone. We therefore compared the He II and [Ar IV] line
profiles in the starburst spectrum, shown in the left panel of Figure 5. The
He II line profile is well fit by a single Gaussian component with standard
deviation $\sigma_{tot}=53$ ${\rm km~{}s}^{-1}$. After correcting the
linewidth for instrumental brodening, $\sigma_{inst}=35$ ${\rm km~{}s}^{-1}$,
the standard deviation is $\sigma_{corr}=40$ ${\rm km~{}s}^{-1}$,
corresponding to 94 ${\rm km~{}s}^{-1}$ FWHM. The [Ar IV] $\lambda 4740$
emission line profile is also well fit by a single Gauassian component,
$\sigma_{corr}=26$ ${\rm km~{}s}^{-1}$. The [Ar IV] linewidth is narrower than
He II and between the expectations for thermal broadening and turbulent
broadening. The larger atomic weight of argon, nine times that of helium,
reduces the thermal linewidth of argon lines by a factor of three relative to
helium lines. Assuming identical turbulent broadening of He II and [Ar IV], we
solved a system of two equations for the thermal velocity dispersion,

$\displaystyle\sigma_{th,HeII}=\sqrt{9/8(\sigma_{tot,HeII}^{2}-\sigma_{tot,ArIV}^{2})}=31{\rm~{}\mbox{${\rm
km~{}s}^{-1}$}},$ (1)

and the turbulent velocity dispersion,

$\displaystyle\sigma_{turb}=\sqrt{(\sigma_{tot,HeII}^{2}-\sigma_{inst}^{2}-\sigma_{th,HeII}^{2}}=25{\rm~{}\mbox{${\rm
km~{}s}^{-1}$}}.$ (2)

This interpretation of the broader He II linewidth yields a gas temperature,
$T=2.4\times 10^{5}{\rm~{}K}(\sigma/31{\rm~{}km/s})^{2}(m/4{\rm~{}amu})$.
Since this implied temperature is much higher than the measured temperature in
the high ionization zone, $T_{e}\approx 2\times 10^{4}$ (Peng et al., 2023),
we consider an alternative explanation, namely that the He II profile includes
a stellar component which does not emit [Ar IV]. Strong stellar winds emit
broad He II lines when the helium is doubly ionized. This He II component is
seen in spectra of young super star clusters in dwarf galaxies (Smith et al.,
2016; Leitherer et al., 2018; Wofford et al., 2023) as well as in classical
Wolf-Rayet stars. Those super star clusters include very massive main-sequence
stars rather than classical Wolf-Rayet stars which are evolved, helium-core
burning stars (Gräfener & Hamann, 2008).

Figure 5: Various Gaussian models fit to the He II line profile. (a) The He II
line profile is well fit with a single Gaussian component ($v=-31\pm 1$ ${\rm
km~{}s}^{-1}$, $125\pm 2$ ${\rm km~{}s}^{-1}$), but the line width is broader
than other nearby high-ionization lines. (b) Joint fit to to He II, [Ar IV],
and He I $\lambda 4713.14$ lines with a single linewidth. The positive
residual in the line wings of He II demonstrate its broader width. In the best
fit to the four lines, the He II component accounts for 96% of the total line
flux defined by the fit in Panel a, and the error bars allow the range from
from 80 to 100%. The broader width of He II may indicate that not all of the
line flux comes from He+2 recombination.

Panel b of Figure 5 illustrates the larger velocity spread of the He II
emission relative to [Ar IV] and He I. Subtraction of a joint fit to the He
II, [Ar IV], and He I $\lambda 4713$ lines from the data reveals positive
residuals extending up to $\pm 200$ ${\rm km~{}s}^{-1}$ around the He II line
center. The He II linewidths of super star clusters in NGC 3125 are 1000 to
1400 ${\rm km~{}s}^{-1}$ FWHM (Wofford et al., 2023), much larger than our
proposed second component in J1044+0353. One reason for the difference could
be the lower gas-phase metallicity of J1044+0353, which is just one-seventh
the LMC-like metallicity of NGC 3125. As stellar metallicity decreases,
stellar rotation is expected to increase, thereby decreasing the effective
escape speed and increasing the effective Eddington ratio (Meynet & Maeder,
2002). At stellar metallicities of $0.1\mbox{${\rm~{}Z_{\odot}}$}$, i.e.
slightly higher than the gas-phase metallicity of J1044+0353, theoretical
models predict He II line widths of just 490 ${\rm km~{}s}^{-1}$ FWHM for WNh
stars, a type of VMS (Gräfener & Vink, 2015). Considering the very young age
of the J1044+0353 starburst, we conclude that VMS are the most likely source
of the second, broader component that we identify in the He II line profile.

Gräfener & Vink (2015) argue that heii emission from low-metallicity stellar
winds could be mistaken for nebular emission due to the low wind velocities,
but this is not the case in J1044+0353. The He II luminosity from J1044+0353
is dominated by nebular emission. This very-high ionization nebula subtends a
larger solid angle than the string of compact clusters in Figure 4, so the He
II emission cannot be dominated by stellar emission. The detection of
collisionally excited [Ar IV] lines, which have no stellar counterpart,
further demonstrate the presence of a very-high ionization nebula. The
amplitudes of the He II and [Ar IV] $\lambda 4711$ line are equal across most
of the nebula; the only exception is a spectrum extracted from a seeing
limited aperture centered on Region 303, where the He II line has a larger
amplitude than the [Ar IV] $\lambda 4711$ line.

Photoionization by an AGN, essentially a point source, is not consistent with
the overall elongation of the nebula in Figure 4 or the optical emission-line
ratios. The high production efficiency of He+-ionizing photons appears to be a
property of very young, low-metallicity stars, a population which likely
includes VMS in J1044+0353. The presence of this stellar He II component
lowers the inferred photoionization rate, Q(He+), slightly. We estimate this
correction by noting that the residual flux in Panel b accounts for
$4^{+16}_{-4}\%$ of the total line flux in Panel a. Hence
$\mbox{Q(He${}^{+}$)}_{corr}=0.96^{+0.04}_{-0.16}\mbox{Q(He${}^{+}$)}$, where
Q(He+) is our measured value below.

####  3.1.2 Hydrogen Recombination

We present a global view of the hydrogen recombination-line nebula in Figure
3. The faint shells contribute very little of the total H$\beta$ luminosity.
Most of the recombinations occur near the compact starburst. There, the
H$\beta$ isophotes are round out to 1.2 kpc (4$\farcs$3). It is therefore
reasonable to describe the nebula photoionized by the starburst as a galactic-
scale Strömgren sphere.

On larger scales, the H$\beta$ surface brightness profile depends on position
angle. In projection on the sky, the starburst to shell separation reaches
$\mbox{R${}_{X}$}=3.3$ kpc (12$\farcs$0). The faint shell protruding directly
north of SB extends to $R=1.99$ kpc (7$\farcs$15). To the west of SB, low
surface brightness emission extends to at least 2.5 kpc. In the east, the
H$\beta$ contours and nearly circular around FEN. The presence of stars hot
enough to ionize hydrogen indicates ongoing star formation in FEN. The Balmer
absorption in the FEN line profiles, however, indicate the ratio of Type A
stars to Type O stars is larger than its ratio in SB. Peng et al. (2023)
modeled the star formation history as a burst 20 Myr ago plus a constant
component. In contrast to the Strömgren spheres powered by FEN and SB, the
peanut shape of the H$\beta$ contours draws attention to the minimum in the
H$\beta$ emission between them. The faintest H$\beta$ contours follow the
extended [O III] filaments visible in Figure 3. These loops are shells of gas
which the outflow has swept out of its path.

###  3.2 Photometry

The compact SB is tiny compared to the nebula it ionizes, a geometry similar
to the H II regions around massive stars, but on a global scale. The centroid
of the He II emission lies in Region 303, and we measured integrated emission-
line fluxes in a series of concentric circular apertures around the centroid
of the He II emission In each aperture, we measured the flux in the
recombination-line images, and we extracted a spectrum from which we directly
measured the integrated line profile. These methods gave consistent fluxes.
The photometric growth curve for the He II $\lambda 4686$ emission reaches a
well-defined asymptotic level, so the total flux is well defined.

Figure 6 presents the radial surface brightness profiles. The ionization
structure follows the expectation for a centrally concentrated source; the
very-high ionization zone is smaller than the other zones in Figure 6. The
seeing-limited PSF shapes the profile at radii less than 1$\farcs$0\. The
surface brightness profile falls steeply between 1$\farcs$0 and 2$\farcs$0\.
The He I surface brightness declines at a rate similar to the He II emission
out to 2$\farcs$0, beyond which the nebula centered on the FEN region flattens
the profile. The edge of the He+2 Strömgren sphere is therefore sharper than
the boundary of the He+ nebulae. The inner H$\beta$ profile can be described
by the same exponential scale length as He II and He I, but contribution from
the FEN region flatten the profile at larger radii. The aperture photometry
describes the nebular structure out to radii where the isophotes can be
described as round but does not capture the filamentary structures at larger
radii.

Figure 6: Surface brightness profiles. The post-starburst region is detected
at radius 5$\farcs$1 in the He I $\lambda 4471$. The H$\beta$ emission is
detected at $10\sigma$ at radius 7$\farcs$14, the largst circular aperture
with complete coverage. Along the long axis of the mosaic, the H$\beta$
emission is detected to radius 12$\farcs$0\.

At many position angles, the H$\beta$ emission extends to the edge of the
mosaic cube. The profile in Figure 6 is truncated at the largest circular
aperture that fits within the rectangular mosaic. The slope of the H$\beta$
growth curve is very close to zero at these distances, so the filamentary
emission would not increase the total observed H$\beta$ flux significantly. In
Table 3, we list two radii determined from visual examination of the binned
recombination-line images. The smaller one, $R_{c}(X)$ where $X$ identifies
the emission line, describes the radius where the shape of the intensity
contours transitions from largely circular, to mainly filamentary, structures.
The scale of the more filamentary emission is position angle dependent and
limited by the mosaic field-of-view for H$\beta$. We interactively examined
the spectrum across the nebula to determine the maximum extent of the H$\beta$
emission. We denote the maximum projected radius relative to the starburst as
$R_{x}(X)$. These radii are marked in Figure 6. These visual classifications
correspond to S/N ratios of 5 and 3 roughly. We adopt the $R_{c}$ values as
the best direct measure the Strömgren radii. The He II emission is detected
within $\mbox{R${}_{C}$}=610$ pc (2$\farcs$2) of the compact starburst at all
position angles and extends to $\mbox{R${}_{X}$}=730$ pc to the northwest and
the southeast.

###  3.3 Ionizing Photon Luminosities

In this section, we assume the nebula is photoionized. Then the recombination
rates determine the ionization rate of H and He+ in the galaxy. We calculate
Q(H) and Q(He+), the photon luminosities at $h\nu\geq\ 13.6$ eV and $h\nu\geq\
54.4$ eV respectively, from the recombination-line luminosities measured from
the KCWI data cubes. These integrated luminosities are listed in Table 3 and
are larger than the values measured from fixed apertures such as the SDSS
fiber spectrum. We discuss modifications for the possible stellar contribution
to the He II luminosity in Section 6.

Following recombination of an electron with a hydrogen or helium ion, the
electron cascades downward to lower energy levels producing the recombination
line emission. In photoionization equilibrium, the H I and He II
recombination-line luminosities determine the photoionization rates of H atoms
and He+ ions, respectively. The luminosities of H$\beta$ and He II $\lambda
4686$ therefore measure the photon luminosity absorbed by the ISM at energies
exceeding 13.6 eV and 4 Rydberg (54.4 eV), respectively. Correction for the
leakage of ionizing radiation from the galaxy then yields the intrinsic EUV
continuum level at these energies, Q(H) and Q(He+) respectively.

In photo-ionization equilibrium, the hydrogen ionizing photon luminosity
follows directly from the luminosity of a recombination line,

$\displaystyle\mbox{Q(H)}=\frac{\alpha_{B}(H,T)}{\alpha^{eff}_{H\beta}h\nu_{H\beta}}\frac{L_{H\beta}}{(1-\mbox{$f_{esc}$})},$
(3)

where $(1-\mbox{$f_{esc}$})$ denotes the fraction of the ionizing photon
luminosity absorbed by the ISM. The on-the-spot approximation for the diffuse
radiation field (Osterbrock & Ferland, 2006) may break down along some
channels through the ISM, but the recombination emission will come from gas
that presents a significant cross section to the ionizing radiation. We
therefore adopt Case B recombination coefficients (Pequignot et al., 1991;
Storey & Hummer, 1995) at $T_{e}=2.0\times 10^{4}$ K, the electron temperature
measured in the highly-ionized O+2 zone (Peng et al., 2023). Inserting the
recombinatino rates in Eqn. 3, the ionizing photon luminosity becomes

$\displaystyle\mbox{Q(H)}=4.35\times
10^{52}{\rm~{}s}^{-1}\frac{L_{H\beta}}{2\times
10^{40}{\rm~{}ergs~{}s}^{-1}(1-\mbox{$f_{esc}$})},.$ (4)

which evaluates to $\mbox{Q(H)}=5.04\pm 0.48\times
10^{52}(1-\mbox{$f_{esc}$})^{-1}{\rm~{}s}^{-1}$ for the extinction-corrected
luminosity of J1044+0353.

The He+-ionizing continuum at energies $\geq\ 54.4$ eV can be directly
measured from the luminosity of the He II $\lambda 4686$ recombination line.
In photoionization equilibrium, we have

$\displaystyle\mbox{Q(He${}^{+}$)}=\frac{\alpha_{B}(\mbox{He${}^{+}$},T)}{\alpha^{eff}_{\lambda
4686}h\nu_{4686}}\frac{L_{\lambda 4686}}{(1-\mbox{$f_{esc,54}$})}$ (5)
$\displaystyle=4.92\times 10^{50}{\rm s}^{-1}\frac{L_{\lambda 4686}}{4\times
10^{38}(1-\mbox{$f_{esc,54}$})},$

giving an ionizing photon luminosity $\mbox{Q(He${}^{+}$)}=4.89\pm 0.47\times
10^{50}(1-\mbox{$f_{esc,54}$})^{-1}$ s-1 for J1044+0353. The result is
sensitive to the electron temperature, and we note that reducing the electron
temperature from $2.0\times 10^{4}$ K to $T_{e}=1.0\times 10^{4}$ K yields the
coefficient in Equation 3.3 1.23 times higher. We introduce $f_{esc,54}$ to
distinguish the escape fraction at $h\nu\geq 4$ Rydberg from $f_{esc}$ at the
Lyman limit.

The spectral bandpass also includes several He I recombination lines, the
strongest of which is He I $\lambda 4471$. The size of the He+ zone is
sensitive to the temperature of the stars because hydrogen and helium present
similar cross sections at the helium ionization edge of 24.6 eV (Osterbrock &
Ferland, 2006). The more abundant H atoms make up for the larger cross-section
of He+ per ion. It follows that the photons produced by recombinations to He0,
however, do not necessarily photoionize another helium atom. No simple
relation exists to convert He I line luminosities to the ionizing photon
luminosity at 24.6 eV (Schaerer & Vacca, 1998).

###  3.4 Root-mean-square Electron Density

The overall size of the H$\beta$ and He II nebulae provides important
constraints on the structure of the ISM. The ionization front driven by the
young star clusters stalls when the ionized volume becomes large enough that
the recombination rate balances the ionizing photon luminosity. Even when the
emitting clouds have high electron density $n_{e}$, nebulae can grow to large
sizes because the clouds typically fill a tiny fraction of the geometrical
volume.

The ionized volume is the geometrical volume multiplied by the volume filling
factor, $\epsilon$, of the emitting clouds (Strömgren, 1948). For a spherical
nebula that volume is $\frac{4\pi}{3}\mbox{R${}_{\rm
S}$(H${}^{+}$)}^{3}\mbox{$\epsilon$}$, and we write the Strömgren radius as

$\displaystyle\mbox{R${}_{\rm
S}$(H${}^{+}$)}=\bigg{(}\frac{3\mbox{Q(H)}(1-\mbox{$f_{esc}$})}{4\pi\alpha_{B}({\rm
H^{0}},T)n_{e}n_{p}}\bigg{)}^{1/3}\mbox{$\epsilon$}^{-1/3},$ (6)

in order to show the inverse dependence on the filling factor. We assume
$\epsilon$ is constant throughout the nebula, then the quantity
$\mbox{$\epsilon$}n_{e}^{2}$ defines the mean-square electron density,

$\mbox{$\langle n_{e}^{2}\rangle$}\equiv\mbox{$\epsilon$}n_{e}^{2}.$ (7)

Using Eqn. 4, we can express Q(H) directly in terms of the observed H$\beta$
luminosity, L(H$\beta$). Solving Equation 6 for the root-mean-square (RMS)
density, we then obtain the RMS density in terms of the measured size and
luminosity of the nebula,

$\displaystyle\mbox{$\langle
n_{e}^{2}\rangle^{1/2}$}=1.22{\rm~{}cm}^{-3}\bigg{(}\frac{1{\rm~{}kpc}}{\mbox{R${}_{\rm
S}$(H${}^{+}$)}}\bigg{)}^{3/2}\bigg{(}\frac{L(\mbox{H$\beta$})}{10^{40}{\rm~{}ergs~{}s}^{-1}}\bigg{)}^{1/2}\times~{}~{}~{}$
(8) $\displaystyle\bigg{(}\frac{1.610\times
10^{-14}{\rm~{}ergs~{}cm^{3}~{}s^{-1}}}{\alpha_{H\beta}^{eff}({\rm
H^{0}},T)}\bigg{)}^{1/2}\bigg{(}\frac{n_{e}/n_{p}}{12/10}\bigg{)}^{1/2}.$

We adopted an electron temperature $T_{e}=2.0\times 10^{4}$ K, appropriate to
the O+2 zone in J1044+0353; cooler temperatures of $1.5\times 10^{4}$ K (or
$1.0\times 10^{4}$ K) would lower the coefficient to 1.07 (or 0.89,
respectively).

The estimated RMS density in the galaxy depends on the aperture. For the H+
Strömgren sphere, Table 3 lists $L(\mbox{H$\beta$})=2.54\times 10^{40}$ ergs
s-1 and RC(H+) = 1990 pc. Substitution of these values into Eqn. 3.4 gives a
global RMS density of 0.69 cm-3. The J1044+0353 H$\beta$ emission is strongly
concentrated near SB, and roughly 82% of that H$\beta$ luminosity is emitted
within the much smaller aperture describing the He+2 Strömgren sphere,
$\mbox{R${}_{C}$}(\mbox{He${}^{+2}$})=610$ pc. Eqn. 3.4 therefore indicates a
higher RMS density, $\mbox{$\langle n_{e}^{2}\rangle^{1/2}$}=3.7$ cm-3, for
the H$\beta$ emission from the very-high ionization zone.

In close analogy to Eqn. 6, the photon luminosity at energies $\geq\ 54.4$ eV
determines the volume of very-highly ionized clouds in photoionization
equilibrium.

$\displaystyle\frac{4\pi}{3}\mbox{R${}_{\rm
S}$(He${}^{+2}$)}^{3}\mbox{$\epsilon_{h}$}=\frac{\mbox{Q(He${}^{+}$)}(1-\mbox{$f_{esc}$}54)}{\alpha_{B}({\rm
He^{+}},T)n_{e}n_{He+2}}.$ (9)

We adopt a fiducial ratio of H-to-He atoms of 10:1 by number, consistent with
the 12 electrons per 10 protons adopted in Eqn. 3.4. This helium abundance is
equivalent to a helium mass fraction $Y=0.286$, a bit higher than the lower
bound set by the primordial fraction, $Y=0.2470$ (Planck Collaboration et al.,
2020). We assume that helium is fully ionized in the very-high ionization
zone. A photon luminosity Q(He+) then ionizes a volume of radius

$\displaystyle\mbox{R${}_{\rm
S}$(He${}^{+2}$)}=208{\rm~{}pc}\bigg{(}\frac{\mbox{Q(He${}^{+}$)}}{10^{50}{\rm~{}s}^{-1}}\bigg{)}^{1/3}\bigg{(}\frac{1{\rm~{}cm}^{-3}}{\mbox{$\langle
n_{e}^{2}\rangle^{1/2}$}}\bigg{)}^{2/3}\times\ ~{}~{}~{}$ (10)
$\displaystyle\bigg{(}\frac{1/10}{n(He^{+2})/n_{e}}\bigg{)}^{1/3}\bigg{(}\frac{9.03\times
10^{-13}{\rm~{}ergs~{}cm}^{3}~{}{\rm
s}^{-1}}{\alpha_{B}(He^{+},T)}\bigg{)}^{1/3}\mbox{$\epsilon_{h}$}^{-1/3},$

where we have evaluated the recombination coefficient at a fiducial
temperature of $T_{e}=2.0\times 10^{4}$ K.444 The radius would be 93% (83%) of
this value at temperatures of $1.5\times 10^{4}$ K or $1\times 10^{4}$ K,
respectively. The filling factor of the very highly ionized gas,
$\epsilon_{h}$, is then defined by the ratio of the RMS density derived from
the He II emission and the local, physical density of the clouds.

Solving for the RMS density in the He+2 zone and substituting for Q(He+) using
Eqn. 3.3, we find

$\displaystyle\mbox{$\langle
n_{e}^{2}\rangle^{1/2}$}=0.25{\rm~{}cm}^{-3}\bigg{(}\frac{600{\rm~{}pc}}{\mbox{R${}_{\rm
S}$(He${}^{+2}$)}}\bigg{)}^{3/2}\bigg{(}\frac{L_{\lambda
4686}}{10^{38}{\rm~{}ergs~{}s}^{-1}}\bigg{)}^{1/2}\times\ ~{}~{}~{}~{}$ (11)
$\displaystyle\bigg{(}\frac{n(H)/n(He)}{10}\bigg{)}^{1/2}\bigg{(}\frac{10/12}{n_{p}/n_{e}}\bigg{)}^{1/2}.$

For a He II luminosity $4\times 10^{38}$ ergs s-1 and RC(He+2) = 610 pc, this
expression yields an RMS density $\langle n_{e}^{2}\rangle^{1/2}$= 0.50 cm-3
for the very-high ionization zone.

In summary, our analysis of the integrated spectrum extracted from the fixed
aperture defined by RC(He+2) finds different RMS densities for the He II\- and
H$\beta$-emitting regions. This result has implications for the physical
structure of the H II region. In the next section, we measure the local
electron density in the different ionization zones. Then, in Sec. 4.1.3, we
compare these physical densities and the RMS densities estimated here, thereby
measuring the filling factor in each ionization zone via Eqn. 7.

##  4 Global Ionization and Density Structure

The mosaic cube detects [O III] $\lambda\lambda 4959,5007$ emission across the
entire field of view shown in Figure 3. The emission-line nebula therefore
extends further beyond the starlight than a casual comparison to Figure 1
suggests. Large spatial gradients in the intensity ratio of the [O III]
doublet relative to the [O II] $\lambda\lambda 3726,29$ doublet have two
practical implications: (1) the area over which the [O II] doublet can be used
to measure the electron density is restricted; and (2) we can map out the
pathways of Lyman continuum leakage. In this Section, we use the collisionally
excited lines to describe the physical properties of the ionized gas: electron
density and volume filling factor in Section 4.1, and the ionization parameter
and LyC leakage in Section 4.2.1.

###  4.1 Density Structure

Densities derived from ions with very different ionization potentials provide
insight into the structure of the ISM. In this section, we derive the gas
density in low and high ionization zones using density-sensitive ratios of
emission lines. We then argue that the volume filling factor of these clouds
is quite low, a conclusion that follows from the observed sizes of the
spatially resolved Strömgren radii.

####  4.1.1 Density Gradients with Ionization Zone

The intensity ratio of lines collisionally excited from the same ground state
depends only on the relative statistical weights of their upper levels when
the gas density is well below the critical densities of the transitions. In
the density range where collisional de-excitation rates are important in one
of the two lines, the line ratio directly measures the electron density. The
bandpass covered by our KCWI observations covers several density-sensitive
emission-line doublets. The [O II] doublet, hereafter ${\rm O2DR~{}}\equiv
F(\lambda 3729)/F(\lambda 3726)$, probes the low-to-intermediate ionization
zones where ions with ionization potentials between 13.6 - 35.1 eV coexist. We
applied PyNeb (Luridiana et al., 2015) to derive densities from line ratio
measurements..

In order to compare the density in different ionization states, we will
discuss a spectrum extracted from a circular aperture of radius RC(He+2). To
place these densities in the broader context of the galaxy, we note that the
SB region has higher $n_{e}$ than FEN, NEN, or MEN. In addition, spectra
extracted from larger apertures have a higher O2DR, equivalent to a lower
$n_{e}$. Within RC(He+2) aperture, the [O II] density is $185^{+9}_{-11}$
cm-3. In the larger aperture defined by RC(H+), but still centered on SB, the
density is $161^{+7}_{-10}$ cm-3.

Figure 2 shows the strong nebular [Ar IV] $\lambda\lambda 4711,40$ emission.
The [Ar IV] doublet ratio $F(\lambda 4711)/F(\lambda 4740)$ is sensitive to
the electron density in regions of the ISM where ions with ionization
potentials between 40.7 - 59.8 eV are found; a region including both the high-
ionization zone where [O III] is emitted (35.1 - 54.9 eV) and the very-high
ionization zone emitting He II. The J1044+0353 KCWI cube detects the [Ar IV]
doublet at high S/N ratio across the entire He II\- emitting zone, as defined
by RC(He+2) in Table 3.

We measured the [Ar IV] line strengths by fitting Gaussian profiles. Table 3
lists the results. Applying PyNEB to these line ratios, we computed an [Ar IV]
density of $1440^{+44}_{-128}$ cm-3. The [O II] density was lower,
$185^{+9}_{-11}$ cm-3, in the same spectrum. This result demonstrates that the
density in the very-high ionization zone is larger than the density in the
low-to-intermediate ionization zone; the large difference provides insight
into the structure of the ISM. We do not detect significant spatial variations
in $n_{e}(\mbox{[Ar IV]})$ within RC(He+2).

Intermediate ionization zone (24.4 - 47.9 eV) densities measured via the [C
III] $\lambda 1907$ / C III] $\lambda 1909$ ratio typically follow a similar
trend where the C III] densities are higher than [O II] densities (James et
al., 2014). The KCWI spectra also cover the [Cl III] $\lambda\lambda 5517,37$
doublet, a density diagnostic probing the intermediate-to-high ionization zone
(23.8 - 39.6 eV). In our KCWI spectra, the [Cl III] doublet lands in a
spectral region attenuated by the dichroic beam splitter. As a result, the S/N
ratio is much lower than we obtained for the other density measurements, and
we do not discuss it further here. In the next section, we provide further
insight into the structure of the low-ionization and very-high ionization gas
clouds by computing their volume filling factors.

For completeness, we note that recombination of oxygen ions can lead to [O II]
emission from very dense regions (Osterbrock & Ferland, 2006), but the gas
density in J1044+0353 is not high enough for recombination to enhance [O II]
emission.

####  4.1.2 Density Gradients with Distance from SB

The density gradient across J1044+0353 may be influenced by several factors
including gas flows and hydrostatic equilibrium. To gain further insight about
the structure of the ISM, we mapped the electron density across this low mass
galaxy. Toward the northwest, west, and southwest of the starburst, the
density declines with radius. The density gradient is shallowest to the east
(PA $\approx 90\deg$), where it remains constant over more than a kiloparsec,
and steepest to the west (PA $\approx 270\deg$). The galactic center of mass
is not well constrained, but the high densities suggest it lies somewhere in
the region between SB and MEN.

We measured the O2DR for a set of spectra extracted from apertures defined by
an adaptively binned map of the [O II] emission. Figure 7 plots the electron
density derived from these O2DR measurements as a function of the projected
distance from SB. The dispersion in the $n_{e}$ values increases with
separation from SB. We have labeled the points by position angle in order to
contrast the density gradients in different directions.

Figure 7: Radial density profile. Projected distance was computed relative to
the compact starburst. Electron density calculated from [O II] $\lambda\lambda
3726,29$ doublet ratio at $T_{e}=1.5\times 10^{4}$ K. The gray band shows the
much higher density derived from [Ar IV] in the very-high ionization zone.
Median values are computed separately for PA = 10∘ to -170∘ counter-clockwise
(large black squares), and position angles from 10∘ to 190∘ clockwise (large
black circles). The standard deviation (black error bars) illustrate the large
spread in density. The density is roughly constant in the direction of the
post-starburst region. In contrast, the density gradient is more negative in
the hemisphere toward the west. The point-spread function flattens the density
profile $n_{e}(r)$ within 250 pc ($2\times$ FWHM) of SB.

####  4.1.3 Volume Filling Factor

The 3D spectroscopy made it possible to resolve the spatial extent of the
ionized nebula and its inner region where helium is doubly ionized. These
measurements determine the RMS electron densities over the associated
macroscopic volume. The RMS densities are much lower than the local density in
the clouds because the emitting clouds do not fill the Strömgren sphere. We
used density-sensitive line ratios to measure the physical density in the
ionized gas clouds. Comparison to the RMS density then defines their volume
filling factor, $\mbox{$\epsilon$}\equiv\frac{\mbox{$\langle
n_{e}^{2}\rangle$}}{n_{e}^{2}}$. Table 1 summarizes the density measurements,
and inspection shows all the physical densities are higher than any of the RMS
densities.

Volume filling factors for H II regions are typically derived from [O II] or
[S II] densities because high-ionization lines are often weak and/or in the
rest-UV spectrum. Over the full Strömgren sphere, RC(H+), our measurements
yield $\epsilon$$\approx(0.69/161)^{2}\approx 1.8\times 10^{-5}$. These
filling factors are very low compared to Galactic H II regions, where the
filling factors are a few percent (Osterbrock & Ferland, 2006). Using the KCWI
cube, spectral extraction from smaller apertures leads to higher filling
factors because the gradient in the physical density is shallower than the
surface brightness profile of the recombination lines. Within the He II
Strömgren sphere of J1044+0353, for example, we estimate
$\epsilon$$\approx(3.7/185)^{2}\approx 4\times 10^{-4}$.

The strong He II and [Ar IV] emission of extreme emission-line galaxies
(EELGs) make it possible to measure another filling factor, $\epsilon_{h}$.
For this very-high ionization gas, we find
$\mbox{$\epsilon_{h}$}\approx(0.5/1440)^{2}\approx 1.2\times 10^{-7}$. Within
the RC(He+2) aperture, we then have
$\mbox{$\epsilon_{h}$}/\mbox{$\epsilon$}\sim 3\times 10^{-4}$. The very-high
ionization emission comes from 3% the photoionized gas by volume.

###  4.2 Ionization Structure

In the canonical picture of a spherical H II region, the ionization structure
is layered in concentric shells, the higher ionization zones lying interior to
the lower ionization zones. When this model is applied to the global structure
of a galaxy, emission from the low ionization-state gas comes from a zone
located at larger radii than the very-high and high ionization zones.
Photoionization models with four zones can fit most of the emission lines (not
He II) in the J1044+0353 spectrum (Berg et al., 2021). Here, however, we argue
that this onion-like structure does not provide an adequate description of the
global ionization structure in J1044+0353. Specifically, the [O II] emission
does not come from a shell. If it did, then we would see the same shell at all
projected radii, but we find a density gradient. The emission from very-high
ionization states is centrally concentrated, but, in contrast to the canonical
picture, we find lower density, less ionized gas co-exists with the highly-
ionized, high-density gas. Neither of these components of ionized gas fill the
majority of the volume, so the observed ionization structure requires a
multiphase medium. Consistent with this picture, we note that radiation
pressure dominated H II regions have a wide range of gas densities (and
therefore ionization parameters) within individual clouds (Stern et al.,
2016).

A population of clouds surrounded by much hotter (and lower density) gas
offers a more realistic schematic. The dynamical importance of both the hot
phase and radiation pressure have been discussed in the context of quasar
outflows (Stern et al., 2016). Propagation of an ionization front through the
clouds stratifies the ionization structure on the spatial scale of individual
clouds. The ionization state decreases with increasing depth into a cloud.
With this picture in mind, the He II emission comes from the cloud surfaces
facing the compact SB, and this very-high ionization zone is ionization-
bounded within a typical cloud. The [O II] emitting volume of a typical cloud
is higher than the He II emitting volume, even within RC(He+2). This structure
requires very dense gas clouds in addition to a hard ionizing spectrum. We
discuss the Stern et al. (2016) models in the context of J1044+0353 further in
Sec. 6. In this section we use traditional photoionization models to interpret
emission-line ratios. We select plane-parallel, rather than spherical
geometry, and picture the resulting ionization gradients as describing
individual clouds. Varying the ionization parameter of the cloud then
describes how radial gradients in cloud density and radiation intensity change
emission-line ratios.

####  4.2.1 Ionization Parameter

The value of the ionization parameter, the number of ionizing photons per
hydrogen atom, causes significant variations in emission-line ratios among H
II regions having the same gas-phase metallicity and ionizing spectrum. The
strongest feature in rest-optical spectra is the strength of collisionally-
excited lines from O+2 relative to O+. At metallicities below solar, the [O
III] / [O II] flux ratio, or O32, increases linearly with the ionization
parameter (Kewley & Dopita, 2002). Previous work identified J1044+0353 as a
candidate Lyman continuum leaker based on the high O32 ratio of the starburst
spectrum (Berg et al., 2022). Mapping the projected ionization structure
provides new insight about the ionization structure of the nebula and the
anisotropic escape of LyC radiation. In this paper, we define the O32 ratio by
the strength of the [O III] $\lambda\lambda 4959,5007$ doublet relative to the
[O II] $\lambda\lambda 3726,3729$ doublet. For comparison to papers that
exclude [O III] $\lambda 4959$ from the numerator of O32, we multiply their
ratios by four-thirds (an addition of 0.125 dex) before comparing to our O32
measurements.

Figure 8 presents our map of the O32 ratio. The mean O32 ratio within
0$\farcs$9 of SB is 18. In the post-starburst region around FEN, the O32 ratio
is much lower, roughly 3. The O32 values are highly position-angle dependent
in the annulus between 0$\farcs$9 and 1$\farcs$5 surrounding SB. The minimum
ratio is O32 $=10$ at PA=225∘. The high O32 ratio to the north north-east of
SB (PA = 44 to -11) defines an ionization cone. The O32 values between
position angles of 0 and $25^{\circ}$ remain between 13 and 18 out to
projected radii of at least 4$\farcs$2 (1.2 kpc). The solid angle subtended by
this cone intersects the bipolar outflow identified by Peng et al. (2023). Its
size is roughly 20% of $4\pi$ steradian near SB. At larger radii, the O32
measurements are truncated because no [O II] emission is detected. Using the
[O III] contours to extend the ionization cone, we measure a lower limit,
${\rm O32}>18$, at the northern edge of the mosaic. This definition of the
ionization cone is extremely conservative. At nearly every position angle, we
detect the [O III] doublet and H$\beta$ lines at the edge of the mosaic.

Figure 8: The O32 line ratio map. The [O III] doublet is detected all the way
to the edge of the map; but the map is black where no [O II] emission is
detected. The red and black circles mark the location of the starbust and FEN,
respectively. The negative radial gradient to the east and west of the
starburst are typical of ionization-bounded galaxies. In contrast, north-
northeast of the starburst, the O32 ratio remains high to the edge of the
frame, a sign of a density bounded nebula. The O32 ratio exceeds 18 in the
region between the two outermost contours to the north-northeast, and we call
this region the ionization cone. The contours denote the [O III] $\lambda
5007$ surface brightness (by factors of two from $1.5\times 10^{-17}$ to
$4.8\times 10^{-16}$ ergs s-1 cm-2 per square arcsecond. See Section 4.2.1 for
further details.

In a plane-parallel geometry, the ionization parameter describes the ratio of
the ionizing flux at the inner surface of a cloud relative to the number
density of hydrogen atoms. To assign an ionization parameter to the O32
values, we adopt the $Z=0.05\mbox{${\rm~{}Z_{\odot}}$}$ relation in Figure 1
of Kewley & Dopita (2002), corresponding to $12+\log(O/H)=7.6$, which is the
closest model to the SB metallicity. The ionization parameter, $q$, describes
the maximum speed at which the radiation field can drive an ionization front
through a cloud complex. O32 ratios of $15,10,5$ indicate $q=1300,1000,{\rm
and~{}}600$ ${\rm km~{}s}^{-1}$, respectively. It follows that the timescale
for ionizing gas clouds, $\tau<1.5\times
10^{5}{\rm~{}yr~{}}(r/100{\rm~{}pc})(q/600{\rm~{}km/s})^{-1}$ on the scale of
the compact star clusters, is shorter than the cluster ages. This result is
consistent with the ionization front having radii of several kpc.

Comparison to other models often requires the dimensionless ionization
parameter, $U\equiv q/c$, where our O32 values of $15,10,5$ correspond to
$\log U=-2.36,-2.48,{\rm and}-2.70$. One systematic uncertainty affecting the
interpretation of O32 ratios is model geometry. Spherical models, for example,
give ionization parameters $\approx 0.2$ dex higher than the plane-parallel
geometry adopted.555 We read the commonly used dimensionless form of the
ionization parameter $U$ off Figure 5 in Berg et al. (2016), where the model
geometry is spherical and compared to Kewley & Dopita (2002). An additional
systematic correction upwards is needed in the He II\- emitting region where
some oxygen has been further ionized to O+3. Relative to three-zone models, a
four-zone model systematically yields higher $U$ by up to 0.5 dex (Berg et
al., 2021). Since the ionization cone extends well above the region where the
gas-phase O/H ratio has been directly measured, we also considered the
potential impact of metallicity gradients. Metals in the central region of
J1044+0353 have likely been diluted by the accretion which fueled the
starburst, and the outflowing gas may have been enriched by cooling of a hot
wind (Peng et al., 2023). At fixed O32, adopting the curve for a higher
metallicity in Figure 1 of Kewley & Dopita (2002) would further increase the
estimated ionization parameter in the ionization cone. Our argument that the
O32 map reveals an ionization cone therefore appears to be robust to these
systematic uncertainties.

####  4.2.2 Evidence for Anisotropic Lyman Continuum Escape

The ionization cone defined by the anisotropic O32 emission in Figure 8
provides new insight into the controversy over the relationship between
measured O32 ratios and $f_{esc}$(Faisst, 2016; Naidu et al., 2018). Previous
studies, see Figure 8a of Izotov et al. (2021) and Figure 10a of Bassett et
al. (2019) for example, demonstrate that galaxies with high $f_{esc}$ also
have high O32 ratios, but that galaxies with high O32 ratios do not
necessarily have directly detectable Lyman continuum. The scatter in the
relationship between the O32 ratio and the directly measured Lyman continuum
leakage has been interpreted as evidence that O32 is a mediocre predicter of
$f_{esc}$ for any individual galaxy (Izotov et al., 2020). While their
conclusion accurately describes leakage directly along our sightline to a
galaxy, we interpret the anisotropic O32 emission in Fig. 8 as evidence that
the global LyC leakage can be highly anisotropic. Our data therefore support
the sketch in Figure 11 of Bassett et al. (2019) where a galaxy with an
ionization cone pointed away (toward) the observer produces high O32
($f_{esc}$) but low $f_{esc}$ (O32). Galaxy orientation may therefore
determine which strong W([O III]) emitters have directly detected Lyman
continuum emission (Tang et al., 2021). If anisotropic escape is common, we
argue that O32 could be a more robust measure of the global $f_{esc}$ than the
residual FUV continuum.

Knowledge of the orientation and covering factor of density-bounded channels,
in addition to measurements of gas-phase metallicity and ionization parameter,
will likely be required in order to accurately calculate the relationship
between the O32 ratio and $f_{esc}$ via photoionization modeling. Density
bounded models produce large O32 ratios because the outer, low-ionization zone
emitting [O II] is truncated; an observed O32 ratio can be produced with a
lower ionization parameter than an ionization-bounded model would require
(Pellegrini et al., 2012; Jaskot & Oey, 2013). An an example, we applied the
density-bounded models shown in Figures 4 and 16 of Plat et al. (2019) to the
ionization cone in Figure 8. Where ${\rm O32}\approx 15$, corresponding to
$F(\lambda 5007)/F(3727)\approx 11$, those density bounded models indicate
$f_{esc}$ of 80% and $\log{\rm U(R_{S})}\approx-3.35$. Making such an
inference from a homogeneous model, however, could be very misleading. For
example, Bassett et al. (2019) fit the empirical relation between $f_{esc}$
and O32 (Faisst, 2016; Izotov et al., 2018) with a series of density-bounded,
photoioniztion models. Although they found a good fit by varying the H I
optical depth, the required metallicity and ionization parameter did not match
the properties of the galaxy sample. Since homogeneous, density bounded models
also predict very weak [O I] emission, the high [O I] $\lambda 6300$ / [O III]
$\lambda 5007$ flux ratio (O1O3) of Lyman continuum leakers presents another,
perhaps related, puzzle (Plat et al., 2019). Two-zone models, where one one
zone is density bounded, offer a possible solution (Ramambason et al., 2020)
which can be tested by resolving the escape channels.

The Ly$\alpha$ profile supports the idea that LyC escape is highly anisotropic
from J1044+0353. Radiative transfer models predict that LyC and Ly$\alpha$
photons escape through the same holes (Kakiichi & Gronke, 2021). The Hubble
COS spectrum of SB reveals a double-peaked Ly$\alpha$ emission-line profile;
however, as is often the case when the spectroscopi aperture is much smaller
than the solid angle of the Ly$\alpha$ nebula, the emission-line profile sits
in a broad, dark H I absorption trough (Hu et al., 2023). Profiles of this
type require a cloudy or multiphase ISM. Hu et al. (2023) fit the H I trough
intensity along with H I damping wings, obtaining an H I column density of
$\log N_{\mbox{H I}}({\rm cm}^{-2})=21.84\pm 0.03$ over $91\pm 3$ percent of
starburst solid angle. In other words, directly along our viewing angle, the
LyC leakage is $\leq\ 9\pm 3\%$. When $\mbox{$f_{esc}$}\leq 0.10$, Ly$\alpha$
peak separations are typically greater than 300 ${\rm km~{}s}^{-1}$ (Izotov et
al. (2018), Figure 10; Gazagnes et al. (2020)). The peak separation in the COS
spectrum, 425 ${\rm km~{}s}^{-1}$, is therefore consistent with little LyC
leakage along the direct sightline, and the asymmetry of the red peak is also
low, $A_{f}=1.84\pm 0.36$ (Hu et al., 2023). The Ly$\alpha$ profile therefore
supports the anisotropic geometry indicated in Figure 8.

Based on the empirical relation (Izotov et al., 2021), we estimate that the
Lyman continuum escape fraction exceeds 5% where the O32 ratio (scaled to our
definition in Section 4.2.1) exceeds 10 in Fig. 8. Perhaps most importantly,
however, our observations strongly suggest that the direction of the bipolar
outflow powered by FEN dictates the direction of the ionization cone powered
by SB, thus underlining the importance of the star formation history. In
future work, we will present a larger sample of O32 maps, addressing the
question of whether J1044+0353 is representative of [O III]-selected, dwarf
galaxies, and then introduce O1O3 maps in order to test two-zone models and
potentially measure $f_{esc}$ indirectly.

##  5 Gas Kinematics in the Starburst Region

In this section we describe the gas kinematics in the starburst region from
multiple perspectives: the ionized gas shells in Sec. 5.1, the blueshifted
resonance absorption lines in Sec. 5.2, and finally the broad emission-line
wings in Sec. 5.3 . Viewed together, these data determine the launch radius
and constrain the mass flux of a nascent wind. We discuss the required
supernova rate in Section 6.2 and the radiative losses from the nascent wind
in Section 6.3.2.

###  5.1 Starburst-driven Superbubbles

In this section, we use the properties of the prominent H$\alpha$ shells in
Figure 4 to estimate the mechanical feedback from the young star clusters. We
predict the expansion velocties of the shells and then use the shell radii and
velocities to quantify their power requirements. We compare our result to the
mechanical feedback predicted by population synthesis models in Section 6.2.

Although there is more energy in the radiation field, supernova are the
dominant feedback mechanism after the shells have emerged from the cloud cores
and the Strömgren radius has propagated through the shells. Thermalization of
the energy from clustered supernovae creates hot bubbles whose thermal
pressure drives shock fronts into the ISM. The shells in Figure 4 consist of
interstellar gas swept up by these shocks. We describe the growth rate of
these shells using the Weaver et al. (1977) model, replacing the stellar wind
by a series of supernova explosions. In a roughly coeval stellar population,
such as a compact star cluster, the supernova rate will be nearly constant for
$\sim 40$ Myr following the first supernova explosion because the slope of the
initial mass function is almost perfectly offset by the age – lifetime
relation of massive stars (Leitherer et al., 1999).

The shells are likely about 2 Myr old. This estimate is the difference between
the upper limit on the starburst age, $4.01\pm 1.13$ Myr (Olivier et al.,
2022), and the 3 Myr delay for the first supernova explosion.666 We note that
this burst age was derived from population synthesis modeling based on BPASS
binary-star models. Starburst 99 models (single stars) indicate a younger age,
$1.04\pm 2.75$ Myr (Olivier et al., 2022), which is marginally too young to
have produced supernova explosions. For a constant rate of mechanical energy
injection, the shell velocity is

$\displaystyle V=100~{}\mbox{${\rm km~{}s}^{-1}$}\
\bigg{(}\frac{R}{340{\rm~{}pc}}\bigg{)}\bigg{(}\frac{2{\rm~{}Myr}}{t}\bigg{)}.$
(12)

In Figure 4, the northwest shell extends 395 pc (1$\farcs$42) from the star
clusters in Region 303, and the southwest shell has a radius of 341 pc
(1$\farcs$12) with respect to the cluster in Region 304. At $t=2$ Myr, the
Weaver et al. (1977) model predicts expansion speeds of 116 ${\rm
km~{}s}^{-1}$ and 100 ${\rm km~{}s}^{-1}$, respectively, for the northwest and
southwest shells.

Combining the equations for the evolution of the shell radius and velocity
(Weaver et al., 1977) gives us an expression for the power requirement:

$\displaystyle L_{w}=1.85\times
10^{41}\mbox{~{}ergs~{}s${}^{-1}$}\bigg{(}\frac{R}{340{\rm~{}pc}}\bigg{)}^{2}\times\
$ (13) $\displaystyle\bigg{(}\frac{V}{100~{}{\rm
km~{}s}^{-1}}\bigg{)}^{3}\bigg{(}\frac{n}{4{\rm~{}cm}^{-3}}\bigg{)}\bigg{(}\frac{\Omega}{4\pi}\bigg{)},$

where $n$ is the hydrogen number density in the ambient medium. We take the
RMS density in the very-high ionization zone, see Section 3.4, as a fiducial
value but discuss this choice further in Sec. 6.2. Following Geen & de Koter
(2022), we include the solid angle $\Omega$ of the shell because the NW and SW
shells emerge from narrow chimineys emanating from the Regions 303 and 304.
Inserting a solid angle $\Omega<4\pi$ in Equation 5.1 reduces the power
requirement linearly. We measure opening angles of approximately $45\deg$ and
$75\deg$ for the northwest and southwest shells, respectively, corresponding
to solid angles of 0.478 and 1.298 steradians. We conclude that the growth of
the northwest and southwest shells requires a mechanical power of
$L_{w}=1.48\times 10^{40}$ ergs s-1 and $L_{w}=1.92\times 10^{40}$ ergs s-1,
respectively, the equivalent of 900 and 1200 supernova explosions over the 2
Myr period.

On galactic scales, however, superbubbles elongate along the direction of the
steepest density gradient (Mac Low & McCray, 1988). Our density measurements,
see Figure 7, indicate that the steepest gradient is to the west of the
starburst region, consistent with superbubble breakthrough in that direction.
In the above argument, we ignored density gradients in the ambient medium
because the seeing-limited KCWI data cube does not resolve the ambient density
profile $n(r)$ on the scale of the superbubble shells resolved by the Hubble
imaging. A shell driven into a medium with an inverse square density profile
does not decelerate, so the predicted expansion speeds would increase to 192
and 166 ${\rm km~{}s}^{-1}$ at 2 Myr. Direct measurements of the expansion
velocities would clearly be valuable. Expansion faster than the predicted
value of 100 to 116 ${\rm km~{}s}^{-1}$ would indicate a steep density
gradient in the ambient medium. Whereas very low shell velocities, or order 10
${\rm km~{}s}^{-1}$, would longer timescales for star formation in the
starburst region.

In the following two subsections, we describe unequivocal kinematic evidence
for strong outflows from the starburst. In the future, mapping infrared
recombination lines with JWST would be the cleanest way to directly measure
shell exansion speeds, as NIRSpec could provide the high spatial and spectral
resolution required.

###  5.2 Ultraviolet Absorption-line Measurements of Outflow Velocity

Absorption in resonance lines is a powerful diagnostic of gaseous outflows and
inflows (Martin et al., 2012; Heckman et al., 2015; Chisholm et al., 2017). We
inspected the Hubble COS UV spectrum of J1044+0353 (Berg et al., 2022; James
et al., 2022). The large red circle in Figure 4 shows that the 2$\farcs$5 (695
pc) diameter aperture encompasses the H$\alpha$ shells. Xu et al. (2022) fit
the Si II $\lambda 1190,1193$, and 1304 profiles with an interstellar
component fixed at $v=0$ and a blueshifted component at $v_{0}=-52\pm 12$
${\rm km~{}s}^{-1}$. The width of the latter is 123 ${\rm km~{}s}^{-1}$ FWHM,
so the outflowing gas reaches speeds of $v_{0}+0.5{\rm FWHM}=114$${\rm
km~{}s}^{-1}$. The maximum speeds therefore match the predicted expansion
speed of the northwest and southwest shells very well. In detail, however,
this comparison is not straightforward because the spatial sampling of the
absorption-line kinematics is determined by the distribution of UV light
behind the shells. In the simplest scenario, projection of the shell expansion
velocity onto our sightline could explain why the line center appears at 52
${\rm km~{}s}^{-1}$. Previously, Xu et al. (2022) attributed the blueshifted
absorption component to a galactic wind. Our analysis demonstrates that some
of these winds are actually expanding shells (Weaver et al., 1977) rather than
the free-flowing winds described by (Chevalier & Clegg, 1985).777 To avoid
confusion, we emphasize that J1044+0353 does have a global wind, which we map
out in Peng et al. (2023). This galactic outflow, however, lies entirely
outside the COS spectrosopic aperture and is driven by an older stellar
population. Consistent with this conclusion, no high-ionization absorption is
detected in Si IV $\lambda\lambda 1393,1402$.

The distinction between superbubbles and a galactic outflow matters because
not all superbubbles blow out of the ISM, a topic we return to in Section
6.3.2; and, by associating the UV absorption with the shells, we determine the
location of the absorbing gas along the line-of-sight. The estimated mass flux
scales linearly with radius, so the detected shells imply mass fluxes orders
of magnitude lower than circumgalactic shells would require.

The UV spectrum constrains the shell column density. We conservatively assume
that Si+ is the dominant ionization state of silicon atoms and ignore
depletion onto dust grains. Then, using the gas-phase metallicity of
0.06${\rm~{}Z_{\odot}}$, we obtain $N({\rm H})\approx 5.2\times
10^{5}(Z/0.06\mbox{${\rm~{}Z_{\odot}}$})^{-1}N({\rm Si~{}II})$, which is
$N(H)\approx 3.6\times 10^{20}$ cm-2 for the Si II column density measured by
Xu et al. (2022). The O I $\lambda 1302$ and C II $\lambda 1334$ line profiles
are similar in shape to the blueshifted Si II absorption troughs, and this
demonstrates that the shell contains some neutral gas. The damped Ly$\alpha$
absorption discussed in Section 4.2.2 has no significant Doppler shift, but
(Hu et al., 2023) identified a second H I component which is blueshifted and
has a lower column density of $\log N_{HI}({\rm cm}^{-2})=20.08\pm 0.44$. We
attribute this second H I component to the expanding shells and emphasize that
it has a very low covering fraction. This may reflect the small area of the
shells in Figure 4 relative to the entire COS aperture.

For spherical shells, the mass flux is well defined at the shell radius:

$\displaystyle\dot{M}(r)=1.8\mbox{${\rm~{}M_{\odot}}$~{}yr${}^{-1}$}\
\bigg{(}\frac{\Omega}{4\pi}\bigg{)}\bigg{(}\frac{N_{H}}{3.6\times
10^{20}{\rm~{}cm}^{-2}}\bigg{)}\times\ ~{}~{}~{}$ (14)
$\displaystyle\bigg{(}\frac{R}{340{\rm~{}pc}}\bigg{)}\bigg{(}\frac{V}{100{\rm~{}km/s}}\bigg{)},$

where $N_{H}$ represents the shell column density (ionized and neutral gas).
If shells like the ones visible in Figure 4 were breaking through the parent
clouds in every direction, then the implied mass flux is $\approx
1.8$${\rm~{}M_{\odot}}$ yr-1. The outflowing mass flux of warm-ionized gas
would be 12 times the aperture SFR from Table 1. Making the correction for the
restricted solid angle of each shell, $\Omega\approx 0.48$ steradian and
$\approx 1.298$ steradian, the two shells combined transport $\dot{M}\approx
0.28$${\rm~{}M_{\odot}}$ yr-1, a value that is roughly twice the (COS)
aperture SFR of 0.15${\rm~{}M_{\odot}}$ yr-1. These mass fluxes demonstrate
significant transport of gas out of the starburst region.

###  5.3 Two Origins of Broad Line Wings

The starburst spectrum in Figure 2 reveals broad wings on the [O III] and
H$\beta$ line profiles. Further inspection of these line profiles across the
KCWI mosaic reveals that this high-velocity [O III] and H$\beta$ emission
comes from a spatially extended region centered on the starburst. In this
section, we characterize the spatial extent, kinematics, and luminosity of
these broad line wings

To map these broad wings, we adaptively binned our [O III] $\lambda 4959$ line
image to SNR $\approx 30$. We extracted a spectrum from each bin and first fit
a single Gaussian profile. To determine whether a second component was
detected, we then fit a double Gaussian profile and applied an F-test, as
defined by Equation 1 in Xu et al. (2022). If the calculated F value was
larger than the theoretical value, a significance level
$\alpha=0.003(3\sigma)$, then the extra component was considered necessary.
The top-right panel of Figure 9 maps this high-velocity gas across J1044+0353.
We detected broad wings in and around the starburst region; these broad wings
are absent in region FEN and the post-starburst region more generally.

Figure 9: Superbubble blowout in the J1044+0353 starburst. The radius of the
red circle is 1$\farcs$53\. Top Left: Hubble H$\alpha$ image. Superbubbles
breakthrough the central gas concentration, reaching 370 pc northwest and 340
pc southwest of the central star cluster. The open filaments beyond these
closed loops are the signature of a superbubble blowout. Top Right: KCWI map
of the [O III] $\lambda 4959$ bins that pass the F-test for a second
component. The surface brightness is the sum of the fitted broad and very-
broad components. Bottom Left: Starburst spectrum of [O III] $\lambda 4959$. A
three-component fit is required: narrow component (97 ${\rm km~{}s}^{-1}$
FWHM), broad component (213 ${\rm km~{}s}^{-1}$ FWHM, $3.17\times 10^{-15}$
ergs s-1 cm-2), and very-broad component (752 ${\rm km~{}s}^{-1}$ FWHM,
$4.5\times 10^{-16}$ ergs s-1 cm-2). The broad and very-broad components are
5.7% and 0.8%, respectively, of the total line flux. Bottom Right: The [O III]
$\lambda 4959$ spectrum extracted along the western filament requires a two-
component fit. The narrow component has a width of 93 ${\rm km~{}s}^{-1}$ FWHM
again. The width of the second component is 698 ${\rm km~{}s}^{-1}$ FWHM, so
we identify it as a physical extension of the very-broad component. The flux
of the very-broad component, $4.1\times 10^{-17}$ ergs s-1 cm-2, is relatively
stronger in the filament, where it accounts for 10% of the [O III] $\lambda
4959$ emission.

The line profiles in the starburst region require three Gaussian components,
all at the systemic velocity. The bottom-left panel of Figure 9 shows the [O
III] fit. Most of the line flux comes from the narrow component; this emission
comes from photoionized clouds in the ISM. A broader component, 213 ${\rm
km~{}s}^{-1}$ FWHM, accounts for almost 6% of the line flux toward the
starburst. This broad component is likely emitted by the expanding shells. The
predicted shell velocities, 100 to 116 ${\rm km~{}s}^{-1}$ in Section 5.1,
correspond to linewidths of 200 to 232 ${\rm km~{}s}^{-1}$ since we would see
emission from their near (approaching) and far (receding) sides.

A two Gaussian fit to the starburst spectrum misses the very-broad wings.
Their fitted width is $752$ ${\rm km~{}s}^{-1}$ FWHM, and we detect this
component at Doppler shifts up to $\pm 1000$ ${\rm km~{}s}^{-1}$. This very-
broad component is not emitted by the expanding shells; its velocity is too
large. In the upper right panel of Figure 9, the spatial mapping of the line
profiles offers important insight about the physical origin of the very-broad
component. To the north, south, and southeast, the very-broad wings are seen
out to 2$\farcs$0 (560 pc), 3$\farcs$3 (920 pc), and 3$\farcs$2 (880 pc),
respectively. This component extends further to the west, up to 4$\farcs$3
(1.20 kpc) from SB. Evidently, one superbubble shell has already accelerated
and fragmented, thereby launching a hot wind (Mac Low & McCray, 1988; Mac Low
et al., 1989). In support of this conclusion, the direction of the western arm
coincides with the steepest density gradient and the elongation of the
superbubble shells. Following Mac Low et al. (1989), who showed that blowout
occurs at roughly three times the fiducial timescale $t_{D}\equiv
H^{5/3}\rho_{0}^{1/3}L_{w}^{-1/3}$, we estimate a minimum age of

$\displaystyle t_{blow}\approx 3t_{D}\approx
3.7~{}{\rm~{}Myr}\bigg{(}\frac{H}{200{\rm~{}pc}}\bigg{)}^{5/3}\times\
~{}~{}~{}$ (15)
$\displaystyle\bigg{(}\frac{10^{40}{\rm~{}ergs}^{-1}}{L_{w}}\bigg{)}^{1/3}\bigg{(}\frac{n}{4{\rm~{}cm}^{-3}}\bigg{)}^{1/3}.$

Seeing-limited data does not resolve the density field on the spatial scale of
the parent gas clouds. We have inserted a fiducial value of 200 pc for the
e-folding scale, $H$, of the density gradient; this value represents the
distance between Region 304 and Region 301. This argument predicts that the
superbubble shells are not only breaking through their parent gas clouds but
are also close to blowout.

The top left panel of Figure 9 draws attention to two fainter H$\alpha$
filaments. These largely radial filaments extend beyond the superbubble
shells, and we suggest they are remnants of ruptured superbubble shells,
relics of superbubble blowout. The KCWI data cube provides kinematic
information on this larger spatial scale, and we show the region where the [O
III] emission line requires a second velocity component in the top right panel
of Figure 9. The odd extension to the west of SB falls within the two radial
H$\alpha$ filaments. Whereas the starburst spectrum, shown in the lower left
panel, has three velocity components, the broad component is absent bewteen in
this western arm; yet the very-broad component persists. The bottom-right
panel of Figure 9 shows our fit to a spectrum extracted along this western
arm. Whereas the very-broad component comprises just 0.8% of the total line
flux in the starburst region, it grows to 10% of the total [O III] flux in the
western arm. We therefore suggest that superbubble blowout, a nascent wind,
produces the very-broad component. The absence of the broad component in the
western arm is consistent with the superbubble shells being outside the
spectroscopic aperture. We estimate that the solid angle of the outflow is at
least 10 % of $4\pi$ steradian. We discuss the implications for launching a
galactic wind further in Section 6.2.

##  6 Discussion

In this section we compare our measurements of the radiative and mechanical
feedback from the J1044+0353 starburst to expectations from stellar population
synthesis modeling. We model the temporal evolution of the ionizing continuum
and the mechanical feedback using the Binary Population and Spectral Synthesis
v2.1 (BPASS) models (Eldridge et al., 2017). In these models mass transfer
between binary stars mixes a larger fraction of the zero-age stellar envelope
into the core, producing longer lived and hotter stars than do traditional
models. Including the evolution of massive-star binaries in evolutionary
synthesis models therefore boosts the He+-ionizing luminosity relative to
population synthesis models based on single stars (Eldridge & Stanway, 2009,
2012). We adopt these models because they produce the hardest ionizing
continua among currently available population synthesis models; but the
systematic uncertainties are substantial, and the soft X-ray emission is not
included.

In Section 6.1, we revisit the spectral hardness problem in J1044+0353 based
on our new measurements of ionizing luminosity from Section 3.3 and discuss
whether stellar populations can produce the ionizing photons. We discuss the
implications of the supernova-driven shells for the starburst properties in
Section 6.2. In Section 6.3.2, we discuss radiative cooling of the hot wind
where the superbubble has blown out of the starburst region.

###  6.1 Ionizing Photon Production Efficiency

The ionizing photon production efficiency, $\xi_{ion}$, compares the
luminosity of the ionizing continuum to the non-ionizing UV continuum. We
expect high $\xi_{ion}$ in J1044+0353 because $\xi_{ion}$ shows a strong
positive correlation with increasing W([O III]) (Chevallard et al., 2018). In
Section 3.1, we measured the ionizing photon luminosity at two energies, 13.6
eV and 54.4 eV. This absorbed photon luminosity is a lower limit on the total
luminosity because a fraction, $f_{esc}$, of the ionizing radiation presumably
escapes from the galaxy.

Dividing the ionizing luminosity by the extinction-corrected, non-ionizing UV
luminosity defines the production efficiency for ionizing photons,
$\mbox{$\xi_{ion}$}(X)\equiv Q(X)/L^{UV}_{\nu},$ where X denotes either the
H-ionizing or the He+-ionizing continuum. This definition of $\xi_{ion}$
implicitly includes the contribution of nebular continuum as well as massive
stars. Defined this way, $\xi_{ion}$ is the parameter needed to determine the
ionization rate of intergalactic gas from the star formation rate density of
the universe and knowledge of $f_{esc}$ (Bouwens et al., 2016). We measured UV
magnitudes for both the Starburst and the entire galaxy in Sec. 2.2, In terms
of the absolute UV magnitude, this expression becomes

$\displaystyle\mbox{$\xi_{ion}$}(X)=2.31\times
10^{24}{\rm~{}Hz~{}erg}^{-1}\bigg{(}\frac{{\rm
Q(X)}}{10^{53}{\rm~{}s}^{-1}}\bigg{)}\times\ ~{}~{}~{}$ (16) $\displaystyle
10^{0.4[M_{UV}+20]}.$

Applying Eqn. 6.1 to the entire galaxy, $M_{UV}=-16.68$ from Table 3, we find
$\log\mbox{$\xi_{ion}$}({\rm Hz~{}erg}^{-1})=25.5$. If we limit measurement to
the UV luminosity of the SB, the ionizing photon production efficiencies
increase to $\log\mbox{$\xi_{ion}$}({\rm Hz~{}erg}^{-1})=25.7$ and
$\log\mbox{$\xi_{ion}({\rm He}^{+})$}({\rm Hz~{}erg}^{-1})=23.6$. For the
harder photons that can ionize He+, Eqn. 6.1 yields $\log\mbox{$\xi_{ion}({\rm
He}^{+})$}({\rm Hz~{}erg}^{-1})=23.5$. In contrast to the growth of W([O III])
with increasing aperture size, the $\xi_{ion}$ and $\xi_{ion}({\rm He}^{+})$
measurements are 1.25 to 1.5 times larger in the Starburst aperture than for
the galaxy as a whole.

When working with population synthesis models, the ionizing photon
luminosities are normalized directly by the stellar continuum. We call this
quantity $\mbox{$\xi_{ion}$}_{*}$ here to distinguish it from $\xi_{ion}$.
Olivier et al. (2022) estimated $\log\mbox{$\xi_{ion}$}_{*}({\rm
Hz~{}erg}^{-1})=25.8\pm 0.9$ by modeling the FUV spectrum of J1044+0353, so
the difference is just $\mbox{$\xi_{ion}$}_{*}-\mbox{$\xi_{ion}$}\approx 0.1$
dex for J1044+0353. We computed $\xi_{ion}$∗ and $\xi_{ion}({\rm He}^{+})$∗
for the BPASS models by integrating the number of photons at energies
exceeding 13.6 eV and 54.4 eV, respectively, and normalizing by the flux
density, $F_{\nu}$ at 1500 Å in the model spectra. In these models, the
production efficiency of H ionizing photons decreases continuously as the
stellar population ages, but the temporal evolution of the He+ ionizing
spectrum depends on stellar metallicity. The tracks for the lowest metallicity
models in Figure 10 show $\xi_{ion}$(He+) rising for several Myr before
beginning a slow decline. The brown track, model bin_z001, corresponds to 12 +
log(O/H) = 7.61 based on Table 2 of Xiao et al. (2018). This track is the
closest match to the gas-phase metallicity of J1044+0353,
$12+\log(O/H)=7.45\pm 0.03$ (Berg et al., 2022).

Figure 10: Ionizing photon production efficiency in J1044+0353 (black symbols)
compared BPASS tracks for a coeval population of binary stars (Eldridge et
al., 2017; Xiao et al., 2018, BPASSv2.1). The stellar metallicities shown
range from $12+\log(O/H)=$ 5.6 (zem5) to 8.52 (z008). Age increases from right
($\log t(yr)=6.0$) to left ($\log t(yr)=7.0$) in steps of 0.1 dex. The model
at the most appropriate metallicity, bin_z001, reaches a maximum
$\xi_{ion}$(He+), $2.2\times 10^{23}$ Hz erg-1, near our measured value but at
a younger age of 2.5 Myr. At the measured age of 4 Myr, the model is full dex
below our measured value. This spectral hardness problem is alleviated by the
models with extremely metal-poor stellar populations.

Figure 10 compares our measurements of the ionizing photon efficiencies to the
BPASS tracks. Model bin_z001 is the nearest match to J1044+0353. Our estimated
$\xi_{ion}$ for J1044+0353 falls within the model range between ages of 3 to 6
Myr. Corrections for escaping Lyman continuum shift our measurements to higher
$\xi_{ion}$, and therefore require younger ages; but the correction is less
than 0.1 dex in $\xi_{ion}$ for $\mbox{$f_{esc}$}\leq 0.21$. The production
efficiency of He+ ionizing photons peaks at an age of 2.5 Myr, where
$\mbox{$\xi_{ion}$}_{*}(\mbox{He${}^{+}$})$, just 70% of our measurement. At
the favored age of 4 Myr, the models are 1.1 dex lower than our measurement.
Correction for the stellar contribution to the He II $\lambda 4686$
luminosity, $4^{+16}_{-4}\%$ in Section 3.1.1, reduces $\xi_{ion}({\rm
He}^{+})$ but leaves a significant excess of hard photons relative to model
bin_z001.

We show two lower metallicity BPASS models in Figure 10 because previous work
suggests that the the excess of He+-ionizing photons in some starburst
galaxies is related to their low metallicity (Chevallard et al., 2018;
Schaerer et al., 2019; Senchyna et al., 2020; Marques-Chaves et al., 2022;
Sander, 2022). The bin_zem5 model with $12+\log(O/H)=5.60$ meets
$\xi_{ion}$(He+) requirement at $t=2.5,5.0,{\rm~{}and~{}}6.3$ Myr. At an age
of 4 Myr there is no significant spectral hardness problem, and a LyC escape
fraction of 16% would be enough to increase $\xi_{ion}$ to the model’s value.
However, such extreme stellar metallicity is not consistent with the stellar
wind and photospheric features previously detected in the FUV spectra of
J1044+0353; they indicate a higher stellar metallicity $Z_{*}=0.14\pm
0.04$${\rm~{}Z_{\odot}}$(Olivier et al., 2022). Super-solar $\alpha/{\rm Fe}$
abundance ratios have also been proposed to explain the hard radiation fields
of young galaxies because iron largely determines the far-UV opacity of stars
(Steidel et al., 2018; Topping et al., 2020; Hirtenstein et al., 2021; Sander,
2022). Interestingly, however, although the photospheric S V $\lambda 1502$
absorption in the J1044+0353 FUV spectrum is weaker than in a comparison
sample, the Fe V absorption in the same spectrum is actually stronger than the
S V absorption (Olivier et al., 2022). It follows that the J1044+0353
starburst does not show the hypothesized super-solar $\alpha/{\rm Fe}$
abuandance ratio. We conclude that BPASS models cannot explain the high
$\xi_{ion}$(He+) in J1044+0353 and confirm the previously recognized spectral
hardness problem (Olivier et al., 2022). In Section 3.1.1, we presented
evidence for VMS in J1044+0353. Including VMS with low-metallicity in stellar
population synthesis models is expected to reduce this tension; whether it
eliminates the discprepancy, however, remains unknown at this time.

###  6.2 Supernova Feedback

The superbubble shells presented in Section 5 demonstrate strong mechanical
feedback from the starburst. In this Section, we discuss the significance of
very young stellar populations producing supernovae and then compare the
number of required supernova to stellar population synthesis models.

Figure 11 illustrates the temporal evolution of the mechanical feedback from
several widely-applied stellar population synthesis models. The delay for the
first supernova is two to three Myr, reflecting a weak dependence on the upper
mass limit. Starburst 99 models fitted to the UV spectrum require a starburst
age of only $1.04\pm 2.75$ Myr (Olivier et al., 2022), and these models
produce the first supernova at an age of 3.3 Myr. Figure 11 demonstrates that
the BPASS binary-star models alleviate this tension. With the same upper mass
limit, the evolutionary tracks used by BPASS produce supernova slightly
sooner, at 2.6 Myr; but, the fitted starburst age increases to $4.03\pm 0.85$
Myr (Olivier et al., 2022). In Section 5.1, we adopted a fiducial starburst
age of 5 Myr and a delay of 3 Myr for the first explosion to obtain 2 Myr for
shell growth. Population synthesis models make assumptions about which massive
stars produce supernovae. The minimum mass for core collapse determines a low-
mass cutoff around 8-10${\rm~{}M_{\odot}}$. Modeling when stellar core
collapse produces a supernova explosion continues to be a challenging problem,
however (Fryer et al., 2023; Heger et al., 2023), as core collapse in stars
with initial masses above 40${\rm~{}M_{\odot}}$ is not always expected to
produce a supernova (Smartt, 2015; Sukhbold & Adams, 2020). Hence, it is not
clear how the presence of very massive stars would change the supernova rate
per unit stellar mass. Feedback observations can therefore help address these
substantial uncertainties regarding the supernova rate from young, metal-poor
stellar populations.

Figure 11:  Mechanical feedback from population synthesis modeling of a
single-age stellar population. Supernova feedback and total feedback including
stellar winds are indicated by black and gray shading, respectively; the
shaded range represents stellar metallicities ranging from roughly 5% to 20%
${\rm~{}Z_{\odot}}$. The ages of the starburst (SB) and the far-eastern
clusters (FEN) are indicated by the vertical shading. The SB region in
J1044+0353 has just begun to produce supernovae, whereas the clusters near the
base of the bipolar outflow have already produced many supernovae.

Starburst 99 population synthesis models run with a Kroupa initial mass
function yield a fiducial mechanical power of $L_{SN}=1.7\pm 0.1\times
10^{40}$ ergs s-1 per $10^{6}\mbox{${\rm~{}M_{\odot}}$}\ $ of stars. Using the
starburst masses from either Table 2 or Table 1, the available power is
$L_{SN}=2.3\times 10^{40}$ ergs s-1 or $L_{SN}=2.6\times 10^{40}$ ergs s-1,
respectively, where we take the difference as an indications of the systematic
uncertainty in the stellar mass estimate. Then using the stellar masses from
Table 2, we estimate that Region 303, at the base of the northwest shell,
produced $L_{SN}=5.9\times 10^{39}$ergs s-1; and Region 304 produced
$L_{SN}=5.8\times\times 10^{39}$ ergs s-1. These values are only 39% and 30%
of our estimated power requirements for the NW and SW shells, respectively, in
Section 5.1. While a stellar population that produces more mechanical power
per unit mass can eliminate this discrepancy, the uncertainties in the
required power allow other solutions.

Those estimated power requirements scale linearly with the ambient density,
for example; see Equation 5.1. In Section 5.1, we adopted the RMS density of
the warm gas as the fiducial ambient density. Below, in Section 6.3.2, we
argue for emergence of a hot phase, a three-phase ISM. The ambient density,
and hence the required mechanical power, would be lower if the hot phase is
the volume-filling phase. Conductive evaporation from the superbubble shells
produces an interior temperture (Weaver et al., 1977)

$\displaystyle T_{b}=1.30\times 10^{7}{\rm~{}K}\bigg{(}\frac{L_{w}}{5\times
10^{39}{\rm~{}erg~{}s^{-1}}}\bigg{)}\bigg{(}\frac{n_{0}}{4{\rm~{}cm^{-3}}}\bigg{)}\times\
~{}~{}~{}$ (17)
$\displaystyle\bigg{(}\frac{2{\rm~{}Myr}}{t}\bigg{)}\kappa_{0}^{-2/7},$

where $\kappa_{0}\approx 1$. A hot phase at a fiducial temperture of
$T_{h}=1\times 10^{7}$ K has of density of $n_{0,h}\approx
0.37{\rm~{}cm}^{-3}(10^{7}{\rm~{}K}/T_{h})$ in pressure equilibrium with in
the warm ($2\times 10^{4}$ K) clouds, where the electron density derived from
the [O II]\- doublet ratio is 185 cm-3. An ambient density of 0.4 cm-3 would
reduce the required mechanical power by a factor of ten, and a normal IMF for
Regions 303 and 304 would produce 4 and 3 times the required mechanical power.
A low ambient density would also reduce the blowout timescale given by
Equation 5.3, just 1.7 Myr after the first supernova for a density $n_{0}=0.4$
cm-3.

###  6.3 Origin of the Very Broad [O III] Wings

The [O III] luminosity of the broad velocity component exceeds the rate of
mechanical energy deposition in J1044+0353. Our schematic places the origin of
this broad velocity component in the superbubble shells, consistent with
photoionization heating this component. The superbubble blowout identified in
Figure 9 is noteworthy because it connects the very-broad velocity component
to the process of launching a galactic wind. The [O III] luminosity of this
very-broad component, which is much lower than that of the broad component,
plausibly comes from cold, outflowing clouds which exchange mass, momentum,
and energy with a hot wind. In Section 6.3.1, we constrain the mass-loading of
that hot wind and discuss whether it cools catastrophically. In Section 6.3.2
we further discuss the luminosity of the very-broad component.

####  6.3.1 Mass Loading of the Hot Wind

In this section we model the linewidth of the very-broad component using a
single-phase wind model that includes radiative cooling. Thompson et al.
(2016) present a picture where entrained clouds seed a cooling instability.
Adding mass to the hot wind lowers the wind temperature. As a hot wind cools
through the temperature range where the cooling coefficient increases, roughly
$2\times 10^{7}>T({\rm K})>2\times 10^{5}$ (chapter 34 Draine, 2011), runaway
radiative cooling can stall the wind.

The mass-loading parameter $\beta$ describes the mass flux of the hot wind in
units of the SFR, $\beta\equiv\dot{M}_{hot}/\dot{M}_{*}$. The Thompson et al.
(2016) model defines the minimum amount of mass loading, $\beta_{cool}$,
required for the cooling radius, where the cooling time drops below the
expansion time, to move inward. Mass loading sufficient to move the cooling
radius inwards all the way to the starburst radius defines a critical mass
loading parameter, $\beta_{crit}$. Writing Equation 10 of Thompson et al.
(2016) in units appropriate to J1044+0353, we have

$\displaystyle\beta_{crit}(r_{cool}=R_{SB})\approx
3.1\alpha^{0.730}\bigg{[}\frac{\zeta}{0.587}\bigg{]}^{0.730}\times\ ~{}~{}~{}$
(18)
$\displaystyle\bigg{[}\frac{\Omega}{4\pi}\frac{R_{SB}}{100{\rm~{}pc}}\frac{0.15\mbox{${\rm~{}M_{\odot}}$~{}yr${}^{-1}$}}{\dot{M}_{*}}\bigg{]}^{0.270},$

where we have included a parameter $\zeta$ to make the aperture SFR give the
same mechanical power as the aperture mass.888 Thompson et al. (2016) adopt a
constant supernova rate equivalent to 1 supernova per unit SFR every 100
years. For a constant SFR, however, the supernova rate rises steadily for the
first 40 Myr, so we cannot apply their fiducial scaling to a very young
starburst. In addition, for consistency with the rest of this paper, where we
have parameterized the mechanical feedback in terms of the burst mass, we
adopted a mechanical power, $L_{SN}=1.76\times
10^{40}{\rm~{}erg~{}s}^{-1}~{}(\zeta/0.587)(\dot{M}_{*}/1.0\mbox{${\rm~{}M_{\odot}}$~{}yr${}^{-1}$})$,
where $\dot{M}_{*}$ is the SFR.

The adiabatic wind model originally described by Chevalier & Clegg (1985)
predicts a maximum wind velocity

$\displaystyle v_{term}=980~{}\mbox{${\rm km~{}s}^{-1}$}\
\alpha^{0.5}\beta^{-0.5},$ (19)

where the parameter $\alpha$ describes the energy deposition rate relative to
the supernova rate. Energy loss through radiative cooling will decelerats the
hot wind. Equation 14 of Thompson et al. (2016) predicts that the terminal
velocity falls to

$\displaystyle v_{term}=440~{}\mbox{${\rm km~{}s}^{-1}$}\
\bigg{(}\frac{\zeta}{0.587}\bigg{)}^{0.135}\bigg{(}\frac{\alpha\xi}{\Omega_{4\pi}}\bigg{)}^{0.135}\times\
~{}~{}~{}$ (20)
$\displaystyle\bigg{(}\frac{\dot{M}_{*}}{0.15~{}\mbox{${\rm~{}M_{\odot}}$~{}yr${}^{-1}$}}\bigg{)}^{0.135}\bigg{(}\frac{100~{}{\rm
pc}}{R_{SB}}\bigg{)}^{0.135},$

in the limit $\beta=\beta_{crit}$. The wind fails if $\beta>\beta_{crit}$, so
this velocity describes the minimum outflow velocity. The metallicity
parameter $\xi\approx 1$ describes a solar metallicity wind, an appropriate
value in a dwarf galaxy (Martin et al., 2002).

In a scenario where the ram pressure of the wind accelerates interstellar gas
clouds, the velocity of the clouds cannot exceed the speed of the hot wind
(Murray et al., 2005). It follows that the terminal velocity of the hot wind
limits the maximum Doppler shift of the [O III] line wings. Hence we can use
observational constraints on $\beta$ and the terminal wind speed in order to
discuss whether strong cooling causes the nascent wind to fail.

In Section 5.2, we estimated the mass flux in the two H$\alpha$ shells. When
these shells undergo blowout, warm gas mass will be entrained in the hot wind
at a rate of 0.28 ${\rm~{}M_{\odot}}$ yr-1. Normalized by the SFR, we estimate
$\beta=0.7-1.9$.999 We obtained an average SFR by dividing the stellar masses
in Table 2 by the corresponding cluster age. Relative to regions 303 and 304,
we obtain a mass loading parameter $\beta\sim 1.1$. Alternatively, normalizing
by the SFR of the entire starburst region gives $\beta\sim 0.7-0.8$. We find a
larger value, $\beta\approx 1.9$, when we use the aperture SFR from Table 1.
The differences between these normalizations define the uncertainty in our
estimate of $\beta$.  These modest mass-loading factors are near unity, the
fiducial value in the adiabatic wind model. In addition, the faint wings of
the [O III] emission line reach $\pm 1000$ ${\rm km~{}s}^{-1}$, a good match
to the terminal wind velocity predicted by Equation 19. We therefore draw two
conclusions: (1) the wind has a hot phase, and (2) some portion of that wind
is nearly adiabatic.

The typical cloud speed is significantly lower than 1000 ${\rm km~{}s}^{-1}$.
The half-width at half maximum intensity of the very-broad component suggests
a speed closer to 380 ${\rm km~{}s}^{-1}$. We combined Equations 6.3.1 and
6.3.1 and found that reducing $v_{term}$ to 380 ${\rm km~{}s}^{-1}$ requires
$\beta_{crit}=3,~{}5,{\rm~{}and~{}}10$ for $\alpha=0.5,1,{\rm~{}and~{}}2$,
respectively. These estimates for the critical mass loading parameter are
independent of the solid angle of the outflow. Since they are significantly
larger than the estimated value of $\beta$, we argue that the nascent wind
will not undergo catastrophic cooling.

Figure 12 illustrates how $\beta_{crit}$ and $\beta_{cool}$ fall with
increasing SFR surface density. We mark the SFR surface density of the
J1044+0353 starburst.101010 The starburst region of J1044+0353 has
$\dot{\Sigma}_{*}\approx 2.39^{+0.32}_{-0.48}$ ${\rm~{}M_{\odot}}$ yr-1 kpc-2
based on the aperture SFR and half-light radius from Berg et al. (2022) and Xu
et al. (2022), respectively. We estimate a lower value value,
$\dot{\Sigma}_{*}\approx 1.0$ ${\rm~{}M_{\odot}}$ yr-1 kpc-2, from the
starburst model in Table 2 and the area of the red circle in Figure 1.  The
shaded region illustrates the range of $\beta$ where our fiducial wind model
cools radiatively but without catastrophic energy losses. The wide horizontal
bar denotes our estimate of the mass loading parameter. The overlap between
the viable models, the horizontal bar, and the starburst SFR/Area selects the
locus of viable models. Catastrophic cooling, in contrast, would require
$\dot{\Sigma}_{*}>30$ ${\rm~{}M_{\odot}}$ yr-1 kpc-2 for the chosen
parameters. Raising $\alpha$ from 1 to 2, which is equivalent to doubling the
number of supernova per unit SFR, would allow SFR surface densities up to
$\dot{\Sigma}_{*}<200$ ${\rm~{}M_{\odot}}$ yr-1 kpc-2 to generate viable
winds.

Figure 12: The minimum mass-loading required to get cooling anywhere in the
wind (solid lines) and the maximum mass-loading beyond which the entire wind
cools (dashed lines), as given by Equations 7 and 10 in Thompson et al.
(2016). Models that lie below the $\beta_{crit}$ line identify viable models
for launching a wind; winds with $\beta>\beta_{crit}$ cool catastrophically.
The negative slope indicates that at fixed $\beta=\dot{M}_{hot}/\dot{M}_{*}$
there exists a critial SFR surface density above which the gas cools too
strongly at blowout to form a steady wind. Our fiducial wind model adopts a
starburst radius $R_{*}=100$ pc, solid angle $\Omega=4\pi$, and
$\alpha\equiv\dot{E}_{hot}/L_{SN}=1$. The intersection of this model with the
SFR surface density of J1044+0353 is shaded orange. Our $\beta$ measurement,
shown by the horizontal bar, has significant overlap with the orange region
where the radiative energy losses are not catastrophic. The thick dashed
(solid) line illustrates how raising (lowering) the supernova rate by settting
$\alpha=2$ ($\alpha=0.5$) raises (lowers) $\beta_{crit}$ ($\beta_{ccool}$).

####  6.3.2 Radiative Cooling of the Hot Wind

Radiative cooling plays a key role in cloud survival, and perhaps growth, as
hydrodynamical instabilities otherwise shred the entrained clouds (Gronke &
Oh, 2018). The single-phase model discussed in Section 6.3.1 describe the wind
with a single temperature at a given radius, so it does not include the mixing
layers between clouds and the hot wind. These mixing layers likely dominate
the total cooling (Fielding & Bryan, 2022). The [O III] luminosity of the
very-broad component provides insight about the cooling mechanism.

The very-broad component has a total, reddening-corrected luminosity
$L(\mbox{[O III]}\ \lambda 4959,5007)\approx 1.1\times 10^{39}$ ergs s-1. We
estimated $L_{SN}\approx 2.64\times 10^{40}$ ergs s-1 from the aperture mass
in Table 1. Since energy is put into the hot wind at a rate
$\dot{E}_{hot}=\alpha L_{SN}$, it follows that the luminosity emitted in the
very-broad wings of the [O III] doublet is $0.04\alpha^{-1}\dot{E}_{hot}$. As
a single-phase wind must cool from an initial temperature $T\sim few\times
10^{7}$ K, first free-free radiation and then emission from highly ionized
metals including O VII, and O VI lines, carry away energy. The gas temperature
must decrease to a $few\times 10^{5}$ K before [O III] becomes an important
coolant, so the the thermal energy has been reduced by a factor of order 100
before [O III] becomes an important coolant. We can confirm this back-of-the-
envelope argument by applying Equation 36 of Thompson et al. (2016). The
result suggests that the energy radiated (primarily in UV and optical lines)
is $L_{UV/O}\sim 7.0\times
10^{38}(\dot{M}/1\mbox{${\rm~{}M_{\odot}}$~{}yr${}^{-1}$})$ ergs s-1 as the
wind cools from $10^{5}$ to $10^{4}$ K. Taking
$\dot{M}=\beta\dot{M}_{*}<5\dot{M}_{*}$, the model predicts a maximum
luminosity $L_{cool}(10^{4}-10^{5}~{}{\rm K})\sim 5.3\times 10^{38}$ ergs s-1,
roughly 2% of $L_{SN}$ as anticipated.

Only a portion of this $L_{cool}$ comes out in the [O III] lines, so the model
requires $L_{\mbox{[O III]}})<L_{cool}$. We examined the (Ploeckinger &
Schaye, 2020) models to gain insight into the fraction of $L_{cool}$ emitted
in [O III] $\lambda 5007$ line. In their models the cooling gas is exposed to
an interstellar radiation field. As the gas cools from $10^{5.5}$ K to
$10^{4}$ K, about 1% of the cooling luminosity is emitted in [O III] $\lambda
5007$. Yet the measured [O III] luminosity of the very-broad component is
comparable to $L_{cool}$ from the single-phase wind model. In future work, we
will further described this luminosity problem using a larger sample of
galaxies with very-broad emission lines. For now, we suggest that the
discrepancy draws attention to the limitations of single-phase wind models and
underlines the importance of multiphase wind modeling. Our results motivate
calculating the fraction $L_{SN}$ that emerges as $L_{\mbox{[O III]}}$ under
realistic conditions in outflows.

Some models for cooling winds predict strong He II emission (Danehkar et al.,
2021), so we revisted our interpretation of the $\pm 400$ ${\rm km~{}s}^{-1}$
He II line wings in Section 3.1.1. The He II line has comparable S/N ratio to
[O III] $\lambda 4363$, so we examined the broad component of [O III] $\lambda
4363$ in order to understand the detectability of broad He II components. The
auroral [O III] line has a broad base. A joint fit with the nebular [O III]
returns the broad component described in Section 5.3. It is therefore
remarkable that a joint fit including He II $\lambda 4686$ produces negative
residuals where a broad He II component might be expected. We therefore
conclude that the gas which emits the broad H$\beta$ and [O III] $\lambda
4959$ components does not contribute to the He II luminosity.

##  7 Summary & Conclusions

In this paper, we presented integral field spectroscopy of the reionization-
era analog, J1044+0353, an [O III]\- selected galaxy with strong He II line
emission. We combined several KCWI pointings into a mosaic cube designed to
cover the galactic Strömgren sphere. At several thousand independent
locations, we analyzed the intensity and profile shape of emission lines in
the 3700 - 5600 Å bandpass and described the global gas kinematics, density
gradients, filling factor, and ionization structure. We leveraged high-
resolution of Hubble images to better understand the morphology of massive
star clusters and superbubbles in the starburst region. We compared the
supernova rates and the ionizing photon productions rates, $\xi_{ion}({\rm
H})$ and $\xi_{ion}({\rm He}^{+})$, predicted by stellar population synthesis
models to our measurements of these quantities in J1044+0353. We summarize the
main results here and draw attention to what they teach us about galaxies
widely believed to have reionized the universe.

##### The Spectral Hardness Problem

Previous studies have concluded that the hardness of the ionizing spectrum in
starburst galaxies increases as metallicity decreases (Senchyna et al., 2017;
Chevallard et al., 2018; Senchyna et al., 2020; Marques-Chaves et al., 2022;
Schaerer et al., 2022; Flury et al., 2022), but whether low-metallicity stars
can supply all the hard ionizing photons remains an active topic of discussion
(Senchyna et al., 2017; Berg et al., 2018; Ravindranath et al., 2020; Saxena
et al., 2020; Olivier et al., 2022; Harish et al., 2023). Our KCWI observation
revealed that the He II $\lambda 4686$ emission from J1044+0353 comes from an
ellipsoid with semi-minor and semi-major axes of 610 and 730 pc, respectively,
centered on a pair of compact star clusters, Region 303. Most of the He II
luminosity has a nebular origin, He+2 recombination, based on the narrow core
of the He II line profile and the spatial coincidence of the He II and [Ar IV]
emission regions. Based on the spatial elongation of this very-high ionization
zone, which follows the locus defined by five massive stars clusters, we
suggest that the source of the He+-ionizing photons is not a single source,
such as an AGN. Normalized by the UV luminosity of the starburst, the ionizing
photon luminosity, Q(He+), corresponds to an ionizing photon efficiency
$\mbox{$\xi_{ion}({\rm He}^{+})$}=23.53$ Hz erg-1, which exceeds the values
produced by BPASS models of an appropriate metallicity and age. Through a
close examination of the He II line profile, which is broader than the [Ar IV]
line profile, we estimated that $4^{+16}_{-4}\%$ of the He II line luminosity
comes from a second component, not of nebular origin. The most likely source
of this broader component is stellar wind emission. Very massive ($>100$
${\rm~{}M_{\odot}}$), hydrogen-core burning stars are the most likely source
of He II wind emission from a stellar population only a few Myr old (Gräfener
& Hamann, 2008; Wofford et al., 2023). The relatively slow stellar wind, i.e.
the width of the second component is just $\pm 200$ ${\rm km~{}s}^{-1}$ from
line center, then indicates the presence of stars with significant rotation,
and thus likely stellar metallicity as low the measured gas-phase metallicity
(Meynet & Maeder, 2002; Gräfener & Vink, 2015). This small stellar
contribution to the He II luminosity does not significantly mitigate the
spectral hardness problem, but it does underline the pressing need for
accurate theoretical modeling of the He+-ionizing continuum from low-
metallicity star clusters.