# Haro 11: The Spatially Resolved Lyman Continuum Sources

Lena Komarova Astronomy Department, University of Michigan, 1085 South
University Avenue, Ann Arbor, MI 48109, USA  M. S. Oey Astronomy Department,
University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA
Svea Hernandez AURA for ESA, Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Angela Adamo Department of Astronomy,
Oskar Klein Centre, Stockholm University, AlbaNova University Centre, SE-106
91 Stockholm, Sweden  Mattia Sirressi Department of Astronomy, Oskar Klein
Centre, Stockholm University, AlbaNova University Centre, SE-106 91 Stockholm,
Sweden  Claus Leitherer AURA for ESA, Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  J. M. Mas-Hesse Centro de
Astrobiología – Dept. de Astrofísica (CSIC-INTA), 28850 Madrid, Spain  Göran
Östlin Department of Astronomy, Oskar Klein Centre, Stockholm University,
AlbaNova University Centre, SE-106 91 Stockholm, Sweden  Edmund Hodges-Kluck
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Arjan Bik
Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova
University Centre, SE-106 91 Stockholm, Sweden  Matthew J. Hayes Department of
Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova University
Centre, SE-106 91 Stockholm, Sweden  Anne E. Jaskot Astronomy Department,
Williams College, Williamstown, MA 01267, USA  Daniel Kunth Institut
d’Astrophysique, Paris, 98 bis Boulevard Arago, F-75014 Paris, France  Peter
Laursen Cosmic Dawn Center (DAWN)  Niels Bohr Institute, University of
Copenhagen, Jagtvej 128, DK-2200, Copenhagen N, Denmark  Jens Melinder
Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova
University Centre, SE-106 91 Stockholm, Sweden  T. Emil Rivera-Thorsen
Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova
University Centre, SE-106 91 Stockholm, Sweden

(Accepted March 31, 2024)

###### Abstract

As the nearest confirmed Lyman continuum (LyC) emitter, Haro 11 is an
exceptional laboratory for studying LyC escape processes crucial to cosmic
reionization. Our new HST/COS G130M/1055 observations of its three star-
forming knots now reveal that the observed LyC originates in Knots B and C,
with $903-912$ Å luminosities of $1.9\pm 1.5\times 10^{40}~{}\rm erg~{}s^{-1}$
and $0.9\pm 0.7\times 10^{40}~{}\rm erg~{}s^{-1}$, respectively. We derive
local escape fractions $f_{\rm{esc,912}}$ $=3.4\pm 2.9$% and $5.1\pm 4.3$% for
Knots B and C, respectively. Our Starburst99 modeling shows dominant
populations on the order of $\sim 1-4$ Myr and $1-2\times 10^{7}$ $\rm
M_{\odot}$ in each knot, with the youngest population in Knot B. Thus, the
knot with the strongest LyC detection has the highest LyC production. However,
LyC escape is likely less efficient in Knot B than in Knot C due to higher
neutral gas covering. Our results therefore stress the importance of the
intrinsic ionizing luminosity, and not just the escape fraction, for LyC
detection. Similarly, the Ly$\alpha$ escape fraction does not consistently
correlate with LyC flux, nor do narrow Ly$\alpha$ red peaks. High observed
Ly$\alpha$ luminosity and low Ly$\alpha$ peak velocity separation, however, do
correlate with higher LyC escape. Another insight comes from the undetected
Knot A, which drives the Green Pea properties of Haro 11. Its density-bounded
conditions suggest highly anisotropic LyC escape. Finally, both of the LyC-
leaking Knots, B and C, host ultra-luminous X-ray sources (ULXs). While stars
strongly dominate over the ULXs in LyC emission, this intriguing coincidence
underscores the importance of unveiling the role of accretors in LyC escape
and reionization.

Lyman-alpha galaxies (978) — Starburst galaxies (1570) — Massive stars (732) —
Young massive clusters (2049) — Stellar feedback (1602) — H II regions (694) —
Dwarf irregular galaxies (417) — Radiative transfer (1335) — Intergalactic
medium (813) — Ultraluminous X-ray sources (2164)

††facilities: HST(COS), VLT(MUSE)††software: Starburst99 (Leitherer et al.,
1999, 2014), Astropy (Astropy Collaboration et al., 2022), lmfit (Newville et
al., 2016), Matplotlib (Hunter, 2007)

##  1 Introduction

The ionizing sources and physical mechanisms responsible for cosmic
reionization at $\rm z>6$ remain a critical unsolved problem in cosmology. The
major contenders for providing the required LyC are active galactic nuclei
(AGN) and massive stars in starbursts, with their relative contributions to
reionization still uncertain. Accreting sources other than AGN, e.g. ultra
luminous X-ray sources (ULXs), may also play a significant role in ionizing
the IGM (Madau & Fragos, 2017; Ross et al., 2017; Sazonov & Khabibullin,
2018). Some studies show that AGN could produce sufficient energy to reionize
the universe (e.g., Madau & Haardt, 2015; Giallongo et al., 2015); while
others suggest that the AGN number density and ionizing emissivity were too
low in the early universe (Shankar & Mathur, 2007; Fontanot et al., 2012;
Hassan et al., 2018; Faucher-Giguère, 2020), pointing to star-forming galaxies
as the dominant source. On the other hand, while dwarf starbursts seem to be
promising candidates for reionization agents (Bouwens et al., 2012; Sharma et
al., 2017; Yeh et al., 2023), they may not have sufficiently high escape
fractions (e.g., Fontanot et al., 2014). JWST is now revealing the blue UV
slopes of the earliest, $z=8-16$, galaxies, pointing to young and dust-poor
stellar populations (Cullen et al., 2023a, b; Topping et al., 2023; Morales et
al., 2023), in further support of galaxy-driven reionization. JWST is moreover
uncovering an abundance of high-redshift starbursts with high ionizing photon
production efficiencies, implying that galaxies could have reionized the
universe with somewhat lower escape fractions than previously assumed (Matthee
et al., 2023; Atek et al., 2024).

Knowing the relevant UV sources is only half of the problem, however. The
other key question is how LyC escapes the immediate environment of the source
without being absorbed by the local interstellar medium (ISM). The standard
paradigm for stellar-driven reionization is that supernovae and stellar winds
clear pathways in the ISM that become optically thin to LyC (e.g., Clarke &
Oey, 2002; Fujita et al., 2003; Ma et al., 2016). There is now evidence of a
radiation-driven feedback mode that also may enable LyC escape in the most
extreme metal-poor starbursts such as extreme Green Peas (GPs) (Jaskot et al.,
2017; Komarova et al., 2021; Flury et al., 2022a, b). Similarly, accretion-
driven sources such as AGN or X-ray binaries may create optically thin
channels through winds and jets (Smith et al., 2020).

Haro 11 is an extreme dwarf starburst galaxy with dozens of young massive
clusters (Adamo et al., 2010; Sirressi et al., 2022) and a gas consumption
timescale of 50 Myr (Östlin et al., 2021). It is one of the most important
local ($z=0.021$) galaxies for advancing our understanding of cosmic
reionization, and a wealth of data across the electromagnetic spectrum exist
for this object. Haro 11 is the first local Lyman continuum emitter (LCE) to
be observationally confirmed (Bergvall et al., 2006; Leitet et al., 2011), and
is moreover the closest one, lying at a distance of only 88.5 Mpc (assuming
$H_{0}=73~{}\rm km~{}s^{-1}~{}Mpc^{-1}$; Sirressi et al., 2022, hereafter
S22). With its intense star formation triggered by a dwarf galaxy merger
(e.g., Östlin et al., 2001, 2015), Haro 11 is dominated by three starburst
knots: A, B, and C (Kunth et al., 2003), indicated in Figure 1. While Knot A
is a purely star-forming region, with properties consistent with those of
Green Peas (Keenan et al., 2017), Knots B and C host ULXs (Prestwich et al.,
2015), including a possible low-luminosity AGN (LLAGN) (Gross et al., 2021).
Thus, Haro 11 is a prime laboratory for pinpointing the nature of UV sources
and LyC escape processes crucial to reionization.

However, the initial LyC detection (Bergvall et al., 2006) did not resolve
which of the three knots is/are responsible for the LyC emission, since the
$30\arcsec\times 30\arcsec$ FUSE aperture encompassed the entirety of the
galaxy. Unveiling the exact source(s) of the Haro 11 LyC leakage is necesssary
to clarify the relative roles of compact objects and massive stars in this
prototypical object. The three knots vary drastically in Ly$\alpha$ emission
properties, extinction, stellar populations, and gas properties. So, putting
the LyC detection in context of knot properties reveals the connection between
environment and LyC escape. In this paper, we present HST/COS observations of
the three knots at wavelengths below the Lyman limit, revealing the LyC-
emitting sources. Section 2 contains the details of observations and data
analysis. We present our results in Section 3, discuss cosmological
implications in Section 4, and summarize our main conclusions in Section 5.

Figure 1: HST color composite image of Haro 11, indicating Knots A, B, and C,
as well as the COS apertures (red = 5500 Å continuum, F550M; green = 3360 Å
continuum, F336W; blue = Ly$\alpha$, F122M).

##  2 Observations and Analysis

We obtained HST/COS FUV spectra of Knots A, B, and C in Haro 11 (Cycle 28,
program ID 16260, PI: Oey). We used the medium-resolution G130M grating in
TIME-TAG mode, coupled with the $2.5\arcsec$ Primary Science Aperture (PSA),
with all four grating offset positions (FP-POS = ALL). The observations were
taken with the grating centered on 1055 Å, and resolving power $R\sim 10,000$.
A total of 13 Visits over the course of 24 orbits were obtained, in 2021
January, May, October, November, and December; and 2022 May and July. For each
Visit on each knot, we obtained 4 sub-exposures at different focal plane
offset positions in the wavelength range 900 – 1200 Å. Some sub-exposures had
unusable data due to acquisition failures and resulted in additional Visits.
All sub-exposures used for the LyC measurements are listed in Table 1.

The raw COS observations were calibrated with CalCOS version 3.4.0 (Soderblom,
2021). Given that the observations for the individual knots were obtained
throughout several different Visits, with more than one association file, we
combined their individual calibrated x1d files using the IDL software by
Danforth et al. (2010, 2016). We combined the spectra, weighting by exposure
times and taking into account the data quality (DQ) flags set by the standard
CalCOS pipeline.

A critical step in the calibration of Lyman continuum observations is an
accurate background correction. The standard CalCOS pipeline estimates the
background contribution in the region of the detector where the science
spectrum is located, by computing the average counts in two predefined regions
external to the science extraction region. For the COS FUV detector, these two
predefined regions are typically located below and above the science
extraction box, depending on the lifetime position (LP) used at the time of
the observations. The COS team has reported that the background levels are
correlated with solar activity (Dashtamirova et al., 2019), detector gain and
high voltage (HV). Additionally, they have found that the background levels
vary throughout the detector and with time, with slightly higher levels
towards the edge of the detector.

Table 1: COS G130M/1055 Observations

Target | UT Start Date | Exposure time aaExposure times used in all Visits incorporated in our LyC measurements. (s)
---|---|---
| 2021-05-06 | 2,672
| 2021-05-07 | 3,406
Knot A | 2021-10-13 | 2,306
| 2021-12-01 | 2,306
| 2021-12-01 | 3,862
| 2022-05-20 | 2,306
|  | Total: 16,858
| 2022-05-20 | 4,924
Knot B | 2022-07-05 | 4,712
| 2022-07-07 | 4,713
|  | Total: 14,349
| 2021-01-27 | 4,841
Knot C | 2021-05-12 | 3,500
| 2021-01-25 | 4,847
| 2021-12-01 | 2,306
|  | Total: 11,994

Notes.

To investigate whether the background correction was accurately accounting for
the expected detector contribution in the science extraction regions, we
collapsed the 2-dimensional flt images in the dispersion direction. For the
G130M/1055 configuration observed at Lifetime Position 2, the background
region boundaries are located at pixels $y\sim 448$ and $y\sim 728$, with
predefined widths of 51 pixels. Inspection of the collapsed profiles showed
that in all exposures the background region centered on pixel $y\sim 728$
showed a slightly higher background level towards the edge of the detector
than that observed closer to the science extraction region. This in turn
caused a slight oversubtraction of the science spectra. To improve the
background correction, we modified the location of the background regions in
the extraction tab reference file (XTRACTAB) to be closer to the science
extraction region (centered at $y\sim$ 588), centering the predefined regions
around pixels $y\sim$ 505 and $y\sim$ 663\. The background estimates using
these new regions were $\sim 3-6$% lower than those calculated at the original
background locations.

Given the extended nature of the Haro 11 targets, we adopted a BOXCAR
extraction technique, available in CalCOS. As detailed in James et al. (2022),
for extended targets, a broader extraction box may be necessary to collect the
full extent of their flux. To select an optimal extraction height value for
the science region, we explored varying the nominal size by $\pm$8 pixels, in
steps of 2 pixels, aiming to improve the signal-to-noise of the extracted
spectrum. We adopted extraction heights for the science regions in COS segment
B of 63 (standard), 71, and 71 pixels for Knots A, B, and C, respectively.

We sum the high-quality spectra from each Visit on each knot, with total
effective exposure times for Knots A, B, and C of 4.7, 4.0, and 3.3 hours,
respectively (Table 1). To maximize the signal-to-noise ratio, we moreover bin
the spectra by 6 G130M/1055 resolution elements, and thus 60 pixels $\approx
0.6$ Å.

Our reduced, median-combined COS LyC spectra for the three knots are shown in
Figure 2. The LyC region we observe is a 9 Å window below the redshifted Haro
11 Lyman limit of 931.4 Å, as indicated in Figure 2. We compute the mean LyC
flux densities $F_{912}$ in the three knots, excluding 1.1 Å-wide geocoronal H
Lyman series line regions at 930.75 Å, 923.15 Å, and 926.25 Å. The flux
measurement remains the same within the error, when including these regions.
There are no other known geocoronal features within the measured region. Our
observed LyC fluxes and luminosities are shown in Table 2. We detect the
strongest LyC emission in Knot B, which has about 2/3 of the total flux,
leaving 1/3 in Knot C and none in Knot A. Our combined LyC flux density from
Knots B and C is $F_{912}=(3.5\pm 2)\times 10^{-15}~{}\rm
erg~{}s^{-1}~{}cm^{-2}$ Å-1, consistent with the $(4.0\pm 0.9)\times
10^{-15}~{}\rm erg~{}s^{-1}~{}cm^{-2}$ Å-1 measured by Leitet et al. (2011) in
the $30\arcsec\times 30\arcsec$ FUSE aperture.

Figure 2: New HST/COS Lyman continuum observations of Knots A, B, and C in
Haro 11. We measure the LyC in a 9 Å window between 922.4 Å (the left edge of
the figure) and the redshifted Lyman limit at 931.4 Å, shown as a black dashed
line. The gray regions represent intervals that were excluded from LyC
measurements due to geocoronal emission lines. Figure 3: Top: Combined new COS
G130M/1055 + Sirressi et al. (2022) G130M/1300 + G160M/1600 + MUSE
observations (black) of Haro 11 Knot A, and our model consisting of 3 SSPs,
indicated by the colored lines as shown. The light green regions represent
intervals that were used for fitting the models to the data. The rest was
excluded to mask ISM and nebular lines, geocoronal emission, and detector
gaps. Bottom: Zoom of age-sensitive O VI, N V, C IV, and Si IV P-Cygni profile
fits. Figure 4: Same as Figure 3 but for Knot B. The model assumes two single
stellar populations and an episode of continuous star formation. Figure 5:
Same as Figure 3 but for Knot C. The model assumes one episode of continuous
star formation plus an older background population (see text). Table 2: LyC
Fluxes and Observationally Derived Properties of Knots A, B, and C.

Parameter | Knot A | Knot B | Knot C | Reference
---|---|---|---|---
$F_{912}$ ($10^{-15}$ erg $\rm s^{-1}cm^{-2}$ Å-1) | $<3.35$ | 2.31 $\pm$ 1.8 | 1.20 $\pm$ 0.94 | 1
$L_{\rm 912,obs}$ ($10^{40}$ erg $\rm s^{-1}$) | $<2.82$ | 1.95 $\pm$ 1.5 | 0.92 $\pm$ 0.72 | 1
12 + log(O/H) | 8.42 $\pm$ 0.02 | 8.64 $\pm$ 0.02 | 7.64 $\pm$ 0.01 | 2
$\rm O_{32}$ | 9 | 8 | 1 | 7
$Q_{0}(\rm H\alpha)\tablenotemark{\rm\scriptsize a}~{}(10^{53}~{}\rm s^{-1})$ | $2.3\pm 0.1$ | $8.2\pm 0.4$ | $2.3\pm 0.1$ | 4
$f_{\rm c,Si~{}II}$ | 0.95 $\pm$ 0.05 | 0.96 $\pm$ 0.04 | 0.55 $\pm$ 0.05 | 6
$f_{\rm c,HI}$ | $\sim 1$ | $\sim 1$ | $\sim 0.8$ | 6
$\log(N_{\rm Si~{}II}/\rm cm^{-2})$ | $15.2\pm 0.3\tablenotemark{\rm\scriptsize b}$ | $15.3\pm 0.3$ | $14.7\pm 0.3$ | 6
$\log(N_{\rm HI}/\rm cm^{-2})$ | $20.7\pm 0.1\tablenotemark{\rm\scriptsize b}$ | $21.0\pm 0.1$ | – | 6
$L_{\rm Ly\alpha}$ ($10^{40}$ erg $\rm s^{-1}$) | 3.2 $\pm$ 0.3 | 6.1 $\pm$ 0.6 | 14.0 $\pm$ 1.4 | 6
$f_{\rm esc,Ly\alpha}$ | 1.2 $\pm$ 0.1% | 0.65 $\pm$ 0.06% | 5.9 $\pm$ 0.6% | 6
$v_{\rm sep,Ly\alpha}$ (km $\rm s^{-1}$) | 532 $\pm$ 20 | 409 $\pm$ 20 | 400 $\pm$ 20 | 6
$L_{\rm X}$ [0.3-8 keV] ($10^{40}$ erg $\rm s^{-1}$) | – | 9.5 $\pm$ 1 | 4.5 $\pm$ 1 | 5

References. — (1) This work; (2) Menacho et al. (2021); (3) James et al.
(2013); (4) Sirressi et al. (2022); (5) Gross et al. (2021); (6) Östlin et al.
(2021); (7) Keenan et al. (2017); (8) Menacho et al. (2019)

##  3 Stellar Populations and LyC escape

To quantify the LyC escape efficiency in each region, we compute the local
escape fraction as:

$f_{\rm{esc,912}}=\frac{L_{\rm 912,obs}}{L_{\rm 912,int}}\quad,$ (1)

where $L_{\rm 912,obs}$ is the observed integrated luminosity in our 9 Å
window of 903 – 912 Å at rest (or 922 – 931 Å in the Haro 11 frame), and
$L_{\rm 912,int}$ is the intrinsic luminosity in this range, similar to Leitet
et al. (2011). Thus, $f_{\rm{esc,912}}$ is an approximation for the fraction
of the total produced ionizing power that escapes the region. We assume that
the ionizing luminosity is dominated by FUV radiation from massive stars, and
we estimate $L_{\rm 912,int}$ by modeling the stellar population in each knot
as described below.

Sirressi et al. (2022) have previously constrained the stellar populations of
Haro 11 using multi-band HST photometry in the F140LP, F220W, F275W, F336W,
F435W, F550M, F555W, F665N, and F814W filters (Östlin et al., 2009; Adamo et
al., 2010). Fitting the resulting SEDs, they estimate individual cluster
parameters in each knot. They additionally use FUV $1150-1800$ Å COS
G130M/1300 + G160M/1600 spectra (Östlin et al., 2021, programs 15352, 13017;
PIs Östlin, Heckman) and optical $4650-7000$ Å MUSE spectra (Menacho et al.,
2019, program 096.B-0923(A); PI Östlin). The G130M/1300 spectra were corrected
for broad Ly$\alpha$ absorption as described in Sirressi et al. (2022). As
detailed in Menacho et al. (2019), the MUSE observations of each knot are the
result of 16 dithered integrations at 4 different position angles. The final
MUSE spectra were extracted with apertures of the same size as for the COS
spectra, 2.5″ in diameter, and corrected for vignetting, with the resulting
FUV and optical continua well matched across the wavelength gap.

We combine these existing spectroscopic observations with our newly obtained
COS G130M/1055 spectra ($900-1200$ Å) to model each knot’s stellar content. To
increase the signal-to-noise in the UV, we further bin all COS spectra by a
factor of 2 in wavelength and median-combine the overlap region. Since the
G130M/1055 spectra have different initial extracted spectral sampling from the
G130M/1300 and G160M/1600 data (0.6 and 0.4 Å $\rm px^{-1}$, respectively),
the final combined UV spectrum has sampling of 1.2 and 0.8 Å $\rm px^{-1}$,
over the respective ranges, while the optical MUSE data sampling is 1.2 Å $\rm
px^{-1}$. The combined spectra of the three knots are shown in Figures 3 – 5.
Below we describe our stellar population modeling assumptions, spectral
features of interest, and fitting procedure.

We use Starburst99 (Leitherer et al., 1999, 2014) model spectra of varying
cluster ages ($1-100$ Myr), masses ($10^{5}-10^{9}~{}\rm M_{\odot}$), and
extinctions ($E(B-V)=0$ to 1). We also consider different star formation
histories, specifically, constant star formation (CSF) or single stellar
populations (SSPs). We fit either CSF or a maximum of 3 SSPs to each knot’s
combined COS G130M/1055 + G130M/1300 + G160M/1600 + MUSE spectrum. While the
S22 photometric analysis revealed at least 7 individual star clusters in each
knot, our goal is to approximate the $1-3$ dominant stellar populations within
each region. We evaluate a number of stellar model assumptions by using a
variety of evolutionary tracks. With respect to S22, we extend the parameter
space to include stellar rotation, testing Geneva 2012 tracks with $v=0$,
$v=0.4\times v_{\rm breakup}$ (Ekström et al., 2012a). We also consider
standard or high mass loss: Ekström et al. (2012a) or Meynet et al. (1994)
tracks. For all models, a Salpeter initial mass function (IMF) with stellar
mass range $0.1-120~{}\rm M_{\odot}$ is assumed, similarly to S22. The stellar
mass normalizations depend on the low-mass end of the IMF, and we estimate
that our resulting masses would be $1.6\times$ lower if a Kroupa IMF had been
assumed. We re-sample the model spectra onto our observed wavelength grid,
matching the variable sampling, and convolve with a Gaussian of width $\sim 1$
Å, equal to the stellar broadening we measure from photospheric C III 1247.

Following S22, we perform a linear interpolation of the model spectra between
the discrete metallicity values available in Starburst99
($Z=0.001,0.008,0.02$), in order to estimate models for metallicities
appropriate to Haro 11. Previous Haro 11 metallicity measurements have yielded
discrepant results spanning a factor of two in each knot (Guseva et al., 2012;
James et al., 2013; Menacho et al., 2021). Menacho et al. (2021) provide a
comprehensive summary and discussion of the previous metallicity measurements.
While these measurements relied on the same, direct method, the discrepancies
are attributed to differences in ionization correction factors and aperture
sizes, as well as some possible real variations between stellar and ionized
gas metallicities probed by different apertures. We therefore use the most
spatially resolved metallicity measurements for Haro 11 by Menacho et al.
(2021), adopting their measured central values for Knots A and B. For Knot C,
which shows a much larger age spread, $1-100$ Myr, we treat metallicity as a
free parameter and obtain values consistent with Menacho et al. (2021).
Lastly, we apply extinction to the models, using the SMC law from Prevot et
al. (1984). If the Calzetti et al. (2000) law is assumed, it again results in
unreasonably high extinctions compared to previous measurements, and thus also
stellar masses.

Since our objective is to constrain the youngest, LyC-emitting population, we
focus on fitting the age-sensitive O VI $\lambda$1038, N V $\lambda$1243, C IV
$\lambda$1548, and Si IV $\lambda$1400 P Cygni features, which trace stellar
winds from young stars, as well as the optical and UV continuum, which helps
constrain the extinctions and mass normalizations. While we do not weigh the P
Cygni profiles differently in fitting the spectrum, we reject models that fit
the continuum but not the P Cygni features. We mask out the narrow
interstellar absorption components before fitting, taking special care in the
absorption regions of P Cygni lines. We further isolate feature-free sections
of the UV and optical continuum, and mask out nebular emission lines, ISM
features, and detector gaps before fitting. The spectral regions we use for
fitting are shown in the green sections in Figures 3 – 5. We note that O VI
$\lambda$1038 has contamination from Ly$\beta$ and moreover, the stellar wind
feature in the model atmospheres is especially uncertain due to combined
effects of wind inhomogeneities and X-rays, as well as low signal-to-noise.
For each knot’s spectrum, discrete age combinations from 1 to 100 Myr are fit
in steps of 1 Myr, determining the best-fitting masses and extinctions using
the differential evolution minimizer in the Python package lmfit. The best
model is then chosen as the age, stellar mass $M$, and extinction $E(B-V)$
combination that results in the lowest $\chi^{2}$.

Testing the effects of stellar rotation and assumed mass loss rates, we find
that, for all three knots, the evolutionary tracks that produce the highest-
quality fits are high mass-loss, non-rotating 1994 Geneva models (Meynet et
al., 1994). In these tracks, the mass-loss rate was doubled from the
“standard” case (Schaerer et al., 1993), to match observations. Since 2012, a
lower, theoretical mass-loss prescription (Vink et al., 2001) has been adopted
by the Geneva group accounting for wind inhomogeneities, based on clumping-
corrected mass-loss rates (Ekström et al., 2012b; Georgy et al., 2013). Models
assuming the opposite cases of fast rotation and low mass-loss produce
reasonable fits for Haro 11, but with consistently higher $\chi^{2}$ values.
The preference for non-rotating tracks in Haro 11 suggests the rotation rates
to be at the slower end of the spectrum defined by the two available rates of
$v=0$ and $v=0.4\times v_{\rm breakup}$. Noting that mass loss and rotation
rates require further study in low-metallicity systems, we present the stellar
population fits assuming high mass-loss and zero rotation, similarly to S22.

For each knot’s stellar population fit, we compute $L_{\rm 912,int}$ from the
model spectrum and obtain the escape fraction from equation 1. For the escape
fraction uncertainty, we combine the error in the detected LyC luminosity with
that in the modelled $L_{\rm 912,int}$, which we obtain by Monte Carlo
sampling our best-fit stellar masses and ages $10^{4}$ times. We also
calculate the predicted $Q_{0}$, the emission rate of H-ionizing photons, to
compare to that inferred from H$\alpha$ observations, $Q_{0}(\rm H\alpha)$.
The latter was measured for all knots by Sirressi et al. (2022) from the MUSE
spectra, extracted with the same apertures as in the FUV COS observations.
HST/WFC3 F665N imaging suggests that $Q_{0}(\rm H\alpha)$ is underestimated by
$25\%$, $15\%$, and $15\%$ for Knots A, B, and C, respectively. We also
estimate the predicted thermal FIR luminosity of each knot that arises from
dust-processed optical and UV radiation, by integrating the difference between
the intrinsic and the reddened model spectra in our fitting ranges $912-1750$
Å and $4500-7000$ Å. To compare to observations, we derive the
$1-1000~{}\mu\rm m$ luminosity from the galaxy-integrated IRAS $F_{60\mu\rm
m}$, $F_{100\mu\rm m}$ fluxes and the Helou et al. (1988) prescription,
obtaining $L_{\rm FIR}=2.7\pm 1.3\times 10^{44}~{}\rm erg~{}s^{-1}$. So, the
dust content in Haro 11 generates a substantial FIR luminosity consistent with
luminous infrared galaxies (LIRGs). This value is reasonably consistent with
the sum of our estimated FIR luminosities from the individual knots obtained
below, which is $1.7\times 10^{44}~{}\rm erg~{}s^{-1}$.

Detailed stellar population fits for each knot are presented below, with
resulting parameters shown in Table 3, together with a discussion of the LyC
escape efficiency and mechanisms. We present our results for the three knots
in order of decreasing detected LyC flux.

Table 3: Parameters Fitted from Spectroscopic Stellar Population Models

Region | ModelaaModel for star formation: single stellar populations (SSPs) or continuous star formation (CSF). | Age | $M$ | $E(B-V)$ |  $Z$bbFrom Menacho et al. (2021). | $L_{\rm FIR}$ | $Q_{0}$ | $\log(\xi_{\rm ion})$ | $f_{\rm esc,912}$ |
---|---|---|---|---|---|---|---|---|---|---
|  | Myr | $10^{6}~{}\rm M_{\odot}$ |  |  | $10^{43}~{}\rm erg~{}s^{-1}$ | $10^{53}~{}\rm s^{-1}$ | $\rm Hz~{}erg^{-1}$ |  |
|  | 3$\pm 0.5$ | 7.4$\pm 0.6$ | 0.05$\pm 0.01$ |  | 0.8 | 2.3 |  |  |
Knot A | 3 SSPs | 4$\pm 0.5$ | 13.4$\pm 4.1$ | 0.46$\pm 0.2$ | 0.011 | 1.9 | 2.2 | 25.23 | $<0.10$ |
|  | 13$\pm 3$ | 4.5$\pm 1.2$ | 0.16$\pm 0.05$ |  | 0.1 | $3\times 10^{-3}$ |  |  |
| CSF | 3$\pm 0.3$ | 1.9$\pm 0.4$ | 0.04$\pm 0.02$ |  | 0.3 | 0.7 |  |  |
Knot B | SSP | 1$\pm 0.2$ | 22$\pm 1.5$ | 0.5$\pm 0.1$ | 0.018 | 7.7 | 9.3 | 25.32 | 0.034$\pm 0.029$ |
| SSP | 13$\pm 2$ | 66$\pm 10$ | 0.21$\pm 0.06$ |  | 4.7 | 0.02 |  |  |
Knot C | CSF | 15$\pm 3$ | 26$\pm 4$ | 0.04$\pm 0.01$ | 0.004 | 1.1 | 4.0 | 25.22 | 0.051$\pm 0.043$ |
| SSP |  $>30$ccFits based on both Starburst99 and Yggdrasil (see text). Ages up to $\sim$ 2 Gyr can be fit. | 360$\pm 80$ | 0.47$\pm 0.25$ | 0.001 | 0.7 | $2\times 10^{-4}$ |  |  |

Notes.

Table 4: S22 Photometric Clusters

Region | Age | $M$ | $E(B-V)$ |  $N_{\rm cl}$aa$N_{\rm cl}$ is the number of clusters in each age bin.
---|---|---|---|---
| Myr | $10^{6}~{}\rm M_{\odot}$ |  |
| 4 | 1.4 | 0.06 | 1
Knot A | 5 | 1.7 | 0.11 | 3
| 6 | 3.3 | 0.06 | 3
| 14 | 4.7 | 0.02 | 1
| 15 | 10.8 | 0.07 | 3
| 1 | 8.7 | 0.5 | 2
Knot B | 2 | 7.2 | 0.4 | 1
| 4 | 13.7 | 0.44 | 2
| 5 | 44.1 | 0.4 | 4
| 14 | 1.6 | 0.3 | 1
Knot C | 6, $8-10$ | 2.2 | 0.06 | 4
| 15 | 44.3 | 0.05 | 3

Notes. ††footnotetext: S22 cluster parameters combined in age bins, where
masses are summed for the bin, and extinctions are averaged. The typical
uncertainties on ages, masses, and extinctions are 1 Myr, $0.5\times
10^{6}~{}\rm M_{\odot}$, and 0.02, respectively. The metallicity fixed for the
cluster analysis is $Z=0.004$.

###  3.1 Knot B

We find that the region of Haro 11 with the brightest measured LyC is Knot B.
Using the distance of 88.5 Mpc to Haro 11, we obtain $L_{\rm 912,obs}$ =
$(2~{}\pm~{}1.5)\times 10^{40}~{}\rm erg~{}s^{-1}$ for Knot B (Table 2). This
region by far dominates the H$\alpha$ emission of Haro 11 and thus has the
highest intrinsic ionizing power among the three knots, with $Q_{0}(\rm
H\alpha)=8\times 10^{53}~{}\rm s^{-1}$ (Sirressi et al., 2022). Its young
massive stars should therefore dominate the production of ionizing photons in
Haro 11. However, this object shows large amounts of gas and dust, and it
therefore has been somewhat overlooked in the literature when considering the
origin of the LyC from this galaxy.

Our stellar population synthesis confirms that Knot B produces the most
ionizing radiation in Haro 11. The best-fitting model for the combined COS
G130M/1055 + G130M/1300 + G160M/1600 + MUSE spectrum of Knot B is shown in
Figure 4. A strong match for the spectrum (reduced $\chi^{2}=1.7$), with age-
sensitive P-Cygni wind features well fitted, consists of a continuous star
formation episode for the last 3 Myr, as well as 1-Myr and 13-Myr single
stellar populations. The 1 Myr component is the youngest age observed in Haro
11, as can be seen in Table 3, where detailed stellar population model
parameters are shown for all knots. Moreover, Knot B has the most massive LyC-
bright population, where the stellar mass in ages $<5$ Myr is $2.4\times
10^{7}$ $\rm M_{\odot}$, which is $3\times$ higher than that found in the
second-brightest LyC-leaking region, Knot C. As shown in Table 3, Knot B thus
has the highest ionizing photon production efficiency among the knots,
$\log(\xi_{\rm ion})=25.3$, compared to 25.2 in Knots A and C. Notably, the
dominant, 1-Myr component in Knot B is heavily obscured, showing $E(B-V)=0.5$
and contributing marginally to the observed FUV flux. The age-sensitive P
Cygni features are thus accounted for by the 3 Myr CSF component alone, while
the dusty 1 Myr component contributes $\sim 30\%$ of the optical luminosity.
But with a mass of $\sim 2\times 10^{7}~{}\rm M_{\odot}$, the 1 Myr population
successfully reproduces the high photoionization rate inferred from H$\alpha$.
As can be seen in Tables 2 and 3, our model gives a total $Q_{0}=1\times
10^{54}~{}\rm s^{-1}$, while that inferred from H$\alpha$ is $Q_{0}(\rm
H\alpha)$ = $0.8\times 10^{54}~{}\rm s^{-1}$, which is in good agreement. Our
predicted FIR luminosity, dominated by the 1-Myr component and tracing dust-
processed radiation, is $L_{\rm FIR}=1.3\times 10^{44}~{}\rm ergs^{-1}$. Our
modeling thus confirms the existence of a massive, obscured young population
that dominates the ionizing power and FIR luminosity of Knot B.

Our results can be compared to the S22 photometrically derived cluster
parameters given in Table 4. S22 photometric modelling pointed to 10 clusters
within the COS aperture with ages 1–5 Myr of varying masses and extinctions,
as well as a 14 Myr cluster. Since we fix the number of separate population
components in our spectroscopic analysis to three, our model for Knot B cannot
account for all the S22 cluster ages. But our model is consistent with the two
youngest S22 clusters, 1 Myr and 2 Myr, where our total mass in these ages is
within 40% of the corresponding mass found by S22, and we find similarly high
extinctions of $E(B-V)\sim 0.5$. Notably, while Wolf-Rayet (WR) emission
features (e.g., He I $\lambda$1640, and the $\lambda$5608 Å and $\lambda$4650
Å bumps) are not included in Starburst99 models, these features observed in
the spectrum of Knot B are signatures of a $\gtrsim$3 Myr population,
consistent with the presence of photometrically derived 4 Myr old clusters
(Table 4). Our 3-Myr CSF component thus accounts both for the presence of the
youngest stars, as well as the more evolved WR components. Lastly, the 13 Myr
component we find is $40\times$ more massive than the 14 Myr cluster, which
may be attributed to the spectroscopic aperture also capturing the extended,
diffuse populations.

Given the large cluster masses and young ages seen in Knot B, it could
feasibly host Very Massive Stars (VMS, $\rm M>120~{}M_{\odot}$). We examine
the spectrum of Knot B for VMS diagnostics, such as O V $\lambda$1371
absorption, He II $\lambda$1640 emission with equivalent width (EW) $>3~{}$Å,
an absent or weak double-peaked red bump, and a blue bump without WR lines
(Kunth & Sargent, 1981; Martins et al., 2023; Wofford et al., 2023). We do not
detect O V $\lambda$1371 in Knot B, and we measure its He II $\lambda$1640 EW
to be $2.3$ Å. The red bump detected in Knot B is broad and smooth, while the
blue bump shows WR C IV $\lambda$4658\. The above are all consistent with
classical WR stars, and not VMS.

Although Knot B is the brightest LyC source in Haro 11, its local escape
fraction is low. Based on the intrinsic LyC luminosity of our modelled stellar
components, we obtain $f_{\rm{esc,912}}$= $0.034$$\pm 0.029$ (Table 3). The
uncertainty is dominated by observational error of $77\%$ on the detected
flux, while the uncertainty in the intrinsic LyC luminosity is $21\%$. To
compare to the escape fraction implied by the component cluster parameters, we
model the clusters’ UV spectra in Starburst99, based on their ages, masses,
and metallicities. Summing their intrinsic LyC luminosities in the $903-912$ Å
range, we obtain $f_{\rm{esc,912}}$= $0.02\pm 0.01$ for the clusters,
consistent with our spectroscopic model. Knot B also shows the lowest
Ly$\alpha$ escape fraction among the three knots, $f_{\rm esc,Ly\alpha}=$ 0.65
% (Östlin et al., 2021), although the observed Ly$\alpha$ peak separation of
400 $\rm km\ s^{-1}$ is similar to that observed in Knot C and other weak LyC
leakers with $f_{\rm esc,LyC}\sim 0.03$ (Izotov et al., 2018; Flury et al.,
2022b). The inefficiency of LyC and Ly$\alpha$ escape in Knot B can be
explained by its significant gas and dust reservoir. It exhibits the highest H
I column densities among the three knots, $\log(\rm N_{HI}/cm^{-2})\approx 21$
(Östlin et al., 2021), and the highest molecular gas mass, $M_{\rm H_{2}}$ = 2
x $10^{9}$ $\rm M_{\odot}$ (Gao et al., 2022). Moreover, it has a neutral
covering fraction in Si II, $f_{\rm c,Si~{}II}=0.96$ (Östlin et al., 2021) and
thus $f_{\rm c,HI}\geq 0.96$ in H I (Chisholm et al., 2018), consistent with
the escape fraction we obtain. With filamentary dust clouds across the region
further obscuring it, Knot B shows significant dust extinction of $E(B-V)=0.5$
(Table 3).

Ionization-parameter mapping (IPM), on the other hand, points to a complex
picture of optical depth. Knot B shows a significant [S II]
$\lambda\lambda$6717,6731 deficiency, $\Delta[\rm S~{}II]=-0.16$ (Östlin et
al., in prep), where $\Delta[\rm S~{}II]$ is the displacement in $\log([\rm
S~{}II]/H\alpha)$ from typical star-forming galaxies (Wang et al., 2021),
which may indicate density-bounded conditions (Pellegrini et al., 2012; Wang
et al., 2021). Values of [O II] $\lambda$3727/H$\alpha<0.1$ indicate low
optical depth, and for Knot B, [O II] $\lambda$3727/H$\alpha\sim 0.05$ in the
central line of sight (Keenan et al., 2017), suggests LyC escape. Knot B also
shows elevated $\rm O_{32}$ = [O III] $\lambda$5007/[O II] $\lambda 3727>8$
overall (Keenan et al., 2017), roughly twice the mean value observed in local
unresolved LCEs (Flury et al., 2022a). Higher $\rm O_{32}$ values point to
density-bounded conditions and correlate with the LyC escape fraction (Izotov
et al., 2016a, b; Flury et al., 2022b). However, the extended region shows a
confined, ionization-bounded morphology in transverse directions (Keenan et
al., 2017). Therefore, Knot B must be leaking LyC through a narrow ionization
cone in the line of sight. Menacho et al. (2019) report evidence of narrow
highly ionized channels with Knot B at the base. They also find 1000 km/s
outflows, likely driven by stellar feedback. We measure
$\rm[S~{}II]\lambda\lambda 6716,6731/H\alpha=0.14$ in our spectrum of Knot B,
too low to signal shock heating by supernovae. The S22 photometrically derived
cluster parameters in Knot B, the youngest of which we capture in our
spectroscopic model, nevertheless imply a total power of $2\times 10^{41}$
erg/s in stellar winds and supernovae.

Since the LyC-leaking Knot B hosts an ultra-luminous X-ray source (ULX), its
potential contribution to the LyC emission needs to be evaluated. This is an
unusually bright, hard ULX (Prestwich et al., 2015). The X-ray emission has
been revealed to originate from at least two objects (Gross et al., 2021).
Based on the X-ray hardness, Gross et al. (2021) suggest that one or both of
the sources is a black hole binary in a low-accretion, hard state, with the
high X-ray luminosity suggesting the presence of an intermediate-mass black
hole (IMBH) of mass $\rm M_{\bullet}>7600~{}M_{\odot}$ in the region.
Alternatively, it may be a low-luminosity AGN (LLAGN), whose signatures are
obscured by the dense gas and dust observed in Knot B or diluted by the
intense star formation. It is thus possible for Knot B to be a unique merger
site of two IMBH’s or LLAGN (Gross et al., 2021).

However, we find that the ULX is unlikely to be sufficiently bright to
contribute to the LyC leakage from Knot B. A generous upper limit to the LyC
luminosity of the ULX can be estimated by extrapolating the X-ray power-law
into the UV. Using the observed slope of $\Gamma=1.7$ (Gross et al., 2021), we
estimate the intrinsic $903-912$ Å luminosity of the ULX to be $L_{912,\rm
ULX}\lesssim 6\times 10^{39}\rm~{}erg~{}s^{-1}$. This is fainter than the
observed LyC of Knot B, $L_{912,\rm obs}=2\times 10^{40}\rm~{}erg~{}s^{-1}$,
although it agrees within the observational uncertainty of $\pm 1.5\times
10^{40}\rm~{}erg~{}s^{-1}$. Besides, more realistic ULX SED models predict
even fainter LyC luminosities than we estimate, by 1-2 orders of magnitude
(e.g., Fernández-Ontiveros et al., 2012; Gierliński et al., 2009). Regardless,
the intrinsic stellar LyC emission we constrain from our population synthesis,
$L_{912,\rm int}=6\times 10^{41}\rm~{}erg~{}s^{-1}$, is two orders of
magnitude higher than that estimated for the ULX. So, the stellar population
alone can fully account for the observed LyC leakage from Knot B, with a
relatively insignificant contribution from the ULX.

###  3.2 Knot C

Knot C is the region with the second-strongest LyC flux of Haro 11, with a
detected LyC luminosity $L_{\rm 912,obs}$ = $(0.9\pm 0.7)\times 10^{40}~{}\rm
erg~{}s^{-1}$ (Table 2), or nominally about half of the Knot B luminosity. We
note, however, that the large uncertainties on the LyC fluxes preclude a
conclusive claim on the relative strengths of the emergent LyC from the knots.
Knot C has been the prime candidate for LyC emission from Haro 11 based on its
highest Ly$\alpha$ escape among the three knots, with $f_{\rm esc,Ly\alpha}=$
6 % for Knot C, and $f_{\rm esc,Ly\alpha}\lesssim 1\%$ for A and B (Östlin et
al., 2021). Why does LyC appear fainter than that in Knot B, and what is the
corresponding local escape fraction?

Our best stellar population fit to the combined COS G130M/1055 + G130M/1300 +
G160M/1600 + MUSE spectrum of Knot C is shown in Figure 5, with detailed
parameters in Table 3. It consists of two components: a population
continuously formed for 15 Myr at $\rm SFR=1.7~{}M_{\odot}~{}yr^{-1}$ and a
single older population with age 100 Myr. The model fits the observed spectrum
well (reduced $\chi^{2}=2.4$), with the P Cygni profiles of O VI, N V, C IV,
and Si IV fit well, excluding ISM absorption, and clearly indicating the
presence of $1-5$ Myr-old stars. As can be seen in Tables 2 and 3, our stellar
population fit gives a total photoionization rate $Q_{0}=4\times 10^{53}~{}\rm
s^{-1}$, which is within a factor of two of the value inferred from H$\alpha$.
As Knot C appears to be a nuclear star cluster, a continuous star formation
history is reasonable (Adamo et al., 2010).

The current episode of constant star formation is superimposed on an older,
background population, which is likely the diffuse, extended, bulge-like
component. We allow the metallicities of the components to vary between
observed values, $Z_{\rm obs}=0.004$ from Menacho et al. (2021), and
$Z=0.001$. We note that the old component we find is a generic background
population whose age is not well determined above $\sim 30$ Myr. Because its
contribution to the spectrum is only in the optical continuum, it can be
fitted with a wide range of age–extinction combinations. The high-resolution
Starburst99 UV models are only available up to ages of 100 Myr. We have
therefore also fit low-resolution Yggdrasil SSP models (Zackrisson et al.,
2011) to the spectrum of Knot C and found a good fit with ages up to 2 Gyr and
extinctions $E(B-V)<0.8$, consistent with previously reported values (Sirressi
et al., 2022). Our predicted FIR luminosity from dust processing is $L_{\rm
FIR}=2\times 10^{43}~{}\rm erg~{}s^{-1}$.

Comparing our results to the S22 photometrically derived cluster parameters in
Tables 3 and 4, we see that S22 find Knot C to be dominated by one massive
cluster of age 15 Myr, which is likely the nuclear cluster. The cluster
properties we obtain from the spectrum are consistent with S22 within the
errors. Thus Knot C, undergoing continuous star formation for the last 15 Myr,
has the oldest mass-weighted age in Haro 11, in contrast to the 2 Myr-
dominated Knot B. While there are $1-5$ Myr stars present, their total mass in
our model is $\sim 9\times 10^{6}~{}\rm M_{\odot}$, which is $3\times$ lower
than that in Knot B. Knot C therefore has a lower intrinsic ionizing
luminosity than Knot B.

Our stellar population model gives a local escape fraction from Knot C of
$f_{\rm{esc,912}}$= $0.051$$\pm 0.043$. Here, the uncertainty is dominated by
observational error of $78\%$, while the model uncertainty is $30\%$. The
escape fraction of Knot C may thus be higher than that of Knot B, $0.034\pm
0.029$ (Table 3). So, despite appearing fainter in LyC than Knot B, Knot C may
leak LyC more efficiently. The observed properties of Knot C are also
consistent with it having the highest LyC escape among the knots. In addition
to the highest Ly$\alpha$ escape fraction, Knot C exhibits the lowest neutral
covering fraction $f_{\rm c,Si~{}II}<0.5$, corresponding to $f_{\rm c,HI}\sim
0.8$ (Chisholm et al., 2018) and the lowest neutral column density, as
measured with the apparent optical depth method (Savage & Sembach, 1991)
applied to Si II, $\log(N_{\rm Si~{}II}/\rm cm^{-2})=14.7$ (Table 2) (Östlin
et al., 2021) .

On the other hand, its extinction is comparable to that of Knot B, suggesting
a weak relation between dust extinction and LyC escape. Ionization-parameter
mapping results are likewise ambiguous. On one hand, Knot C shows a strong [S
II] $\lambda\lambda$6717,6731 deficiency of $\Delta[\rm S~{}II]=-0.14$ (Östlin
et al., in prep), suggestive of LyC escape (Wang et al., 2021; Pellegrini et
al., 2012). On the other hand, the knot appears to be in a low ionization
state, based on low $\rm O_{32}\leq 3$ and [O III] $\lambda$5007/H$\alpha\leq
0.5$ and high [O II] $\lambda$3727/H$\alpha\sim 0.3$ (Keenan et al., 2017). As
shown in Figure 5 of Keenan et al. (2017), it appears to have a high-
ionization region extending to the east, but overall, Knot C appears to have
low ionization relative to the rest of the galaxy. It moreover shows a
confined morphology for the high-ionization region, with the $\rm O_{32}$
ratio transitioning smoothly and quickly to lower values into optically thick
envelopes.

Despite some signatures of high optical depth, Knot C has the advantage of
extensive stellar feedback. With a continuous star-formation history of 15 Myr
(Table 3), Knot C has had a significant supernova history. S22 suggest that
the mechanical feedback has been taking place even longer, over the last 40
Myr (Sirressi et al., 2022). The S22 feedback model can account for the
energetics of the observed soft diffuse X-ray emission in Haro 11 reported by
Grimes et al. (2007). We measure $\rm[S~{}II]\lambda\lambda
6716,6731/H\alpha=0.26$ in our spectrum of Knot C, consistent with stellar
photoionization. But Menacho et al. (2019) find that the combination of
$\rm[O~{}I]\lambda 6300/H\alpha$ and $\rm[O~{}III]\lambda 5007/H\alpha$ ratios
on the outskirts of Haro 11 indicates $200-600$ km/s shocks. They also find a
$\sim 2$ kpc, high-ionization structure with Knot C at the center, that is
likely a superbubble. Knot C also shows 1000 km $\rm s^{-1}$ gas (Menacho et
al., 2019), which is difficult to explain with supernovae or stellar winds,
and may instead be a signature of radiation-driven outflows (Komarova et al.,
2021). Overall, the significant stellar feedback in Knot C, whether radiation-
or mechanically dominated, may have cleared optically thin channels in its
ISM, through which LyC photons can escape.

Similar to Knot B, Knot C contains a ULX that, based on its $0.3-8.0$ keV
spectrum and luminosity, might be one of the most luminous soft ULXs known,
with $L_{\rm X}=4.5\times 10^{40}~{}\rm erg\ s^{-1}$ (Prestwich et al., 2015;
Gross et al., 2021). The X-ray emission likewise originates in two point
sources, where the secondary object shows $L_{\rm X}\sim 2\times 10^{40}~{}\rm
erg\ s^{-1}$ (Gross et al., 2021). Prestwich et al. (2015) and Gross et al.
(2021) explain its high luminosity by an IMBH of mass $\rm
M_{\bullet}>20~{}M_{\odot}$, undergoing super-Eddington accretion (Swartz et
al., 2011; Kaaret et al., 2017), although it is also possible to be a neutron
star binary, since compact object masses are poorly constrained. In turn, the
blowout of inner-disk material in this intense accretion phase can result in a
disk wind (e.g., Middleton et al., 2015). If the super-Eddington accretion
drives an outflow with a mechanical luminosity similar to its X-ray luminosity
(e.g., Justham & Schawinski, 2012), the outflow power would be comparable to
the stellar wind power of Knot C, estimated to be $6\times 10^{40}~{}\rm erg\
s^{-1}$ (Sirressi et al., 2022). So, the ULX feedback may contribute
significantly to gas clearing, promoting LyC escape.

To understand the role of the ULX in the LyC leakage from the knot, we
estimate the ULX LyC output by extrapolating its X-ray power-law into the UV.
For the observed spectral index $\Gamma=2.1$, we estimate $L_{912,\rm ULX}\leq
2\times 10^{40}\rm~{}erg~{}s^{-1}$. This is a generous upper limit, as
Vinokurov et al. (2013) and Kaaret & Corbel (2009) show that more realistic
model SEDs of supercritical accretion disks are even fainter in the UV. Our
estimated ULX LyC luminosity is twice as bright as the observed LyC luminosity
from Knot C, and thus the ULX can plausibly contribute LyC. However, our
modelled stellar LyC of $2\times 10^{41}~{}\rm erg\ s^{-1}$ exceeds that of
the ULX by at least an order of magnitude. Thus, if the ULX is responsible for
some of the LyC emission from the region, its contribution may be unimportant
compared to that of the stellar population. However, if its mechanical
feedback dominates the gas clearing for the observed LyC, its LyC contribution
may still be significant.

So, the stellar population of Knot C can alone account for the observed LyC
leakage from this region. The soft, luminous ULX observed in this knot may be
sufficiently bright to contribute to the LyC emission, but its intrinsic
production is $<10$% of the stellar emission. The ULX may be able to aid in
the escape of Lyman radiation through mechanical feedback. Further
investigation of the mechanical and radiative feedback of the ULXs is needed
to conclusively establish their roles in the LyC leakage from Knot C.

###  3.3 Knot A

Knot A has also been predicted to be the LyC-leaking knot, based on
ionization-parameter mapping, which Keenan et al. (2017) use to show that Knot
A is responsible for a large, $\sim$kpc-sized region with high O32. They find
a central $\rm O_{32}\sim 9$, consistent with the most extreme Green Peas,
which are the largest class of local LyC leakers (e.g., Flury et al., 2022a).
However, we do not detect LyC in Knot A. The $2\sigma$ upper limit of the LyC
flux density from Knot A is $F_{\rm 912}<3.3\times 10^{-15}~{}\rm erg\
s^{-1}\rm cm^{-2}$ Å-1 in the range $903-912$ Å, or $L_{\rm 912,obs}<2.8\times
10^{40}~{}\rm erg\ s^{-1}$.

To estimate the local escape fraction upper limit, we use our population
synthesis results. Our stellar model for the spectrum of Knot A consists of
three SSPs of ages 3 Myr, 4 Myr, and 13 Myr, as shown in Figure 3 and detailed
in Table 3. The model fits the spectrum reasonably well (reduced
$\chi^{2}=3.7$), accounting for both the continuum and age-sensitive features.
The P Cygni profiles of O VI, N V, C IV, Si IV clearly indicate the presence
of 3 Myr-old stars. Although the Wolf-Rayet He I $\lambda$1640 and 5608 Å and
4650 Å bumps are not included in Starburst99 models, these features are
consistent with a $\sim$4 Myr population. With $1\times 10^{7}~{}\rm
M_{\odot}$, this 4 Myr component dominates the stellar mass of Knot A, but is
highly obscured, with $E(B-V)=0.5$. It therefore has no FUV contribution,
including in the age-sensitive P Cygni features, and the 3 Myr component can
alone account for these profiles, while the dusty population contributes
$25\%$ of the optical luminosity. The 3 Myr and 4 Myr populations reproduce
the observed photoionization rate from H$\alpha$ within a factor of two
(Tables 2 and 3). However, given that Knot A has the least dust among the
knots, with previously reported extinction $E(B-V)\sim 0.2$ (Menacho et al.,
2021), this discrepancy may indicate a modest mass overestimate in the 4-Myr
component. Comparing with the S22 clusters in Table 4, our 4-Myr spectroscopic
component is $10\times$ more massive than the 4-Myr S22 cluster, but has an
$8\times$ higher extinction. Similar to the 14 Myr cluster found by S22, we
find a 13 Myr component, with a matching mass but $2\times$ higher extinction.
Lastly, the total FIR luminosity predicted by our model for the knot is
$2.8\times 10^{43}~{}\rm erg\ s^{-1}$.

The massive young populations we uncover in Knot A may possibly host VMS. This
was previously suggested as an explanation for Knot A’s broad blue bump
(Keenan et al., 2017). Similar to Knot B, we evaluate the spectrum of Knot A
for VMS signatures. We do not detect the blue-shifted O V $\lambda$1371
absorption common in VMS, and we measure the He II $\lambda$1640 EW = 2 Å,
while VMS show $>3$ Å (Martins et al., 2023). The observed red bump is broad
and smooth, while the blue bump shows WR C IV $\lambda$4658\. So, as in Knot
B, we do not find evidence of VMS in Knot A. The above observations are
instead consistent with classical WR stars.

From our modelled intrinsic LyC luminosity of the stellar populations, we
estimate a 2$\sigma$ upper limit to the escape fraction
$f_{\rm{esc,912}}$$\leq 0.10$ (Table 3). Our non-detection thus does not
conclusively rule out significant LyC escape in Knot A. Indeed, multiple lines
of evidence point to some degree of LyC escape in Knot A. First, despite
showing the lowest observed Ly$\alpha$ luminosity among the knots, $L_{\rm
Ly\alpha}=3\times 10^{40}~{}\rm erg\ s^{-1}$, the Ly$\alpha$ escape fraction
of Knot A is $f_{\rm esc,Ly\alpha}=1.2\pm 0.12$% (Östlin et al., 2021), which
is $2\times$ that of Knot B (Table 2). Nevertheless, this value is lower than
that suggested by its reddening, implying that Ly$\alpha$ is strongly
attenuated by dust scattering (Östlin et al., 2021). The neutral covering
fraction is similar to that of Knot B, $f_{c}=0.95\pm 0.05$ (Östlin et al.,
2021), and so is the column density, log($\rm N_{HI}/cm^{-2}$) = 20.7 (Östlin
et al., 2021).

Moreover, Knot A shows a number of signatures of low optical depth in
ionization parameter mapping. Its [S II] deficiency is most extreme among the
knots, $\Delta[\rm S~{}II]=-0.19$ (Östlin et al., in prep), predicting the
highest $f_{\rm esc,LyC}$ (Wang et al., 2021). Most importantly, Keenan et al.
(2017) report high-$\rm O_{32}$ gas originating at Knot A and spreading to
distances of $>$ 2 kpc, and Menacho et al. (2019) observe it to $>$ 4 kpc as
high [O III]/H$\alpha$, signaling a transparent medium in the transverse
directions. In fact, the $\rm O_{32}$ ratio does not transition smoothly into
lower values at the edges of this region, strongly implying low optical depth
in the plane of the sky. Menacho et al. (2019) show that this structure
exhibits the highest ionization, $\rm[O~{}III]/H\alpha\gtrsim 3$, in velocity
bins from $-300~{}\rm km\ s^{-1}$ to $-150~{}\rm km\ s^{-1}$, while the
central $\rm[O~{}III]/H\alpha\sim 1$. This suggests that the high-ionization
structure is likely an optically thin outflow driven by LyC from the knot, but
its axis is not coincident with our line of sight.

The main implication of our results in light of IPM observations is therefore
that LyC escape must be highly anisotropic. While we do not detect LyC in Knot
A, it is likely that the leakage, if any, occurs through a channel not
coincident with our line of sight (e.g., Zastrow et al., 2011). This is
consistent with the Ly$\alpha$ peak velocity separation, which is highest in
Knot A, $v_{\rm sep,Ly\alpha}\sim 500$ km $\rm s^{-1}$, compared to the other
two knots. Such large $v_{\rm sep,Ly\alpha}$ is consistent with escape
fractions $<1\%$ (Flury et al., 2022b).

##  4 Discussion

We find that the LyC-emitting regions in Haro 11 are Knots B ($L_{\rm
912,obs}$ = $2.3\pm 1.8\times 10^{40}~{}\rm erg~{}s^{-1}$) and C ($L_{\rm
912,obs}$ = $0.9\pm 0.7\times 10^{40}~{}\rm erg~{}s^{-1}$). We determine their
respective LyC escape fractions to be $f_{\rm{esc,912}}$$=3.4\pm 2.9\%$ and
$5.1\pm 4.3$% (Table 3). The total LyC-luminosity-weighted escape fraction of
Haro 11 is $f_{\rm esc,912}=3.9\pm 3.4\%$, consistent with $3.3\pm 0.7$%
obtained by Leitet et al. (2011). Knot B appears to dominate in LyC
luminosity, with the caveat that the low signal-to-noise in our LyC
observations prevents a conclusive determination of which knot dominates. At
face value, Knot B is responsible for $2/3$ of the total observed flux, and it
has $\sim 3\times$ higher mass in 1 – 5 Myr-old stars than Knot C (Table 3),
which is the age of peak LyC production. Thus Knot B strongly dominates the
ionizing photon production. However, it has almost full neutral covering
(Östlin et al., 2021) and high extinction (Table 3). In comparison, Knot C is
a more evolved region of constant star formation for the last 15 Myr, with a
commensurate history of supernova feedback, possible LLAGN feedback, and low
neutral covering fraction (Sirressi et al., 2022; Östlin et al., 2021). Our
findings highlight the sensitivity of LyC escape to the star-formation rate,
age, and optical depth.

Thus, Knot B both has more young stars and emits more strongly in LyC, but has
a slightly lower escape fraction than Knot C. Our results underscore the fact
that LyC escape fraction and escaping LyC luminosity are separate quantities.
In Haro 11, Knot B seems to produce the greatest LyC emission because it
strongly dominates in LyC production. But its large gas and dust content
apparently keeps its local escape fraction lower than the more evolved, and
cleared, Knot C. On the other hand, Knot A shows similar age and extinction to
Knot B, but we do not detect it in LyC. It is thus the interplay of star
formation intensity, age, gas clearing, and/or line-of-sight orientation that
determines the efficiency and detection of LyC escape. Tracers of
$f_{\rm{esc,912}}$ alone are insufficient indicators of escaping LyC
luminosity along the line of sight.

Flury et al. (2022b) connect LyC properties to star formation density,
Ly$\alpha$ properties, and $\rm O_{32}$ in 66 local star-forming galaxies,
including Green Peas. They find possible evidence for two modes of LyC escape
in the most extreme starbursts, dictated by whether stellar feedback is wind-
dominated or radiation-dominated. One population shows younger ages, higher
$\rm O_{32}$, and low metallicity, suggesting strong radiation-dominated
feedback, which may be linked to optically thin, radiation-driven winds
(Komarova et al., 2021). On the other hand, the somewhat older starbursts with
lower $\rm O_{32}$ and higher metallicity likely leak LyC with the help of
superwinds (e.g., Heckman et al., 2011; Zastrow et al., 2013), and these show
lower $f_{\rm esc,LyC}$. If Knot A is indeed an anisotropic LyC emitter as
evidenced by ionization-parameter mapping (Keenan et al., 2017), then it falls
into the radiation-dominated category, consistent with its very high
ionization parameter, and other radiation-dominated features (Section 3.3;
Keenan et al., 2017). Knot B may likewise be radiation-dominated based on its
high $\rm O_{32}$ and young age, with the caveat that its metallicity is close
to solar. Finally, Knot C has low $\rm O_{32}$ and is apparently dominated by
mechanical feedback, given its extensive supernova history (Sirressi et al.,
2022). Since we observe multiple stellar generations in each region, it is
likely that both modes are at play, with their relative importance to be
established.

Haro 11 is a local Green Pea analog, showing $\rm O_{32}$ and other radiation-
dominated properties characteristic of these objects (Micheva et al., 2017;
Keenan et al., 2017). These properties are linked to Knot A, which turns out
to not show direct detection in LyC, although IPM strongly implies that the
region is optically thin in other sightlines. This additionally stresses the
dependence of the escape fraction on the geometrical distribution of dust and
neutral gas along the line of sight. The fact that the detected LyC from Haro
11 originates from regions other than Knot A further demonstrates that LyC
emission from Green Peas may be more complex than consideration of a single
starburst (e.g., Micheva et al., 2018). Notably, the global $\rm O_{32}$ for
Haro 11 is only 2.5 (James et al., 2013), which is on the low side for a GP.
This $\rm O_{32}$ is consistent with that of unresolved GPs in the Low-
redshift Lyman Continuum Survey (LzLCS), containing the largest sample of
local LCEs to date. The LzLCS GPs show $f_{\rm esc,LyC}=1-4\%$ for similar
$\rm O_{32}$. In addition, the ionizing photon production efficiencies we
estimate are $\log(\xi_{\rm ion})=25.32$ for Knot B and 25.22 for Knots A and
C. These values are lower than that of reionization-era galaxies, $25.4-25.8$
(Simmonds et al., 2023; Saxena et al., 2023; Atek et al., 2024), or local
strong LyC leakers, $25.6-26$ (Schaerer et al., 2016). But they are consistent
with the standard values assumed in reionization models (Robertson et al.,
2015).

Our spatially resolved study thus uncovers how some of the global properties
we observe at higher redshifts may arise from multiple star-forming regions of
widely differing properties. In particular, the regions dominating the
ionization may not be the primary LyC sources in our line of sight, and more
evolved regions’ contribution should not be discounted.

###  4.1 Implications for $Ly\alpha$ as a Diagnostic of LyC

As LyC emission is difficult to observe directly, indirect tracers are
required to identify LyC leakers. Ly$\alpha$ is expected to correlate with LyC
escape, as it is also sensitive to the hydrogen column density. Radiative
transfer simulations (e.g., Verhamme et al., 2015) and observations of several
LCEs (Verhamme et al., 2017) confirm a tight relationship between Ly$\alpha$
and LyC escape. The LzLCS survey shows that Ly$\alpha$ width and peak
separation are some of the strongest indirect predictors of LyC escape
fractions (Flury et al., 2022b). With their new, larger sample, the authors
reproduce the anti-correlation between Ly$\alpha$ peak velocity separation
$v_{\rm sep,Ly\alpha}$ and LyC escape fraction, $f_{\rm esc,LyC}$, first
established by Izotov et al. (2018).

Our findings in Haro 11 are consistent with these Ly$\alpha$ predictions. As
shown in Table 2, the Ly$\alpha$ peak separations observed in Knots A, B, and
C are 530, 409, and 400 $\rm km~{}s^{-1}$, respectively (Östlin et al., 2021),
decreasing with increasing $f_{\rm esc,LyC}$. Knots B and C are consistent
with the Flury et al. (2022b) $v_{\rm sep,Ly\alpha}-f_{\rm esc,LyC}$ relation,
which predicts $f_{\rm esc,LyC}\sim 0.03$ for $v_{\rm sep,Ly\alpha}=400~{}\rm
km~{}s^{-1}$. The two LCE knots are also consistent with a $f_{\rm
esc,Ly\alpha}-f_{\rm esc,LyC}$ correlation, where Knot C has a higher $f_{\rm
esc,Ly\alpha}$ and a slightly higher $f_{\rm esc,912}$. We calculate the
luminosity-weighted average Ly$\alpha$ escape fraction for the three knots to
be $f_{\rm esc,Ly\alpha}=4.2\pm 0.6$%. We also estimate the global Ly$\alpha$
peak velocity separation by combining the Ly$\alpha$ profiles of the three
knots, obtaining $v_{\rm sep,Ly\alpha}=410\pm 70~{}\rm km\ s^{-1}$. Our
luminosity-weighted LyC escape fraction $f_{\rm esc,912}=3.9\pm 3.4$% is
consistent with the averaged $v_{\rm sep,Ly\alpha}$ according to the Flury et
al. (2022b) relation. These values account only for Ly$\alpha$ emission from
the knots, and not diffuse Ly$\alpha$ observed outside of them (Östlin et al.,
2009).

The efficiency of Ly$\alpha$ escape, or its escape fraction, can provide
insight into LyC radiative transfer. Knot C has both the highest Ly$\alpha$
luminosity and $f_{\rm esc,Ly\alpha}$ among the knots, and its $f_{\rm
esc,912}$ is indeed likely the highest despite the LyC luminosity being only
half of that in Knot B. So, while the LyC and Ly$\alpha$ escape fractions
correlate, the emerging LyC and Ly$\alpha$ luminosities do not necessarily do
so.

As for the shape of the Ly$\alpha$ profiles, Knot A and B clearly show broader
Ly$\alpha$ red peaks than Knot C, with their respective widths $302~{}\rm km\
s^{-1}$, $338~{}\rm km\ s^{-1}$, and $195~{}\rm km\ s^{-1}$, suggesting higher
optical depth and low LyC escape (Östlin et al., 2021). The Ly$\alpha$ red
peak width specifically points to Knot B as the weakest leaker, inconsistent
with our results. Since both Knots A and B show neutral covering fractions
close to unity (Table 2; Östlin et al., 2021), it appears that the correlation
of Ly$\alpha$ red peak width with $f_{\rm esc,LyC}$ is not as strong when
considering individual star-forming regions instead of integrated galaxy
properties. The Ly$\alpha$ profile of Knot C, on the other hand, is narrow,
but with a higher red peak asymmetry than in Knots A and B (Rivera-Thorsen et
al., 2017; Östlin et al., 2021). Kakiichi & Gronke (2021) show that higher red
peak asymmetry around its center points to LyC escape through a hole-ridden
ISM, or a picket fence structure, as the asymmetry is seen to correlate with
ISM porosity in their radiation-hydrodynamic simulations. This arises because
the asymmetry of the red peak traces the presence of both optically thin and
thick channels, while a symmetric red peak indicates isotropic leakage.

One inconsistency we see in Ly$\alpha$ predictions is in the non-detected Knot
A, where its peak velocity separation is highest among the knots, but its
Ly$\alpha$ escape fraction is twice that of Knot B. The peak velocity
separation thus implies the lowest LyC $f_{\rm esc,LyC}$ among the knots,
while the Ly$\alpha$ escape fraction suggests it should be higher than that in
Knot B. A likely explanation for this may be that the LyC is emerging in
directions transverse to our line of sight, as implied by ionization-parameter
mapping (Section 3.3), while Ly$\alpha$ can be enhanced by scattering into our
line of sight. The observed Ly$\alpha$ luminosity of Knot A is notably the
lowest among the knots, and the intrinsic LyC we estimate is $3\times$ lower
than that of Knot B (Table 3), while exhibiting a similar covering fraction,
column density, and thus similarly broad red Ly$\alpha$ peak. Lastly, it is
important to note that the Ly$\alpha$ relations described above were
established for unresolved galaxies, which likely also consist of multiple
star-forming regions with varying properties. In our spatially resolved study,
we connect Ly$\alpha$ profiles from smaller, $\sim$1-kpc apertures to
individual knot properties, providing a view of separate components that may
make up the integrated Ly$\alpha$ observations.

Although most of our results agree well with Ly$\alpha$ predictions, there are
still significant differences in the escape conditions for LyC vs.
Ly$\rm\alpha$. The maximum column density at which Ly$\rm\alpha$ can escape is
log($N_{\rm HI}$/$\rm cm^{-2}$) $<$ 13, while the threshold for LyC is
log($N_{\rm HI}$/$\rm cm^{-2}$) $<$ 17\. So, the four orders of magnitude of
difference provide a parameter space where gas can be optically thin in LyC
but not in Ly$\rm\alpha$. Another major difference in the escape mechanisms is
scattering: Ly$\rm\alpha$ scatters strongly, modulating its escape path and
additionally promoting dust absorption relative to LyC. Dijkstra et al. (2016)
simulate the radiative transfer of Lyman radiation in multiphase ISM to
investigate the relationship between LyC and Ly$\rm\alpha$ escape fractions.
They find a positive correlation, as expected, but with significant scatter at
higher Ly$\rm\alpha$ escape fractions that is driven by gas covering fraction.
The corresponding LyC escape fractions in this region are lower than expected
from the correlation, decreasing with higher covering fractions. The scatter
extends at least two orders of magnitude, consistent with observations (Flury
et al., 2022b), showing that the ISM porosity introduces appreciable
stochasticity to the relationship between LyC and Ly$\rm\alpha$ radiation.
This can also provide a context for the Ly$\alpha$ vs LyC observations of Knot
A.

Thus, Ly$\alpha$ properties can provide clues to LyC escape conditions, though
not without additional independent tracers. Our Haro 11 study shows that
Ly$\alpha$ peak velocity separation is consistent with it tracing $f_{\rm
esc,LyC}$, where the knots fall on the observed relation within scatter. But
the Ly$\alpha$ luminosity, red peak width, and escape fraction do not
correlate directly with LyC escape in regions of varying gas optical depth and
covering. Our results underscore the significant distinctions in Ly$\alpha$
and LyC radiative transfer, and that further study into the effects of ISM
morphology and anisotropy of LyC escape is needed.

###  4.2 IPM and Anisotropy of LyC Escape

Ionization-parameter mapping relies on nebular emission line ratios with
different ionization potentials to find regions of low optical depth, where
high-ionization species dominate. This serves as an indirect tracer for LyC
escape (Pellegrini et al., 2012).

Our resolved observations of the LyC-emitting regions of Haro 11 are not fully
consistent with predictions from IPM, if isotropic escape is assumed. Most
importantly, in Knot A, IPM demonstrates LyC escape conditions transverse to
the line of sight (Keenan et al., 2017, ; Section 3.3). Yet we find no LyC
detection in the COS aperture, which probes the line of sight. As noted above,
this suggests that LyC escape may be extremely anisotropic. This is also
implied from our observations of Knot B, which leaks LyC despite being almost
fully covered in neutral gas, implying a narrow escape path. Knot C, on the
other hand, is seen to be in a low-ionization state, implying high optical
depth in LyC. Yet we find it to have the highest LyC escape fraction. Thus the
IPM predictions appear to be sensitive to line-of-sight effects.

Hydrodynamic simulations of starbursts at $z=4-6$ by Cen & Kimm (2015) confirm
that the LyC escape fraction depends strongly on the viewing angle. They find
that only highly ionized, evacuated channels with small solid angles allow
significant LyC propagation. Indeed, the observed scatter in LyC escape
fractions with respect to correlated starburst properties at $z\sim 3$ and
$z<0.4$ points to a line of sight effect, where partial ISM clearing results
in limited transparent paths (Flury et al., 2022b; Nestor et al., 2011).
Several nearby starbursts with active stellar feedback likewise exhibit narrow
ionization cones (Zastrow et al., 2011, 2013). The line-of-sight bias is thus
crucial to account for in reionization studies at all redshifts. Anisotropy of
LyC escape, dictated by the gas-clearing mechanisms such as winds and
supernovae, as well as the optical depth and ionization structure, need to be
accounted for when using IPM as a predictor for LyC escape. This may be even
more important if accretion-driven feedback is responsible for the necessary
gas clearing.

Thus, our observations of Haro 11 offer important data on what information IPM
does and does not provide on LyC escape. While it gives insight into the
ionization structure in the plane of the sky, additional tracers probing line-
of-sight conditions are required to identify objects in which LyC can be
directly detected.

###  4.3 X-ray Sources and LyC Escape

The two LyC-leaking Knots B and C both host ULXs, while the undetected Knot A
is purely star-forming. This interesting coincidence raises the question of
the role of accretors in LyC escape, which has been a major problem in cosmic
reionization (e.g, Volonteri & Gnedin, 2009; Madau & Haardt, 2015). Although
the ULX LyC contribution is likely negligible in Knot B, the ULX in Knot C,
which may be an LLAGN, may contribute to the LyC leakage from the region at
the $\lesssim 10$% level. On the other hand, the focused mechanical feedback
that can be expected from these X-ray sources might contribute to clearing
optically thin channels for LyC escape in both Knots B and C. Low-accretion,
hard X-ray sources as that seen in Knot B may produce jets of substantial
mechanical power (e.g., Merloni & Heinz, 2013; May et al., 2018). Likewise,
the soft, super-Eddington sources, such as the one hosted by Knot C, can also
drive radiation-driven disk winds that may be important (e.g., Middleton et
al., 2015). These feedback mechanisms may significantly enhance ULX
contributions to LyC escape.

The question of the role of accretors in cosmic reionization is all the more
important in light of the observed excess X-ray emission in early-universe
galaxy analogs. For instance, Kaaret et al. (2011) observe an increased X-ray
luminosity per star formation rate $L_{\rm X}/\rm SFR$ in blue compact dwarfs
(BCDs), and Brorby et al. (2014) find the BCD X-ray luminosity function to
have $10\times$ the normalization observed in solar-metallicity galaxies.
Also, Douna et al. (2015) find $10\times$ more high-mass X-ray binaries
(HMXBs) per $\rm SFR$ bin in $<0.2~{}Z_{\odot}$ galaxies than in solar-
metallicity objects. Extremely metal poor galaxies (XMPG,
$Z<0.05~{}Z_{\odot}$) show more ULXs than higher-metallicity galaxies
(Prestwich et al., 2013). Moreover, Basu-Zych et al. (2013) find that a sample
of $z<0.1$ Lyman break analogs (LBAs), including Haro 11, likewise exhibit a
higher $L_{\rm X}/\rm SFR$ than solar-metallicity galaxies. Their
interpretation is that lower metallicity results in more luminous HMXBs, as
weaker stellar winds lead to more massive compact objects. Brorby et al.
(2016) quantify this $L_{\rm X}-\rm SFR-Z$ relation, where $L_{\rm X}/\rm SFR$
increases with lower $\log(\rm O/H)$. Finally, Dittenber et al. (2020) find
that the majority of local Ly$\alpha$ emitters, and thus candidates for LyC
escape, may be driven by ULXs, finding a connection between Ly$\alpha$ escape
and HMXBs and/or LLAGN. Thus, metal-poor starbursts in the early universe
likely formed an overabundance of X-ray binaries, which may have contributed
to the process of reionization.

From our LyC study of Haro 11, we see that LyC-emitting regions may coincide
with ULX sites, but the role, if any, of accretors in ionizing radiation
production and mechanical feedback remains to be clarified.

##  5 Summary and Conclusions

Haro 11 is a key object to understanding cosmic reionization as it is the
closest, and first, confirmed local LyC emitter. It is dominated by three
star-forming Knots A, B, and C of widely varying properties, and since the
original, spatially unresolved detection, it has remained unclear which of
these three regions is responsible for the observed LyC emission. We therefore
obtained new HST/COS G130M/1055 observations of each of the knots in the range
$900-1200$ Å, which reveal that Knots B and C are the LyC emitters toward our
line of sight. Their respective $903-912$ Å luminosities are $1.9\pm 1.5\times
10^{40}~{}\rm erg~{}s^{-1}$ and $0.9\pm 0.7\times 10^{40}~{}\rm erg~{}s^{-1}$.
So, Knot B seems to dominate the leaking LyC luminosity of Haro 11, accounting
for 66 $\pm$ 50$\%$ of the detected flux in $903-912$ Å, and Knot C accounts
for 35 $\pm$ 25$\%$. The total Haro 11 LyC luminosity is $2.9\pm 1.0\times
10^{40}~{}\rm erg~{}s^{-1}$.

We perform stellar population synthesis to constrain the stellar parameters
and thus local LyC escape fractions of each knot. For this, we combine our new
COS G130M/1055 with Sirressi et al. (2022) COS G130M/1300 + G160M/1600, as
well as Menacho et al. (2019) MUSE spectra. Fitting Starburst99 models to each
knot’s spectrum, we find that Knot B, the brightest LyC emitter, is dominated
by a $\sim 1$ Myr, $2\times 10^{7}$ $\rm M_{\odot}$ heavily obscured
component. Knot C is best fit with a continuous star formation history for
$\sim 15$ Myr with stellar mass of $3\times 10^{7}$ $\rm M_{\odot}$ and low
obscuration. Knot A is dominated by a $\sim 4$ Myr, $1\times 10^{7}$ $\rm
M_{\odot}$ highly obscured population. Our modeling thus uncovers massive,
young, obscured stellar populations in Knot B and Knot A, which dominate the
ionizing photon production in their respective regions but have minimal FUV
imprint and contribute $<40\%$ in the optical. We do not see conclusive
evidence of Very Massive Stars in any of the knots. Instead, we find clear
signatures of classical WR stars in Knots A and B.

The primary tracer of the young populations is the large H$\alpha$ luminosity,
which we reproduce with our models within a factor of two or better. Moreover,
Haro 11 shows a large FIR luminosity, $L_{\rm FIR}=2.7\times 10^{44}~{}\rm
erg~{}s^{-1}$, that qualifies it as a LIRG. Our dusty stellar population
models for the three knots combined account for $1.7\times 10^{44}~{}\rm
erg~{}s^{-1}$, which is within 60% of the observed integrated FIR radiation.
Thus the model estimates and observed values are in reasonable agreement,
especially considering that the IRAS aperture includes the entire galaxy.
Sirressi et al. (2022) photometrically detect $7-11$ clusters in each knot,
with ages $1-15$ Myr and masses up to $5\times 10^{7}$ $\rm M_{\odot}$ (Table
3). Our 3-component spectroscopic models capture the aggregate young, UV-
dominant populations, as well as the $\sim 15$ Myr generation that does not
contribute significant LyC. The spectrum of Knot C additionally shows an old
background population with an age up to $\sim 2$ Gyr.

The corresponding LyC escape fractions in the range $903-912$ Å are
$f_{\rm{esc,912}}$ $=3.4\pm 2.9\%$ for Knot B and $5.1\pm 4.3$% for Knot C. In
the case of Knot A, we place a $2\sigma$ upper limit to the escape fraction of
$f_{\rm{esc,912}}$$\leq 10$%. The luminosity-weighted escape fraction for the
entirety of Haro 11 is $f_{\rm esc,912}=3.9\pm 3.4$%, consistent with Leitet
et al. (2011).

Our results underscore that the LyC escape fraction ($f_{\rm{esc,912}}$) and
escaping LyC luminosity ($L_{\rm 912,obs}$) are distinct fundamental
parameters for characterizing LyC escape. Although we find that the majority
of Haro 11’s LyC flux likely originates from Knot B, the values above
demonstrate that its local escape fraction appears to be lower than that of
the LyC-fainter Knot C. The reason is that Knot B has by far the highest
ionizing photon production; but it exhibits the largest amount of neutral gas
among the knots, and in particular, it has a covering fraction close to unity
(Gao et al., 2022; Östlin et al., 2021). This results in a potentially lower
escape fraction. On the other hand, Knot C is intrinsically fainter in LyC due
to an older mass-weighted age, but it is significantly less obscured. Some of
its gas has likely been evacuated by feedback, as its HI covering fraction is
$f_{\rm c,HI}\sim 0.8$ and its neutral column density is the lowest among the
knots (Table 2; Östlin et al., 2021). Characterizing LyC emission by only
$f_{\rm{esc,912}}$ would prioritize Knot C over B, whereas Knot B is in fact
the more luminous LyC emitter. Thus, although the escape fraction is often
emphasized in the literature, the relevant parameter is the convolution of the
escape fraction and the intrinsic LyC luminosity, since both drive the
observable LyC luminosity and the number of ionizing photons leaked into the
IGM.

We also use our new observations to test Ly$\alpha$ as a LyC escape tracer by
comparing our results to Ly$\alpha$-based predictions. We see a correlation of
$f_{\rm{esc,912}}$ with the observed Ly$\alpha$ luminosity, and an inverse
relation with Ly$\alpha$ peak velocity separation $v_{\rm sep,Ly\alpha}$, as
expected. But the Ly$\alpha$ escape fraction $f_{\rm esc,Ly\alpha}$ does not
consistently trace LyC escape. While Knot C has both the highest Ly$\alpha$
and LyC escape fractions, the lowest $f_{\rm esc,Ly\alpha}$ in Knot B
incorrectly predicts it to be the weakest LyC leaker. The Ly$\alpha$ red peak
width similarly points to Knot C as the strongest leaker and Knot B as the
weakest, with the latter not consistent with our observations. So, the
observed Ly$\alpha$ luminosity and peak velocity separation appear to more
consistently correlate with LyC escape, while $f_{\rm esc,Ly\alpha}$ and
Ly$\alpha$ red peak width are less consistent tracers. However, we stress that
these results are based on only three star-forming knots in this galaxy.

There are important implications from the fact that Knot A is not detected in
the LyC, despite it driving the Green Pea properties of Haro 11. Green Peas
are the largest class of local LyC emitters, and their radiation properties
strongly correlate with $f_{\rm esc,LyC}$ (Flury et al., 2022a, b). Knot A has
therefore been predicted to be a leaking knot based on its ionization
parameter, which is the highest of the three knots (Keenan et al., 2017).
First, its non-detection highlights the potential importance of multiple star-
forming regions driving the LyC escape in Green Peas. In Haro 11, we see that
while the knot dominating the GP properties is not directly detected, other
knots are LyC emitters. Many GPs, like Haro 11, are mergers hosting multiple
star-forming knots. So, while GP properties and signatures of radiation-
dominated feedback are alone insufficient to predict LyC emission , the
multiplicity of starbursts may be important. Micheva et al. (2018) noted the
potential role of two-stage starbursts in LyC emitters. On the other hand,
Knot A shows clear signatures of density-bounded conditions in the plane of
the sky (Keenan et al., 2017). It is therefore likely that Knot A is leaking
LyC transverse to our line of sight, as seen in, e.g., NGC 5253 (Zastrow et
al., 2011). This implies that LyC escape must be highly anisotropic, as
predicted by simulations (Cen & Kimm, 2015) and inferred from observations
(Nestor et al., 2011; Flury et al., 2022b).

Lastly, we note the intriguing coincidence that the two LyC-leaking knots are
the hosts of the only two ULXs in Haro 11 (Prestwich et al., 2015; Gross et
al., 2021). Neither of the ULXs appear to contribute significantly to the LyC
emission, especially considering that the stellar populations dominate the UV
light by $1-2$ orders of magnitude. Nevertheless, the X-ray sources may
promote LyC escape through accretion-dominated mechanical feedback, where
powerful disk winds and jets may clear optically thin channels. A multitude of
studies show that reionization-era analogs, such as LBAs and BCDs, have an
overabundance of ULXs (e.g., Dittenber et al., 2020; Brorby et al., 2014;
Basu-Zych et al., 2013; Kaaret et al., 2011). Further investigation is thus
required into both the ultraviolet emission and mechanical feedback of ULXs,
in order to determine their role in LyC escape.

We thank Paul Draghis and Daniel Schaerer for useful discussions. This work
was supported by NASA grant to M.S.O HST-GO-16260. A. Adamo acknowledges the
support of the Swedish Research Council, Vetenskapsrådet grant 2021-05559.
JMMH is funded by Spanish MCIN/AEI/10.13039/501100011033 grant
PID2019-107061GB-C61. D.K. is supported by the Centre National d’Études
Spatiales (CNES)/Centre National de la Recherche Scientifique (CNRS);
convention no 230400. The Cosmic Dawn Center (DAWN) is funded by the Danish
National Research Foundation under grant DNRF140. Lastly, we thank the
anonymous referee for helpful suggestions that improved the manuscript. This
research is based on observations made with the NASA/ESA Hubble Space
Telescope obtained from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5–26555. These observations are associated with
programs 16260, 15352, and 13017. This research is also based on observations
collected at the European Southern Observatory under ESO programme
096.B-0923(A). . The HST data used in this paper can be found in MAST:
http://dx.doi.org/10.17909/p40f-yy93 (catalog 10.17909/p40f-yy93)

## References

  * Adamo et al. (2010) Adamo, A., Östlin, G., Zackrisson, E., et al. 2010, MNRAS, 407, 870, doi: 10.1111/j.1365-2966.2010.16983.x
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Atek et al. (2024) Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Nature, 626, 975, doi: 10.1038/s41586-024-07043-6
  * Basu-Zych et al. (2013) Basu-Zych, A. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2013, ApJ, 774, 152, doi: 10.1088/0004-637X/774/2/152
  * Bergvall et al. (2006) Bergvall, N., Zackrisson, E., Andersson, B. G., et al. 2006, A&A, 448, 513, doi: 10.1051/0004-6361:20053788
  * Bouwens et al. (2012) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 752, L5, doi: 10.1088/2041-8205/752/1/L5
  * Brorby et al. (2014) Brorby, M., Kaaret, P., & Prestwich, A. 2014, MNRAS, 441, 2346, doi: 10.1093/mnras/stu736
  * Brorby et al. (2016) Brorby, M., Kaaret, P., Prestwich, A., & Mirabel, I. F. 2016, MNRAS, 457, 4081, doi: 10.1093/mnras/stw284
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Cen & Kimm (2015) Cen, R., & Kimm, T. 2015, ApJ, 801, L25, doi: 10.1088/2041-8205/801/2/L25
  * Chisholm et al. (2018) Chisholm, J., Gazagnes, S., Schaerer, D., et al. 2018, A&A, 616, A30, doi: 10.1051/0004-6361/201832758
  * Clarke & Oey (2002) Clarke, C., & Oey, M. S. 2002, MNRAS, 337, 1299, doi: 10.1046/j.1365-8711.2002.05976.x
  * Cullen et al. (2023a) Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023a, MNRAS, 520, 14, doi: 10.1093/mnras/stad073
  * Cullen et al. (2023b) Cullen, F., McLeod, D. J., McLure, R. J., et al. 2023b, arXiv e-prints, arXiv:2311.06209, doi: 10.48550/arXiv.2311.06209
  * Danforth et al. (2010) Danforth, C. W., Keeney, B. A., Stocke, J. T., Shull, J. M., & Yao, Y. 2010, ApJ, 720, 976, doi: 10.1088/0004-637X/720/1/976
  * Danforth et al. (2016) Danforth, C. W., Keeney, B. A., Tilton, E. M., et al. 2016, ApJ, 817, 111, doi: 10.3847/0004-637X/817/2/111
  * Dashtamirova et al. (2019) Dashtamirova, D., White, J., & Sahnow, D. 2019, Changes in the COS/FUV Dark Rate: Impact on the Monitoring Program and Background Extraction Regions, Instrument Science Report COS 2019-11, 17 pages
  * Dijkstra et al. (2016) Dijkstra, M., Gronke, M., & Venkatesan, A. 2016, ApJ, 828, 71, doi: 10.3847/0004-637X/828/2/71
  * Dittenber et al. (2020) Dittenber, B., Oey, M. S., Hodges-Kluck, E., et al. 2020, ApJ, 890, L12, doi: 10.3847/2041-8213/ab723e
  * Douna et al. (2015) Douna, V. M., Pellizza, L. J., Mirabel, I. F., & Pedrosa, S. E. 2015, A&A, 579, A44, doi: 10.1051/0004-6361/201525617
  * Ekström et al. (2012a) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012a, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * Ekström et al. (2012b) —. 2012b, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614, doi: 10.1093/mnras/staa302
  * Fernández-Ontiveros et al. (2012) Fernández-Ontiveros, J. A., Prieto, M. A., Acosta-Pulido, J. A., & Montes, M. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012006, doi: 10.1088/1742-6596/372/1/012006
  * Flury et al. (2022a) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022a, ApJS, 260, 1, doi: 10.3847/1538-4365/ac5331
  * Flury et al. (2022b) —. 2022b, ApJ, 930, 126, doi: 10.3847/1538-4357/ac61e4
  * Fontanot et al. (2014) Fontanot, F., Cristiani, S., Pfrommer, C., Cupani, G., & Vanzella, E. 2014, MNRAS, 438, 2097, doi: 10.1093/mnras/stt2332
  * Fontanot et al. (2012) Fontanot, F., Cristiani, S., & Vanzella, E. 2012, MNRAS, 425, 1413, doi: 10.1111/j.1365-2966.2012.21594.x
  * Fujita et al. (2003) Fujita, A., Martin, C. L., Mac Low, M.-M., & Abel, T. 2003, ApJ, 599, 50, doi: 10.1086/379276
  * Gao et al. (2022) Gao, Y., Gu, Q., Shi, Y., et al. 2022, A&A, 661, A136, doi: 10.1051/0004-6361/202142309
  * Georgy et al. (2013) Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103, doi: 10.1051/0004-6361/201322178
  * Giallongo et al. (2015) Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83, doi: 10.1051/0004-6361/201425334
  * Gierliński et al. (2009) Gierliński, M., Done, C., & Page, K. 2009, MNRAS, 392, 1106, doi: 10.1111/j.1365-2966.2008.14166.x
  * Grimes et al. (2007) Grimes, J. P., Heckman, T., Strickland, D., et al. 2007, ApJ, 668, 891, doi: 10.1086/521353
  * Gross et al. (2021) Gross, A. C., Prestwich, A., & Kaaret, P. 2021, MNRAS, 505, 610, doi: 10.1093/mnras/stab1331
  * Guseva et al. (2012) Guseva, N. G., Izotov, Y. I., Fricke, K. J., & Henkel, C. 2012, A&A, 541, A115, doi: 10.1051/0004-6361/201118742
  * Hassan et al. (2018) Hassan, S., Davé, R., Mitra, S., et al. 2018, MNRAS, 473, 227, doi: 10.1093/mnras/stx2194
  * Heckman et al. (2011) Heckman, T. M., Borthakur, S., Overzier, R., et al. 2011, ApJ, 730, 5, doi: 10.1088/0004-637X/730/1/5
  * Helou et al. (1988) Helou, G., Khan, I. R., Malek, L., & Boehmer, L. 1988, ApJS, 68, 151, doi: 10.1086/191285
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Izotov et al. (2016a) Izotov, Y. I., Orlitová, I., Schaerer, D., et al. 2016a, Nature, 529, 178, doi: 10.1038/nature16456
  * Izotov et al. (2016b) Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016b, MNRAS, 461, 3683, doi: 10.1093/mnras/stw1205
  * Izotov et al. (2018) Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018, MNRAS, 478, 4851, doi: 10.1093/mnras/sty1378
  * James et al. (2013) James, B. L., Tsamis, Y. G., Walsh, J. R., Barlow, M. J., & Westmoquette, M. S. 2013, MNRAS, 430, 2097, doi: 10.1093/mnras/stt034
  * James et al. (2022) James, B. L., Berg, D. A., King, T., et al. 2022, arXiv e-prints, arXiv:2206.01224.  https://arxiv.org/abs/2206.01224
  * Jaskot et al. (2017) Jaskot, A. E., Oey, M. S., Scarlata, C., & Dowd, T. 2017, ApJ, 851, L9, doi: 10.3847/2041-8213/aa9d83
  * Justham & Schawinski (2012) Justham, S., & Schawinski, K. 2012, MNRAS, 423, 1641, doi: 10.1111/j.1365-2966.2012.20985.x
  * Kaaret & Corbel (2009) Kaaret, P., & Corbel, S. 2009, ApJ, 697, 950, doi: 10.1088/0004-637X/697/1/950
  * Kaaret et al. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, Annual Review of Astronomy and Astrophysics, 55, 303, doi: 10.1146/annurev-astro-091916-055259
  * Kaaret et al. (2011) Kaaret, P., Schmitt, J., & Gorski, M. 2011, ApJ, 741, 10, doi: 10.1088/0004-637X/741/1/10
  * Kakiichi & Gronke (2021) Kakiichi, K., & Gronke, M. 2021, ApJ, 908, 30, doi: 10.3847/1538-4357/abc2d9
  * Keenan et al. (2017) Keenan, R. P., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 848, 12, doi: 10.3847/1538-4357/aa8b77
  * Komarova et al. (2021) Komarova, L., Oey, M. S., Krumholz, M. R., et al. 2021, ApJ, 920, L46, doi: 10.3847/2041-8213/ac2c09
  * Kunth et al. (2003) Kunth, D., Leitherer, C., Mas-Hesse, J. M., Östlin, G., & Petrosian, A. 2003, ApJ, 597, 263, doi: 10.1086/378396
  * Kunth & Sargent (1981) Kunth, D., & Sargent, W. L. W. 1981, A&A, 101, L5
  * Leitet et al. (2011) Leitet, E., Bergvall, N., Piskunov, N., & Andersson, B.-G. 2011, Astronomy & Astrophysics, 532, A107, doi: 10.1051/0004-6361/201015654
  * Leitherer et al. (2014) Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14, doi: 10.1088/0067-0049/212/1/14
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Kasen, D., et al. 2016, MNRAS, 459, 3614, doi: 10.1093/mnras/stw941
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Madau & Haardt (2015) Madau, P., & Haardt, F. 2015, ApJ, 813, L8, doi: 10.1088/2041-8205/813/1/L8
  * Martins et al. (2023) Martins, F., Schaerer, D., Marques-Chaves, R., & Upadhyaya, A. 2023, A&A, 678, A159, doi: 10.1051/0004-6361/202346732
  * Matthee et al. (2023) Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2023, ApJ, 950, 67, doi: 10.3847/1538-4357/acc846
  * May et al. (2018) May, D., Rodríguez-Ardila, A., Prieto, M. A., et al. 2018, Monthly Notices of the Royal Astronomical Society: Letters, 481, L105, doi: 10.1093/mnrasl/sly155
  * Menacho et al. (2019) Menacho, V., Östlin, G., Bik, A., et al. 2019, MNRAS, 487, 3183, doi: 10.1093/mnras/stz1414
  * Menacho et al. (2021) —. 2021, MNRAS, 506, 1777, doi: 10.1093/mnras/stab1491
  * Merloni & Heinz (2013) Merloni, A., & Heinz, S. 2013, in Planets, Stars and Stellar Systems. Volume 6: Extragalactic Astronomy and Cosmology, ed. T. D. Oswalt & W. C. Keel, Vol. 6 (Springer Science), 503, doi: 10.1007/978-94-007-5609-0_11
  * Meynet et al. (1994) Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103, 97
  * Micheva et al. (2017) Micheva, G., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 845, 165, doi: 10.3847/1538-4357/aa830b
  * Micheva et al. (2018) Micheva, G., Oey, M. S., Keenan, R. P., Jaskot, A. E., & James, B. L. 2018, ApJ, 867, 2, doi: 10.3847/1538-4357/aae372
  * Middleton et al. (2015) Middleton, M. J., Heil, L., Pintore, F., Walton, D. J., & Roberts, T. P. 2015, Monthly Notices of the Royal Astronomical Society, 447, 3243, doi: 10.1093/mnras/stu2644
  * Morales et al. (2023) Morales, A. M., Finkelstein, S. L., Leung, G. C. K., et al. 2023, arXiv e-prints, arXiv:2311.04294, doi: 10.48550/arXiv.2311.04294
  * Nestor et al. (2011) Nestor, D. B., Shapley, A. E., Steidel, C. C., & Siana, B. 2011, ApJ, 736, 18, doi: 10.1088/0004-637X/736/1/18
  * Newville et al. (2016) Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Astrophysics Source Code Library, record ascl:1606.014.  http://ascl.net/1606.014
  * Östlin et al. (2001) Östlin, G., Amram, P., Bergvall, N., et al. 2001, A&A, 374, 800, doi: 10.1051/0004-6361:20010832
  * Östlin et al. (2009) Östlin, G., Hayes, M., Kunth, D., et al. 2009, AJ, 138, 923, doi: 10.1088/0004-6256/138/3/923
  * Östlin et al. (2015) Östlin, G., Marquart, T., Cumming, R. J., et al. 2015, A&A, 583, A55, doi: 10.1051/0004-6361/201323233
  * Östlin et al. (2021) Östlin, G., Rivera-Thorsen, T. E., Menacho, V., et al. 2021, ApJ, 912, 155, doi: 10.3847/1538-4357/abf1e8
  * Pellegrini et al. (2012) Pellegrini, E. W., Oey, M. S., Winkler, P. F., et al. 2012, ApJ, 755, 40, doi: 10.1088/0004-637X/755/1/40
  * Prestwich et al. (2015) Prestwich, A. H., Jackson, F., Kaaret, P., et al. 2015, ApJ, 812, 166, doi: 10.1088/0004-637X/812/2/166
  * Prestwich et al. (2013) Prestwich, A. H., Tsantaki, M., Zezas, A., et al. 2013, ApJ, 769, 92, doi: 10.1088/0004-637X/769/2/92
  * Prevot et al. (1984) Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389
  * Rivera-Thorsen et al. (2017) Rivera-Thorsen, T. E., Östlin, G., Hayes, M., & Puschnig, J. 2017, ApJ, 837, 29, doi: 10.3847/1538-4357/aa5d0a
  * Robertson et al. (2015) Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19, doi: 10.1088/2041-8205/802/2/L19
  * Ross et al. (2017) Ross, H. E., Dixon, K. L., Iliev, I. T., & Mellema, G. 2017, MNRAS, 468, 3785, doi: 10.1093/mnras/stx649
  * Savage & Sembach (1991) Savage, B. D., & Sembach, K. R. 1991, ApJ, 379, 245, doi: 10.1086/170498
  * Saxena et al. (2023) Saxena, A., Bunker, A. J., Jones, G. C., et al. 2023, arXiv e-prints, arXiv:2306.04536, doi: 10.48550/arXiv.2306.04536
  * Sazonov & Khabibullin (2018) Sazonov, S., & Khabibullin, I. 2018, Monthly Notices of the Royal Astronomical Society, 476, 2530, doi: 10.1093/mnras/sty442
  * Schaerer et al. (1993) Schaerer, D., Charbonnel, C., Meynet, G., Maeder, A., & Schaller, G. 1993, A&AS, 102, 339
  * Schaerer et al. (2016) Schaerer, D., Izotov, Y. I., Verhamme, A., et al. 2016, A&A, 591, L8, doi: 10.1051/0004-6361/201628943
  * Shankar & Mathur (2007) Shankar, F., & Mathur, S. 2007, ApJ, 660, 1051, doi: 10.1086/512845
  * Sharma et al. (2017) Sharma, M., Theuns, T., Frenk, C., et al. 2017, MNRAS, 468, 2176, doi: 10.1093/mnras/stx578
  * Simmonds et al. (2023) Simmonds, C., Tacchella, S., Maseda, M., et al. 2023, MNRAS, 523, 5468, doi: 10.1093/mnras/stad1749
  * Sirressi et al. (2022) Sirressi, M., Adamo, A., Hayes, M., et al. 2022, MNRAS, 510, 4819, doi: 10.1093/mnras/stab3774
  * Smith et al. (2020) Smith, B. M., Windhorst, R. A., Cohen, S. H., et al. 2020, ApJ, 897, 41, doi: 10.3847/1538-4357/ab8811
  * Soderblom (2021) Soderblom, D. R. 2021, in COS Data Handbook v. 5.0, Vol. 5 (STScI), 5
  * Swartz et al. (2011) Swartz, D. A., Soria, R., Tennant, A. F., & Yukita, M. 2011, The Astrophysical Journal, 741, 49, doi: 10.1088/0004-637x/741/1/49
  * Topping et al. (2023) Topping, M. W., Stark, D. P., Endsley, R., et al. 2023, arXiv e-prints, arXiv:2307.08835, doi: 10.48550/arXiv.2307.08835
  * Verhamme et al. (2015) Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7, doi: 10.1051/0004-6361/201423978
  * Verhamme et al. (2017) Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13, doi: 10.1051/0004-6361/201629264
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Vinokurov et al. (2013) Vinokurov, A., Fabrika, S., & Atapin, K. 2013, Astrophysical Bulletin, 68, 139, doi: 10.1134/S1990341313020028
  * Volonteri & Gnedin (2009) Volonteri, M., & Gnedin, N. Y. 2009, ApJ, 703, 2113, doi: 10.1088/0004-637X/703/2/2113
  * Wang et al. (2021) Wang, B., Heckman, T. M., Amorín, R., et al. 2021, ApJ, 916, 3, doi: 10.3847/1538-4357/ac0434
  * Wofford et al. (2023) Wofford, A., Sixtos, A., Charlot, S., et al. 2023, MNRAS, 523, 3949, doi: 10.1093/mnras/stad1622
  * Yeh et al. (2023) Yeh, J. Y. C., Smith, A., Kannan, R., et al. 2023, MNRAS, 520, 2757, doi: 10.1093/mnras/stad210
  * Zackrisson et al. (2011) Zackrisson, E., Rydberg, C.-E., Schaerer, D., Östlin, G., & Tuli, M. 2011, ApJ, 740, 13, doi: 10.1088/0004-637X/740/1/13
  * Zastrow et al. (2013) Zastrow, J., Oey, M. S., Veilleux, S., & McDonald, M. 2013, ApJ, 779, 76, doi: 10.1088/0004-637X/779/1/76
  * Zastrow et al. (2011) Zastrow, J., Oey, M. S., Veilleux, S., McDonald, M., & Martin, C. L. 2011, ApJ, 741, L17, doi: 10.1088/2041-8205/741/1/L17