Resolving the Mechanical and Radiative Feedback in J1044+0353 with Keck Cosmic Web Imager Spectral Mapping

We present integral field spectroscopy toward and around J1044+0353, a rapidly growing, low-metallicity galaxy that 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 λ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.


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 M ⊙ produce Lyman continuum (LyC) photons at a significant rate, so LyC leakage from starforming 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 re-Corresponding author: Crystal Martin cmartin@ucsb.eduliable 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 densitybounded 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-ofsight 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 λ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 × 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 λ1640 (Berg et al. 2016(Berg et al. , 2022)).Recombination of He +2 produces He II λ1640 and optical He II λ4686 emission but requires photons with energies ≥ 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 λ1640 emission from J1044+0353 (Berg et al. 2019(Berg et al. , 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. 2022Mingozzi et al. , 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α 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β 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.
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, ≈ 8−10 M ⊙ , 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 massmetallicity 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 lowdensity 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 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. ′′ 5 (≈ 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  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 star-burst 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 A V = 0.20 ± 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 interface1 .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×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 ZDIST = 0.01363 ± 0.00068 (Willick et al. 1997).We adopt a flat ΛCDM cosmology with Ω Λ = 0.7, Ω m = 0.3, and H 0 = 70 km s −1 Mpc −1 .The luminosity distance is then 59.0 Mpc, and the angular scale is 278 pc/ ′′ .

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 M r = −18 and with NSA stellar masses less than 10 9 M ⊙ .We searched this sample for restframe [O III] λ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. 2014Smit et al. , 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 W([O III] λ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 λ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.

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. ′′ 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. ′′ 35 wide, and the 24 slices subtend an 8. ′′ 4 × 20.′′ 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 ≈ 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.′′ 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. ′′ 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] λλ4959, 5007 over the entire field of view.

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 starfree 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 λ3650 to λ5625 bandpass.The extracted, flux-calibrated spectrum of Feige 34 is consistent with its reference spectrum to within ±3% across this bandpass.We applied this inverse sensitivity function to each cube.The measure-ments reported in this paper use this absolute flux calibration.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 package3 .This procedure produces a mosaic cube with 0. ′′ 146 × 0. ′′ 146 × 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 ≈ 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] λλ4959, 5007 and Hβ emission lines in slits that intersect the spatially unresolved core of the starburst.No slits show saturated [O III] λ4959 or Hβ emission in the short 300 s exposures.Whenever a measurement requires integration over a masked spaxel, we simply measure the Hβ and [O III] λ4959 line flux from the matched mosaic of 300 s cubes.Measurements of the [O III] λ5007 line in the 300 s mosaic are assigned a flux equal to three times the [O III] λ4959 flux, consistent with the respective radiative decay rates of the two transitions.

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. ′′ 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. ′′ 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.
We measured W([O III]) across a set of circular apertures centered on the starburst.In a circular aperture of radius 1. ′′ 5, we measure a rest-frame W(λ5007) of nearly 1400 Å.The measured EW grows rapdily with increasing aperture radius out to a radius of 2. ′′ 2, where it reachs 1780 Å, and then slowly rises to 1930 Å within an aperture of radius of 7. ′′ 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. ′′ 2 slitlet (6-10 kpc in width) would capture most of the nebular emission.However, spectra obtained through narrow 0. ′′ 2 NIRSpec shutters (1 to 1.6 kpc) would yield lower equivalent widths and miss much of the nebular emission.

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 narrowband 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 λ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.

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. ′′ 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 Z ⊙ , Z = 0.05 Z ⊙ , 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×10 5 M ⊙ and 3.2×10 5 M ⊙ , 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.face 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. ′′ 2 by 2. ′′ 56 box around the starburst region, and a 8. ′′ 42 by 4. ′′ 74 box that included most of the UV light from the galaxy.Uncorrected for extinction, the starburst and total magnitudes were m 150 = 18.70 and 18.34, respectively.After correcting for Galactic and internal extinction, see Table 1, the absolute magnitudes are M SB F U V = −16.32± 0.17 and M T ot F U V = −16.68± 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 ±15 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.′′ 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 • 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 • .This edge-on orientation suggests that the radial outflow speed is significantly higher than the line-of-sight outflow speeds, which are 40 -50 km s −1 .

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-meansquare (RMS) electron density.

Recombination Line Images
Prominent emission lines from [Ar IV] λλ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.

He +2 Recombination
Nebular He II λ1640 and λ4686 emission is produced by the recombination of He +2 , and thus demonstrates the presence of hard ionizing photons with energies hν > 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 emissionline 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 M ⊙ (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 λ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 σ tot = 53 km s −1 .After correcting the linewidth for instrumental brodening, σ inst = 35 km s −1 , the standard deviation is σ corr = 40 km s −1 , corresponding to 94 km s −1 FWHM.The [Ar IV] λ4740 emission line profile is also well fit by a single Gauassian component, σ corr = 26 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, and the turbulent velocity dispersion, This interpretation of the broader He II linewidth yields a gas temperature, T = 2.4 × 10 5 K(σ/31 km/s) 2 (m/4 amu).Since this implied temperature is much higher than the measured temperature in the high ionization zone, T e ≈ 2 × 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).
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 λ4713 lines from the data reveals positive residuals extending up to ±200 km s −1 around the He II line center.The He II linewidths of super star clusters in NGC 3125 are 1000 to 1400 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 Z ⊙ , i.e. slightly higher than the gas-phase metallicity of J1044+0353, theoretical models predict He II line widths of just 490 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] λ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] λ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, lowmetallicity 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 Q(He + ) corr = 0.96 +0.04 −0.16 Q(He + ), where Q(He + ) is our measured value below.

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β luminosity.Most of the recombinations occur near the compact starburst.There, the Hβ isophotes are round out to 1.2 kpc (4.′′ 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β surface brightness profile depends on position angle.In projection on the sky, the starburst to shell separation reaches R X = 3.3 kpc (12.′′ 0).The faint shell protruding directly north of SB extends to R = 1.99 kpc (7.′′ 15).To the west of SB, low surface brightness emission extends to at least 2.5 kpc.In the east, the Hβ contours and nearly circular around FEN.The presence of stars hot enough to ionize hy- These loops are shells of gas which the outflow has swept out of its path.

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 λ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. ′′ 0. The surface brightness profile falls steeply between 1. ′′ 0 and 2. ′′ 0. The He I surface brightness declines at a rate similar to the He II emission out to 2. ′′ 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β 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. ′′ 1 in the He I λ4471.The Hβ emission is detected at 10σ at radius 7. ′′ 14, the largst circular aperture with complete coverage.Along the long axis of the mosaic, the Hβ emission is detected to radius 12. ′′ 0.
At many position angles, the Hβ 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β growth curve is very close to zero at these distances, so the filamentary emission would not increase the total observed Hβ 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β.We interactively examined the spectrum across the nebula to determine the maximum extent of the Hβ 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 R C = 610 pc (2.′′ 2) of the compact starburst at all position angles and extends to R X = 730 pc to the northwest and the southeast.

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ν ≥ 13.6 eV and hν ≥ 54.4 eV respectively, from the recombinationline 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β and He II λ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, where (1 − 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 × 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 which evaluates to Q(H) = 5.04 ± 0.48 × 10 52 (1 − f esc ) −1 s −1 for the extinction-corrected luminosity of J1044+0353.
The He + -ionizing continuum at energies ≥ 54.4 eV can be directly measured from the luminosity of the He II λ4686 recombination line.In photoionization equilibrium, we have giving an ionizing photon luminosity Q(He + ) = 4.89 ± 0.47×10 50 (1−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 × 10 4 K to T e = 1.0 × 10 4 K yields the coefficient in Equation 61.23 times higher.We introduce f esc,54 to distinguish the escape fraction at hν ≥ 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 λ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 He 0 , 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).

Root-mean-square Electron Density
The overall size of the Hβ 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, ϵ, of the emitting clouds (Strömgren 1948).For a spherical nebula that volume is 4π 3 R S (H + ) 3 ϵ, and we write the Strömgren radius as in order to show the inverse dependence on the filling factor.We assume ϵ is constant throughout the nebula, then the quantity ϵn 2 e defines the mean-square electron density, ⟨n 2 e ⟩ ≡ ϵn 2 e .
Using Eqn. 4, we can express Q(H) directly in terms of the observed Hβ luminosity, L(Hβ).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, Hβ (H 0 , T ) .We adopted an electron temperature T e = 2.0 × 104 K, appropriate to the O +2 zone in J1044+0353; cooler temperatures of 1.5 × 10 4 K (or 1.0 × 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(Hβ) = 2.54 × 10 40 ergs s −1 and R C (H + ) = 1990 pc. Substitution of these values into Eqn.9 gives a global RMS density of 0.69 cm −3 .The J1044+0353 Hβ emission is strongly concentrated near SB, and roughly 82% of that Hβ luminosity is emitted within the much smaller aperture describing the He +2 Strömgren sphere, R C (He +2 ) = 610 pc.Eqn. 9 therefore indicates a higher RMS density, ⟨n 2 e ⟩ 1/2 = 3.7 cm −3 , for the Hβ emission from the very-high ionization zone.
In close analogy to Eqn. 6, the photon luminosity at energies ≥ 54.4 eV determines the volume of very-highly ionized clouds in photoionization equilibrium.
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. 9.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 where we have evaluated the recombination coefficient at a fiducial temperature of T e = 2.0×10 4 K. 4 The filling factor of the very highly ionized gas, ϵ 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.6, we find n(H)/n(He) 10 1/2 10/12 .
For a He II luminosity 4×10 38 ergs s −1 and R C (He +2 ) = 610 pc, this expression yields an RMS density ⟨n 2 e ⟩ 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 R C (He +2 ) finds different RMS densities for the He II-and Hβ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.

GLOBAL IONIZATION AND DENSITY STRUCTURE
The mosaic cube detects [O III] λλ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  (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.

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.

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 O2DR ≡ F (λ3729)/F (λ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 R C (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 R C (He +2 ) aperture, the [O II] density is 185 +9 −11 cm −3 .In the larger aperture defined by R C (H + ), but still centered on SB, the density is 161 +7 −10 cm −3 .Figure 2 shows the strong nebular [Ar IV] λλ4711, 40 emission.
The [Ar IV] doublet ratio F (λ4711)/F (λ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 R C (He +2 ) in Table 3.
We measured the [Ar IV] line strengths by fitting Gaussian profiles.−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 ([Ar IV]) within R C (He +2 ).
Intermediate ionization zone (24.4 -47.9 eV) densities measured via the [C III] λ1907 / C III] λ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] λλ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.

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 ≈ 90 deg), where it remains constant over more than a kiloparsec, and steepest to the west (PA ≈ 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.

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, ϵ ≡ ⟨n 2 e ⟩ n 2 e . 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 highionization lines are often weak and/or in the rest-UV spectrum.Over the full Strömgren sphere, R C (H + ), our measurements yield ϵ≈ (0.69/161) 2 ≈ 1.8 × 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 ϵ≈ (3.7/185) 2 ≈ 4 × 10 −4 .

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 ionizationstate 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 R C (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.

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] λλ4959, 5007 doublet relative to the [O II] λλ3726, 3729 doublet.For comparison to papers that exclude [O III] λ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. ′′ 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. ′′ 9 and 1. ′′ 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 • remain between 13 and 18 out to projected radii of at least 4. ′′ 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π 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, 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β lines at the edge of the mosaic.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 Z ⊙ 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, and 600 km s −1 , respectively.It follows that the timescale for ionizing gas clouds, τ < 1.5 × 105 yr (r/100 pc)(q/600 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 ≡ q/c, where our O32 values of 15, 10, 5 correspond to log U = −2.36,−2.48, and − 2.70.One systematic uncertainty affecting the interpretation of O32 ratios is model geometry.Spherical models, for example, give ionization parameters ≈ 0.2 dex higher than the plane-parallel geometry adopted. 5An 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.

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 O32 ≈ 15, corresponding to F (λ5007)/F (3727) ≈ 11, those density bounded models indicate f esc of 80% and log U(R S ) ≈ −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 densitybounded, 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] λ6300 / [O III] λ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α profile supports the idea that LyC escape is highly anisotropic from J1044+0353.Radiative transfer models predict that LyC and Lyα photons escape through the same holes (Kakiichi & Gronke 2021).The Hubble COS spectrum of SB reveals a double-peaked Lyα emission-line profile; however, as is often the case when the spectroscopi aperture is much smaller than the solid angle of the Lyα 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 H I (cm −2 ) = 21.84±0.03over 91±3 percent of starburst solid angle.In other words, directly along our viewing angle, the LyC leakage is ≤ 9 ± 3%.When f esc ≤ 0.10, Lyα peak separations are typically greater than 300 km s −1 (Izotov et al. (2018), Figure 10; Gazagnes et al. (2020)).The peak separation in the COS spectrum, 425 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 ± 0.36 (Hu et al. 2023).The Lyα 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.

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.

Starburst-driven Superbubbles
In this section, we use the properties of the prominent Hα 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 ∼ 40 Myr following the first supernova explosion because the slope of the initial mass function is almost perfectly offset by the agelifetime 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 ± 1.13 Myr (Olivier et al. 2022), and the 3 Myr delay for the first supernova explosion.6For a constant rate of mechanical energy injection, the shell velocity is In Figure 4, the northwest shell extends 395 pc (1.′′ 42) from the star clusters in Region 303, and the southwest shell has a radius of 341 pc (1.′′ 12) with respect to the cluster in Region 304.At t = 2 Myr, the Weaver et al. (1977) model predicts expansion speeds of 116 km s −1 and 100 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: 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 Ω of the shell because the NW and SW shells emerge from narrow chimineys emanating from the Regions 303 and 304.Inserting a solid angle Ω < 4π in Equation 14 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×10 40 ergs s −1 and L w = 1.92 × 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 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 km s −1 would indicate a steep density gradient in the ambient medium.Whereas very low shell velocities, or order 10 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.

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. ′′ 5 (695 pc) diameter aperture encompasses the Hα shells.Xu et al. (2022) fit the Si II λ1190, 1193, and 1304 profiles with an interstellar component fixed at v = 0 and a blueshifted component at v 0 = −52 ± 12 km s −1 .The width of the latter is 123 km s −1 FWHM, so the outflowing gas reaches speeds of v 0 + 0.5FWHM = 114km 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 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). 7he 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 lineof-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 Z ⊙ , we obtain N (H) ≈ 5.2 × 10 5 (Z/0.06Z ⊙ ) −1 N (Si II), which is N (H) ≈ 3.6 × 10 20 cm −2 for the Si II column density measured by Xu et al. (2022).The O I λ1302 and C II λ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α 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 (cm −2 ) = 20.08 ± 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: 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 ≈ 1.8 M ⊙ 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, Ω ≈ 0.48 steradian and ≈ 1.298 steradian, the two shells combined transport Ṁ ≈ 0.28 M ⊙ yr −1 , a value that is roughly twice the (COS) aperture SFR of 0.15 M ⊙ yr −1 .These mass fluxes demonstrate significant transport of gas out of the starburst region.

Two Origins of Broad Line Wings
The starburst spectrum in Figure 2 reveals broad wings on the [O III] and Hβ line profiles.Further inspection of these line profiles across the KCWI mosaic reveals that this high-velocity [O III] and Hβ 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] λ4959 line image to SNR ≈ 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 α = 0.003(3σ), then the extra component was considered necessary.The topright 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.
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 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 km s −1 in Section 5.1, correspond to linewidths of 200 to 232 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 km s −1 FWHM, and we detect this component at Doppler shifts up to ±1000 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. ′′ 0 (560 pc), 3. ′′ 3 (920 pc), and 3. ′′ 2 (880 pc), respectively.This component extends further to the west, up to 4. ′′ 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 ≡ H 5/3 ρ 1/3 0 L −1/3 w , we estimate a minimum age of 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α 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α 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π steradian.We discuss the implications for launching a galactic wind further in Section 6.2.A three-component fit is required: narrow component (97 km s −1 FWHM), broad component (213 km s −1 FWHM, 3.17 × 10 −15 ergs s −1 cm −2 ), and very-broad component (752 km s −1 FWHM, 4.5 × 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] λ4959 spectrum extracted along the western filament requires a two-component fit.The narrow component has a width of 93 km s −1 FWHM again.The width of the second component is 698 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 × 10 −17 ergs s −1 cm −2 , is relatively stronger in the filament, where it accounts for 10% of the [O III] λ4959 emission.
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.

Ionizing Photon Production Efficiency
The ionizing photon production efficiency, ξ ion , compares the luminosity of the ionizing continuum to the non-ionizing UV continuum.We expect high ξ ion in J1044+0353 because ξ 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 extinctioncorrected, non-ionizing UV luminosity defines the production efficiency for ionizing photons, ξ ion (X) ≡ Q(X)/L U V ν , where X denotes either the H-ionizing or the He + -ionizing continuum.This definition of ξ ion implicitly includes the contribution of nebular continuum as well as massive stars.Defined this way, ξ 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 Applying Eqn. 17 to the entire galaxy, M U V = −16.68from Table 3, we find log ξ ion (Hz erg −1 ) = 25.5.
If we limit measurement to the UV luminosity of the SB, the ionizing photon production efficiencies increase to log ξ ion (Hz erg −1 ) = 25.7 and log ξ ion (He + )(Hz erg −1 ) = 23.6.
For the harder photons that can ionize He + , Eqn. 17 yields log ξ ion (He + )(Hz erg −1 ) = 23.5.In contrast to the growth of W([O III]) with increasing aperture size, the ξ ion and ξ ion (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 ξ ion * here to distinguish it from ξ ion .Olivier et al. (2022) estimated log ξ ion * (Hz erg −1 ) = 25.8 ± 0.9 by modeling the FUV spectrum of J1044+0353, so the difference is just ξ ion * − ξ ion ≈ 0.1 dex for J1044+0353.We computed ξ ion * and ξ ion (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 ν 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 ξ 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 ± 0.03 (Berg et al. 2022).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 ξion(He + ), 2.2 × 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 ξ 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 ξ ion , and therefore require younger ages; but the correction is less than 0.1 dex in ξ ion for f esc ≤ 0.21.The production efficiency of He + ionizing photons peaks at an age of 2.5 Myr, where ξ ion * (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 λ4686 luminosity, 4 +16 −4 % in Section 3.1.1,reduces ξ ion (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 ξ ion (He + ) requirement at t = 2.5, 5.0, 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 ξ 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 ± 0.04 Z ⊙ (Olivier et al. 2022).Supersolar α/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 λ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 α/Fe abuandance ratio.We conclude that BPASS models cannot explain the high ξ 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 lowmetallicity in stellar population synthesis models is expected to reduce this tension; whether it eliminates the discprepancy, however, remains unknown at this time.

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 com-pare 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 ± 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 ± 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 M ⊙ .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 M ⊙ 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, metalpoor stellar populations.
Starburst 99 population synthesis models run with a Kroupa initial mass function yield a fiducial mechanical power of L SN = 1.7 ± 0.1 × 10 40 ergs s −1 per 10 6 M ⊙ of stars.Using the starburst masses from either Table 2 or Table 1, the available power is L SN = 2.3×10 40 ergs s −1 or L SN = 2.6 × 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 × 10 39 ergs s −1 ; and Region 304 produced L SN = 5.8 × ×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 14.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 evapora-tion from the superbubble shells produces an interior temperture (Weaver et al. 1977) , where κ 0 ≈ 1.A hot phase at a fiducial temperture of T h = 1 × 10 7 K has of density of n 0,h ≈ 0.37 cm −3 (10 7 K/T h ) in pressure equilibrium with in the warm (2 × 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 16, just 1.7 Myr after the first supernova for a density n 0 = 0.4 cm −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.

Mass Loading of the Hot Wind
In this section we model the linewidth of the verybroad 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 × 10 7 > T (K) > 2 × 10 5 (chapter 34 Draine 2011), runaway radiative cooling can stall the wind.
The mass-loading parameter β describes the mass flux of the hot wind in units of the SFR, β ≡ Ṁhot / Ṁ * .The Thompson et al. (2016) model defines the minimum amount of mass loading, β 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, β crit .Writing Equation 10 of Thompson et al. (2016) in units appropriate to J1044+0353, we have , where we have included a parameter ζ to make the aperture SFR give the same mechanical power as the aperture mass.8 The adiabatic wind model originally described by Chevalier & Clegg (1985) predicts a maximum wind velocity where the parameter α 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 , in the limit β = β crit .The wind fails if β > β crit , so this velocity describes the minimum outflow velocity.The metallicity parameter ξ ≈ 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 β 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α shells.When these shells undergo blowout, warm gas mass will be entrained in the hot wind at a rate of 0.28 M ⊙ yr −1 .Normalized by the SFR, we estimate β = 0.7 − 1.9. 9These 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 ±1000 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 km s −1 .The half-width at half maximum intensity of the very-broad component suggests a speed closer to 380 km s −1 .We combined Equations 21 and 19 and found that reducing v term to 380 km s −1 requires β crit = 3, 5, and 10 for α = 0.5, 1, 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 β, we argue that the nascent wind will not undergo catastrophic cooling.
Figure 12 illustrates how β crit and β cool fall with increasing SFR surface density.We mark the SFR surface density of the J1044+0353 starburst. 10The shaded region illustrates the range of β 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 Σ * > 30 M ⊙ yr −1 kpc −2 for the chosen parameters.Raising α from 1 to 2, which is equivalent to doubling the number of supernova per unit SFR, would allow SFR surface densities up to Σ * < 200 M ⊙ yr −1 kpc −2 to generate viable winds.

Radiative Cooling of the Hot Wind
9 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 β ∼ 1.1.Alternatively, normalizing by the SFR of the entire starburst region gives β ∼ 0.7 − 0.8.We find a larger value, β ≈ 1.9, when we use the aperture SFR from Table 1.The differences between these normalizations define the uncertainty in our estimate of β. 10 The starburst region of J1044+0353 has Σ * ≈ 2.39 +0.32 −0.48 M ⊙ 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, Σ * ≈ 1.0 M ⊙ yr −1 kpc −2 , from the starburst model in Table 2 and the area of the red circle in Figure 1. .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 βcrit line identify viable models for launching a wind; winds with β > βcrit cool catastrophically.The negative slope indicates that at fixed β = Ṁhot / Ṁ * 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 Ω = 4π, and α ≡ Ėhot /LSN = 1.The intersection of this model with the SFR surface density of J1044+0353 is shaded orange.Our β 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 α = 2 (α = 0.5) raises (lowers) βcrit (β ccool ).
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, reddeningcorrected luminosity L([O III] λ4959, 5007) ≈ 1.1 × 10 39 ergs s −1 .We estimated L SN ≈ 2.64 × 10 40 ergs s −1 from the aperture mass in Table 1.Since energy is put into the hot wind at a rate Ėhot = αL SN , it follows that the luminosity emitted in the very-broad wings of the [O III] doublet is 0.04α −1 Ėhot .As a single-phase wind must cool from an initial temperature T ∼ f ew × 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 f ew × 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 U V /O ∼ 7.0 × 10 38 ( Ṁ /1 M ⊙ yr −1 ) ergs s −1 as the wind cools from 10 5 to 10 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] λ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 [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 ±400 km s −1 He II line wings in Section 3.1.1.The He II line has comparable S/N ratio to [O III] λ4363, so we examined the broad component of [O III] λ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 λ4686 produces negative residuals where a broad He II component might be expected.We therefore conclude that the gas which emits the broad Hβ and [O III] λ4959 components does not contribute to the He II luminosity.

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 de-signed 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, ξ ion (H) and ξ ion (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 λ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 ξ ion (He + ) = 23.53Hz 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 M ⊙ ), 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 ±200 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.
The Supernova Delay Problem -Simultaneous fits of several line ratios are typically used to determine f esc ; this information breaks degeneracies between ionizationbounded and density-bounded models (Plat et al. 2019) or between density-bounded channels and densitybounded galaxies (Ramambason et al. 2020).We previously mapped the gas-phase oxygen abundance and interstellar reddening in J1044+0353 (Peng et al. 2023).In this paper, we identified a density-bounded ionization cone based on the anisotropic structure seen an F([O III] λλ4959, 5007)/F ([O II]λλ3726, 29) flux ratio map.The O32 ratio is high in the starburst region, declines with radius at most position angles, but remains high throughout a density-bounded ionization cone extending north-northeast.A key part of this argument relied on [O III] detections at large separations, 1.5 -2.5 kpc, from the starburst region where no [O II] emission was detected, yielding a lower limit on O32.
Along the line of sight to J1044+0353, only a small fraction of the starburst is covered by low column density H I (Hu et al. 2023).Hence, we argue that the LyC leakage from J1044+0353 is primarily in the direction of the ionization cone.This anisotropic escape results from the ionizing luminosity, produced in the starburst region, encountering the cavity excavated by a bipolar outflow launched from region FEN.The older stellar age of FEN makes this region a strong source of mechanical feedback even though its current production rate of ionizing photons is negligible compared to the compact starburst.At least in J1044+0353, the star formation history solves the supernova delay problem.
The Catastrophic Cooling Problem -In numerical simulations, galactic winds form through superbubble breakthrough and blowout (Mac Low & McCray 1988;Mac Low et al. 1989).Remarkably, we show that superbubble breakthrough and blowout imprint distinct signatures on the starburst spectrum.Toward the starburst region, the [O III] and Hβ line profiles show a broad (213 km s −1 FWHM) and a very broad (752 km s −1 FWHM) components.We mapped these components individually.The broad component has the higher luminosity, is spatially extended, and is centered on the starburst region.The two most prominent superbubbles within this region extend 350 and 400 pc, respectively, from the compact star clusters labeled Region 304 and 303.Growing these bubbles requires shell velocities v SW ≈ 100(t/2 Myr) −1 km s −1 and v N W ≈ 120(t/2 Myr) −1 km s −1 .The starburst age from Olivier et al. (2022) therefore produces a good match between the shell velocity and the blueshift of UV absorption lines (Xu et al. 2022).We therefore suggest that the outflowing gas identified by Xu et al. (2022) resides in a superbubble shell rather than a multiphase galactic wind; the actual bipolar wind form J1044+0353 is outside their spectroscopic aperture.The superbubble dynamics constrain the supernova rate in star clusters believed to have formed some very massive stars, a population where empirical constraints on supernova rates are needed (Fryer et al. 2023).A filamentary arm extending 1 kpc to the west of the starburst appears only in the lower intensity, very-broad component.Close examination of the Hubble Hα image reveals faint Hα emission on each side of this very-high velocity gas.These features are consistent with superbubble blowout, and the velocity of the very-broad component matches expectations for the velocity of a free-flowing, hot wind (Chevalier & Clegg 1985).These observations suggest a picture where the young starburst is starting to launch a galactic wind.A comparison to the Thompson et al. (2016) wind model demonstrated that the radiative cooling from this nascent wind is not catastrophic.The [O III] luminosity of the wind is higher than single-phase wind models can explain, however; and the discrepancy may provide a useful target for multi-phase wind models.These results aid the interpretation of the emission-line wings detected in 25-40% of redshift 3 < z < 9 galaxies (Carniani et al. 2023;Xu et al. 2023).
The galaxy population believed to drive cosmic reionization had low masses, very low abundances of heavy elements, and high star-formation efficiencies (Ouchi et al. 2009;Wise & Cen 2009;Bouwens et al. 2015b).JWST is finding a large population of such galaxies via their strong [O III] emission (Oesch et al. 2023;Tang et al. 2023;Matthee et al. 2023).In this paper, we have demonstrated that the star formation history across a local analog, J1044+0353, plays a critical role in forming pathways for LyC escape.The LyC photons generated by the very young starburst encounter an ISM whose structure has already been modified by a galactic wind launched 10-20 Myr prior from a post-starburst region.These results support the idea that galaxy orientation determines which strong W([O III]) have directly detected LyC (Tang et al. 2021).We suggest that due to anisotropic LyC escape, spatial mapping of the O32 ratio may turn out to be a more accurate indicator of LyC escape than direct detection of the Lyman cotinuum, which requires an optically thin channel to be aligned with our sightline.
Our schematic picture of J1044+0353 is qualitatively similar to the scenario seen in full-physics simulations of reionization-era galaxies; the SFR is highly variable in time, and LyC escape is anisotropic (Martin-Alvarez et al. 2023;Witten et al. 2024).Bursty star formation histories at redshift z ∼ 10 − 13 also offer an explanation (Muñoz et al. 2023;Sun et al. 2023b,a) for the surprising luminosity excess observed (Naidu et al. 2022;Finkelstein et al. 2023;Casey et al. 2023).The rapid build up of galaxies during the reionization era produces many galaxy -galaxy mergers (Witten et al. 2024).In contrast, it is not obvious that local analogs selected to resemble high-redshift galaxies based on their spectral properties would share this environmental property.Indeed, the stellar mass of J1044+0353 is well below the limits of pair surveys (Stierwalt et al. 2017;Besla et al. 2018), so the probability of a dwarf -dwarf merger in the M * ∼ 10 6 M ⊙ to 10 7 M ⊙ range is not constrained observationally.J1044+0353 has neither a luminous companion nor any known tidal tails, i.e. features that would confirm an interaction or merger hypothesis.We note, however, the unexplained presence of a band of high extinction, elongated perpendicular to a line drawn from SB to FEN (Peng et al. 2023); the coincident string of older star clusters there resembles Mrk 709 (Kimbro et al. 2021).Hence we suggest that a dwarf -dwarf merger triggered the recent star formation activity in J1044+0353.Mergers of dwarf galaxies have been observed to produce more widely distributed star formation than mergers of massive galaxies, which funnel gas into a central starburst (Privon et al. 2017), and we would expect the first passage to produce a starburst in both dwarfs (Kimbro et al. 2021).A merger scenario predicts a positive radial age gradient in the starburst region, similar to the one measured in 30 Doradus (Brandl et al. 1996;Cignoni et al. 2015) and some interstellar gas flung out to large distances (Pearson et al. 2016(Pearson et al. , 2018) ) and possibly detectable via 21-cm mapping.a Scaled to the distance adopted in this paper gives MF U V = −16.23MR = −15.78.
b Assumes a luminosity distance of 55 Mpc.Scale by 1.15 to match the distance adopted throughout this paper.

Figure 5 .
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 ± 1 km s −1 , 125 ± 2 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 λ4713.14lines 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.drogen 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β contours draws attention to the minimum in the Hβ emission between them.The faintest Hβ 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.
suggests.Large spatial gradients in the intensity ratio of the [O III] doublet relative to the [O II] λλ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

Figure 7 .
Figure 7. Radial density profile.Projected distance was computed relative to the compact starburst.Electron density calculated from [O II] λλ3726, 29 doublet ratio at Te = 1.5 × 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 ne(r) within 250 pc (2× FWHM) of SB.

Figure 8 .
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] λ5007 surface brightness (by factors of two from 1.5 × 10 −17 to 4.8 × 10 −16 ergs s −1 cm −2 per square arcsecond.See Section 4.2.1 for further details.

Figure 9 .
Figure 9. Superbubble blowout in the J1044+0353 starburst.The radius of the red circle is 1. ′′ 53.Top Left: Hubble Hα image.Superbubbles breakthrough the central gas concentration, reaching 370 pc northwest and 340 pc 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] λ4959 bins that pass the F-test for a second component.The surface brightness is the sum of the fitted broad and verybroad components.Bottom Left: Starburst spectrum of [O III] λ4959.A three-component fit is required: narrow component (97 km s −1 FWHM), broad component (213 km s −1 FWHM, 3.17 × 10 −15 ergs s −1 cm −2 ), and very-broad component (752 km s −1 FWHM, 4.5 × 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] λ4959 spectrum extracted along the western filament requires a two-component fit.The narrow component has a width of 93 km s −1 FWHM again.The width of the second component is 698 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 × 10 −17 ergs s −1 cm −2 , is relatively stronger in the filament, where it accounts for 10% of the [O III] λ4959 emission.
Figure10.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 ξion(He + ), 2.2 × 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.
Figure12.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 inThompson et al. (2016).Models that lie below the βcrit line identify viable models for launching a wind; winds with β > βcrit cool catastrophically.The negative slope indicates that at fixed β = Ṁhot / Ṁ * 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 Ω = 4π, and α ≡ Ėhot /LSN = 1.The intersection of this model with the SFR surface density of J1044+0353 is shaded orange.Our β 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 α = 2 (α = 0.5) raises (lowers) βcrit (β ccool ).
4 K. Taking Ṁ = β Ṁ * < 5 Ṁ * , the model predicts a maximum luminosity L cool (10 4 − 10 5 K) ∼ 5.3 × 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 [O III] ) < L cool .We examined the (Ploeckinger & Schaye 2020) models to gain insight into the fraction of L cool emitted in [O III] λ5007 line.

Table 3 .
KCWI Measurements of Radiative Feedback