# Observations of Extremely Metal-Poor O Stars: Weak Winds and Constraints for
Evolution Models

O. Grace Telford Carnegie-Princeton Fellow  Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA;
<EMAIL_ADDRESS>The Observatories of the Carnegie Institution for
Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA  Rutgers
University, Department of Physics and Astronomy, 136 Frelinghuysen Road,
Piscataway, NJ 08854, USA  John Chisholm Department of Astronomy, The
University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX
78712-1205, USA  Andreas A. C. Sander Zentrum für Astronomie der Universität
Heidelberg, Astronomisches Rechen-Institut, Mönchhofstr. 12-14, 69120
Heidelberg, Germany  Varsha Ramachandran Zentrum für Astronomie der
Universität Heidelberg, Astronomisches Rechen-Institut, Mönchhofstr. 12-14,
69120 Heidelberg, Germany  Kristen B. W. McQuinn Rutgers University,
Department of Physics and Astronomy, 136 Frelinghuysen Road, Piscataway, NJ
08854, USA  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Danielle A. Berg Department of Astronomy, The
University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX
78712-1205, USA

###### Abstract

Metal-poor massive stars drive the evolution of low-mass galaxies, both
locally and at high redshift. However, quantifying the feedback they impart to
their local surroundings remains uncertain because models of stellar
evolution, mass loss, and ionizing spectra are unconstrained by observations
below 20% solar metallicity ($Z_{\odot}$). We present new Keck Cosmic Web
Imager optical spectroscopy of three O stars in the nearby dwarf galaxies Leo
P, Sextans A, and WLM, which have gas-phase oxygen abundances of 3–14%
$Z_{\odot}$. To characterize their fundamental stellar properties and
radiation-driven winds, we fit PoWR atmosphere models to the optical spectra
simultaneously with Hubble Space Telescope far-ultraviolet (FUV) spectra and
multi-wavelength photometry. We find that all three stars have effective
temperatures consistent with their spectral types and surface gravities
typical of main-sequence dwarf stars. Yet, the combination of those inferred
parameters and luminosity for the two lower-$Z$ stars is not reproduced by
stellar evolution models, even those that include rotation or binary
interactions. The scenario of multiple-star systems is difficult to reconcile
with all available data, suggesting that these observations pose a challenge
to current evolution models. We highlight the importance of validating the
relationship between stellar mass, temperature, and luminosity at very low $Z$
for accurate estimates of ionizing photon production and spectral hardness.
Finally, all three stars’ FUV wind profiles reveal low mass-loss rates and
terminal wind velocities in tension with expectations from widely adopted
radiation-driven wind models. These results provide empirical benchmarks for
future development of mass-loss and evolution models for metal-poor stellar
populations.

††facilities: Keck:II (KCWI), HST (COS, ACS, WFC3, WFPC2)††software: PoWR
(Gräfener et al., 2002; Hamann & Gräfener, 2003; Sander et al., 2015), Astropy
(Astropy Collaboration et al., 2013, 2018, 2022), CWITools (O’Sullivan & Chen,
2020), KCWI DRP (Morrissey et al., 2018), iPython (Perez & Granger, 2007),
Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011; Harris et al.,
2020), SAOImageDS9 (Joye & Mandel, 2003), SciPy (Virtanen et al., 2020)

##  1 Introduction

The recently launched James Webb Space Telescope (JWST) is now successfully
detecting large numbers of galaxies at very high redshift ($z$; e.g., Roberts-
Borsani et al. 2022; Arrabal Haro et al. 2023; Curtis-Lake et al. 2023), well
into the epoch when the intergalactic medium (IGM) was reionized at $z\gtrsim
6$ (e.g., Becker et al., 2001; Fan et al., 2006). At these early times (less
than 1 Gyr after the Big Bang), metal deposition by stars and supernovae had
not yet enriched galaxies’ star-forming gas to the metallicities ($Z$)
observed in galaxies of comparable mass in the present-day universe. Indeed,
JWST studies are now routinely requiring very low stellar $Z$ ($\lesssim$ 10%
of the solar metallicity, $Z_{\odot}$) to explain both the nebular emission
and spectral energy distributions (SEDs) of $z>6$ galaxies (e.g., Cameron et
al. 2023; Endsley et al. 2023; Topping et al. 2024). The hottest and most
massive O-type stars in these metal-poor galaxies dominate the rest-frame UV
and optical light that we observe; produce the ionizing radiation that powers
their nebular emission; and are the progenitors of supernova explosions that
regulate star formation and facilitate the escape of ionizing photons into the
IGM (e.g., Byler et al., 2017; Trebitsch et al., 2017; Dayal & Ferrara, 2018).

To understand the nature of these early galaxies, stellar population synthesis
(SPS) models are required to infer their physical properties from observed
integrated light (e.g., Tinsley, 1980; Conroy, 2013). Models of stellar
evolution and spectra, both key ingredients in SPS codes, are sensitive to the
treatment of various physical processes in massive-star atmospheres, including
mixing (convection, rotation), metal opacities, and mass-loss via radiation-
driven winds (e.g., Maeder & Meynet, 2000; Eldridge & Stanway, 2022; Vink,
2022). In the very low-$Z$ regime ($\lesssim$ 20% $Z_{\odot}$), SPS models
remain fundamentally limited due to a lack of empirical data that can assess
the validity of the theoretical approaches that are typically adopted.

Both optical and far-ultraviolet (FUV) spectra of individual massive stars are
required to measure key stellar and wind properties for comparison to stellar-
model predictions. Large datasets of suitable quality and wavelength coverage
have been assembled over decades for O and B stars in the Milky Way
($Z_{\odot}$; e.g., Prinja et al. 1990; Przybilla et al. 2010; Holgado et al.
2020), and in the lower-$Z$ environments of the Large and Small Magellanic
Clouds (LMC and SMC; 50% and 20% $Z_{\odot}$, respectively; Hunter et al.
2007; Trundle et al. 2007) with the latter recently boosted by the ULLYSES and
X-Shooting ULLYSES programs (Roman-Duval et al., 2020; Vink et al., 2023).
Below 20% $Z_{\odot}$, however, few spectra of massive stars exist (Bresolin
et al., 2006, 2007; Evans et al., 2007, 2019; Garcia & Herrero, 2013; Lorenzo
et al., 2022; Gull et al., 2022) because the necessary low-metallicity
environments are almost exclusively found in more distant dwarf galaxies at
the edge of or beyond the Local Group, which are faint and expensive to
observe. Additional observations in this low-$Z$ regime are required to
constrain how the properties of massive stars and their evolutionary pathways
change at the metallicities relevant in the early universe.

Massive stars lose mass via radiation-driven winds, which impacts the
evolutionary path and endpoint for a star of a given initial stellar mass
($M_{\star}$) and $Z$, as well as the ionizing photon production of that star
over its lifetime (Lamers & Cassinelli, 1999; Kudritzki & Puls, 2000; Martins
& Palacios, 2021). These winds, accelerated to high terminal velocities
($v_{\infty}$), also deposit substantial energetic feedback into their
surroundings and contribute to regulating the star-formation process
(Leitherer et al., 1992; Rosen, 2022). As hot-star winds are largely driven by
momentum transfer via absorption and scattering of photons in line transitions
of ionized metals (Castor et al., 1975; Pauldrach et al., 1986; Vink et al.,
1999), theory predicts that mass-loss rates ($\dot{M}_{\star}$) decrease
toward lower $Z$. However, different theoretical treatments make remarkably
different predictions of how $\dot{M}_{\star}$ and the $v_{\infty}$ scale with
$Z$ and other stellar properties (particularly the bolometric luminosity,
$L_{\star}$, and the effective temperature, $T_{\mathrm{eff}}$; e.g., Vink et
al. 2001; Krtička & Kubát 2018; Björklund et al. 2021).

Observational studies are therefore essential to clarify the $Z$ dependence of
mass loss from OB stars and avoid incorrect extrapolations from higher to
lower metallicity. Mokiem et al. (2007) constructed an empirical scaling
relation between $\dot{M}_{\star}$ and $L_{\star}$ at Milky Way, LMC, and SMC
metallicities, finding both a normalization and $Z$ dependence consistent with
the modeling of Vink et al. (2001). However, those $\dot{M}_{\star}$ drawn
from the literature were inhomogeneous in terms of data quality, wavelength
coverage, and analysis technique. The Vink et al. (2001) mass-loss
prescription has since been found to over-predict $\dot{M}_{\star}$ for O
stars at 20% $Z_{\odot}$ (e.g., Bouret et al., 2013), and multiple
observational analyses have found evidence for a steeper scaling between
$\dot{M}_{\star}$ and $Z$ (Ramachandran et al., 2019; Marcolino et al., 2022;
Rickard et al., 2022).

Some efforts have been made to observationally constrain the stellar and wind
properties of OB stars in more distant dwarf galaxies ($\sim$0.7–1.5 Mpc).
Tramper et al. (2011, 2014) analyzed optical spectra of several stars in the
galaxies IC 1613, NGC 3109, and WLM, finding higher $\dot{M}_{\star}$ than
theoretical predictions based on metallicity estimates from the gas-phase
oxygen abundances of those galaxies (9–12% $Z_{\odot}$; Marble et al. 2010).
Bouret et al. (2015) re-analyzed three of those stars by fitting atmosphere
models to FUV spectra, and concluded that the $\dot{M}_{\star}$ measured from
H$\alpha$ emission alone were overestimated and that the stars’ iron
abundances may be SMC-like. Though the sample size was small, the inferred
$\dot{M}_{\star}$ and stellar parameters did not follow a well-defined power-
law scaling. For eight OB stars in IC 1613, Garcia et al. (2014) studied the
observed $v_{\infty}$ via UV spectroscopy, finding no clear differences
compared to similar SMC and LMC stars. Contrary to expectations from measured
gas-phase and B supergiant oxygen abundances (Marble et al., 2010; Bresolin et
al., 2007), the iron abundance of the young stars in IC 1613 seems to be
similar to the SMC value. In summary, while these previously studied metal-
poor galaxies are interesting sub-solar $\alpha$/Fe environments, few massive
stars with sub-SMC iron abundances have been observed.

Table 1: Basic properties of the three target stars and their host galaxies.  Star | R.A. | Decl. | $m_{V}$ | $A_{V}$ | Spectral | Host | 12+$\log$(O/H) | $Z_{\mathrm{gas}}/Z_{\odot}$ | Distance
---|---|---|---|---|---|---|---|---|---
Name | (J2000) | (J2000) | (mag) | (mag) | Type | Galaxy | (nebular) | (nebular) | (Mpc)
LP26 | 10:21:45.1217 | +18:05:16.93 | 21.51 | 0.04 | O7–8 V | Leo P | $7.17\pm 0.04$ | 3% | $1.62\pm 0.15$
S3 | 10:10:58.1866 |  $-$04:43:18.45 | 20.80 | 0.09 | O9 V | Sextans A | $7.49\pm 0.04$ | 6% | $1.38\pm 0.05$
A15 | 00:02:00.5333 |  $-$15:29:52.41 | 20.25 | 0.04 | O7 V((f)) | WLM | $7.83\pm 0.06$ | 14% | $0.98\pm 0.04$

Telford et al. (2021) increased the small dataset that can constrain wind
properties below 20% $Z_{\odot}$ by presenting FUV spectra of three O-dwarf
(luminosity class V) stars in nearby, metal-poor dwarf galaxies (with gas-
phase oxygen abundances of 3–14% $Z_{\odot}$; see Table 1) observed with the
Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST). Those
authors characterized the observed photospheric and wind features and
estimated the fundamental stellar parameters from the observed stellar SEDs,
but lacked the optical spectroscopy necessary for quantitative analysis and
$\dot{M}_{\star}$ inference. Here, we build upon that work by presenting new
optical spectra of the three stars observed with the Keck Cosmic Web Imager
(KCWI), as well as multi-wavelength measurements of the stellar and wind
properties obtained from detailed atmosphere modeling of the optical and FUV
spectra.

The remainder of the paper is organized as follows. Section 2 presents the
observations used in our analysis: FUV and optical spectroscopy, as well as
stellar photometry measured from archival HST imaging. Section 3 describes our
stellar atmosphere modeling technique, and the best-fit models and inferred
stellar and wind parameters are presented in Section 4. In Section 5, we
compare our results to the predictions of widely used stellar evolution
models, mass-loss prescriptions, and grids of ionizing spectra predicted by
atmosphere models. Finally, we summarize our conclusions in Section 6.

Throughout, we adopt the Asplund et al. (2009) solar abundances as our
reference, with bulk metal mass fraction $Z_{\odot}=0.0142$. Mass fractions of
individual elements are denoted as, e.g., $X_{\mathrm{C}}$ for the carbon mass
fraction. In general, when discussing the metallicities of galaxies, we mean
the gas-phase oxygen abundance (12+$\log$(O/H)), while stellar metallicities
refer to the bulk metal mass fraction $Z$. As mentioned above, evidence exists
for non-solar O/Fe abundance ratios in metal-poor dwarf galaxies, with
important implications for the opacities and $\dot{M}_{\star}$ of massive
stars that are primarily set by the Fe abundance. Thus, we use the
$X_{\mathrm{Fe}}/X_{\mathrm{{Fe,}\odot}}$ of observed stars to determine the
stellar model $Z$ against which to compare. Wherever model predictions are
discussed, we rescale the bulk metal mass fraction $Z$ from the solar
reference adopted by each set of models to express all metallicities relative
to the Asplund et al. (2009) solar abundance scale. Finally, we use the terms
extremely or very metal-poor to mean $\leq 10$% $Z_{\odot}$, following Kunth &
Östlin (2000).

##  2 Observations and Data Reduction

###  2.1 Optical Spectroscopy

Table 2: Details of the Keck/KCWI observations.  Star | R.A. | Decl. | Position Angle | Exposure Time | Typical Seeing | Continuum SNR | Flux
---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (∘ E of N) | (s) | (′′) | (4500–4700 Å) | Standard
LP26 | 10:21:44.91 | +18:05:20.6 | 150 | 12000 [$10\times 1200$] | 1.9 | 32 | Feige34
S3 | 10:10:58.42 |  $-$04:43:25.8 | 110 | 9600 [$8\times 1200$] | 1.3 | 50 | Feige34, GD108
A15 | 00:02:00.05 |  $-$15:29:52.1 | 90 | 7200 [$6\times 1200$] | 1.0 | 46 | LTT949
Figure 1: KCWI fields observed in the three dwarf galaxies. From left to
right, white-light images of Leo P, Sextans A, and WLM are shown. Each image
is constructed by integrating the flux in the co-added KCWI observations of
each galaxy over the full wavelength range. The three stars analyzed here are
labeled and encircled in blue. The 1D spectra used in this analysis are the
summed spectra in the spaxels falling within those blue contours, which
enclose 50% of the stellar flux. White bars in the bottom-left corners show 10
pc at the distance of each galaxy for scale.

Spectroscopy at visible wavelengths is required to determine the fundamental
properties of O-type stars. Specifically, the relative strengths of He i and
He ii photospheric absorption lines is the primary $T_{\mathrm{eff}}$
diagnostic for mid-late O stars, while the Balmer line profiles encode the
surface gravity, ($\log(g)$; e.g., Simón-Díaz 2020). We therefore observed
optical spectra of the three metal-poor O stars in our sample using the Keck
Cosmic Web Imager (KCWI) on the 10-m Keck II Telescope. This integral field
unit (IFU) spectrograph is highly sensitive at blue wavelengths, ideal for
measuring the key diagnostic absorption features between 4000–5000 Å.

All three dwarf galaxy targets were observed by NASA Keck program N194111The
KCWI data presented in this paper can be downloaded from the Keck Observatory
Archive: https://koa.ipac.caltech.edu/. Leo P was observed on 2021 January 13
(2020B semester, PI: J. Chisholm), followed by WLM on 2021 August 10 and
Sextans A on 2022 January 5 (2021B semester, PI: O. G. Telford). We used the
same instrument configuration for all observations, selecting the medium
slicer, BM grating, and a central wavelength (cenwave) of 4500 Å. This setup
yields a spectral resolution of $R$$\sim$4000, a $16\farcs 5\times 20\farcs 0$
field of view (FoV), spatial pixels (or spaxels) of angular size $0\farcs
29\times 0\farcs 68$, and covers $\sim$3970–4930 Å. Table 2 summarizes the
exposure times and typical conditions during the observations of the three
galaxies. We interleaved our science exposures with observations of off-galaxy
fields for both Sextans A and WLM to measure the sky background. In the case
of Leo P, its bright stellar and nebular emission covers a small enough area
that we were able to obtain a background observation by nodding the FoV up to
$5\farcs 4$ to the northeast (with the position angle fixed), keeping the
target star LP26 in the FoV throughout.

Leo P was also observed by the same program for 30 minutes on 2022 January 5
at a different cenwave of 4700 Å to cover the [O iii]$\lambda\lambda$ 4959,
5007 nebular lines, and the emission-line spectra from both Leo P observations
were presented in Telford et al. (2023). However, we do not include the 2022
data in this analysis because the continuum signal-to-noise ratio (SNR) of
$\sim$20 is lower than achieved in the earlier observations, and we verified
that combining both datasets does not improve the final data quality. We
prefer to use only the deeper 2021 data at the same cenwave observed for all
three stars in our sample to maintain consistent methodology.

We reduce the KCWI data using the IDL-based KCWI data reduction
pipeline222https://github.com/Keck-DataReductionPipelines/KcwiDRP (DRP;
version 1.2.1; Morrissey et al. 2018). We follow the same process described in
Section 2.1 of Telford et al. (2023) to calibrate the images using bias, dark,
and internal/dome flat frames; reject cosmic rays; obtain a wavelength
calibration using arc-lamp frames; and transform the data from 2D images to 3D
data cubes. The offset field observations (or off-galaxy regions of the data
cubes in the case of Leo P) were used to subtract the sky background. We use
the interactive tool in the KCWI DRP to derive inverse sensitivity curves from
the flux standard spectra obtained during all observing runs (Table 2) and
transform units of counts to flux. Residuals in the fits to the standard-star
continua were $\lesssim$ 5%.

We combine the reduced and flux-calibrated individual exposures using CWITools
(O’Sullivan & Chen, 2020). A WCS solution is obtained for each exposure, given
the known coordinates of the target star, from a Gaussian fit to the star’s
light profile. The exposures are then aligned and interpolated onto a common
spatial and wavelength grid with square pixels 0$\farcs$29 on a side and 0.5-Å
wavelength sampling (or $\sim$37.5 km s-1). Finally, the data cubes
constructed from individual exposures are averaged, weighted by exposure time.
Figure 1 shows white-light images constructed by integrating each spaxel in
the coadded data cubes over all wavelengths. The target star in each image is
encircled in blue and labeled, where the circle size corresponds to the
inferred seeing (i.e., the full width at half maximum, or FWHM, of a 2D
Gaussian fit to the star’s light profile; see next paragraph). Stars observed
under better conditions appear elongated along the horizontal direction in the
coadded images because the KCWI detector pixels are rectangular, and have been
resampled to a uniform grid of square pixels. Each target star is nicely
separated from any comparably bright OB stars in high-resolution HST imaging
(0.05′′ pixels) of the same fields.

Figure 2: KCWI optical spectra of the three O stars. Continuum-normalized KCWI
spectra of LP26 (top, purple), S3 (middle, blue), and A15 (bottom, red)
plotted as a function of rest-frame wavelength. An integer offset is added for
clarity, and the continuum level for each star is indicated by a horizontal
black line with gray shading to show $\pm$ 5%. The key diagnostics of O-star
$T_{\mathrm{eff}}$ (He i and He ii lines) and $\log(g)$ (Balmer lines) are
labeled at the top of the plot. All three spectra show narrow nebular emission
filling in the centers of the Balmer lines, and the He i transitions are also
affected for S3 (primarily the 4471 Å line) and LP26 (all lines). These
normalized spectra are used in our PoWR modeling analysis in Sections 3.3 and
4 below.

From these coadded data cubes, we extract the 1D spectrum for each star by
fitting a 2D Gaussian model to its spatial profile in the white-light image
and summing the spectra in all spaxels within the FWHM, which captures 50% of
the total flux from each star. Adopting a larger radius (out to two times the
FWHM, which encloses 94% of the flux) does not increase the final SNR
substantially, but does increase the contribution of nebular emission to the
final spectra. We allow for different FWHM in the two spatial dimensions to
account for the rectangular spaxels in the raw KCWI images. In the case of
LP26, the 2D Gaussian was fit to the profile of the H ii region in which the
star is embedded, which is larger than the star itself. We tested instead
using the FWHM of the flux-standard star to define the aperture within which
the 1D spectrum of LP26 is constructed, but found that the FWHM of the
standard was similar to that of the H ii region.

Because observing conditions inevitably vary throughout the night, we
considered extracting a 1D spectrum within the FWHM of the center of each
target star’s light profile in the individual exposures, then combining those
1D spectra. However, we found that this does not improve the SNR of the final
1D stellar spectrum, and more importantly, does not change the equivalent
width or shape of any of the stellar lines used in our atmosphere modeling.
Thus, we prefer to use the CWITools coadds for consistency with our previous
analysis of KCWI observations of Leo P (Telford et al., 2023). We also tried
excluding the exposures taken under the poorest seeing during each observing
run, but found that including all available data maximizes the continuum SNR,
and this choice does not change the stellar absorption line shapes.

For comparison to the optical spectra predicted by stellar atmosphere models
(see Sections 3.3 and 4 below), observed optical spectra are typically
normalized to a continuum level of 1. For each star, we model the continuum as
a fifth-degree polynomial, masking all absorption and emission lines, and
divide the observed flux by that model continuum to obtain the normalized
spectrum used in this analysis. All three spectra contained a residual
sinusoidal feature over $\sim$4170–4300 Å with amplitude $\lesssim$ 5%, which
was also visible in the standard-star spectra and is thus an instrumental
artifact. To remove this unphysical feature, we divide out a spline fit to the
continuum in that region only, masking the he ii 4200 Å absorption line, but
this choice does not affect the conclusions drawn from our modeling below. The
SNR reported in Table 2 is calculated from the normalized continuum rms
between 4500–4700 Å, with the absorption lines in that range masked. Finally,
we correct the observed wavelengths to the heliocentric reference frame, with
small velocity shifts of $\sim$20 km s-1 (i.e., less than one pixel at this
spectral resolution).

The final, normalized spectra used in our analysis below are shown in Figure
2, vertically offset for clarity, in order of increasing host-galaxy
metallicity from top to bottom: LP26 (purple), S3 (blue), and A15 (red). Black
lines indicate the continuum level for each star, with gray shading to
demonstrate that the normalized data fall within $\pm$ 5% of the flat
continuum level. All three spectra are affected by narrow nebular emission,
but this is cleanly separable from the Balmer-line wings that are critical to
constrain $\log(g)$. The He i photospheric absorption lines are completely
filled in by nebular emission for LP26, which is surrounded by a bright H ii
region, but are clearly detected in the spectra of S3 and A15. All three
spectra contain clear He ii absorption, confirming their classification as
O-type stars. Based on the ratio of the equivalent widths of the observed he
ii 4540 Å and he i 4471 Å lines (Martins, 2018), we confirm the spectral types
of S3 and A15 as O9 and O7, respectively. We cannot quantitatively assign a
spectral type to LP26 due to the infilling of all optical he i lines by
nebular emission. All three stars have broad Balmer wings, consistent with
luminosity class V.

###  2.2 Far-Ultraviolet Spectroscopy

Most important diagnostics of stellar mass loss and metal abundances for hot
stars are highly ionized metal transitions found in the FUV part of the
spectrum. Thus, we use the HST/COS FUV spectra of the three target stars
presented in Telford et al. (2021) in this analysis. A comprehensive
discussion of the data and reduction can be found in Section 2.3 of that
paper, but we briefly summarize key aspects here.

The COS spectra were obtained using the moderate-resolution G130M and G160M
gratings during April–June of 2020 (GO-15967; PI: J. Chisholm). For each star,
the CalCOS pipeline-reduced spectra333The HST/COS data analyzed in this paper
can be downloaded from MAST: http://dx.doi.org/10.17909/hs6f-a370 (catalog
10.17909/hs6f-a370) from the two gratings were downloaded from the Mikulski
Archive for Space Telescopes (MAST), resampled onto a common wavelength grid
at the native sampling of the G160M grating, and coadded. Based on the
velocity FWHM of Milky Way interstellar medium (ISM) absorption features in
the spectra, the spectral resolution is better than 50 km s-1. We bin the
spectra by 12 pixels to maximize SNR while ensuring 1.5–2 samples per
resolution element. Given our goal of modeling the weak photospheric
absorption features in these metal-poor stars, achieving the highest SNR
possible is essential.

The FUV data will be compared to model spectra in Section 4 below, providing
both line diagnostics and a constraint on the overall SED shape ($\pm 5$%
flux-calibration uncertainties444https://hst-docs.stsci.edu/cosihb). The model
spectra account for dust extinction, both in Milky Way and internal to the
host dwarf galaxies, so we do not apply any extinction correction to the
observed spectra. For the star LP26 in Leo P only, we correct the FUV
observations for the contribution of nebular continuum emission from its
surrounding H ii region, which is not included in the model spectra. Telford
et al. (2021) found that a substantial contribution to the FUV flux from
nebular continuum (up to a maximum of 19% at the reddest wavelengths in the
G160M data) was required to explain the shape of LP26’s SED. Thus, we divide
the observed spectrum by the best-fit nebular continuum in that SED modeling
to ensure a fair comparison between the observed COS spectrum of LP26 and the
PoWR model spectra. No significant nebular continuum contribution was detected
for S3 or A15, consistent with their environments outside of bright H ii
regions.

###  2.3 Near-Ultraviolet, Optical, and Near-Infrared Photometry from Archival
Imaging

Finally, the stellar SED constructed from photometry across a wide wavelength
baseline is required to constrain $L_{\star}$, $M_{\star}$, and reddening
$E(B-V)$ due to dust along the line of sight to the star. We employ the high-
quality stellar photometry measured from archival near-ultraviolet (NUV),
optical, and near-infrared (NIR) HST imaging presented in Table 3 of Telford
et al. (2021). Details on the process of measuring point-spread function (PSF)
photometry with dolphot (Dolphin, 2000, 2016), including parameter choices and
quality metrics, are given in Section 2.2 of that paper. These measurements
typically have quite high SNR of $>$ 100 in the NUV and optical, and
$\sim$10–30 in the NIR. We use the observed magnitudes and fluxes, uncorrected
for Milky Way foreground dust along the line of sight to the targets, and
include dust extinction in the atmosphere models to which we compare below.

Accurate distances to the host galaxies are essential to translate observed
into absolute magnitudes. Thus, we use distances determined from the tip of
the red-giant branch and, for Leo P, the horizontal branch in optical color-
magnitude diagrams (Table 1). These are among the most precise distance-
measurement techniques for nearby galaxies (e.g., Beaton et al., 2018), and we
account for distance uncertainties in the $L_{\star}$ measurements reported
below.

##  3 PoWR Stellar Atmosphere Models

###  3.1 The PoWR Code

To analyze the observations, we employ the Potsdam Wolf-Rayet PoWR stellar-
atmosphere code (Gräfener et al., 2002; Hamann & Gräfener, 2003; Sander et
al., 2015). PoWR is a 1D state-of-the-art code solving the statistical
equations together with the radiative transfer in an expanding, stationary
outflow outside of local thermodynamic equilibrium (i.e., in a non-LTE
environment). The radiative transfer is performed in the co-moving frame to
properly account for the stellar wind and the quasi-hydrostatic photosphere.
After a model atmosphere has converged, the emergent spectrum is obtained by
computing the formal integral along rays in the observer’s frame. Fundamental
input parameters are $T_{\star}$, $L_{\star}$, $M_{\star}$ (or alternatively,
$\log(g)$), and the chemical composition. Formally, $T_{\star}$ denotes the
effective temperature corresponding to the inner boundary of the model, at a
Rosseland optical depth ($\tau_{\mathrm{Ross}}$) of 20. For a better
comparison, we will list the more common $T_{\mathrm{eff}}$ defined at
$\tau_{\mathrm{Ross}}=2/3$ throughout this work, which is an output parameter
of PoWR. All listed radii and surface gravity values are also referring to the
same optical depth.

The velocity field in the supersonic wind is prescribed by a $\beta$-type
velocity law:

$v(r)=v_{\infty}\left(1-\frac{R_{\star}}{r}\right)^{\beta},$ (1)

where we adopt $\beta=0.8$ and adjust the free parameter $v_{\infty}$ as
necessary to reproduce the wind features in the UV spectra. The velocity and
density stratification in the subsonic part of the atmosphere is obtained from
a consistent solution of the hydrostatic equation, accounting for the full
radiative force from continua and lines (Sander et al., 2015), and is smoothly
connected to the $\beta$-law supersonic wind regime.

The $\dot{M}_{\star}$ of the wind is a free parameter that is adjusted to
match to the spectral appearance of wind-sensitive features, particularly
metal resonance lines in the FUV spectra. Wind inhomogenieties are treated in
the microclumping approach: the wind is assumed to be composed of small clumps
with density $D\rho(r)$, where $D$ is the clumping factor and $\rho(r)$ is the
density profile of a homogenous wind with the same $\dot{M}_{\star}$ (Hamann &
Koesterke, 1998). The interclump space is void, so the wind volume-filling
factor is $D^{-1}$. We adopt a depth-dependent $D$, increasing from unity to
$D_{\infty}=10$ between the sonic point and $10\,R_{\star}$. This is a
commonly adopted value (e.g., Hainich et al. 2019) and is supported by a
recent analysis of similar O V stars in the 20% $Z_{\odot}$ SMC by Rickard et
al. (2022), who found $5\leq D_{\infty}\leq 20$. Because
$\dot{M}_{\star}$$\sqrt{D}$ is constant for microclumping, we account for this
factor of 2 uncertainty in $D$ in our $\dot{M}_{\star}$ measurements reported
below.

The PoWR code offers the option to include X-rays in the wind (Baum et al.,
1992), but our analysis revealed that there was no need to include X-rays to
match the wind-sensitive metal resonance lines in the observed FUV spectra.
While X-ray emission from the stellar wind appears ubiquitous for O stars in
the Milky Way and LMC (e.g., Rauw et al., 2015; Crowther et al., 2022), no
observations exist that can constrain the X-ray fluxes of our target stars (or
even of massive stars in the SMC). We thus use the least complex model that
can explain the observations to avoid introducing additional free parameters
and do not include X-rays in our presented models (following, e.g., Rickard et
al. 2022).

For the microtubulent velocity ($\xi$) in the formal integral, we adopt a
depth-dependent approach (Shenar et al., 2015) where $\xi(r)=0.3\,v(r)$, with
a minimum value of $\xi_{\mathrm{min}}$. The resulting element-dependent
Doppler velocity $v_{\text{Dop}}$ used in the observer’s frame calculation is
then calculated via:

$v_{\text{Dop}}(r)=\sqrt{v_{\text{th}}(r)^{2}+\xi(r)^{2}},$ (2)

where $v_{\text{th}}$ is the thermal velocity of the considered element of the
stellar atmosphere.

Detailed model atoms of H, He, C, N, O, Mg, Si, P, and S are included in the
non-LTE calculations, while the vastly more numerous and complex transitions
of Fe-group elements (atomic numbers 21–28) are treated in a superlevel
approach (described in Gräfener et al. 2002). Though our data cannot directly
constrain the abundances of all considered elements, they are important
opacity sources in the stellar atmospheres that must be accounted for to
enable accurate inference of fundamental stellar properties from their
observed spectra (e.g., Sander et al., 2017).

###  3.2 Choice of Initial Parameters for PoWR Models

For each of the three low-$Z$ stars in this analysis, we aim to find a PoWR
atmosphere model that simultaneously matches the HST/COS FUV spectrum,
Keck/KCWI optical spectrum, and HST photometry (Section 2). Initial guesses
for each parameter of interest in our PoWR models are chosen based on existing
information about the target stars to minimize the number of computationally
expensive model iterations needed to find a good fit to the observations. We
adjust these initial parameters, following the procedure in Section 3.3 below,
to find the parameter set that best matches all available data. Each iteration
requires human intervention and is not automated, so the resulting best-fit
models are not sensitive to the initial-guess parameters. We explain our
choice of initial parameters here, and provide the values in Table A1 in
Appendix A.

We initially adopt the stellar properties ($T_{\mathrm{eff}}$, $\log(g)$, and
$L_{\star}$) inferred by Telford et al. (2021) from model SEDs fit to the FUV
spectra and HST photometry (see their Table 7). For the wind properties, we
pick the nearest point in the publicly available SMC-metallicity PoWR
grid555http://www.astro.physik.uni-potsdam.de/~wrh/PoWR/SMC-OB-III/ (the OB-
III grid with “low” mass-loss and wind strength parameter $\log(Q)=-14.0$;
Hainich et al. 2019) and scale their $\dot{M}_{\star}$ and $v_{\infty}$ to the
presumed iron abundance ($X_{\mathrm{Fe}}$) for each of our O-star targets.
For this, we assume the power-law scalings $\dot{M}_{\star}\propto Z^{0.69}$
(Vink et al., 2001) and $v_{\infty}\propto Z^{0.13}$ (Leitherer et al., 1992).

Initial mass fractions of most metals are scaled to 3%, 6%, and 20% of the
solar values for models of LP26 in Leo P, S3 in Sextans A, and A15 in WLM,
respectively, assuming Asplund et al. (2009) solar abundances (except Fe-group
elements, for which we adopt the updated solar abundances in Scott et al.
2015). These choices are informed by the gas-phase oxygen abundances measured
in the host galaxies and by the best-fit SED models from Telford et al.
(2021). The higher $Z$ of A15 compared to WLM’s gas-phase oxygen abundance
(14% $Z_{\odot}$; Lee et al. 2005) is also consistent with the atmosphere
modeling of an O-supergiant star in the same galaxy by Bouret et al. (2015).
The initial He mass fraction ($Y$) is interpolated between primordial
($Y_{\mathrm{p}}=0.2453$; Aver et al. 2022) and solar ($Y_{\odot}=0.2485$)
values assuming a linear scaling with the total metal mass fraction $Z$.

We adopt initial mass fractions that deviate from the solar pattern for the
CNO elements, as galaxies’ gas-phase C/O and N/O ratios are known to correlate
with O/H (though, with substantial scatter). The typical observed abundance
ratios (by number) in metal-poor dwarf galaxies of $\log\mathrm{(C/O)}=-0.7$
and $\log\mathrm{(N/O)}=-1.4$, respectively (Berg et al., 2019), are
substantially lower than the solar values of
$\log\mathrm{(C/O)}_{\odot}=-0.26$ and $\log\mathrm{(N/O)}_{\odot}=-0.86$
(Asplund et al., 2009). We adopt in our initial models the observed gas-phase
abundances of N and O in the three host galaxies. The observed
$\log\mathrm{(N/O)}$ in Leo P, Sextans A, and WLM are $-1.35$, $-1.54$, and
$-1.49$, respectively (Skillman et al., 2013; Kniazev et al., 2005; Lee et
al., 2005), which are all close to the typical value for metal-poor galaxies
quoted above. Since we do not have C abundance measurements for the three host
galaxies, we adopt initial C mass fractions in our PoWR models such that
$\log\mathrm{(C/O)}=-0.7$. For the model of LP26 only, we also adopt the
abundances of S and He measured in its surrounding H ii region (Skillman et
al., 2013; Aver et al., 2022).

We convolve each model spectrum to the instrumental resolution (of COS in the
FUV and of KCWI in the visible part of the spectrum) and apply a rotational
broadening kernel, taking as an initial guess for $v\,\mathrm{sin}(i)$ the
values reported in Telford et al. (2021) from modeling the c iii 1176 Å
feature in the COS spectra. The resolution and SNR of our spectra preclude
using a Fourier-transform technique (e.g., iacob-broad; Simón-Díaz & Herrero,
2014) to measure $v\,\mathrm{sin}(i)$, so we instead estimate this quantity
from the observed profiles of the photospheric He and metal lines. We adopt
$\xi_{\text{min}}=10\,\mathrm{km}\,\mathrm{s}^{-1}$, a typical value for SMC O
dwarfs (e.g., Bouret et al., 2003), as an initial guess for the minimum
microturbulent velocity.

Finally, we apply dust extinction to each model. We model the extinction due
to Milky-Way foreground dust toward each star with a fixed reddening $E(B-V)$
from the Green et al. (2015) dust map (see Table 1) and a Fitzpatrick (1999)
extinction law. We also include the extinction internal to each host galaxy
with $E(B-V)_{\mathrm{host}}$ calculated from the best-fit
$A_{V}^{\mathrm{host}}$ from the SED fitting in Telford et al. (2021) as an
initial guess, and adopting the preferred extinction law from that modeling
(Fitzpatrick 1999 for LP26 and S3, and Gordon et al. 2003 for A15).

Table 3: Stellar and wind parameters in the best-fit PoWR models. Parameter |  LP26 |  S3 |  A15
---|---|---|---
$T_{\mathrm{eff}}$ (kK) |  $34\pm 1$ |  $31\pm 1$ |  $36\pm 1$
$\log(g/\mathrm{cm\,s}^{-2})$ |  $4.15\pm 0.1$ |  $3.9\pm 0.1$ |  $4.0\pm 0.05$
$\log(L_{\star}/L_{\odot})$ |  $5.0\pm 0.09$ |  $5.15\pm 0.05$ |  $5.15\pm 0.06$
$M_{\star}^{\mathrm{spec}}$ ($M_{\odot}$) |  $43^{+32}_{-18}$ |  $49^{+30}_{-19}$ |  $34^{+15}_{-10}$
$\log(\dot{M}_{\star}/M_{\odot}\,\mathrm{yr}^{-1})$ |  $-$$8.8\pm 0.35$ |  $-$$9.1\pm 0.35$ |  $-$$8.4\pm 0.35$
$v_{\infty}$ (km s-1) |  $500\pm 200$ |  $300\pm 200$ |  $1100\pm 300$
$v\,\mathrm{sin}(i)$ (km s-1) |  $200\pm 80$ |  $200\pm 50$ |  $45\pm 20$
$\xi_{\mathrm{min}}$ (km s-1) |  $5\pm 5$ |  $5\pm 5$ |  $5\pm 5$
$X_{\mathrm{He}}$ |  ($0.25\pm 0.1$) |  ($0.25\pm 0.05$) |  $0.35\pm 0.05$
$\log(X_{\mathrm{C}})$ |  ($-4.58\pm 0.3$) |  ($-4.29\pm 0.3$) |  ($-3.92\pm 0.3$)
$\log(X_{\mathrm{N}})$ |  ($-5.17\pm 0.7$) |  ($-5.06\pm 0.7$) |  $-3.64\pm 0.3$
$\log(X_{\mathrm{O}})$ |  ($-3.76\pm 0.7$) |  ($-3.46\pm 0.7$) |  ($-3.09\pm 0.3$)
$\log(X_{\mathrm{Fe}})$ |  ($-4.44\pm 0.3$) |  ($-4.14\pm 0.3$) |  ($-3.62\pm 0.3$)
$E(B-V)$host (mag) |  $0.000\pm 0.01$ |  $0.056\pm 0.01$ |  $0.022\pm 0.01$

###  3.3 Procedure for Fitting PoWR Models to Observations of the Metal-Poor O
Stars

Starting from the initial model for each star described above, we assess by
eye the agreement of the calculated model spectrum with the observations, and
adjust the model parameters iteratively until the model spectrum acceptably
reproduces diagnostic features in the data. We first find the
$T_{\mathrm{eff}}$ and $\log(g)$ that produce a good match to the he i and he
ii line strengths and Balmer line wings in the optical spectrum. Metal and he
ii lines in the FUV spectra provide secondary diagnostics of these fundamental
stellar parameters, especially $T_{\mathrm{eff}}$. Because these lines are
also sensitive to He and metal abundances, we simultaneously adjust mass
fractions of individual elements as required to match all observed
photospheric transitions of that element. If there are no detected lines of a
given element, or there is no clear need to adjust the abundance of an element
to be consistent with the data, we adopt our initial guess for its abundance
in our final model. Given the lines detected in our observations, we
ultimately only vary the abundances of He, C, N, O, and Fe in our PoWR
modeling described in Section 3.3 below.

We also adjust $v\,\mathrm{sin}(i)$ to reproduce the observed profiles of the
photospheric He and metal lines. Rotational broadening strongly affects these
line shapes and the depths of the Balmer absorption lines, but importantly,
does not affect the wings of the Balmer lines that we use to determine
$\log(g)$. We then adjust $\xi_{\text{min}}$ to match the observed profiles of
optical he i lines (only in the case of A15, which has low
$v\,\mathrm{sin}(i)$ resulting in line profiles narrow enough to constrain
microturbulence) and observed strengths of metal absorption lines in the FUV
(for S3 and LP26; discussed further in Section 4 below).

Next, we adjust two parameters that affect the shape and normalization of the
model: $E(B-V)_{\mathrm{host}}$ and $L_{\star}$. The observed FUV spectrum is
normalized by the emergent model continuum. The high density of metal
absorption lines in the FUV makes it challenging to determine the continuum
level directly from the data. The quality of the FUV normalization and its
agreement with the normalized model spectrum is therefore a sensitive
diagnostic of $L_{\star}$ and the total (foreground + host galaxy) dust
extinction. We adjust $E(B-V)_{\mathrm{host}}$ and $L_{\star}$ iteratively to
achieve the best possible match to the observed FUV continuum shape and HST
NUV, optical, and NIR photometry.

Finally, we find the wind parameters required to match the observed metal
resonance line profiles in the FUV. The terminal velocity $v_{\infty}$
determines the blueward extent of the broad absorption due to the wind, while
$\dot{M}_{\star}$ sets the strength of the absorption or (weak, unsaturated)
P-Cygni profiles. The n iii 4634–4640 Å emission, which is only detected in
the spectrum of A15, could potentially provide an additional constraint for
our models, but this emission is highly sensitive to the details of line
blanketing and non-LTE effects (Rivero González et al., 2011) rather than
providing a direct wind diagnostic.

Once the set of model parameters that best matches the full suite of
observations for each star is identified, we run additional PoWR models to
estimate uncertainties. We adjust one model parameter at a time, changing each
higher and lower from the preferred value until the resulting model is
inconsistent with the observations. Though this method does not capture the
degeneracies among the various parameters, it is sufficient to approximate the
range of allowed stellar properties given the fundamental assumption that the
stars do not have companions and the choices made in the PoWR modeling.

##  4 Results: PoWR Models of the Metal-Poor O Stars

Here we present the PoWR models that best match the observed FUV and optical
spectra and HST photometry of the three target stars. We discuss each model in
turn and report the corresponding best-fit stellar and wind parameters in
Table 3. Figures 3, 4, and 5 compare our preferred PoWR models (red lines) for
LP26, S3, and A15, respectively, to the observations (black points and lines)
of each star. In all three figures, the FUV spectrum is shown in the three
panels of the left column; the stellar SED is shown in the top-right panel;
and the optical spectrum is shown in the middle-right and bottom-right panels.
The FUV spectra are included in the stellar SED plots because they extend the
wavelength range sampled by the HST photometry.

###  4.1 LP26 in Leo P

Figure 3: PoWR model of LP26 in Leo P. The best PoWR model (red) compared to
the observations of LP26 (black). Left column: The FUV spectra are plotted
with the observed wavelength range divided across the three panels, normalized
by the emergent model continuum. Prominent ISM absorption, nebular emission,
and airglow (1215, 1302–1306 Å) features are apparent, but only Ly$\alpha$
absorption is included in the PoWR model. Right column: The top panel shows
the absolute flux of the stellar SED plotted as a function of wavelength (both
log-scaled), including both the observed FUV spectrum (black line) and HST
photometry (black circles). Error bars on the photometry are smaller than the
point size. The middle and bottom panels show the normalized optical spectra,
with the observed wavelength range divided across the two panels. Again,
narrow nebular emission is obvious in the Balmer lines, and for this star
only, fills in all stellar He i absorption lines.

LP26 in Leo P is the most distant and faintest of the three stars; resides in
the lowest-$Z$ host galaxy; and powers a bright H ii region. The continuum SNR
of $\sim$10 in the FUV spectrum and $\sim$30 in the optical spectrum combined
with high $v\,\mathrm{sin}(i)$ prohibit detection of many weak diagnostic
lines in this very metal-poor star, and strong nebular emission lines
unfortunately fill in most of the He i transitions in the optical. Despite
these challenges, we are able to constrain many key stellar properties thanks
to clearly detected FUV absorption features (especially c iii 1176 Å, c iv
1550 Å, and he ii 1640 Å) and well-resolved Balmer line wings and He ii
absorption in the KCWI data (Figure 3).

The optical and FUV He ii lines favor a $T_{\mathrm{eff}}$ of 34 kK, within
the uncertainties of the $37.5^{+5.9}_{-5.5}$ kK previously estimated from SED
fitting by Telford et al. (2021). The lower best-fit $T_{\mathrm{eff}}$
reported here may be more consistent with the hardness of the stellar ionizing
spectrum implied by the relative strengths of nebular H$\beta$ and he i 4471 Å
emission from the surrounding H ii region (Telford et al., 2023); we return to
this point in Section 5.4.3 below. The observed he ii 4686 Å absorption is
substantially weaker than predicted by our best PoWR model and has an
asymmetric profile, both of which are puzzling. While this is the most
susceptible of the optical He ii transitions to infilling by nebular and/or
stellar wind emission, such effects should also impact the FUV he ii 1640 Å
line, which is well matched by the model. It is also likely that nebular
continuum emission from the H ii region contributes to the optical spectrum
(on the order of 10% of the observed continuum level), so dividing out the
continuum level could create the appearance of weaker optical he ii absorption
than the lines actually formed in the stellar photosphere. Yet, tests with
model spectra showed that this effect cannot account for the combination of
well-matched he ii 4541 Å and weak observed he ii 4686 Å absorption.

The Balmer line profiles favor a relatively high $\log(g)$ of 4.15, typical
for dwarf stars and consistent with the position of LP26 on the main sequence
in the HST color-magnitude diagram of Leo P (McQuinn et al., 2015). This value
is consistent within the uncertainties with the $\log(g)=4.0\pm 0.2$ reported
by Telford et al. (2021). The H$\beta$ stellar absorption feature is more
dramatically filled in by nebular emission than H$\delta$, but the wings of
the three Balmer lines are all well matched by a single $\log(g)$. We find
that the best overall match to the shapes of the FUV and optical line profiles
is given by a $v\,\mathrm{sin}(i)$ of 200 km s-1, lower than the $370\pm 90$
km s-1 measured from the FUV c iii 1176 Å profile (Telford et al., 2021), but
the shapes of the Balmer wings in the PoWR models are not substantially
different for those two different $v\,\mathrm{sin}(i)$ values. We attribute
the need for a lower $v\,\mathrm{sin}(i)$ here to a combination of a
difference in intrinsic line broadening between the PoWR model and the tlusty
models used in Telford et al. (2021), and the additional information provided
by the He i and He ii line profiles in the KCWI data included in this
analysis. We recommend adopting the updated $v\,\mathrm{sin}(i)$ reported here
because it is informed by additional data compared to the Telford et al.
(2021) value, but emphasize that (1) the estimated uncertainty of $\pm 80$ km
s-1 is large, due to low SNR and broad lines, and (2) the key conclusion that
this star is a fast rotator remains unchanged.

Figure 4: PoWR model of S3 in Sextans A. As in Figure 3, the best model (red)
compared to observations of S3 (black).  Figure 5: PoWR model of A15 in WLM.
As in Figure 3, the best model (red) compared to observations of A15 (black).
The narrow absorption lines in the spectra indicate a lower
$v\,\mathrm{sin}(i)$, making the metal lines in the FUV spectrum more clearly
discernible than in the spectra of the other two stars at the
$\mathrm{SNR}\sim 10$ of the HST/COS data.

We find that $\log(L_{\star})=5.0$ with only foreground Milky Way dust
extinction (Green et al., 2015), and no contribution from dust in the host
galaxy, produces the best match to the shape of LP26’s SED and FUV spectrum.
This is consistent with the $\log(L_{\star}/L_{\odot})=5.1\pm 0.2$ and host-
galaxy extinction $A_{V}^{\mathrm{host}}=0.01^{+0.06}_{-0.01}$ inferred from
modeling the stellar SED (Telford et al., 2021), and with the small dust
content implied by the Balmer decrement of the H ii region emission (Skillman
et al., 2013; Telford et al., 2023).

Few metal lines are clearly detected in the COS and KCWI spectra of LP26,
preventing robust constraints on its metal abundances. For our initially
adopted abundance set, we do not find strong discrepancies between the
observed and modeled metal lines for any ions except for c iii. Given the
best-fit $T_{\mathrm{eff}}$ and $\log(g)$, the PoWR model predicts a
substantially stronger c iii 1176 Å feature than is observed, but we found
that even decreasing the C abundance to the lowest
$\log\mathrm{(C/O)}\simeq-1.1$ observed in metal-poor dwarf galaxies did not
resolve the discrepancy. Instead, decreasing $\xi_{\mathrm{min}}$ from 10 km
s-1 to 5 km s-1 while keeping the initial $\log\mathrm{(C/O)}=-0.7$ was able
to roughly reproduce the observed strength of the c iii 1176 Å line, and did
not affect the appearance of any other spectral features in the model. This
result is consistent with low microturbulent velocities previously measured
for metal-poor O V stars in the Magellanic Bridge (Ramachandran et al., 2021).

The metal resonance lines in the FUV spectrum of LP26 indicate a weak stellar
wind. The n v 1238, 1242 Å and si iv 1393, 1402 Å doublets appear to be purely
photospheric, though the si iv profile is affected by either narrow (non-wind)
emission or a high noise level. The c iv 1548, 1550 Å profile is challenging
to interpret due to the superimposed narrow Milky Way ISM absorption lines,
but it appears to be broadened due to a stellar wind, with a blueward extent
several pixels beyond the bluest narrow ISM absorption line. We find that a
model with $\log(\dot{M}_{\star})=-8.8$ and a quite low $v_{\infty}=500$ km
s-1 produces the best match to the observed c iv 1548, 1550 Å profile. This
best-fit $v_{\infty}$ is substantially lower than expected from extrapolating
observations of the winds of more metal-rich stars to 3% $Z_{\odot}$ (the
initial guess for $v_{\infty}$ was 1720 km s-1 for LP26), as we discuss in
Section 5.3 below.

###  4.2 S3 in Sextans A

The FUV spectrum of the other very metal-poor star in our sample, S3 in
Sextans A, has the lowest SNR of the three, but the strongest photospheric and
wind lines are still clearly detected. On the other hand, the KCWI spectrum of
S3 is the highest-quality of the three, reaching SNR of $\sim$50 in the
continuum. S3 does not reside inside a bright H ii region, but the Balmer
lines clearly show nebular emission from diffuse ionized gas, and slight
nebular infilling is discernible in some He i lines (particularly 4471 Å). The
optical He i and He ii lines imply a $T_{\mathrm{eff}}$ of 31 kK for this
star, consistent with the $32.5^{+2.9}_{-3.0}$ kK inferred from its observed
SED (Telford et al., 2021). This agrees nicely with the prominent si iii
absorption at 1294–1300 Å (additional, redder transitions are obscured by
geocoronal emission lines) and multiple strong c iii lines in the FUV
spectrum. Not all optical He lines are exactly matched by the model, but we
checked that changing neither $T_{\mathrm{eff}}$ nor the He abundance can
improve the model’s agreement with the data. We attribute this to infilling by
nebular He i emission, or possibly wind emission for the he ii 4686 Å line.

The Balmer line wings are best matched by a $\log(g)$ of 3.9, indicative of an
unevolved star still on the main sequence and significantly higher than the
$\log(g)=3.5^{+0.1}_{-0.2}$ estimated from modeling the stellar SED. This
higher $\log(g)$ is further supported by the clearly detected FUV si iii
absorption at 1294–1300 Å, which becomes much weaker in lower-$\log(g)$
models. As for LP26, we find that adopting $v\,\mathrm{sin}(i)$ of 200 km s-1
gives the best match the observed profiles of optical and FUV photospheric
absorption lines. This value is lower than the $290\pm 90$ km s-1 inferred
from the c iii 1176 Å profile (Telford et al., 2021), but consistent within
the uncertainty. The choice of $v\,\mathrm{sin}(i)$ affects the depths of the
Balmer lines more than the shapes of the wings, thus, adopting the higher
$v\,\mathrm{sin}(i)$ would not change our best-fit $\log(g)$. We find that the
SED shape is best matched by $\log(L_{\star})=5.15$ and total $E(B-V)=0.085$
(0.029 Milky Way foreground; Green et al. 2015), assuming a Fitzpatrick (1999)
extinction law in Sextans A, both consistent with previous SED fitting results
(Telford et al., 2021). However, this $L_{\star}$ is unexpected in combination
with the lower $T_{\mathrm{eff}}$ and higher $\log(g)$ constrained by the KCWI
spectrum, as we discuss in Section 5.1 below.

We did not find any need to adjust the metal abundances from the initial
guesses to match the observed metal lines in the FUV spectrum and the He lines
in the optical for S3. Similar to our findings for LP26 (Section 4.1),
adopting a lower $\xi=5$ km s-1 provides a good match to the observed c iii
FUV lines given the adopted $\log\mathrm{(C/O)}=-0.7$ and best-fit
$T_{\mathrm{eff}}$ of 31 kK, and is more physically motivated than
significantly decreasing the C mass fraction. Again, the FUV n v and si iv
lines appear purely photospheric, but the c iv 1548, 1550 Å line is clearly
affected by the stellar wind, and is best matched by a low $v_{\infty}=300$ km
s-1 and $\log(\dot{M}_{\star})=-9.1$.

###  4.3 A15 in WLM

A15 in WLM is the closest and brightest of the three stars. It is also the
only one with low $v\,\mathrm{sin}(i)$ and therefore narrow lines which are
more easily detectable above the noise level in the FUV and optical spectra.
The best-fit $T_{\mathrm{eff}}$ and $\log(g)$ of 36 kK and 4.0, respectively,
agree nicely with the values reported in Telford et al. (2021) from modeling
the stellar SED ($37.5^{+7.2}_{-2.7}$ kK and $4.0^{+0.2}_{-0.1}$), though we
find that a slightly lower $L_{\star}$ of 10${}^{5.15}\,L_{\odot}$ better
matches the data than their reported
$\log(L_{\star}/L_{\odot})=5.3^{+0.2}_{-0.1}$. A $v\,\mathrm{sin}(i)$ of 45 km
s-1 was required to match the observed FUV and optical line profiles, which
again is consistent with, but lower than, the $80\pm 45$ km s-1 measured from
the c iii 1176 Å profile by Telford et al. (2021). The low
$v\,\mathrm{sin}(i)$ of this star enabled us to constrain
$\xi=5\,\mathrm{km\,s}^{-1}$ based on the narrow observed profiles of the
optical He lines and FUV metal lines (c iii 1247 Å, S v 1502 Å). This is the
same lower $\xi$ value that was required to match the strength of the c iii
1176 Å profiles in the COS spectra of LP26 and S3, but determined by a
different method, suggesting that low $\xi$ is common among metal-poor O
dwarfs.

This star’s lower $v\,\mathrm{sin}(i)$ and higher $Z$ compared to the other
two enabled tighter abundance constraints. We found that the N and He
abundances both had to be increased to match simultaneously the observed
strengths of all the FUV and optical transitions of those elements. The
enhanced N and He abundances are consistent with CNO-cycle products exposed at
the stellar surface. This suggests that A15 is either evolved or may be mixed
by strong rotation. The latter is only possible if A15 is rotating
substantially faster than the low observed $v\,\mathrm{sin}(i)$ (i.e., if we
are viewing the star at a low inclination). However, we find no evidence of C
depletion, as would be expected for CNO-cycle products (C is quickly depleted
by the CN cycle, but O destruction occurs more slowly; e.g., Maeder et al.
2014). The factor of 2 uncertainty in our inferred C and N abundances implies
a minimum $X_{\mathrm{C}}/X_{\mathrm{N}}\geq 0.13$, well above the CN-
equilibrium value (0.025; Maeder 1987). Thus, while we do not observe
abundances reflecting CNO-cycle equilibrium in A15, its high $X_{\mathrm{N}}$
fits in with the observed trend of stronger N enhancement in O stars at lower
$Z$ (e.g., Martins et al., 2024). The finding that unevolved O dwarfs in the
SMC, similar to A15, are already N-enriched is consistent with the expectation
of more efficient chemical transport in compact, metal-poor stars (Meynet &
Maeder, 2002; Brott et al., 2011).

$\dot{M}_{\star}$ and $v_{\infty}$ were constrained from the n v 1238, 1242 Å
and c iv 1548, 1550 Å wind profiles in the COS spectrum of A15 (s iv 1393,
1402 Å again appears photospheric). We prioritized matching the blueward
extent and depth of the wind lines, but no combination of $\dot{M}_{\star}$
and $v_{\infty}$ could also reproduce the peaks of the n v emission, given the
fixed parameters adopted in our models. The best-fit wind parameters are
$v_{\infty}=1100$ km s-1 and $\log(\dot{M}_{\star})=-8.4$, both higher than
found for the two more metal-poor stars in our sample, as expected for line-
driven winds. The best-fit stellar and wind parameters do not reproduce the
observed weak n iii 4634–4640 Å emission in the optical, but this is not
concerning because the these lines are highly sensitive to non-LTE effects and
are therefore not a reliable wind diagnostic.

##  5 Discussion: Comparison to Theoretical Models at Low Metallicity

Next, we place the inferred stellar and wind properties of the three stars
into the context of theoretical expectations at very low $Z$. We compare the
best-fit parameters of the PoWR model for each star presented in Section 4
above to stellar evolution models, $\dot{M}_{\star}$ and $v_{\infty}$
scalings, and predicted ionizing spectra.

Figure 6: Comparison to genec stellar evolution models. Left: $L_{\star}$
plotted as a function of $T_{\mathrm{eff}}$ (both on a logarithmic scale).
Stellar evolution model tracks at 3% $Z_{\odot}$ from Groh et al. (2019) are
shown as gray lines for initial stellar masses of 12, 15, 20, 25, and 32
$M_{\odot}$. Solid lines are non-rotating models, while dashed lines are
models initially rotating at $0.4\,v_{\mathrm{crit}}$. The two stars in our
sample with metallicities closer to 3% $Z_{\odot}$, LP26 and S3, are shown as
purple and blue circles, respectively. Error bars indicate the uncertainties
on our best-fit parameters, estimated as the range of values that produces
PoWR models consistent with the observations of each star (see Section 3.3 and
Table 3). Right: Same as the left panel, except that stellar evolution models
at 14% $Z_{\odot}$ from Georgy et al. (2013) are shown. These higher-$Z$
models are shifted $\sim$0.03 dex toward cooler $T_{\mathrm{eff}}$ compared to
the models in the left panel. A15 in WLM is indicated by the red circle with
error bars.

###  5.1 Stellar Evolution Models

In Figure 6, we compare the properties inferred from our PoWR modeling of the
three low-$Z$ stars to stellar models calculated with the Geneva Stellar
Evolution Code (genec) at 3% (left panel; Groh et al. 2019) and 14%
$Z_{\odot}$ (right panel; Georgy et al. 2013). The Geneva models are widely
used in SPS codes (e.g., Starburst99; Leitherer et al. 1999, 2014) and have
been computed for both rotating and non-rotating stars. In the rotating
models, stars are assumed to initially rotate at $v=0.4v_{\mathrm{crit}}$,
where $v_{\mathrm{crit}}$ is the critical velocity at which the star breaks
apart. This yields rotation speeds between 110–220 km s-1 on the main sequence
(Ekström et al., 2012), comparable to the $v\,\mathrm{sin}(i)$ inferred for
the two fast-rotating stars S3 and LP26.

Model tracks of single stars in the Hertzprung-Russell diagram (HRD; i.e.,
$L_{\star}$ vs. $T_{\mathrm{eff}}$) are shown as gray lines in Figure 6, with
solid and dashed lines indicating non-rotating and rotating models,
respectively. The initial mass of each model is labeled at the left of each
panel, near the zero-age main sequence. Stars evolve toward higher $L_{\star}$
(upward) and lower $T_{\mathrm{eff}}$ (to the right), exiting the main
sequence when they are no longer burning H in their cores (sharp rightward
turns in the model tracks). The two very low-$Z$ stars in our sample, LP26 and
S3, are shown as purple and blue circles, respectively, in the left panel
compared to the 3% $Z_{\odot}$ models of Groh et al. (2019). A15 is shown as a
red circle in the right panel for comparison to the models of Georgy et al.
(2013) at 14% $Z_{\odot}$ (their adopted SMC metallicity).

The locations of the three observed stars in the HRD lie near the evolutionary
tracks with initial masses of $\sim$20 $M_{\odot}$ (for LP26 and S3) and
$\sim$25 $M_{\odot}$ (for A15). This is consistent with the $M_{\star}$
reported by Telford et al. (2021) from modeling the observed stellar SEDs,
constrained by the parsec evolution models (Bressan et al., 2012):
22${}^{+7}_{-5}$, 21${}^{+4}_{-2}$, and 29${}^{+13}_{-3}$ $M_{\odot}$ for
LP26, S3, and A15. Yet, only A15 is hot enough to be on the main sequence of
the considered evolutionary tracks in that mass range, as would be expected
from the relatively high $\log(g)$ of all three stars (3.9–4.15; Table 3).
Both LP26 and S3 appear to have evolved off the main sequence based on their
positions in the HRD alone, but both stars have broad Balmer wings that are
inconsistent with the lower $\log(g)\lesssim 3.5$ expected for post-main-
sequence stars. For stars with more efficient mixing processes, the main
sequence can extend down to $\sim$$25\,$kK (Brott et al., 2011; Higgins &
Vink, 2019). We also compared our results to the parsec (Bressan et al., 2012;
Tang et al., 2014) and BoOST (Brott et al., 2011; Szécsi et al., 2022) stellar
evolution models. For both, the end of the main sequence is indeed shifted to
lower $T_{\mathrm{eff}}$ such that LP26 lies near its end, but S3 still
appears more evolved. Yet, even late-main-sequence stars would not necessarily
appear as dwarfs, though this would be in line with recent predictions by
Martins & Palacios (2021) assuming lower $\dot{M}_{\star}$ than the genec
models.

More striking is the discrepancy between the evolutionary
($M_{\star}^{\mathrm{evol}}$) and spectroscopic ($M_{\star}^{\mathrm{spec}}$)
stellar masses, where

$M_{\star}^{\mathrm{spec}}=L_{\star}\left(\frac{5780}{T_{\mathrm{eff}}}\right)^{4}10^{\log(g)-4.44},$
(3)

with $T_{\mathrm{eff}}$ in K, $g$ in cm s-2, and $M_{\star}^{\mathrm{spec}}$,
$L_{\star}$ in solar units (Heap et al., 2006). We find
$M_{\star}^{\mathrm{spec}}$ of 43${}^{+32}_{-18}$, 49${}^{+30}_{-19}$, and
34${}^{+15}_{-10}$ $M_{\odot}$ for LP26, S3, and A15. These values are higher
than the $M_{\star}^{\mathrm{evol}}$ suggested by comparison to the genec
models for all three stars, and the uncertainties in stellar parameters
estimated from our PoWR modeling allow for minimum $M_{\star}^{\mathrm{spec}}$
of 25, 30, and 24 $M_{\odot}$ for LP26, S3, and A15, respectively. The large
uncertainties on $M_{\star}^{\mathrm{spec}}$ are consistent with the implied
$M_{\star}^{\mathrm{evol}}$ for A15, but there is substantial tension in the
two mass estimates for S3 and LP26. This “mass discrepancy” is well known for
massive stars in the Milky Way and Magellanic Clouds (e.g., Herrero et al.,
1992; Serenelli et al., 2021), but it is more common to find
$M_{\star}^{\mathrm{evol}}$$\,>\,$$M_{\star}^{\mathrm{spec}}$, in the opposite
sense of the result here. Yet, some studies have found
$M_{\star}^{\mathrm{evol}}$$\,<\,$$M_{\star}^{\mathrm{spec}}$ (e.g.,
Ramachandran et al., 2018), especially for high-mass stars $\gtrsim 35$
$M_{\odot}$ (Markova et al., 2018; Bestenlehner et al., 2020). The cause of
the mass discrepancy remains unknown.

The combinations of $T_{\mathrm{eff}}$, $L_{\star}$, and
$M_{\star}^{\mathrm{spec}}$ we infer are challenging to explain, especially
for the lower-$Z$ stars LP26 and S3. Their $T_{\mathrm{eff}}$ are lower than
typical expectations for main-sequence stars of their
$M_{\star}^{\mathrm{spec}}$, suggesting that they could have an unusual
evolutionary history. At the same time, the genec evolution models for stars
above 40 $M_{\odot}$ would predict far higher $L_{\star}$ than we observe
($\log(L_{\star}/L_{\odot})\gtrsim 5.5$). While stars that are observed to be
over-luminous can be explained by advanced evolutionary status, or perhaps
recent interaction with a companion star, it is harder to identify a physical
mechanism that can cause a star to appear under-luminous. Therefore, the
discrepancy we find for these low-$Z$ stars may be pointing to an incomplete
treatment of relevant physics in stellar evolution and/or atmosphere models.
For example, lower-$Z$ stars are expected to be more compact due to lower
metal opacities in their atmospheres, but such model predictions have never
been empirically calibrated against detailed observations of very metal-poor
OB stars. Thus, it is possible that the physics implemented in current stellar
evolution codes result in an incorrect relationship between stellar size and
other fundamental properties in the metal-poor regime. Alternatively, the
Balmer-line wings could be impacted by some unaccounted-for physics that
impacts their utility as a $\log(g)$ diagnostic at very low $Z$.

Another possibility, particularly at the large distances of these target
stars, is that they are unresolved multiple-star systems. Recent years have
seen an accumulation of evidence that the majority of massive stars form in
multiple systems (e.g., Sana et al., 2012; Offner et al., 2023). However,
detecting a companion star in single-epoch spectroscopy with modest SNR and
resolution can be challenging. If a system contains two stars with similar
$T_{\mathrm{eff}}$ and brightness, then there may be nothing unusual about the
line profiles or strengths of $T_{\mathrm{eff}}$-sensitive lines, depending on
the orbital phase when the system was observed. In this case, the only hint of
a multiple system would be a higher $L_{\star}$ than expected from single-star
evolution models or the measured properties of similar spectral-type stars
(e.g., Bouret et al., 2013; Rickard et al., 2022).

At first, this seems like a natural solution to the puzzling properties of
LP26 and S3. But inspection of Figure 6 shows that to bring LP26 and S3 onto
the main sequence of lower-$M_{\star}$ tracks would require a reduction in
$L_{\star}$ of more than a factor of two ($>$ 0.3 dex), implying multiple
similar-mass companions (i.e., triple systems). Moreover, it is unlikely to
find a binary system with two identical spectral types (let alone two such
systems), so we would expect to see spectroscopic signatures of
different-$T_{\mathrm{eff}}$ stars if these systems were unresolved binaries,
but no such discrepancies are observed.

We also checked for changes in the stars’ radial velocities or absorption line
profiles across the individual COS visits and KCWI exposures, and also checked
coadds of the first and second half of the KCWI observations for each star (to
increase the SNR over individual 20-minute exposures). But we found no changes
in stellar absorption-line widths or velocities over time suggestive of
companions. S3 is also included in a recent catalog of Sextans A OB stars
(Lorenzo et al., 2022). Those authors report that the star (s029 in their
naming convention) is a single-lined spectroscopic binary, but state that
significant radial velocity shifts were only detected in the very weak he ii
4541 Å line and do not show the multi-epoch spectroscopy. Thus, S3 could be in
a multiple-star system, but we do not have clear observational evidence in
hand. Altogether, while it is tempting to invoke multiplicity to explain the
unusual combination of stellar properties that we infer, this scenario does
not easily explain all of the observations of LP26 and S3.

Many massive-star binaries will interact during their lifetimes (e.g., Sana et
al., 2012; Eldridge & Stanway, 2022), raising the possibility that we are
observing the products of recent mass-transfer events. To assess whether this
can explain the properties of LP26 and S3, we compared their locations in the
HRD to the predictions of BPASS binary evolution models (Eldridge et al.,
2017), similar to the comparison to single-star models shown in Figure 6. The
BPASS models follow the evolution of the primary (donor) star following mass
transfer, but not the secondary (accretor). We selected models with initial
$M_{\star}^{\mathrm{primary}}$ between 15–25 $M_{\odot}$, comparable to the
$M_{\star}^{\mathrm{evol}}$ of LP26 and S3 implied by our comparison to genec
single-star models, and a representative range of mass ratios (i.e.,
$M_{\star}^{\mathrm{secondary}}/M_{\star}^{\mathrm{primary}}$) and orbital
periods. This exercise revealed no BPASS tracks that could simultaneously
match the positions of LP26 and S3 in the HRD and their high $\log(g)$, ruling
out the possibility that they are the inflated, evolved primary stars in post-
interaction systems.

Both LP26 and S3 are fast rotators, consistent with expectations for secondary
stars that were spun up via mass gain from an evolved companion (de Mink et
al., 2013). Elevated N abundances are not required to explain the observations
(though the uncertainties allow for up to a factor of 5 higher
$X_{\mathrm{N}}$ than the host galaxies’ gas-phase abundances; Table 3), so
recent accretion from a companion star’s unevolved envelope is plausible
(Maeder et al., 2009). Though relatively few models following the evolution of
secondary stars following mass transfer exist in the literature, Renzo &
Götberg (2021) showed that the internal structure of accretor stars can differ
markedly from that predicted by single-star models. Given that the impacts of
accretion on stellar structure and evolution remain poorly understood, further
theoretical work in this area will be crucial to constrain possible formation
pathways of stars with the observed properties of LP26 and S3.

The primary star in such a binary system could potentially have lost enough
mass that it is now a hot, compact stripped star that would appear faint in
the FUV and optical compared to a companion O-type star, and therefore be
challenging to detect in our data (e.g., Götberg et al., 2017). It could also
have exploded as a supernova (or directly collapsed to a black hole; e.g.,
Sukhbold et al. 2016). While there are no known supernova remnants (SNRs) near
LP26 or S3, it is entirely possible that supernovae did explode nearby fairly
recently but have now faded from view (given SNR lifetimes of $\sim 10^{5}$
yr). Yet, Telford et al. (2021) found no large offsets between the stellar
radial velocities and the systemic velocities of the host galaxies that would
suggest that the stars were ejected from the natal system by a supernova kick.
We also lack evidence for a disk of accreted material around either star
(e.g., double-peaked Balmer emission).

Alternatively, the stars could be the products of mergers, which are expected
to occur in $\sim$30% of massive-star binaries (Sana et al., 2012). After a
merger, stars are predicted to undergo a luminous phase during which they
rapidly rotate before shedding mass and angular momentum and contracting back
to the main sequence as slowly rotating magnetic stars (Schneider et al.,
2019). LP26 or S3 could potentially be observed in the over-luminous and fast-
rotating post-merger phase, but this cannot account for their high $\log(g)$.
Moreover, it would be very unlikely to find two stars in such a short-lived
evolutionary phase.

To conclude, we have presented a surprising set of inferred properties for
both very metal-poor O stars in our sample, though the star A15 appears
consistent with single-star evolution models within the uncertainties. It is
not clear that either effectively single or interacting binary companions can
explain all of the available observations of LP26 and S3, though we cannot
definitively rule out that possibility. Additional observations of a larger
sample of metal-poor O stars are required to assess whether these two stars
are outliers or typical at very low $Z$. These results should inform future
development of massive-star evolution models below 20% $Z_{\odot}$.

###  5.2 Stellar Mass-Loss Rates

Stellar evolution codes must account for massive stars’ mass loss over their
lifetimes (e.g., Brott et al., 2011; Bressan et al., 2012; Ekström et al.,
2012). These models typically adopt prescriptions for $\dot{M}_{\star}$ as a
function of stellar properties, particularly $L_{\star}$, $T_{\mathrm{eff}}$,
and $Z$, calibrated by detailed theoretical calculations of line driving in
hot-star atmospheres. Empirical $\dot{M}_{\star}$ measured from optical and
FUV spectra of stars in the Milky Way and Magellanic Clouds have been used to
test these models (e.g., Mokiem et al., 2007; Ramachandran et al., 2019;
Marcolino et al., 2022), but few $\dot{M}_{\star}$ measurements for massive
stars in galaxies below the 20% $Z_{\odot}$ of the SMC have been reported
(e.g., Tramper et al., 2014; Bouret et al., 2015). Moreover, the analysis of
three stars in IC 1613 and WLM by Bouret et al. (2015) suggests that those
galaxies may have SMC-like Fe abundances, consistent with our results for the
star A15 in WLM. Thus, substantial uncertainty remains in the dependence of
$\dot{M}_{\star}$ on $Z$, especially in the very low-$Z$ regime.

Figure 7: Comparison to theoretical mass-loss rates. The power-law
relationships between $\dot{M}_{\star}$, $L_{\star}$, and $Z$ reported by
Krtička & Kubát (2018) and Björklund et al. (2021) are shown as solid and
dashed lines, respectively. The Vink et al. (2001) $\dot{M}_{\star}$
prescription depends upon additional parameters ($T_{\mathrm{eff}}$,
$M_{\star}^{\mathrm{spec}}$, $v_{\infty}$/$v_{\mathrm{esc}}$), so we
calculated that model’s predictions for the inferred stellar parameters for
each of the three stars in our sample, shown as squares. Red, blue, and purple
correspond to metallicities of 20%, 6%, and 3% $Z_{\odot}$ for both model
predictions and for the three stars analyzed here (circles with error bars).
For comparison, we also show the measurements of Bouret et al. (2015) for
three O stars in WLM and IC 1613 as light red diamonds (all three had inferred
$X_{\mathrm{Fe}}$ of 20% solar). The inferred $\dot{M}_{\star}$ for our
low-$Z$ stars lie at the low-$\dot{M}_{\star}$ end of the range of model
predictions.

The stellar and wind properties inferred from our PoWR models enable a first
comparison to the predictions of theoretical mass-loss recipes extending to
substantially sub-SMC metallicities. In Figure 7, we plot $\dot{M}_{\star}$ as
a function of $L_{\star}$ for the three stars in our observational sample as
colored circles, where purple, blue, and red correspond to LP26, S3, and A15,
respectively (same as in Figure 6). For comparison, we show predicted power-
law relationships between $\dot{M}_{\star}$ and $L_{\star}$ at 3%, 6%, and 20%
$Z_{\odot}$ (purple, blue, and red, corresponding to our adopted Fe abundances
for the three observed stars) from Björklund et al. (2021) and Krtička & Kubát
(2018) as solid and dashed lines, respectively. We also calculate the
$\dot{M}_{\star}$ predicted by the widely used Vink et al. (2001) prescription
given the best-fit $L_{\star}$, $T_{\mathrm{eff}}$,
$M_{\star}^{\mathrm{spec}}$, and the ratio of $v_{\infty}$ to the surface
escape velocity ($v_{\mathrm{esc}}$) for each star, shown as squares with
colors encoding $Z$ as for the other two model predictions. If we adopt the
approximate $M_{\star}^{\mathrm{evol}}$ of the three stars instead, the Vink
et al. (2001) $\dot{M}_{\star}$ predictions increase by $\sim$0.2-0.5 dex.

Of the three models to which we compare, the empirical $\dot{M}_{\star}$ we
obtain for these three metal-poor O stars agree best with the Björklund et al.
(2021) mass-loss recipe. The best-fit $\dot{M}_{\star}$ for all three stars
lie well below the predictions across 3–20% $Z_{\odot}$ of both the Vink et
al. (2001) and Krtička & Kubát (2018), except perhaps for LP26, which lies
between the Björklund et al. 2021 and Krtička & Kubát 2018 models.
Importantly, this result holds despite uncertainties due to the possibility of
artificially high $L_{\star}$ due to companion stars: if the stars’ true
$L_{\star}$ are lower than our inferred values by $\sim$0.2–0.6 dex (see
Section 5.1), then the inferred $\log(\dot{M}_{\star})$ would also shift lower
by $\Delta\log(\dot{M}_{\star})=0.75\Delta\log(L_{\star})$ to reproduce the
observed FUV wind line profiles. Even if all three stars shifted toward lower
$L_{\star}$ (to the left in Figure 7) in that scenario, the inferred
$\dot{M}_{\star}$ values would remain lower than predicted by either the Vink
et al. (2001) or Krtička & Kubát (2018) models. Thus, uncertainty in these
parameters cannot bring our results into agreement with the
higher-$\dot{M}_{\star}$ models.

Figure 8: Impact of $\bm{\dot{M}_{\star}}$ and $\bm{v_{\infty}}$ on FUV wind
profiles. Normalized flux plotted as a function of rest-frame wavelength for
the regions surrounding the two main FUV wind diagnostic features: N v 1238,
1242 Å (left column) and C iv 1548, 1550 Å (right column). The observations
for the three stars are shown in black for LP26, S3, and A15 from top to
bottom. Three models are shown in each panel: in red, our best-fit PoWR model
for that star (presented in Section 4); in blue, a model with the same
parameters except for $\dot{M}_{\star}$ set equal to the prediction from the
Vink et al. (2001) mass-loss recipe (Section 5.2); and in purple, a model with
the same parameters as the best fit except for $v_{\infty}$ set equal to the
expectation from extrapolating the empirical $v_{\infty}$ vs.
$T_{\mathrm{eff}}$ scaling of Hawcroft et al. (2023) to lower $Z$ (Section
5.3). The Vink et al. (2001) $\dot{M}_{\star}$ yield FUV wind profiles that
are obviously inconsistent with the observations. Similarly, the $v_{\infty}$
expected from the Hawcroft et al. (2023) relation causes the c iv profiles to
extend too far to the blue. The purple models produce nearly photospheric c iv
absorption for S3 and LP26 that are clearly too narrow compared to the
observed profiles, demonstrating that low, but non-negligible,
$\dot{M}_{\star}$ are required.

The need for low $\dot{M}_{\star}$ to explain the observations of these stars
is emphasized by Figure 8, which shows the n v (left column) and c iv (right
column) wind profiles for the three stars. As in Figures 3–5, the COS spectra
are shown in black and our best-fit PoWR models are shown in red, with the
best $\dot{M}_{\star}$ and $v_{\infty}$ values reported in red text. For
comparison, the blue lines show PoWR models computed with the
$\dot{M}_{\star}$ predicted by the Vink et al. (2001) mass-loss recipe,
reported in blue text (and shown as squares in Figure 7). All other parameters
are held fixed, and we checked that setting $D_{\infty}=1$ to ensure a fair
comparison to the Vink et al. (2001) models produced negligible changes in the
resultant wind profiles. The c iv wind profiles in the blue models with high
$\dot{M}_{\star}$ are clearly ruled out by the observations.

The purple lines in Figure 8 show PoWR models with higher $v_{\infty}$. We
discuss these models in Section 5.3 below, but here remark that for both LP26
and S3, this change causes the c iv profiles to appear nearly photospheric
(i.e., not showing any signature of a stellar wind). The photospheric c iv
lines in the purple models do not reproduce the observed c iv profiles, which
have more extended blueshifted absorption. Thus, low but non-negligible
$\dot{M}_{\star}$ are required to explain the FUV observations of these two
very low-$Z$ stars.

The analysis presented here provides empirical evidence that winds of
low-$L_{\star}$ O stars below 20% $Z_{\odot}$ are best described by
theoretical models with lower $\dot{M}_{\star}$ values. The only other
reported $\dot{M}_{\star}$ from FUV spectra of O stars in galaxies more metal-
poor than the SMC are shown for comparison in Figure 7 (Bouret et al., 2015).
Their three higher-$L_{\star}$ targets in WLM and IC 1613 are shown as light
red diamonds, along with their expected $\dot{M}_{\star}$ from the Vink et al.
(2001) prescription. While those authors inferred SMC-like $X_{\mathrm{Fe}}$
of 20% solar for all three of their targets, they too found lower
$\dot{M}_{\star}$ than predicted by the Vink et al. (2001) model.
Interestingly, there is substantial scatter in their $\dot{M}_{\star}$
measurements, such that the highest-$\dot{M}_{\star}$ star is consistent with
the Krtička & Kubát (2018) model – but this is at the high-$L_{\star}$ extreme
where differences among the models are smaller. Taken together, the
observations of all six O stars in low-$Z$ galaxies agree best with the
Björklund et al. (2021) models and favor lower $\dot{M}_{\star}$ than
typically assumed in the literature.

Our low $\dot{M}_{\star}$ are also consistent with recent studies of OB stars
in the Milky Way and Magellanic Clouds that found evidence for a steeper
decline in $\dot{M}_{\star}$ with $Z$ (Ramachandran et al., 2019) and
$L_{\star}$ (Rickard et al., 2022) than expected from either mass-loss models,
or from earlier empirical studies (e.g., Mokiem et al., 2007). Björklund et
al. (2021) also show that their predicted $\dot{M}_{\star}$ scaling yields a
better match than that of Vink et al. (2001) to empirical $\dot{M}_{\star}$
inferred for SMC O stars by Bouret et al. (2013). Both theoretical and
empirical work have also proposed a $Z$-dependent slope of the
$\dot{M}_{\star}$–$L_{\star}$ relationship (e.g., Björklund et al., 2021;
Marcolino et al., 2022), but its strength (and even sense) remains unsettled.
Our sample is too small and diverse in $Z$ to characterize the dependence of
$\dot{M}_{\star}$ on both $Z$ and $L_{\star}$ in the very metal-poor regime,
but this open question motivates future analysis of FUV stellar wind features
for a larger sample of low-$Z$ O stars spanning a wider range of $L_{\star}$.

The reported uncertainties on $\dot{M}_{\star}$ reflect the range of values
that produce wind profiles consistent with the observations, given the best-
fit values of all other free parameters in our PoWR models, plus an additional
0.15 dex to account for the factor of 2 uncertainty in $D_{\infty}$ (see
Section 3.1). Yet, particularly at low $L_{\star}$, some wind material could
be shocked to high enough temperatures to emit X-rays (e.g., Lucy & White,
1980; Lagae et al., 2021). It is therefore possible that the total
$\dot{M}_{\star}$ across all phases of the wind could be higher than the
$\dot{M}_{\star}$ only in the phase traced by the wind profiles in the FUV
observations. So far, no insights about the strength of wind-intrinsic X-rays
exist at very low $Z$, so the FUV-based analysis presented here represents the
state-of-the-art constraint on their wind properties.

In principle, X-rays from a shocked wind could also impact the modeled FUV
wind profiles. For the temperature regimes of our targets, atmosphere models
use additional X-rays mainly to increase the strength of the n v 1238, 1242 Å
line. Yet, this line is already well reproduced in our models, even without
including additional X-rays. Moreover, independent analyses of the same
observed stars with and without including X-rays in atmosphere models have
yielded similar wind parameters (e.g., Bouret et al., 2013; Rickard et al.,
2022). We conclude that while treating X-rays might be able to improve the
PoWR models’ ability to match simultaneously the n v and c iv wind profiles,
such modifications are unlikely to change our inferred $\dot{M}_{\star}$
enough to reconcile our results with the order-of-magnitude higher
$\dot{M}_{\star}$ predictions of some models seen in Figure 7.

###  5.3 Terminal Wind Velocities

Figure 9: Comparison to extrapolation of empirical terminal wind velocities at
higher $\bm{Z}$. Solid lines show the scaling relation between $v_{\infty}$,
$T_{\mathrm{eff}}$, and $Z$ fit to observations in the Milky Way, LMC, and SMC
by Hawcroft et al. (2023), extrapolated to lower $Z$ (color-coded as in Figure
7). Circles with error bars show the measured $v_{\infty}$ as a function of
$T_{\mathrm{eff}}$ for the three stars analyzed here. For comparison, we also
show measurements for O stars in WLM and IC 1613 (all with inferred SMC-like
$X_{\mathrm{Fe}}$) drawn from Garcia et al. (2014) and Bouret et al. (2015) as
light red squares and diamonds, respectively. We find lower $v_{\infty}$ for
all three metal-poor stars studied here than expected from the best-fit
relation to observations at 20–100% $Z_{\odot}$, though our result for A15 in
WLM lies within the scatter of published $v_{\infty}$ of individual stars at
similar $Z$ and $T_{\mathrm{eff}}$ in metal-poor galaxies beyond the SMC.

A complete physical understanding of radiation-driven winds, particularly
their momenta and mechanical feedback deposited into the surrounding ISM,
depends critically on $v_{\infty}$ (e.g., Castor et al., 1975; Weaver et al.,
1977). For Milky Way OB stars, $v_{\infty}$ has long been known to correlate
with both $v_{\mathrm{esc}}$ $\left(\propto\sqrt{M_{\star}/R_{\star}}\right)$
and $T_{\mathrm{eff}}$ (Kudritzki & Puls, 2000). Its predicted weak $Z$
dependence, $v_{\infty}\propto Z^{0.13}$ (Leitherer et al., 1992), is
generally found to be consistent with observations. For example, recent
empirical studies of $v_{\infty}$ measured from the blueward extent of FUV
wind-line profiles of OB stars in the Milky Way, LMC, and SMC reported similar
$v_{\infty}\propto Z^{0.08-0.22}$ (Marcolino et al., 2022; Hawcroft et al.,
2023). Yet, like for $\dot{M}_{\star}$, few $v_{\infty}$ measured from FUV
wind profiles have been reported for OB stars in lower-$Z$ galaxies than the
SMC, leaving open the question of whether the same behavior holds in the
extremely metal-poor regime.

Figure 9 shows the best-fit $v_{\infty}$ as a function of $T_{\mathrm{eff}}$
for our three metal-poor stars as circles (colored by adopted Fe abundance as
in Figures 6 and 7). For comparison, we show the $v_{\infty}$ vs.
$T_{\mathrm{eff}}$ relationships of Hawcroft et al. (2023) at the adopted $Z$
of each star as solid lines, extrapolated from the empirical relationship
between $v_{\infty}$, $Z$, and $T_{\mathrm{eff}}$ fit to observed OB stars at
20–100% $Z_{\odot}$. Shaded regions indicate the uncertainty in their model
parameters, not the typical scatter in individual $v_{\infty}$ measurements
($\pm$500 km s-1), which is larger than the shaded regions. For all three
metal-poor stars, the inferred $v_{\infty}$ lies well below the expectation
for their $T_{\mathrm{eff}}$ and $Z$ from the Hawcroft et al. (2023)
relationship.

Figure 8 illustrates that the $v_{\infty}$ values expected from the empirical
scaling for Milky Way, LMC, and SMC O stars is incompatible with the
observations of our metal-poor stars. The purple lines show PoWR models with
the same parameters as our best-fit models (red lines), except that
$v_{\infty}$ has been set to the value calculated from the Hawcroft et al.
(2023) relation for the star’s $T_{\mathrm{eff}}$ and $Z$. As discussed in
Section 5.2, the higher $v_{\infty}$ at fixed $\dot{M}_{\star}$ changes the c
iv profiles in the LP26 and S3 models such that they are nearly photospheric.
Yet, for S3 (middle row), the purple line is discernibly below the red line
just blueward of the c iv blue edge in the best model. If $\dot{M}_{\star}$ is
increased, that blueward absorption becomes more pronounced, producing a c iv
profile that is clearly affected by the wind, and also inconsistent with the
data (black lines). The same is true for the LP26 model; both extremely
low-$Z$ stars require low $v_{\infty}$ to reproduce the blueward extent of
their observed c iv profiles. For A15, the Hawcroft et al. (2023) $v_{\infty}$
produces broader n v and c iv wind profiles; the former is a comparably good
match to the data as the best-fit model, but the c iv absorption extends
beyond the observed blue edge. Like for the other two stars, $\dot{M}_{\star}$
would need to be increased to compensate for the higher $v_{\infty}$,
producing an even more discrepant c iv profile, so the higher $v_{\infty}$ is
ruled out. Still, our best-fit $v_{\infty}$ of 1100 km s-1 for A15 is
consistent within the scatter of measurements for individual SMC stars in
Hawcroft et al. (2023).

To place our results in context, Figure 9 also shows published $v_{\infty}$ of
11 OB stars in the metal-poor galaxies WLM and IC 1613 inferred from FUV
spectroscopy. Measurements from Garcia et al. (2014) and Bouret et al. (2015)
are shown as light red squares and diamonds, respectively. Both papers report
that all of their targets have SMC-like $X_{\mathrm{Fe}}$ and sample
luminosity classes from dwarfs to supergiants. As expected for that Fe
abundance, the empirical $v_{\infty}$ from the literature scatter about the
best-fit $v_{\infty}$–$T_{\mathrm{eff}}$–$Z$ relationship of Hawcroft et al.
(2023) at 20% $Z_{\odot}$. Our best-fit $v_{\infty}$ for the star A15 in WLM,
with the same inferred $X_{\mathrm{Fe}}$, lies below the empirical
$v_{\infty}$ vs. $T_{\mathrm{eff}}$ for the SMC, but interestingly, lies
within a group of stars in IC 1613 that scatter to lower $v_{\infty}$ than the
overall trend. Though the sample size is small, this cluster of
lower-$v_{\infty}$ stars around $\sim$35–38 kK is striking. Future
observations of a larger sample of OB stars in galaxies with gas-phase O
abundances below 20% $Z_{\odot}$ would help to clarify whether this decrease
in typical $v_{\infty}$ is real, or if the scatter across the full
$T_{\mathrm{eff}}$ range is larger and extends to lower $v_{\infty}$ than
suggested by the currently available data.

Radiation-driven wind theory predicts a linear scaling between $v_{\infty}$
and $v_{\mathrm{esc}}$ set by the force multiplier parameter, $\alpha$ (e.g.,
Kudritzki et al., 1989; Lamers & Cassinelli, 1999). Empirical studies of
Galactic hot stars ($\geq 21$ kK) have found a ratio of
$v_{\infty}$/$v_{\mathrm{esc}}$$\,\simeq\,$2.65, with lower ratios observed
for cooler stars with winds driven by different ions (Lamers et al., 1995;
Kudritzki & Puls, 2000). However, it is unclear whether a constant
$v_{\infty}$/$v_{\mathrm{esc}}$ ratio holds at lower $Z$ (Garcia et al., 2014;
Hawcroft et al., 2023). Given the $v_{\mathrm{esc}}$ in our PoWR models of the
three metal-poor stars (1300, 1230, and 1100 km s-1 for LP26, S3, and A15,
respectively), we find $v_{\infty}$/$v_{\mathrm{esc}}$ of 0.38, 0.24, and 1.0.
These values are well below the canonical ratio of 2.65, and lower for the two
very metal-poor stars. This may hint at an extension of the trend for lower
$v_{\infty}$/$v_{\mathrm{esc}}$ for stars in the SMC compared to the LMC
reported by Hawcroft et al. (2023).

Echoing the earlier discussion in Section 5.2, these $v_{\infty}$ measurements
are for the phase of the wind that can be probed by FUV spectroscopy. It is
possible that these stars have higher $\dot{M}_{\star}$ and/or $v_{\infty}$ if
super-ionization due to shocks is indeed important in weak winds.
Unfortunately, we do not have observational diagnostics of this scenario
(e.g., X-rays, or high-ionization UV lines like o vi 1032, 1038 Å; Puls et al.
2020). The lack of strong n v wind features in the two stars with lower $Z$
and $\dot{M}_{\star}$ suggests that super-ionization is not important in their
winds, but examples of X-ray-emitting Galactic O stars that also have weak n v
features are known (e.g., Martins et al. 2005; Huenemoerder et al. 2012).
$v_{\infty}$ inferred from unsaturated P-Cygni profiles are often regarded as
lower limits (Lamers et al., 1995; Crowther et al., 2016). Yet, encouragingly,
Hawcroft et al. (2023) showed that even weak-wind stars in their sample follow
the same $v_{\infty}$ scaling relationships as the stars with saturated wind
profiles. Moreover, the stars with $v_{\infty}$ drawn from the literature
shown in Figure 9 show stronger P-Cygni profiles than the three metal-poor
stars analyzed here, and still are consistent with our results. Given the
available data, our analysis provides the best constraints to date on the
strength of stellar winds in the very low-$Z$ regime.

###  5.4 Ionizing Spectra

Figure 10: Comparison of model ionizing spectra. Total flux emitted from the
stellar surface plotted as a function of wavelength from 20 Å to 2 $\mu$m,
both on logarithmic scales, for LP26, S3, and A15, ordered from top to bottom.
Vertical gray lines indicate the wavelengths below which photons are energetic
enough to ionize H, He, and He+. The dark red line in each panel shows the
best-fit PoWR model presented in Section 4. For reference, a blackbody
spectrum for the same $T_{\mathrm{eff}}$ and $R_{\star}$ as the best-fit PoWR
model is plotted as the light red line, highlighting that even at very low
$Z$, a blackbody is not an appropriate approximation of an O star’s ionizing
spectrum. The best-fit models to the observed stellar SEDs reported by Telford
et al. (2021) are shown in blue; these ionizing spectra are all harder than
the PoWR models reported here due to their higher $T_{\mathrm{eff}}$ (Table 4;
Section 5.4.2). Finally, the purple models show the theoretical expectations
for stars with similar $T_{\mathrm{eff}}$ and $\log(g)$ to each of the three
observed stars at SMC metallicity, drawn from the Hainich et al. (2019) grid
of PoWR models. The lower $L_{\star}$ of the theoretical SMC models results in
lower ionizing flux compared to the best models of our three metal-poor O
stars.

An important application of stellar evolution and spectral models at very low
$Z$ is estimating the ionizing photon production rate of massive stellar
populations in low-mass galaxies. This impacts our understanding of the role
played by early, chemically unevolved galaxies in the process of cosmic
reionization, as well as feedback and nebular emission in low-$Z$ galaxies at
all cosmic epochs. Yet, because we lack observational constraints on
$\dot{M}_{\star}$ prescriptions, stellar evolution, and the ionizing spectra
of O stars below 20% $Z_{\odot}$, substantial uncertainty remains in the shape
and normalization of massive stellar populations’ ionizing spectra in the very
low-$Z$ regime. The atmosphere models tailored to FUV and optical spectra of
very low-$Z$ O stars presented here enable a first comparison to theoretical
models representative of those adopted in SPS codes to interpret the light
observed from metal-poor galaxies.

####  5.4.1 Comparison to Theoretical Atmosphere Models at SMC Metallicity

Table 4: Comparison to other model SEDs. Parameter | LP26 | S3 | A15
---|---|---|---
PoWR Models (This Work)
$T_{\mathrm{eff}}$ (kK) | 34 | 31 | 36
$\log(g/\mathrm{cm\,s}^{-2})$ | 4.15 | 3.9 | 4.0
$\log(L_{\star}/L_{\odot})$ | 5.0 | 5.15 | 5.15
$M_{\star}$ ($M_{\odot}$) | 43 | 49 | 34
$\log(Q(\mathrm{H)/s}^{-1})$ | 48.26 | 48.13 | 48.63
$\log(Q(\mathrm{He)/s}^{-1})$ | 46.78 | 45.69 | 47.51
$Q$(He)/$Q$(H) | 0.034 | 0.003 | 0.075
SMC OB-III PoWR Models
$T_{\mathrm{eff}}$ (kK) | 34 | 31 | 36
$\log(g/\mathrm{cm\,s}^{-2})$ | 4.2 | 4.0 | 4.0
$\log(L_{\star}/L_{\odot})$ | 4.55 | 4.54 | 4.97
$M_{\star}$ ($M_{\odot}$) | 17 | 15 | 23
$\log(Q(\mathrm{H)/s}^{-1})$ | 47.87 | 47.61 | 48.49
$\log(Q(\mathrm{He)/s}^{-1})$ | 46.38 | 45.43 | 47.39
$Q$(He)/$Q$(H) | 0.033 | 0.007 | 0.080
Telford et al. (2021) SED Models
$T_{\mathrm{eff}}$ (kK) | 37.5 | 32.5 | 37.5
$\log(g/\mathrm{cm\,s}^{-2})$ | 4.0 | 3.5 | 4.0
$\log(L_{\star}/L_{\odot})$ | 5.1 | 5.1 | 5.3
$M_{\star}$ ($M_{\odot}$) | 22 | 21 | 29
$Z$/$Z_{\odot}$, Grid | 2%, W | 3%, V | 20%, S
$\log(Q(\mathrm{H)/s}^{-1})$ | 48.63 | 48.55 | 48.84
$\log(Q(\mathrm{He)/s}^{-1})$ | 47.89 | 47.30 | 48.00
$Q$(He)/$Q$(H) | 0.184 | 0.056 | 0.146

The top section of Table 4 summarizes the key properties that impact ionizing
photon production for our best-fit PoWR models to the three metal-poor O
stars. $T_{\mathrm{eff}}$, $\log(g)$, $\log(L_{\star})$, and $M_{\star}$, all
of which set the shape and normalization of the ionizing SED, are reproduced
from Table 3. We also report the production rate of photons capable of
ionizing H ($Q$(H); $\geq 13.6$ eV) and He ($Q$(He); $\geq 24.6$ eV) in the
models. These models predict $\log(Q(\mathrm{H)/s}^{-1})$ ranging from
48.13–48.63, broadly consistent with previous atmosphere modeling of observed
O7–9 stars in the SMC (Ramachandran et al., 2019; Rickard et al., 2022).

First, we consider how our results compare to expectations from a purely
theoretical grid of model SEDs at low $Z$. The middle section of Table 4
summarizes the properties of models with similar $T_{\mathrm{eff}}$ and
$\log(g)$ to our best-fit models to the three metal-poor stars drawn from the
Hainich et al. (2019) grid of PoWR models at SMC metallicity. Figure 10 also
presents a graphical comparison of the emergent flux integrated over the
stellar surface as a function of wavelength for LP26, S3, and A15, ordered
from top to bottom. Our best-fit PoWR models (the same models shown in the
top-left panels of Figures 3–5) are shown as red lines, and the Hainich et al.
(2019) models with comparable stellar parameters at SMC metallicity are shown
as purple lines. For reference, radiation to the left of the gray vertical
lines is capable of ionizing H, He, and He+ (from right to left).

The most obvious difference between our PoWR models and the SMC-star models in
Figure 10 is the normalization, particularly for LP26 and S3. The PoWR model
grid adopts $L_{\star}$ for a given $T_{\mathrm{eff}}$ and $\log(g)$ from
single-star evolution models (Brott et al., 2011), and these theoretical
$L_{\star}$ are lower than our best-fit values (similar to the discussion in
Section 5.1). $L_{\star}$ strongly affects the normalization of the stellar
SED, and therefore the absolute $Q$(H) and $Q$(He): both are higher in our
PoWR models than in the comparable SMC-star models. The high $L_{\star}$
inferred for our three metal-poor O stars thus imply higher ionizing photon
production rates than would be predicted from purely theoretical models of
stellar evolution and ionizing spectra.

On the other hand, we find good agreement in the hardness of the ionizing
spectra. While small changes in the shapes of the red and purple models are
visually apparent in Figure 10, particularly for He-ionizing photons, the
y-axes span 5 orders of magnitude, so subtle differences where the flux is low
do not substantially change $Q$(He)/$Q$(H). Interestingly, this comparison
suggests little evolution in the hardness of O stars’ ionizing spectra across
a factor of several in $Z$, from 3–20% $Z_{\odot}$.

####  5.4.2 Comparison to SED-Fitting Results

Next, we compare to the ionizing spectra predicted from our earlier modeling
of these stars’ SEDs, published before optical spectra were available to
constrain their $T_{\mathrm{eff}}$ and $\log(g)$. Stellar photometry is
readily available for many metal-poor galaxies in the local universe, whereas
optical spectroscopy of the necessary resolution and SNR is expensive and
therefore rare. Thus, it is important to understand how SED-based constraints
on $Q$(H) and the hardness of ionizing spectra compare to the quantitative
spectroscopic analysis presented here.

The bottom section of Table 4 reproduces the $T_{\mathrm{eff}}$, $\log(g)$,
$\log(L_{\star})$, and $M_{\star}$ estimated for the three metal-poor O stars
by Telford et al. (2021) from fitting their HST/COS FUV spectra and HST NUV-
through-NIR photometry. We also report $Q$(H), $Q$(He), and $Q$(He)/$Q$(H) for
those best-fit SED models, drawn from Telford et al. (2023) for LP26 and newly
calculated for S3 and A15 here by integrating the model fluxes above the
minimum photon energy required to ionize H and He. The updated
$T_{\mathrm{eff}}$ and $\log(L_{\star})$ reported here for our best-fit PoWR
models are consistent within the (large) uncertainties on the SED-based
quantities. However, the SED-based $T_{\mathrm{eff}}$ are also higher for all
three stars than found in our models fit to the Keck/KCWI optical spectra in
particular (and informed by all available observations), which enabled the
Telford et al. (2021) models to account for the high $L_{\star}$ required to
match the normalization of the observed SEDs (because $L_{\star}\propto
T_{\mathrm{eff}}^{4}$) within the constraints of the parsec single-star
evolution models adopted in that analysis.

The impact of systematically higher $T_{\mathrm{eff}}$ in the SED-based models
is apparent in Figure 10. The Telford et al. (2021) models are shown as blue
lines, all of which have higher $Q$(H) and $Q$(He) compared to the best-fit
PoWR models (red lines), though both sets of models reproduce the observed
FUV-through-NIR parts of the SEDs. The SED-based models also predict harder
ionizing spectra (i.e., higher $Q$(He)/$Q$(H)) for all three stars than our
best-fit PoWR models. This tendency toward higher $T_{\mathrm{eff}}$ and
$Q$(He)/$Q$(H) to explain high $L_{\star}$ (as would be observed for
unresolved binaries) should be kept in mind for SED-fitting analyses of
resolved massive-star populations, and underscores the importance of optical
spectroscopy to constrain $T_{\mathrm{eff}}$ for O stars.

####  5.4.3 The Special Case of LP26

Finally, we consider LP26, the one star in our sample for which additional
constraints on the ionizing spectrum are available from its surrounding H ii
region. Analysis of the nebular emission lines in the Keck/KCWI data
constrained $\log(Q(\mathrm{H)/s}^{-1})=48.57^{+0.07}_{-0.09}$ and
$Q(\mathrm{He})/Q(\mathrm{H})=0.12\pm 0.03$ (Telford et al., 2023).
Interestingly, while the SED-based $Q$(H) is consistent with the constraint
from the H ii region emission, the $Q$(H) predicted by our best-fit PoWR model
is lower (Table 4). The H ii-region analysis adopted an escape fraction
$f_{\mathrm{esc}}=0$, so if ionizing photons are escaping the H ii region,
then the true $Q$(H) of the ionizing O star would need to be even higher, and
thus in stronger tension with the results presented here. Moreover, the PoWR
model of LP26 predicts a softer ionizing spectrum than allowed by the observed
H ii region emission. These differences between the best-fit PoWR model (red
line) and tlusty model fit to the stellar SED (blue line) are readily apparent
in the top panel of Figure 10.

Whereas the estimated stellar properties from modeling LP26’s observed SED
were consistent with its H ii-region emission, we now find that the additional
constraints on the stellar $T_{\mathrm{eff}}$ and $\log(g)$ from the stellar
optical spectrum complicate our interpretation. The optical He i and He ii
photospheric absorption lines do not allow for a high enough
$T_{\mathrm{eff}}$ to explain the hardness of the ionizing spectrum powering
the surrounding H ii region. The most plausible solution to this puzzle is
likely a multiple-star system with a history of interaction. For example, a
compact stripped-star companion would be challenging to detect in FUV/optical
light, but could contribute significant ionizing flux (Götberg et al., 2017).
Past mass transfer from that companion star could also have spun up LP26 to
its high $v\,\mathrm{sin}(i)$. However, this scenario does not readily explain
the combination of high $L_{\star}$ and high $\log(g)$; this may require
another luminous star to be part of the system, or some other mechanism that
is not yet understood. Long-term monitoring of LP26 would be useful to
constrain the multiplicity of the system and rule out some possible
explanations for its peculiar properties.

###  5.5 Implications for Modeling Metal-Poor Massive Stellar Populations

Here, we connect our analysis of these three very low-$Z$ O stars to the
broader landscape of massive stellar population modeling in metal-poor
galaxies. Models of stellar evolution, spectra, and SEDs are essential tools
in various areas of astrophysics: they are used to interpret observations of
evolved stars, supernovae, and gravitational-wave events; predict feedback and
nucleosynthesis from massive stars and supernovae; and are fundamental
ingredients in the SPS codes required to infer galaxies’ star-formation
histories and ionizing photon production. Massive-star models informed by
empirical results are thus highly important, but the dearth of observations
has been the limiting factor below 20% $Z_{\odot}$.

This analysis has revealed some interesting tensions between our observations
of very metal-poor O stars and model expectations. The fundamental stellar
properties that we infer for the two stars in galaxies below 10% $Z_{\odot}$
are not simultaneously matched by any publicly available stellar evolution
models, even those that treat rotation or binary interactions (Section 5.1).
Either the stars are under-luminous for their $M_{\star}^{\mathrm{spec}}$,
which is challenging to explain, or they are over-luminous for their apparent
$M_{\star}^{\mathrm{evol}}$. In the latter case, multiple-star systems may be
able to alleviate the tension, but the fine-tuned scenarios that could explain
all observational constraints are unlikely. We also highlight the large impact
of $L_{\star}$ on stellar ionizing photon production (Section 5.4),
emphasizing the need to understand the relationship between $M_{\star}$ and
$L_{\star}$ at very low $Z$ to accurately predict metal-poor galaxies’
ionizing fluxes. Altogether, the inferred properties of LP26 and S3 pose a
challenge to current stellar evolution models in the still-unexplored
extremely low-$Z$ regime.

Our $\dot{M}_{\star}$ and $v_{\infty}$ results (Section 5.2–5.3) paint a
picture of weaker winds than expected at and below 20% $Z_{\odot}$. That we
infer $\dot{M}_{\star}$ more than an order of magnitude lower than the Vink et
al. (2001) recipe, which is often adopted in stellar evolution calculations,
is particularly striking, and implies that metal-poor massive stars’
evolutionary pathways and end products are less affected by mass loss than
current models predict. However, a physical understanding of these weak
stellar winds remains elusive because no theoretical models match all of the
available constraints from observations. For example, while the models of
Björklund et al. (2021) broadly agree with our low inferred $\dot{M}_{\star}$
at low $Z$, they also predict that $v_{\infty}$ increases strongly with
decreasing $\dot{M}_{\star}$, resulting in $\gtrsim 3000$ km s-1 for their SMC
models. This is in tension with the low $v_{\infty}$ we infer for our three
metal-poor stars, as well as empirical results at 20% $Z_{\odot}$ and higher
(e.g., Garcia et al., 2014; Bouret et al., 2015; Hawcroft et al., 2023),
suggesting that further theoretical work is needed to understand the driving
of weak stellar winds at low $Z$. Improved models of radiation-driven winds,
informed by FUV spectra of metal-poor OB stars, are critical to ensure
accurate interpretation of metal-poor galaxies’ evolution, particularly now
that observations of chemically unevolved galaxies both locally and at high
$z$ are becoming increasingly routine.

##  6 Conclusions

We have presented new Keck/KCWI optical spectra of three very metal-poor O
stars, LP26, S3, and A15, in nearby dwarf galaxies (Section 2.1, Figures 1–2).
To determine their stellar and wind properties, we fit PoWR atmosphere models
to these observations combined with published HST/COS FUV spectra and HST NUV,
optical, and NIR photometry. Our key findings are:

  1. 1.

We refine fundamental stellar parameter measurements over previous estimates
from modeling the stellar SEDs (Telford et al., 2021). The $T_{\mathrm{eff}}$,
$\log(g)$, and $L_{\star}$ determined from our PoWR models are all consistent
within the large uncertainties on the SED-based estimates, except that we
require a significantly higher $\log(g)$ to match S3’s broad Balmer wings
($3.9\pm 0.1$, vs. previous $\log(g)=3.5^{+0.1}_{-0.2}$). Yet, we also measure
$T_{\mathrm{eff}}$ systematically lower by 1.5–3.5 kK than the best-fit values
from SED fitting, underscoring the importance of optical spectroscopy for
accurate stellar parameter inference (Section 4; Figures 3–5; Tables 3–4).

  2. 2.

The star A15 has strong FUV and optical He and N lines indicating enhanced
abundances of both elements, but the data rule out the low C/N ratio expected
for pure CNO-cycle products. This moderate N enrichment in a dwarf star is
consistent with efficient mixing at low $Z$. While abundances of the two
lower-$Z$ stars are more uncertain due to their relatively high best-fit
$v\,\mathrm{sin}(i)=200$ km s-1 and lower-SNR FUV spectra, neither shows clear
evidence of elevated N. For all three stars, a low $\xi_{\mathrm{min}}=5$ km
s-1 is required to match the observed c iii 1176 Å (for LP26 and S3) and
optical he i lines (for A15), suggesting that microturbulent velocities may
typically be small for O-dwarf stars at very low $Z$ (Section 4; Table 3).

  3. 3.

The inferred $T_{\mathrm{eff}}$, $\log(g)$, and $L_{\star}$ for all three
stars yield higher $M_{\star}^{\mathrm{spec}}$ than the
$M_{\star}^{\mathrm{evol}}$ implied by single-star evolution models, though
$M_{\star}^{\mathrm{spec}}$ and $M_{\star}^{\mathrm{evol}}$ are consistent
within the uncertainties for A15. No evolution tracks, even those including
binary interactions, are able to simultaneously match all inferred properties
of the extremely metal-poor stars LP26 and S3. Though the possibility of
binaries could alleviate some of the tension, it is challenging to identify a
plausible scenario that can account for both stars’ apparent
$T_{\mathrm{eff}}$, $\log(g)$, and $L_{\star}$, as well as the lack of
evidence for companions in the data. Our observations of these low-$Z$ stars
pose an interesting challenge for current evolution models (Section 5.1;
Figure 6).

  4. 4.

We compare our inferred $\dot{M}_{\star}$ from the observed FUV metal
resonance line profiles to the predictions of radiation-driven wind models,
and find the best overall agreement with the Björklund et al. (2021) scaling
at the low-$\dot{M}_{\star}$ end of the range. For all three stars, the data
constrain $\dot{M}_{\star}$ to more than an order of magnitude lower than that
calculated from the Vink et al. (2001) mass-loss recipe, which is widely
adopted in stellar evolution codes. Our findings agree with previous empirical
results for OB stars in the SMC, and for higher-$L_{\star}$ O stars in more
metal-poor dwarf galaxies. This emerging empirical support for lower
$\dot{M}_{\star}$ at low $Z$ should inform the assumptions in stellar
evolution modeling, as $\dot{M}_{\star}$ strongly affects massive stars’
evolution and end products, which propagate to SPS model predictions (Sections
5.2 and 5.5; Figures 7–8).

  5. 5.

We also compare $v_{\infty}$ inferred from the observed wind profiles,
particularly the c iv 1550 Å feature, to the empirical relationship between
$v_{\infty}$, $T_{\mathrm{eff}}$, and $Z$ at 20–100% $Z_{\odot}$ reported by
Hawcroft et al. (2023), extrapolated to lower $Z$. The best-fit $v_{\infty}$
of all three stars are lower than expected, but A15 lies within the scatter of
published $v_{\infty}$ for 11 other OB stars with similar abundance patterns
in galaxies beyond the SMC. The low empirical $v_{\infty}$ in the SMC and more
metal-poor galaxies, including our results here, are in strong tension with
the very high $v_{\infty}$ predicted by the Björklund et al. (2021) models
that best reproduce the observed $\dot{M}_{\star}$ for the same metal-poor
stars. Improved modeling of radiation-driven winds is needed in the low-$Z$
regime (Sections 5.3 and 5.5; Figures 8–9).

  6. 6.

We compare the ionizing photon production rates, $Q$(H), and hardness of the
ionizing spectra in our best-fit PoWR models to expectations from similar
models at SMC metallicity, and to the Telford et al. (2021) best-fit models to
these stars’ SEDs. The theoretical models for SMC stars, constrained by
single-star evolution model expectations, predicts lower $Q$(H) than our best-
fit PoWR models, emphasizing the need to clarify the relationship between
$L_{\star}$ and $M_{\star}$ at very low $Z$ for accurate predictions ionizing
photon production by metal-poor galaxies. We also show that SED fitting is
subject to biases in $T_{\mathrm{eff}}$, emphasizing the need for optical
spectra to constrain $Q$(H) and ionizing spectral hardness (Section 5.4–5.5,
Figure 10, Table 4).

  7. 7.

Finally, we highlight that our inferred properties of LP26 predict $Q$(H) and
$Q$(He)/$Q$(H) that are inconsistent with previous constraints on the ionizing
source of its surrounding H ii region from the observed nebular emission
(Telford et al., 2023). Possible explanations include one or more companion
stars, past interactions, and/or exotic stellar physics at very low $Z$
(Section 5.4.3).

The combined Keck/KCWI and HST/COS spectra analyzed here have given us a first
look at the stellar and wind properties of very metal-poor O stars in dwarf
galaxies near the edge of the Local Group. Yet, the present sample size is too
small for empirical characterization of stellar wind recipes below 20%
$Z_{\odot}$. The initial results from our PoWR modeling emphasize that we have
much to learn about low-$Z$ stellar astrophysics from these spectra, despite
the limited resolution and SNR that can be achieved for faint and distant O
stars. Additional data are particularly needed to clarify the implications of
our findings for the calculation of star-formation rates and ionizing photon
production in reionization-era galaxies. In future work, we will employ
similar-quality spectroscopy from new COS programs (HST-GO-16767, HST-
GO-17491; PI: O. G. Telford) targeting larger samples of extremely metal-poor
O stars to more rigorously test and calibrate stellar models in this low-$Z$
regime.

OGT thanks Jenny Greene, D. John Hillier, Calum Hawcroft, Jon Sundqvist, Alex
de Koter, Ciarán Furey, Frank Backs, Jan Eldridge, and Carles Badenes for
helpful discussions; Wolf-Rainer Hamann and Helge Todt for assistance with the
PoWR modeling; Rosalie McGurk and Luca Rizzi for support with the Keck/KCWI
observations; and the referee for helpful suggestions that improved the
quality of this paper. JC thanks Joseph Cassinelli whose care and guidance
sparked and inspired the initial drive to understand these massive star winds.
Based on observations with the NASA/ESA Hubble Space Telescope obtained at the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555. Support for Program number GO-15967 was provided through a grant
from the STScI under NASA contract NAS5-26555. This work was supported by a
NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute.
Data presented herein were obtained at the W. M. Keck Observatory from
telescope time allocated to the National Aeronautics and Space Administration
through the agency’s scientific partnership with the California Institute of
Technology and the University of California. The Observatory was made possible
by the generous financial support of the W. M. Keck Foundation. The authors
wish to recognize and acknowledge the very significant cultural role and
reverence that the summit of Maunakea has always had within the indigenous
Hawaiian community. We are fortunate to have the opportunity to conduct
observations from this mountain. OGT acknowledges generous support from a
Carnegie-Princeton Fellowship, through Princeton University and the Carnegie
Observatories, and from the Momental Foundation. This work was performed in
part at the Aspen Center for Physics, which is supported by National Science
Foundation grant PHY-2210452, and was partially supported by a grant from the
Simons Foundation. This research used NASA’s Astrophysics Data System,
adstex666https://github.com/yymao/adstex, and the arXiv preprint server. AACS
and VR acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) in the form of an Emmy Noether Research Group – Project-
ID 445674056 (SA4064/1-1, PI Sander). AACS further acknowledges support by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
138713538 – SFB 881 (“The Milky Way System”, subproject P04). This work was
partially supported by funding provided by the Federal Ministry of Education
and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of
the Excellence Strategy of the German Federal and State Governments. This
publication has benefited from discussions in a team meeting (PI: Oskinova)
sponsored by the International Space Science Institute (ISSI) at Bern,
Switzerland. Several figures in this work were created with WRplot, developed
by W.-R. Hamann.

##  Appendix A Parameters in the Initial PoWR Models

Table A1 reports the values of all PoWR model parameters that were adjusted in
our analysis that we adopted in the starting model for each star. The choice
of initial guesses was informed by previous analysis of these stars;
observations of their host galaxies; and theoretical wind-parameter scaling
relationships, as detailed in Section 3.2. Because we judged the goodness of
fit by eye and adjusted parameters by hand, the specific initial guesses have
little impact on the final result, but starting parameters as close as
possible to the best values are beneficial to reduce computation time. We
report our initial model parameters here to facilitate reproducibility of this
analysis.

Table A1: Initial guesses for parameters in the PoWR models. Parameter |  LP26 |  S3 |  A15
---|---|---|---
$T_{\mathrm{eff}}$ (kK) |  37 |  32 |  37
$\log(g/\mathrm{cm\,s}^{-2})$ |  4.0 |  3.5 |  4.0
$\log(L_{\star}/L_{\odot})$ |  5.1 |  5.1 |  5.3
$\log(\dot{M}_{\star}/M_{\odot}\,\mathrm{yr}^{-1})$ |  $-$8.7 |  $-$8.3 |  $-$8.1
$v_{\infty}$ (km s-1) |  1720 |  1500 |  2160
$v\,\mathrm{sin}(i)$ (km s-1) |  370 |  290 |  80
$\xi_{\mathrm{min}}$ (km s-1) |  10 |  10 |  10
$X_{\mathrm{He}}$ |  0.25 |  0.25 |  0.25
$\log(X_{\mathrm{C}})$ |  $-4.58$ |  $-4.29$ |  $-3.92$
$\log(X_{\mathrm{N}})$ |  $-5.17$ |  $-5.06$ |  $-4.64$
$\log(X_{\mathrm{O}})$ |  $-3.76$ |  $-3.46$ |  $-3.09$
$\log(X_{\mathrm{Fe}})$ |  $-4.44$ |  $-4.14$ |  $-3.62$
$E(B-V)$host (mag) |  0.003 |  0.048 |  0.018

## References

  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, Nature, 622, 707, doi: 10.1038/s41586-023-06521-7
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * 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
  * Aver et al. (2022) Aver, E., Berg, D. A., Hirschauer, A. S., et al. 2022, MNRAS, 510, 373, doi: 10.1093/mnras/stab3226
  * Baum et al. (1992) Baum, E., Hamann, W. R., Koesterke, L., & Wessolowski, U. 1992, A&A, 266, 402
  * Beaton et al. (2018) Beaton, R. L., Bono, G., Braga, V. F., et al. 2018, Space Sci. Rev., 214, 113, doi: 10.1007/s11214-018-0542-1
  * Becker et al. (2001) Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850, doi: 10.1086/324231
  * Berg et al. (2019) Berg, D. A., Erb, D. K., Henry, R. B. C., Skillman, E. D., & McQuinn, K. B. W. 2019, ApJ, 874, 93, doi: 10.3847/1538-4357/ab020a
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918, doi: 10.1093/mnras/staa2801
  * Björklund et al. (2021) Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36, doi: 10.1051/0004-6361/202038384
  * Bouret et al. (2003) Bouret, J. C., Lanz, T., Hillier, D. J., et al. 2003, ApJ, 595, 1182, doi: 10.1086/377368
  * Bouret et al. (2015) —. 2015, MNRAS, 449, 1545, doi: 10.1093/mnras/stv379
  * Bouret et al. (2013) Bouret, J. C., Lanz, T., Martins, F., et al. 2013, A&A, 555, A1, doi: 10.1051/0004-6361/201220798
  * Bresolin et al. (2006) Bresolin, F., Pietrzyński, G., Urbaneja, M. A., et al. 2006, ApJ, 648, 1007, doi: 10.1086/506200
  * Bresolin et al. (2007) Bresolin, F., Urbaneja, M. A., Gieren, W., Pietrzyński, G., & Kudritzki, R.-P. 2007, ApJ, 671, 2028, doi: 10.1086/522571
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115, doi: 10.1051/0004-6361/201016113
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Cameron et al. (2023) Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115, doi: 10.1051/0004-6361/202346107
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157, doi: 10.1086/153315
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Crowther et al. (2022) Crowther, P. A., Broos, P. S., Townsley, L. K., et al. 2022, MNRAS, 515, 4130, doi: 10.1093/mnras/stac1952
  * Crowther et al. (2016) Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624, doi: 10.1093/mnras/stw273
  * Curtis-Lake et al. (2023) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622, doi: 10.1038/s41550-023-01918-w
  * Dalcanton et al. (2009) Dalcanton, J. J., Williams, B. F., Seth, A. C., et al. 2009, ApJS, 183, 67, doi: 10.1088/0067-0049/183/1/67
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * de Mink et al. (2013) de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166, doi: 10.1088/0004-637X/764/2/166
  * Dolphin (2016) Dolphin, A. 2016, DOLPHOT: Stellar photometry, Astrophysics Source Code Library, record ascl:1608.013.  http://ascl.net/1608.013
  * Dolphin (2000) Dolphin, A. E. 2000, PASP, 112, 1383, doi: 10.1086/316630
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * Eldridge & Stanway (2022) Eldridge, J. J., & Stanway, E. R. 2022, ARA&A, 60, 455, doi: 10.1146/annurev-astro-052920-100646
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Endsley et al. (2023) Endsley, R., Stark, D. P., Whitler, L., et al. 2023, MNRAS, 524, 2312, doi: 10.1093/mnras/stad1919
  * Evans et al. (2007) Evans, C. J., Bresolin, F., Urbaneja, M. A., et al. 2007, ApJ, 659, 1198, doi: 10.1086/511382
  * Evans et al. (2019) Evans, C. J., Castro, N., Gonzalez, O. A., et al. 2019, A&A, 622, A129, doi: 10.1051/0004-6361/201834145
  * Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415, doi: 10.1146/annurev.astro.44.051905.092514
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293
  * Garcia & Herrero (2013) Garcia, M., & Herrero, A. 2013, A&A, 551, A74, doi: 10.1051/0004-6361/201219977
  * Garcia et al. (2019) Garcia, M., Herrero, A., Najarro, F., Camacho, I., & Lorenzo, M. 2019, MNRAS, 484, 422, doi: 10.1093/mnras/sty3503
  * Garcia et al. (2014) Garcia, M., Herrero, A., Najarro, F., Lennon, D. J., & Alejandro Urbaneja, M. 2014, ApJ, 788, 64, doi: 10.1088/0004-637X/788/1/64
  * Georgy et al. (2013) Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103, doi: 10.1051/0004-6361/201322178
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Götberg et al. (2017) Götberg, Y., de Mink, S. E., & Groh, J. H. 2017, A&A, 608, A11, doi: 10.1051/0004-6361/201730472
  * Gräfener et al. (2002) Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244, doi: 10.1051/0004-6361:20020269
  * Green et al. (2015) Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25, doi: 10.1088/0004-637X/810/1/25
  * Groh et al. (2019) Groh, J. H., Ekström, S., Georgy, C., et al. 2019, A&A, 627, A24, doi: 10.1051/0004-6361/201833720
  * Gull et al. (2022) Gull, M., Weisz, D. R., Senchyna, P., et al. 2022, ApJ, 941, 206, doi: 10.3847/1538-4357/aca295
  * Hainich et al. (2019) Hainich, R., Ramachandran, V., Shenar, T., et al. 2019, A&A, 621, A85, doi: 10.1051/0004-6361/201833787
  * Hamann & Gräfener (2003) Hamann, W.-R., & Gräfener, G. 2003, A&A, 410, 993, doi: 10.1051/0004-6361:20031308
  * Hamann & Koesterke (1998) Hamann, W. R., & Koesterke, L. 1998, A&A, 335, 1003
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hawcroft et al. (2023) Hawcroft, C., Sana, H., Mahy, L., et al. 2023, arXiv e-prints, arXiv:2303.12165, doi: 10.48550/arXiv.2303.12165
  * Heap et al. (2006) Heap, S. R., Lanz, T., & Hubeny, I. 2006, ApJ, 638, 409, doi: 10.1086/498635
  * Herrero et al. (1992) Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209
  * Higgins & Vink (2019) Higgins, E. R., & Vink, J. S. 2019, A&A, 622, A50, doi: 10.1051/0004-6361/201834123
  * Holgado et al. (2020) Holgado, G., Simón-Díaz, S., Haemmerlé, L., et al. 2020, A&A, 638, A157, doi: 10.1051/0004-6361/202037699
  * Huenemoerder et al. (2012) Huenemoerder, D. P., Oskinova, L. M., Ignace, R., et al. 2012, ApJ, 756, L34, doi: 10.1088/2041-8205/756/2/L34
  * Hunter et al. (2007) Hunter, I., Dufton, P. L., Smartt, S. J., et al. 2007, A&A, 466, 277, doi: 10.1051/0004-6361:20066148
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90