A Panchromatic Study of Massive Stars in the Extremely Metal-Poor Local Group Dwarf Galaxy Leo A

We characterize massive stars (M>8 M_sun) in the nearby (D~0.8 Mpc) extremely metal-poor (Z~5% Z_sun) galaxy Leo A using Hubble Space Telescope ultra-violet (UV), optical, and near-infrared (NIR) imaging along with Keck/LRIS and MMT/Binospec optical spectroscopy for 18 main sequence OB stars. We find that: (a) 12 of our 18 stars show emission lines, despite not being associated with an H II region, suggestive of stellar activity (e.g., mass loss, accretion, binary star interaction), which is consistent with previous predictions of enhanced activity at low metallicity; (b) 6 are Be stars, which are the first to be spectroscopically studied at such low metallicity -- these Be stars have unusual panchromatic SEDs; (c) for stars well-fit by the TLUSTY non-local thermodynamic equilibrium (non-LTE) models, the photometric and spectroscopic values of T_eff and log(g) agree to within ~0.01 dex and ~0.18 dex, respectively, indicating that NUV/optical/NIR imaging can be used to reliably characterize massive (M ~ 8-30 M_sun) main sequence star properties relative to optical spectroscopy; (d) the properties of the most massive stars in H II regions are consistent with constraints from previous nebular emission line studies; and (e) 13 stars with M>8 M_sun are>40 pc from a known star cluster or H II region. Our sample comprises ~50% of all known massive stars at Z<10% Z_sun with derived stellar parameters, high-quality optical spectra, and panchromatic photometry.


Introduction
Metal-poor (Z  10% Z e ) massive stars (M  8 M e ) are central to a wide variety of astrophysics. In low-mass galaxies, their feedback (i.e., stellar winds, supernovae) regulates star formation, enriches the interstellar medium, and may drive the formation of dark matter cores (e.g., Mac Low & Ferrara 1999;Stinson et al. 2007;Governato et al. 2010;El-Badry et al. 2016;Emerick et al. 2018;Eldridge & Stanway 2022). At low metallicity, those processes are expected to change significantly in comparison to higher-metallicity environments (e.g., decreased line-driven winds, properties of supernovae; Heger et al. 2003;Eldridge & Tout 2004;Young et al. 2008;Smartt 2009;Vink & Sander 2021) due to changes in stellar properties (e.g., rotation, opacity, internal mixing; Maeder et al. 2005;Meynet et al. 2008Meynet et al. , 2017Sanyal et al. 2017;Groh et al. 2019). For example, in the early universe, ionizing photons from early generations of massive metal-poor stars may power cosmic reionization, since they experience very high effective temperatures over a sizable fraction of their lives, resulting in an extreme ionizing ultraviolet (UV) radiation field (e.g., Wang 2013; Madau & Dickinson 2014;Stark et al. 2014;Robertson et al. 2015;Robertson 2021). At all redshifts, single and binary massive metal-poor stars are thought to be the progenitors to many explosive transients, including gravitational wave sources (e.g., Smartt 2009; Smith 2014; Spera et al. 2015).
Due to their close proximity, which enables resolved star studies, nearby star-forming dwarf galaxies are among the best targets for studies of individual massive stars and their environments at extremely low metallicities (e.g., Dalcanton et al. 2009;Hunter et al. 2012;Ott et al. 2012). However, despite substantial observational investments in resolving gas and stars in nearby dwarf galaxies, our detailed knowledge of their resolved massive star content has been limited. Many massive stars are hot (15,000 K) meaning their spectral energy distributions peak at UV wavelengths. Massive stars are not always the brightest, most easily identifiable stars in the optical (e.g., Massey 2003), which limits the utility of the existing optical imaging. Older Hubble Space Telescope (HST) UV imaging (e.g., with WFPC2) exists for some local dwarf galaxies, but the data quality is modest relative to the optical, and the resulting studies are generally limited to the hottest and/or brightest objects (e.g., Bianchi et al. 2012). Other UV imaging surveys focus on the ensemble population of starforming regions and young clusters as hosts of massive stars, but are less focused on the properties of the individual stars (e.g., LEGUS; Calzetti et al. 2015).
There is a similar paucity of spectra for massive metal-poor stars. While high-resolution optical and UV spectra of massive stars has improved our understanding of massive stars in the metal-rich environments (e.g., Bianchi & Garcia 2002;Bouret et al. 2003;Heap et al. 2006;Sota et al. 2011Sota et al. , 2014Barbá et al. 2017;Russeil et al. 2017;Li 2021), data of similar quality is less abundant at lower metallicities. Much of our understanding of metal-poor massive stars is based primarily on studies of the Small Magellanic Cloud (SMC, Z ∼1/5Z e ; Martins et al. 2004;Castro et al. 2018;Ramachandran et al. 2019). However, the metallicity of the SMC does not probe the extremely metal-poor regime. One example of the need for extremely low-metallicity spectra is the need to calibrate the relative strengths of He I absorption lines that underlie He I emission lines in nebular spectra. This is an important correction for the determination of the primordial helium abundance and currently relies entirely on simulated spectra (e.g., Aver et al. 2021). Existing stellar models for lower metallicities predict significant differences in evolutionary pathways with implications for time-integrated feedback, final fates, and compact remnants (e.g., Meynet & Maeder 2002;Eldridge et al. 2008;Chen et al. 2015;Costa et al. 2021). Lastly, spectral atmospheric models for massive metal-poor stars are currently calibrated from higher-metallicity environments. However, at lower metallicities, models predict, for example, more emission features due to increased stellar activity and deeper absorption features due to the increased effective temperature for a given spectral type, or broader absorption features at a given spectral type (e.g., Kubátová et al. 2019;Martins & Palacios 2021).
Historically, the study of metal-poor massive stars has been challenging, since most star-forming metal-poor dwarf galaxies are located at distances of 1 Mpc, meaning that even the most luminous massive stars are often faint and crowded. The limiting factors above have typically confined the study of lowmetallicity massive stars to a small number of stars in the closest dwarf galaxies. While providing important insights into chemical evolution and massive stellar physics, the challenge of acquiring these data limited the sample sizes and homogeneity required for systematic progress in understanding metal-poor massive star evolution and how these stars impact their surroundings. Several recent efforts have sought to improve this situation. For example, only recently has there been a large survey published of ∼150 low-resolution spectra (R ∼ 1000, signal-to-noise ratio, S/N, of ∼8-256 with a majority S/N ∼ 30-100; Lorenzo et al. 2022). The study currently only provides spectra and spectral type identifications, but no stellar parameters. This is currently the only study of its kind at such low metallicities. Other than this low resolution, optical spectroscopy is only available for a handful of massive stars in several sub-SMC dwarf galaxies within ∼1.5 Mpc (e.g., Venn et al. 2004;Bresolin et al. 2007;Kudritzki et al. 2008;Hosek et al. 2014;Tramper et al. 2014;Berger et al. 2018;Evans et al. 2019;Garcia et al. 2019b;Telford et al. 2021b). In these same systems, there are also concerted efforts to obtain UV spectroscopy of O-type stars with HST (e.g., García et al. 2015;Telford et al. 2021a;Wofford et al. 2021), including a dozen stars as part of the Hubble UV Legacy Library of Young Stars as Essential Standards program (e.g., Roman-Duval et al. 2020.
Motivated by the need for a panchromatic photometric and high-quality spectroscopic inventory, we use Leo A as a testbed galaxy. We pursue the following aims: (a) to characterize the population of massive stars in Leo A and (b) test the reliability of photometric-only techniques in this metallicity regime. The latter is especially important as there are few extremely metal-poor systems in which individual massive stars can be resolved and high-quality spectra can be acquired. Leo A is an ideal test bed as it is the closest of the star-forming extremely metal-poor dwarf galaxies known in the local universe (D ∼ 760 kpc, Z ∼ 5% Z e ; e.g., van Zee et al. 2006;Nagarajan et al. 2022). Despite its favorable distance, low metallicity, and sizable massive star population (e.g., Cole et al. 2007), there has not yet been a systematic investigation of massive metal-poor star properties in Leo A. Existing studies that focus on Leo A's young stars do so in the context of H II region physics or as dynamical tracers of dark matter (e.g., Skillman et al. 1989;van Zee et al. 2006;Brown et al. 2007;Cole et al. 2007;Ruiz-Escobedo et al. 2018;Leščinskaitė et al. 2022;Ricotti et al. 2022;. Recently, Leo A was observed as part of the HST Local Ultra Violet Infrared Treasury program (LUVIT; GO-15275, PI Gilbert and GO-16162, PI Boyer). This survey has obtained new HST UV, optical, and near-IR (NIR) imaging with HST for 19 nearby dwarf galaxies in GO-16162 and 21 nearby dwarf galaxies in GO-15275 (the overlap between the two programs consists of 19 galaxies). Combined with archival HST optical imaging, it provides broad information about the spectral energy distribution (SED) for a variety metal-poor stars across the H-R diagram.
Using LUVIT imaging for target selection, we acquired optical spectroscopy with the Low-Resolution Imaging Spectrograph (LRIS) on Keck and Binospec on MMT for two dozen of the most UV luminous stars in Leo A. This is one of the largest uniform collections of high-quality optical spectroscopy for extremely low-metallicity massive stars. Between the panchromatic HST imaging and the high-fidelity optical spectra, we undertake here a systematic study of hot, massive, metal-poor stars in Leo A and their relationship to the broader galactic environment. This paper is organized as follows. In Section 2, we describe the multiwavelength HST observations and data reduction, along with target selection, observations, and data reduction for the Keck/LRIS and MMT/Binospec optical spectroscopy. In Section 3, we detail the process of measuring stellar properties from the HST-based SEDs. In Section 4, we expand on how we measure the stellar properties from modeling the optical spectra, and from the conventional approach of measuring equivalent widths (EWs). We present the results of each method in Section 5 and present a detailed discussion of the spectrum for each star along with a comparison of the photometric and spectroscopic parameters we measure. In Section 6, we discuss our results in the context of nebular studies of H II regions, field OB-stars, Be stars, and stars in the extremely metal-poor environments, as well as a comparison of SED fitting to full spectral fitting.

HST Observations and Photometry
We make use of HST UV, optical, and NIR broadband imaging for selecting the spectroscopic sample and for measuring stellar properties based on HST photometry. Deep optical imaging was acquired in the Advanced Camera for Surveys (ACS)/ Wide Field Channel (WFC) F475W and F814W filters as part of GO-10590 (Cole et al. 2007). WFC3/ UVIS F275W and F336W and WFC3/IR F110W and F160W imaging was collected as part of LUVIT. We note that the ACS/WFC imaging has the largest spatial coverage of Leo A, followed by WFC3/UVIS imaging, while WFC3/IR imaging covers the smallest area. We only focus on stars for which we have at least 4-band UV+optical imaging, and include NIR when possible.
The LUVIT data are reduced by utilizing the photometric reduction pipelines that have been developed and improved for HST simultaneous N-band photometry . We refer the readers to the LUVIT survey paper (K. M. Gilbert 2022, in preparation) for the details of data collection and reduction for the entire LUVIT sample. Here we provide some of the Leo A specific details. We first aligned all flc (ACS/WFC, WFC3/UVIS) and flt (WFC3/IR) images to the Gaia DR2 astrometric solution following Bajaj (2017) 19 and combined images in each band weighted by exposure time using AstroDrizzle (Hack et al. 2020) after flagging and masking bad pixels and cosmic rays on input exposures. These aligned images, as well as the individual exposures with updated data quality extensions, are processed via the DOL-PHOT package (Dolphin 2000(Dolphin , 2016 for photometry. We use the combined full-depth F475W image as the astrometric reference image for our final photometry. DOLPHOT reports measurements of point-spread function photometry for every object detected on any exposures within the reference image footprint. These include the flux, flux uncertainty, magnitude in the Vega system, S/N, and various photometry quality metrics. A raw output catalog of all objects identified by DOLPHOT is then subjected to quality cuts in each filter to include objects: (1) with S/N 4; (2) with crowding <1.3 and <2.25 for UV and optical/IR filters, respectively; and (3) with square of sharpness (SHARP 2 ) <0.15 and <0.2 for UV/IR and optical filters, respectively, producing the "GST" catalog. These quality cuts were applied in each band independently. To remove objects that are associated with diffraction spikes after performing the GST cuts, we first identified individual tiny regions with severe diffraction spike features using the median filter technique and then culled objects belonging to those regions (K. M. Gilbert 2022, in preparation). To construct an initial spectroscopic target sample, we further required the S/N in F475W and F814W to be >15 and F475W SHARP 2 + F814W SHARP 2 to be 0.3. Additional criteria for the final spectroscopic target selection are presented in Section 2.2. Even though the photometric depth reaches a couple of magnitudes fainter than the old main-sequence turnoff, we limited our stellar SED modeling to the bright targets of our interest. We note that we did not use the NIR medium-band photometry (F127M, F139M, and F153M) in this study, since it does not provide any further significant constraints in this stellar regime.
Photometric accuracy, precision, and completeness were computed from artificial star tests (ASTs). We generated ∼135,000 dust-extinguished SED models in our UV, optical, and NIR bands that cover the entire H-R diagram, but more heavily populate the observed color and magnitude ranges (Choi et al. 2020). We assigned more input ASTs in higher stellar density regions of Leo A to better capture nonlinear trends in photometric characteristics as a function of source crowding. Photometric measurements for input artificial stars were made in the exact same way as we did for real objects, and then identical quality cuts were applied to determine if the injected artificial stars were recovered. Comparisons of the input and output magnitudes are used to quantify our photometry uncertainty. The 50% completeness limits for F275W, F336W, F475W, F814W, F110W, and F160W are 24.97, 25.89, 29.10, 27.89, 26.58, and 25.27, respectively. Figure 1 shows select UV, optical, and NIR color-magnitude diagrams (CMDs) of Leo A, with the spectroscopic sample highlighted, as described in Section 2.2. The NIR CMD contains a smaller sample because the NIR imaging field of view covers a smaller spatial area of Leo A. The CMDs are close to 50% or above of the completeness limit. We overplot the rotating (v/vcrit = 0.4) MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016;Dotter 2016) for reference. We shift the tracks to the Gaia eDR3-anchored distance modulus of μ = 24.40 mag reported by Nagarajan et al. (2022) and A v = 0.03 mag, consistent with the foreground maps of Schlafly & Finkbeiner (2011). We do not apply any further dust corrections, since Leo A shows little internal redding (Skillman et al. 1989;van Zee et al. 2006).

Spectroscopic Target Selection
We observed Leo A for one night using Keck/LRIS in the 2020A semester. To fill the slit mask, we prioritized the stars brightest in the HST F336W − F475W CMD (Figure 1). We were able to place nine UV bright stars on the LRIS slit mask. This number was limited by crowding, LRIS slit mask alignment constraints, and science goals. The remaining slits were filled with filler stars (i.e., bright stars in ground-based photometric observations of Leo A). To achieve a sufficient S/N (see below), we integrated for a full night on this single slit mask.
We also observed Leo A on MMT/Binospec in the 2020A semester. Targets were selected using the Binospec mask software with following criteria as input: F475W − F814W − 0.1, F275W 20, and −0.8 F275W − F336W − 0.2, thereby selecting the bluest and brightest sources.
The two spectroscopic data sets have stars in common, which allows us to check for systematic changes and gauge consistency. Figure 1 shows the LRIS and Binospec spectroscopic samples overplotted on the HST CMDs. The majority of the brightest UV stars have spectroscopic coverage from at least one of the telescopes. The handful that do not were excluded due to slit mask and/or crowding constraints. We plan to target them in a future observing run. Figure 2 shows zoomed-in CMDs that highlight the spectroscopic sample. We label them based on the origin of the spectra used for analysis (i.e., Keck stars: K, Binospec stars: B) and in decreasing F475W magnitude. In the UV-optical CMD (left panel), we see that the stars appear to have masses of at least 8 M e , since they pool around the 10 M e track. In the optical CMD (right panel), the picture is not as clear as in the UV-optical CMD. Interestingly, a handful of stars (K4, K5, K6, K8, B5, and B8) visually appear to be main-sequence (MS) stars in the UV-optical CMD, but appear on the blue core helium burning (BHeB) sequence in the optical CMD. We discuss these stars in detail in Section 5.3. In total, we have spectroscopic observations for 18 stars. Figure 3 shows the spatial position of the spectroscopic sample overplotted on the HST ACS/F475W image.Observational properties of these stars are listed in Table 1. We now describe the details of each spectroscopic data set.

LRIS and Binospec Observation
We observed Leo A on the night of 2020 February 20 using LRIS (Oke et al. 1995) at the Keck telescope at W.M. Keck Observatory in Hawaii, and 2020 February 23, 24, and 29 using the Binospec instrument, mounted on the 6.5 m MMT telescope. For LRIS, we obtained spectra using a custom made 0 7 slit-mask. The gratings for the LRIS red channel were set to 1200/7500, covering 7224−8855 Å. The grism for the LRIS blue channel was set to 600/4000, covering 3330−5910 Å, and the dichroic was set to D680, yielding a nominal resolution of R ∼ 1800. For Binospec, we observed two nights using the 1000 lpmm grating, providing spectra from 3900-5400 Å at a nominal resolution of R ∼ 3600. The last night, we used the 270 lpmm grating and obtained spectra spanning 3900-9240 Å, with a resolution of R ∼ 1250.
For LRIS, we integrated for a total of 6.3 hr using multiple 1200-1800 s exposures. We reduced the spectra using the Pypelt pipeline (Prochaska et al. 2019) and then coadded each single epoch spectrum. The resulting S/Ns range from ∼30 to ∼83 per Å at ∼5000 Å and ∼5 to ∼17 per Å at ∼8000 Å. For Binospec, we integrated for a total of 2 hr on night 1, 1 hr on night 2, and 1 hr on night 3, in intervals of 1200 s exposures. We observed the same stars all three nights. We reduced the spectra with the Binospec Data Reduction Pipeline (Kansky et al. 2019). We coadded single-exposure spectra for each star from the first two nights. The coadded spectra from the first two nights have S/Ns ranging from ∼9 to There are fewer stars on the NIR CMD because WFC3/IR has a smaller field of view and does not entirely overlap the UV and optical imaging. All spectroscopic stars are massive stars (i.e., M > 8M e ). The NUV CMD shows that the stars in our spectroscopic sample are not very evolved. Figure 2. This plot is identical to Figure 1, only zoomed in on luminous blue stars selected for spectroscopy in the UV and optical CMDs. We label the stars by the names in Table 1, where K stands for stars with Keck and B for stars with Binospec spectra. For reference, we overplot select MIST tracks Interestingly, a handful of stars that appear on the main sequence in the UV-optical CMDs appear in the blue core helium burning region of the F475W-F814W CMD. We discuss this in Section 5.3.
∼31 per Å at ∼5000 Å. Spectra from the third night have S/Ns ranging from ∼11 to ∼40 per Å. Table 1 lists the basic properties of our LRIS and Binospec targets. As noted in Table 1, four stars were observed with both Binospec and LRIS: K1/B1, K2/B2, K3/B3, and K7/B6. For simplicity, we refer to the Keck stars as K1-K8. We note that a second star fell into the same slit K2; however, K2 is unaffected by this. This was the only slit to do so, we do not include the second star in the sample.

Spectroscopic Normalization
We normalize all spectroscopic observations using PySpecKit (Ginsburg & Mirocha 2011). We exclude the parts of the spectrum showing strong absorption or emission features. We normalize all of our spectra using a fifth-order polynomial, in order to avoid over fitting. We show examples of normalized sample spectra with a detailed indication of spectral features in Figure 4.

Measuring Stellar Parameters from Panchromatic Photometry
We use the Bayesian Extinction and Stellar Tool (BEAST; Gordon et al. 2016) to infer physical properties of luminous blue stars in Leo A from the HST broadband multiwavelength imaging (F275W, F336W, F475W, F814W, F110W, and F160W). The BEAST is a probabilistic approach to simultaneously modeling the stellar and line-of-sight dust properties of individual stars in nearby galaxies by carefully considering observational uncertainties arising from photon noise and crowding. Here, we briefly summarize how the BEAST works and its application to stars in Leo A. Detailed descriptions of the BEAST technical underpinnings and its application to a nearby low-metallicity dwarf galaxy can be found in Gordon et al. (2016) and Choi et al. (2020), respectively.
Briefly, the BEAST creates a physics model grid by mapping a stellar evolution library onto a stellar atmospheric model using stellar evolutionary parameters (initial mass M ini , age A, and metallicity Z) and stellar atmospheric parameters (effective temperature T log eff ( ), surface gravity g log( ), and metallicity Z). Lower-temperature stars (i.e., A type and cooler) are paired with local thermodynamic equilibrium models from Castelli & Kurucz (2004), while hotter stars are paired with nonlocal thermodynamic equilibrium (non-LTE) atmospheres from the TLUSTY OSTAR and BSTAR grids (Lanz & Hubeny 2003. This process covers intrinsic properties for each star. Stars are then moved to the distance of the galaxy of interest (or can be left as a free parameter) and then the intrinsic light is extinguished by a dust model. The dust model used by the BEAST is a flexible mixture model that allows for smooth variation between average Milky Way-like and SMC-like attenuation curves, including the curve steepness and changes in the 2175 Å bump strength (Gordon et al. 2022). For each of the resultant dust-extinguished stellar models, an observational noise model is determined based on ASTs.
For our analysis of stars in Leo A, we adopt the following priors on various physical parameters. First, we fix the distance modulus to μ = 24.40 mag, which is the Gaia eDR3-anchored RR Lyrae distance reported by Nagarajan et al. (2022). Second, we use the MIST stellar evolutionary models with rotation (v/ vcrit = 0.4) in preference over the MIST nonrotation (v/ vcrit = 0.0) models and the PARSEC Bressan et al. (2012) models. This choice is primarily made due to the evidence of rotation in the spectra, which is accounted for in the MIST Figure 3. The Leo A HST ACS/F475W image zoomed in on the central region of the galaxy. Stars with spectroscopic observations are color-coded by telescope: Keck/LRIS (blue triangles), MMT/Binospec (red circles), and both (magenta squares). Contours of Hα emission from ground-based imaging (Local Volume Legacy Survey; Kennicutt et al. 2007) are overplotted in green. Note that many of the massive stars are not obviously associated with cataloged H II regions including one O-star (K7).

Table 1
Observing Details of the Spectroscopic Data Taken with LRIS (KECK) and Binospec (MMT)

Star
Alternative Note. The columns are: (1) The name of each star used throughout the paper; (2) alternative ID; (3) R.A.; (4) decl.; (5) observing date; (6) slit width; (7) integration time; (8)-(10) select HST magnitudes and colors; (11) and (12) S/N at specified wavelength. A-Same star, with different name in the observing logs. K-Keck, B-Binospec. rotation model. In Appendix A, we discuss the difference between the two MIST models as well as the PARSEC models. Third, we adopt an SMC-like average extinction curve with R V = 2.74 (Gordon et al. 2003), and allow A V to vary from a minimum of A V = 0.00 mag to a maximum of A V = 0.5 mag, which is the amount of internal dust needed to match the width of the upper MS in nearby dwarf galaxies based on fitting optical CMDs (e.g., Dolphin et al. 2003). We note that some stars can have higher A V values than the typical value set by the width of the MS. Finally, we adopt a metallicity of Z = 0.1 Z e throughout the paper. This is slightly higher than the H II region metallicity of Leo A, but it is: (a) the minimum value for both of the publicly available TLUSTY OSTAR and BSTAR highresolution spectral grids, and (b) may be more appropriate for the massive stars, which have been found to be more metal-rich (from iron absorption lines) than gas-phase metallicities (from oxygen emission lines) in other nearby dwarf galaxies (Garcia et al. 2014;Bouret et al. 2015;Telford et al. 2021b). Adopting this single metallicity ensures uniform analysis of the spectra and photometry. We report the median value (_ p50) and the uncertainty computed from the 68% confidence interval (i.e., 0.5 * (_ p84-_ p16)) of a marginalized 1D posterior probability distribution function (PDF) for each parameter. We show example BEAST fits in Figure 5 and discuss them in Section 5.1.

Measuring Stellar Parameters from Optical Spectroscopy
We use two methods to characterize the stars from spectroscopy. The first, which also derives stellar parameters (e.g., T log eff ( ), g log( ), v i sin , A V , and mass), is through an adaptation of the full spectral fitting code, "The Payne" (Ting et al. 2019) for massive stars. The Payne has primarily been used for cool stars in the Galaxy, but its flexibility makes it readily adaptable to different types of stellar spectra, including massive stars (e.g., Xiang et al. 2022). The second is through the measurement of EWs. EWs are conventionally used in the field of metal-poor massive stars (e.g., Bresolin et al. 2006;Evans et al. 2007;Tramper et al. 2011Tramper et al. , 2014Camacho et al. 2016;Evans et al. 2019;Garcia et al. 2019b), and provide a sanity check on our quantitative fitting. The extreme metal-poor nature of Leo A means there is a lack of empirical low-metallicity templates for comparison, while the low-resolution nature of the spectroscopy means that the EWs cannot always be as cleanly determined compared to higher- resolution spectroscopy. We use the EWs to provide guidance rather than rigorous parameter determination.

The Payne
For quantitative spectral fitting, we have adapted The Payne (Ting et al. 2019) for use with the TLUSTY non-LTE massive star model atmospheres (Lanz & Hubeny 2003). Specifically, we use the public OSTAR and BSTAR grids available on the TLUSTY website. 20 At its core, The Payne is a highly efficient, multidimensional interpolator that can generate accurate models based on sparse precomputed grid points by using a neural net (see Ting et al. 2019 for details).
We train two different neural nets that separately cover O-stars (Lanz & Hubeny 2003) and B-stars (Lanz & Hubeny 2007). Spectral features in these stars are sufficiently different, and the public TLUSTY grids are sufficiently sparse, that training nets based on separate models lead to more accurate synthetic spectra when compared to trying to interpolate across both model grids at once.
The two precomputed TLUSTY model grids are not uniform in their coverage of parameter space. For consistency, we adopt a fixed metallicity of Z/Z e = 0.1 in both grids, similarly to the BEAST models chosen. Lastly, we note that a solar helium abundance is assumed and by number is fixed to He/H = 0.1. We list all other parameters in Table 2.
Beyond T log eff ( ) and g log( ), we also fit for radial velocity, V R , and rotational velocity, v i sin , of each star. For v i sin , we use the Fourier method as outlined by Simón-Díaz & Herrero (2007); however, it is currently not a label we train the neural net on. Furthermore, we note that the accuracy of v i sin is limited by the resolution of the spectra. The resolution limit on v i sin is ∼170 km s −1 for Keck and ∼80 km s −1 for Binospec. For stars where we have both Binospec and Keck spectra, we report the Binospec derived v i sin value. We fix the macroturblent velocity to V macro = 20 km s −1 , which is motivated by various literature studies of hot, massive stars in the SMC (e.g., Mokiem et al. 2006;Penny & Gies 2009;Grassitelli et al. 2016). We further fix the instrumental resolution to R ∼ 1800 for the Keck spectra and R ∼ 3600 for the Binospec data. We note that having fixed resolution, V macro = 20 km s −1 and V micro = 20 km s −1 , all might slightly bias g log( ) and v i sin , since they all influence the line profile.
We use The Payne to find best-fit stellar parameters for our spectra via χ 2 minimization. We only fit LRIS and Binospec spectra blueward of 7500 Å, the red edge of the TLUSTY grid, and mask emission lines, which are not included in TLUSTY, as described in Section 5.3.  Note. The columns are: (1) grid name; (2) number of models in grid used to train the neural net; (3) temperature range in log-space; (4) temperature range in kilokelvin; (5) surface gravity range; and (6) fixed microturbulence value. We use the optimization results to initialize a Markov Chain Monte Carlo (MCMC) sampler that explores the posterior probability distribution for each star. We use a Gaussian loglikelihood function and adopt uniform priors for T log eff ( ), g log( ), v i sin , and V R . The ranges of allowed values for each are the grid edges for T log eff ( ) and g log( ) of the respective model grid and [0, 1000] km s −1 for v i sin and [−50, 150] km s −1 for V R .
We use the MCMC sampler emcee (Foreman-Mackey et al. 2013) to sample the posterior. We first use 512 walkers to first burn in our sample for 200 steps. We then initialize 512 walkers, applying a Gaussian scatter of 0.001 and set a maximum step number of 10 5 steps. However, the sampling stops if one of our convergence criteria is reached: (1) the autocorrelation time, τ, has changed by <1% over the last 100 steps, or (2) the sampler has run for >20τ steps. To remove any residual effects of burn-in or corrected samples, we discard the first 5τ samples from each walker and thin each chain samples by a factor of ∼τ/2.
For each star, we report the median value of the marginalized posterior PDFs along with the 16th and 84th percentiles as a measurement of uncertainty. If the uncertainty is less than 0.01, we round it to 0.01.
Once we have the parameters from the full spectral fitting, we use those to determine the luminosity (log (L)), the evolutionary mass (M evo ), the initial mass (M 0 ), and the spectroscopic mass (M spec ). We note that both M evo and M spec are present-day masses. The luminosity is determined by using both the spectroscopic parameters as well as the photometric data. Based on the spectroscopic parameters, we choose the appropriate model to scale in order to fit to the observed and dust-corrected photometric SED from the HST photometry. We use the ratio between the observed dust-corrected SED, f v , to the intrinsic SED (i.e., the model: F v ) and the distance to Leo A to obtain the spectroscopic radius R sing the spectroscopic radius, T log eff ( ), and the Stefan-Boltzmann constant (σ SB ), we calculate the luminosity s =´Ẃ e use the spectroscopic radius and g log( ) value to determine M spec : Finally, we use the derived log (L) and spectroscopically T log eff ( ) and the rotating MIST stellar models to find the mass track that the star falls on to determine M 0 . The models are spaced at 1 M e , and we determine the mass by visual inspection. For M evo we use the stellar parameters to find the stellar mass of the MIST stellar models by looking at the mass at the current position on the track.

Equivalent Widths
To support our qualitative analysis, we measure the EWs of each spectral line from the 1D spectrum for the stars in our sample. After normalizing the spectrum as described in Section 2.3, we combine the specutils (an extensible spectroscopic analysis toolkit for astronomy; Earl et al. 2019) package with scipy (Virtanen et al. 2020) to determine the EWs. We use the CLOUDY line list (Ferland et al. 2017) to identify the correct wavelengths. We specify the region of the spectral lines of interest and provide initial guesses for the line center, amplitude, and width. We choose to fit a Gaussian profile to all of our lines. In the scenario of a mixed absorption-emission feature, we employ a Gaussian fit to the emission component and subtract it out of the spectra before fitting the absorption feature for the EW. The results are used in Section 5.3 to support the commentary and in Table 3.

Analysis of Emission Line Stars
For a subset of stars (K4, K5, K6, K8, B5, B8), both the spectroscopic method and photometric method described above yield high reduced χ 2 (150) values. These stars also have unusual SEDs such that they are bright in the UV, red-optical, and NIR (when available), i.e., their SEDs are not well approximated by a blackbody as expected for normal stars. Their spectra usually feature strong emission, despite not being in close proximity to an H II region (except B5). They further show no or little extended emission beyond the stellar trace in the 2D spectra. We posit that these stars may have circumstellar disks (e.g., Iqbal & Keller 2013;Lazzarini et al. 2021) causing both the emission features seen in the spectra as well as the IR excess in the photometry leading to unusual SEDs. We use an alternative fitting technique to constrain the properties of the star and disk.
We fit emission lines stars in two different ways. First, we use full spectral fitting, which enables the recovery of v i sin and T log eff ( ). Second, we fit the photometric SEDs, which provides T log eff ( ) and A V . For spectroscopic fitting, we use The Payne but emphasize the He I/He II features by inflating uncertainties on other features, such that they effectively only provide continuum constraints. For Be stars, we mask the Balmer lines owing to strong emission. We recover T log eff ( ) and v i sin but not g log( ), which is derived from the Balmer lines.
Photometrically, we opt to fit a dual model of a star and potential disk to the UV-optical-NIR SEDs. We use pystellibs 21 to interpolate between model SEDs, and pyphot 22 to calculate synthetic photometry. Only stellar models are currently publicly available, hence we opt to use stellar models for the disk as well, since the disk is a nuisance parameter. This will still allow us to capture the essential underlying parameters, temperature and extinction, although we might miss the broader physics (e.g., density, velocity, and chemical composition). For the stellar component, we continue to use the TLUSTY models, similarly to the BEAST, while for the potential disk component we use the BaSeL (Lejeune 2002) library, since we need cooler temperatures than TLUSTY can provide. We fit the following parameters: . We adopt the same extinction curve as used in the BEAST fitting. We keep distance as well as mass constant (M star = 10M e for K4, K5, K6, and K8, and M star = 8M e for B5 and B8). We use the MCMC sampler emcee to determine the median value. We first use 64 walkers to first burn in our sample. We burn in the sample with 50 steps and use 128 walkers; we then set a maximum step number of 1000 steps. We determine convergence by the same criteria as 21 https://github.com/mfouesneau/pystellibs 22 https://mfouesneau.github.io/pyphot/ K1 0.16 ± 0.10 0.13 ± 0.10 0.43 ± 0.10 1.08 ± 0.10 0.61 ± 0.10 0.56 ± 0.10 0.18 ± 0.20 0.44 ± 0.20 0.71 ± 0.10 The Payne fitting. We use a log-Gaussian likelihood function and use a Gaussian prior for each parameter. We allow a broad range of priors so that we can simultaneously check whether the SED could also be explained by a companion star. We allow the parameters to sample over following ranges We note that we use the same extinction value for both components, since we are only interested in the stellar component. Lastly, we initialize the sampler by using the spectroscopic values. The results are used in Section 5.3.3 to support the commentary.

Stellar Parameters from Photometry
For comparison with our spectroscopic sample, we only ran the BEAST on stars with F475W < 22 and F475W − F814W < 0.5 (i.e., M  5M e on the MS). For select science cases (e.g., Section 5.1), we used the BEAST on slightly fainter stars.
To facilitate discussion and comparison, we divide the resulting photometric SED fits into three groups based on the χ 2 values. Fits with χ 2 < 20 are considered good fits, 20 χ 2 < 80 are decent fits, and χ 2 80 are poor fits. These χ 2 divisions are based on qualitative assessment of the fits, following other BEAST papers in the literature (e.g., Van De Putte et al. 2020). They are only meant to be used qualitatively to facilitate discussion. Furthermore we note that the uncertainties are likely underestimated, since they take the models at face value. Figure 5 shows three example SEDs in each of our χ 2 categories. These examples all have spectroscopic data enabling further evaluation. The top panel illustrates the case of a good SED fit for star K1 (χ 2 ∼ 18 for the stellar, dust, and crowding bias fit). All SED points for the median model (combined star, dust, and bias model) are consistent with the observed SED within 4σ. The derived physical parameters ( T log eff ( ), g log( ), A V , initial mass M ini , present-day mass M act , and L log ; ( ) see Section 5.3.1) all appear reasonable given the star's location on the observed CMDs. Overall, they indicate that K1 is a late-type O dwarf. We compare the SED and spectroscopically determined parameters in Section 5.3.
The middle panel of Figure 5 illustrates a decent SED fit using star B10 (χ 2 ∼ 72). The residuals show that the median model SED is within ∼4σ for five bands and ∼12σ of the observed SED for one NIR band. Like the good χ 2 for K1, the stellar parameters for B10 are in reasonable agreement with its position on the CMD and a B0.5 dwarf. The star, however, appears to have quite high extinction (A V = 0.49), which is puzzling. We compare the SED and spectroscopically determined parameters in Section 5.3.
The bottom panel of Figure 5 displays an example of a poor SED fit (χ 2 ∼ 1259) for star K6. In general, the median model values range in agreement from ∼5σ for the two UV SED points to ∼15σ for the F475W SED point. K5 falls into the category of stars that show an unusual SED, bright in the UV and red-optical, possibly due to the presence of a disk (see Section 5.3.3). Hence, it is not surprising that a single star model cannot be well fit to the SED.
Having illustrated our SED fitting with a few examples, we now consider general trends in the SED fits of luminous stars in Leo A. Figure 6 shows observed CMDs in Leo A color-coded by the results of our photometric SED fitting for select stellar parameters ( T log eff ( ), g log( ), M ini , A V ) along with the χ 2 value of the fit. The values assigned to these points reflect the median of the posterior distributions reported by the BEAST.
In general, as expected, the temperature decreases for redder CMD colors (second row of Figure 6). This is particularly evident in the UV-optical CMDs, where the hottest stars are located at the blue edge of the MS. We find that T log eff ( ) only weakly depends on χ 2 ; poor fits seem to have a qualitatively reasonable correspondence between UV-optical color and T log eff ( )(i.e., the redder the stars, the cooler the stars). However, the situation changes slightly when considering optical-only colors. Here, while the majority still follow the expected trend, we find several outliers. Namely, several decently fit stars that appear on the MS in the UV-optical CMDs and have hot SED-based values of T log eff ( ) appear to be cooler stars (i.e., BHeB region stars) in the optical only CMD. In the NIR CMD we see, similarly to the optical-only CMD, that the outliers are primarily the same stars that are displaying odd temperature to color combinations.
For a subset of these color outliers, we have spectra that we use to cross-check the reasonability of the photometrically determined parameters. In many cases, these color outliers appear to have emission features that suggest they are not nebular in origin and instead may indicate these are Be stars. We compare spectroscopic and photometric parameters for this subset in Section 5.3 and discuss the candidate Be stars in Section 6.3. Similar to T log eff ( ), trends in g log( ) (third row of Figure 6) are generally consistent with expectations. The bluest stars should have the highest values of g log( ), which is what we find. For redder stars, we expect and observe a decrease in g log( ). We find that several of the stars causing discrepancies in T log eff ( ) in the optical-only CMD show questionable values for g log( ) as well. In the UV CMD, we observe the same outliers as in the optical CMD. In the NIR CMD, we similarly observe that the outliers in T log eff ( ) are outliers in g log( ). We observe a broad trend between M ini and g log( ), expected from stellar models, with a few exceptions. The most-massive stars are generally the brightest and hottest, while fainter stars are less massive. Stars with poor fits in g log( ) are exceptions. Their values of M ini are higher than one expects from their location on the panchromatic CMDs. Stars in the BHeB region of the optical CMD have M ini  13 M e , but are all poor fits, suggesting we should use caution in taking their masses at face value.
Several studies of Leo A find it has very little dust (e.g., van Zee et al. 2006;Ruiz-Escobedo et al. 2018). This is reflected in the BEAST fitting as well. The MS stars all show little to no dust (A V 0.08), though the BHeB sequence shows higher A V values. This difference can potentially be explained by elevated dust production in supergiants (e.g., Waters 2010; de Wit et al. 2014), but we note that several of the high A V stars also have poor fits, as is especially visible in the NIR CMD. Several of the previously mentioned outliers (i.e., the stars with peculiar colors) may have higher extinction, since the BEAST is trying to fit the excess IR emission created by the disk as we discuss in Section 5.3. Another highly probable explanation for the poor fits of stars appearing on the BHeB sequence in the optical CMD with unphysically high dust is the lack of calibration of stellar models for BHeB stars at the extreme low metallicity (e.g., McQuinn et al. 2012). Figure 6. A panchromatic gallery of CMDs for Leo A ranging from UV-optical on the left to NIR on the right. Points in each row are color-coded by results from SED fitting with the BEAST. From top to bottom, the rows reflect: χ 2 , T log eff ( ), g log( ), M ini , and A V . For χ 2 , we divide the fits into three categories: good (pink), decent (purple), and poor (cyan). Stars that have spectroscopic data are plotted as squares.

Radial Velocities from Spectroscopy
We determined the radial velocity (V R ) of the coadded spectra for each star using two methods. First, we model the full spectrum of each star using The Payne, for which V R is a free parameter as described Section 4. Second, we measure the Doppler displacement of specific spectral lines using a Gaussian-fitting method to obtain the radial velocity. The latter method is more commonly used in the literature (e.g., Sana et al. 2013;Camacho et al. 2016;Evans et al. 2019), and we use it as a check on the full spectral fitting approach. We follow the cited literature studies and fit Gaussian line profiles to a collection of Balmer absorption lines (Hα, Hβ, Hγ, Hδ) and He absorption lines (He I 4387, He I 4713, and He II 4686) that have been identified as only mildly impacted by winds in the temperature regime. We exclude any lines that show strong emission. If the emission is small relative to the absorption feature, we mask out pixels with an emission component. Hα is only measured for stars in the Binospec sample, since it is not in the wavelength range of the Keck spectra. We report the radial velocity of each star as the errorweighted mean of the velocities measured for each line. We apply heliocentric corrections of v helio = −5.08 km s −1 to our Keck measurements and v helio = −6.49 km s −1 to our Binospec measurements.
We find that all stars agree within 1σ between the two methods (see Table 4). This confirms that full spectral fitting can provide accurate V R measurements for S/Ns as low as ∼10 (lowest S/N in sample). V R measurements may be accurately measured at lower S/Ns as well, but will be investigated in future work. We further find that the measurements yielded by the full spectral fitting are more precise (σ Gauss  32 km s −1 , σ Payne  11 km s −1 ). For the remainder of the analysis, we use the V R measurements derived from the full spectral fitting and reported in Table 4. Figure 7 shows the distribution of measured radial velocities with select literature measurements overplotted for comparison. For our eight Keck stars, we find a mean radial velocity of V R = 23.7 ± 16.3 km s −1 . For our 10 Binospec stars, we find a mean radial velocity of V R = 29.4 ± 16.2 km s −1 . They are consistent with each other and agree well with radial velocities of supergiants in Leo A (V R = 22.9 ± 10.2 km s −1 ; Brown et al. 2007) and of the old red giant branch stars (V R = 26.2 ± 9.0 km s −1 ; Kirby et al. 2017). Our combined Keck and Binospec samples for nonoverlapping stars yield a mean V R = 26.9 ± 12.3 km s −1 . The stellar measurements are in line with previous gas velocity studies (Young & Lo 1996;Hunter et al. 2012).
Three of the four stars observed with both telescopes yield consistent velocities. The Keck spectrum of the fourth star (K2) yields V R = −13 ± 2 km s −1 , whereas the Binospec spectrum yields V R = 16 ± 4 km s −1 . There is no obvious issue with the data that would cause this level of discrepancy. One possibility is that this is a binary star, and that the different velocities are the result of observing the system at different points in the orbital phase. We are in the process of obtaining time resolved spectroscopy to better quantify the degree of binarity in the massive stars of Leo A.  (2) spectroscopic mass; (3) initial mass; (4) evolutionary mass; (5) effective temperature; (6) surface gravity; (7) luminosity; (8) rotational velocity; (9) heliocentric corrected radial velocity from full spectral fitting; and (10) heliocentric corrected radial velocity measurements using Gaussian profile fitting. Columns (4), (5), (7), and (8) are the stellar parameters that directly resulted from the full spectral fitting with The Payne. Columns (2) and (6) are derived using the previous stellar parameters in addition to the SED. Column (3) lists the evolutionary mass determined by eye. Figure 7. The distribution of radial velocities for stars in Leo A. Our LRIS (blue) and Binospec (pink) velocities are consistent with the Keck/DEIMOSbased velocities of red giants (RG) measured by Kirby et al. (2017), and the average velocity measured from the Supergiants (SG) by Brown et al. (2007). Each sample is normalized by counts.

Stellar Parameters from Spectroscopic Analysis
We present a discussion of our spectroscopic sample. This section includes an example of each stellar type found in the study; a detailed individual discussion of all stars can be found in the Appendix B. We determine these stellar parameters for the stars using full spectral fitting via The Payne and via photometric SED fitting with the BEAST. We also use EWs to provide spectral types using conventional stellar classification schemes in the literature. We discuss agreement between the methods. The spectrum of each star is shown in Figure 8, and Table 3 provides a full list of lines we detect and measure in each spectrum. For guiding the discussion, we group the stars into broad categories (e.g., O-stars and B(e)-stars) and, for brevity, focus only on the most salient features in the spectra and photometry.

O-type Stars
We classify four stars as O-type stars in our sample (K1/B1, K2/B2, K7/B6, and B13) based on their EW measurements. Overall, the parameters derived from spectroscopy and photometry yield consistent values with the EW classification. We present a detailed example of the analysis of K1 and highlight interesting features of the other stars. K1/B1: K1/B1 is the brightest star in our Leo A sample (F475W VEGA = 20.09 and F336W VEGA = 18.36). Its locations in the optical and UV CMDs are typical of an O-star. Its Figure 8. Keck/LRIS and MMT/Binospec normalized spectra (with an arbitrary offset added for clarity) for all 18 stars in the spectroscopic sample. For stars with both Keck and MMT data, we show the Keck spectra in the 4000-5100 Å region. If there is no LRIS spectrum, we show the Binospec (1000 lpmm) spectrum. The spectra in the 6520-6610 Å region are from the Binospec 270 lpmm grating and are shown for all stars that have Binospec spectra. The spectra are sorted by spectral type from O-type stars at the top to B-type and Be-type stars at the bottom. This is one of the largest set of optical spectra for sub-SMC metallicity massive dwarf stars. spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of a typical late-type O-star: strong absorption of He I and He II. Most of the Balmer series is in absorption, though Hβ appears to exhibit weak emission, likely due to K1ʼs location in an H II region (see Figure 3).
We first estimate the stellar type of K1 using the relevant EWs, which are listed in Table 3. We use both Galactic (Conti & Alschuler 1971;Sota et al. 2011;Martins 2018) and subsolar metallicity classification schemes (Castro et al. 2008), prioritizing the latter whenever possible. From the EW(He I λ4471)/EW(He II λ4542) ratio, we find that K1 is of spectral type O9, or perhaps later, in the traditional spectral classification scheme of Conti & Alschuler (1971). We also consider EW ratios of the (He I λ4471/He II λ4542) and (He I λ4388/ He II λ4542) in the classification scheme of Martins (2018). They indicate a star of type O8.5-9. The He II λ4542 to He I λ4388 comparison yields a spectral type O9.5 in the Sota et al. (2011) classification scheme. Lastly, using the Castro et al. (2008) classification scheme, we find the star to be an O9. K1 could be either a dwarf or a subgiant star based on the comparison with the Martins (2018) classification scheme using (log(EW(He II λ4686))+log(EW(He I λ4388))), since its value falls within both ranges. However, using Castro et al. (2008), we determine the star is likely a dwarf. Our best estimate from spectroscopy is that K1 is an O9V star. We employ the same EW classification scheme for our other O stars.
We fit the full Keck spectrum of K1 using The Payne. The median parameters for the spectral fitting are listed in Table 4. For K1, we find T log eff ( )= 4.49 ± 0.01 dex and g log( ) = 3.69 ± 0.01 dex, which, like the EW ratios, indicates that K1 is a late-type O star. Using Equations (1)-(3), we derive M spec = 17 ± 3M e and M 0 = 19 ± 1M e , and  L L log( )= 4.90 0.09  dex. All of these values agree with the classifications derived from the EW measurements and location in the CMD. The photometric analysis yields T log eff ( )= 4.52 ± 0.02 dex and g log( ) = 3.92 ± 0.05 dex. While T log eff ( ) is consistent, g log( ) shows a 0.23 dex difference. This difference may be explainable by three different reasons: (1) keeping the instrument resolution, and micro-and macroturbulence constant might impact g log( ) since they all affect the line widths; (2) subtraction of the nebular emission during data reduction; if the Balmer lines are slightly over/under-subtracted, then the line width and g log( ) can be affected; and (3) assumptions within stellar atmosphere models give rise to systematic differences (∼0.2 dex; Markova et al. 2018;Dufton et al. 2019). Computing our own systematic uncertainties requires delving into stellar atmosphere physics for very metal-poor stars, and is beyond the scope of this paper. We find that M 0 and M ini are consistent with each other, as seen in Tables 4 and 5.
We show the spectrum of K1 with the models generated from the two different methods in Figure 9; we overplot synthetic spectra of K1 from our spectroscopy (pink) and photometric (blue) parameters. Both synthetic spectra provide a similar match to the He I/He II feature, as expected since our T log eff ( ) is in agreement. We find the difference in g log( ) is best displayed in the Hβ feature, primarily in the wings. The Balmer lines are well fit by the models, which we highlight by showing a zoom in on Hβ. However, while the model appears to match the He I/He II ratio, it struggles with the recovering the depth of the He I lines, in that the model is shallow. This is clear from He I λ4471. A lower temperature may improve the He I fit, but this would make the fit to He II worse. This difference could be due to the fixed resolution of the spectra or, alternatively, could be due to the H to He ratio that is expected to change at lower metallicity (e.g., Martins & Palacios 2021), but is currently fixed at the solar value. Despite some modest disagreements, the model spectrum appears to provide a good overall characterization of K1. Our formal fitting consistently suggests K1 is a late-type O dwarf.
K2/B2: K2/B2 is an O9.7 V type star. We find a small mismatch in g log( ) between spectroscopy and photometry. The emission line contamination of the Balmer lines may affect our spectroscopic determination of g log( ). The discrepancy in g log( ) is likely why M spec is in disagreement with M act . We further note that M 0 < M spec , which is physically impossible.
sample. We suggest that K2 may be a binary star system. This is potentially further supported by the odd present-day-toinitial-mass relation (see Section 17). Currently it is not possible to verify the binary status with our limited time information. We are in the process of acquiring time series spectra with Keck to further investigate this possibility. K7/B6: K7/B6 is an O9.5 V type star and yields consistent parameters from photometry and spectroscopy. Although K7 is classified as an O-type star, it is not near any of the known H II regions in Leo A. We discuss the possibility of it being a runaway star in Section 6.2.
B13: B13 yields an S/N of ∼12 making the classification challenging. The star is likely an O9.7 V type star, but could be consistent with a B0 V type star as well. Its Hα feature is fully in emission, and the remaining Balmer lines are in absorption with some light emission within the feature. Given that the star shows narrow emission features, poor spectroscopic fitting, and extended emission in the 2D spectra, we postulate that the star must have some nebular contamination.

B-type Stars
We classify eight stars as B-type star in our sample (K3/B3, B4, B7, B9, B10, B11, B12, and B14) based on their EW measurements. Overall the parameters derived from spectroscopy and photometry yield consistent values with the EW classification. We present a detailed example of the analysis of K3/B3 and highlight interesting features of the other stars.
K3/B3: K3 is the third brightest star in our sample (F336W VEGA = 18.87 and F475W VEGA = 20.4). Its locations in the optical and UV CMDs are typical of an MS star. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum broadly suggests the star is a B-type star: He I, Mg II, Si II, and Si III are in absorption. The Balmer series shows no Hα, while Hβ is primarily in absorption. For Hβ, there is emission in the red wing, which could either be consistent with the wavelength of two O II features at λ4865Å and λ4872Å and the forbidden [Fe IV] emission at λ4868Å, or the kinematic components of some Hβ emission. The remaining Balmer lines appear to have little to no emission in their red wings. The Hβ in the 2D spectra shows that the emission is spatially extended with a slope toward redder wavelengths, as opposed to perpendicular to the trace. The features appear to spatially extend to ∼5.5 pc. We find Fe II at λ6540Å in emission to the blue side of where Hα should be. We further find Fe II in emission at λ8617Å and [S IV] at λ8575Å. Based on the EW, we use the Castro et al. (2008) and Evans et al. (2015) classification schemes and find that K3 is most likely a B1.5 V star.
From full spectral fitting, we find K3 has median values of T log eff ( ) = 4.42 ± 0.01 dex and g log( ) = 3.75 ± 0.05 dex, which, like the EW ratios, indicate that K3 is an early-type B star. Using Equations (1)-(3), we derive M spec = 20 ± 5M e , M 0 = 13 ± 1M e , and  L L log 4.53 0.12 ( )=  dex. These are all in line with the classification derived from the EW measurements and location in the CMD. The spectroscopic parameters are consistent with the results from SED fitting, which yield T log eff ( ) = 4.41 ± 0.01 dex, g log( ) = 3.66 ± 0.01 dex, and M ini = 12.9 ± 0.17M e .
We postulate that the majority of the emission is due to the nebula southwest of the star, which appears to fall into our slit, based on the HST imaging. The nebula can explain the emission features near Hβ as well as Fe II and other emission at redder wavelengths. The 2D spectra further support this scenario, since the emission primarily extends to only one side of the stellar trace. Furthermore, the nebula shape suggests that we may be observing a bow shock nebula (e.g., Gull & Sofia 1979;van Buren & McCray 1988;Meyer et al. 2014). This is supported by the following: K3 is a runaway candidate (see Figure 9. Comparison between the parameters derived for K1 using full spectral fitting and SED fitting. We show a plot of the observed spectra (black dots), the TLUSTY-based model generated by the median parameters derived from the full spectral fitting (pink), and the TLUSTY-based model generated by the median parameters derived from the SED fitting (blue). While both models reasonably replicate the observed He I 4471 and He II 4541 fits, both struggle to fully fill the absorption feature of the helium lines. This shortcoming is potentially due to the fixed resolution, metallicity, or macroturbulence. We see the effect of the different g log( ) values between the methods most prominently in the Balmer lines (see the Hβ feature). The gray bar shows the part of the spectrum that we masked. Section 6.2), the nebula is the most visible in the F110W and F160W imaging, and the absence of Hα. To capture the intrinsic characteristics of the nebula, its formation, and the star-nebula interaction, follow-up observations are necessary, ideally with an integral field spectrograph.
B4: B4 is a B0 V type star and yields consistent parameters from photometry and spectroscopy. The presence of emission lines coupled with the star not being near an H II region suggest the emission is not gaseous in origin. Instead, this may indicate the presence of a circumstellar nebula or of stellar winds. Given that Hα shows only weak emission and He II 4686 in absorption, the stellar wind scenario seems unlikely (e.g., Lamers & Cassinelli 1999). Alternatively, the star could have a companion that provides the emission components to the spectra.
B7:B7 is a P-Cyg B1 V or a Be1 V type star. The star resembles the traditional P Cyg star (Rivet et al. 2020, and references therein). We do not observe metal lines with P-Cygni profiles. This is not surprising given Leo A's low metallicity. However, the lines are narrower than what would be expected from stellar wind emission lines. Alternatively, the star could be a Be star, whose disk is oriented so that we coincidentally capture the emission resembling the P-Cygni profile. This could also explain the lower v i sin value, since Be stars typically have higher v i sin values (e.g., Iqbal & Keller 2013;Rivinius et al. 2013;Arcos et al. 2018). Potentially this could also explain why, in comparison to the other Be stars in this sample (see Section 5.3.3), the star does not show signs of potential NIR excess in the optical CMD. Extended spectroscopic coverage into the UV could help shed light on wind parameters, or photometric coverage into the NIR could shed light on a potential disk. We discuss a more clear example in Section 5.3.3.
B9: B9 is a B1.5 V type star. Hα appears in emission, and all other Balmer lines are in absorption. B9 yields consistent parameters from photometry and spectroscopy.
B10: B10 is a B0.5/Be0.5 V type star. The photometric parameters and spectroscopic parameters are consistent. Both methods yield decent to high χ 2 values (see Tables 4 and 5) suggesting that potentially a circumstellar disk or nebulae could affect the star. The stellar emission and the high v i sin value would support a Be star scenario. We discuss more clear examples in Section 5.3.3. B11: B11 is a B2 V type star. Most Balmer lines are in absorption with variable amounts of emission. Since the B11 spectrum shows an S/N of <10, we cannot reliably recover the spectroscopic parameters. Given the S/N is ∼10, we refrain from making any further speculations as to the origin of the emission.
B12: B12 is a B1.5 V type star. The star shows signs of stellar emission, but the low S/N (∼12) prohibits detailed analysis. We do note that we only have four bands in the photometry. B12 is one of the few stars that we do not find consistency between the two methods, although both χ 2 values would suggest they are good or decent fits.
B14: B14 is a B3 V type star. Hα is in absorption, while the other Balmer lines show potential emission in their wings. Due to the low S/N (∼15) of the spectrum, the weak emission in the Balmer lines could be artifacts of the data.
T log eff ( ) is consistent between the two methods, while g log( ) shows a discrepancy likely due to the affected Balmer lines.

Classical Be Stars
Our sample contains several objects that appear to be classical Be stars: rapidly rotating, main-sequence B stars with emission lines that originate in a circumstellar accretion disk (Rivinius et al. 2013). The key spectroscopic feature is emission in the Balmer lines (usually strongest in Hα), while the key photometric feature is an IR excess due to the extended circumstellar disk. We provide the analysis of the clearest example (K6) in this section. In total we find six stars (K4, K5, K6, K8, B5, and B8) for which we can identify both key features. They all show strong Balmer emission lines, and their particular SEDs cannot be reproduced by a "normal" single star model.
K6: K6 (F475W VEGA = 20.84 and F336W VEGA = 19.35) is located on the BHeB sequence in the optical CMD, but appears on the MS in the UV-optical CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features. Stars K4, K5, K8, and B5 share similar CMD characteristics. The Keck spectrum of K6 is shown in Figure 8.
The spectral features suggest K6 is a fast rotating early B-type star. K6 shows He I in absorption and emission, and Mg II, S II, and Si IV in absorption. He II appears to potentially be a mix of emission with absorption. We find strong emission in Hβ. We find an emission feature within an absorption feature for Hγ and Hδ. However, Hò appears to be mostly in absorption. All of the features are very shallow, which is typical of fast rotators. The signs of strong rotation and emission mean that detection of EWs might be harder to recover for certain features (e.g., Si, Mg, and He II), since they become so shallow that they are indistinguishable from the noise. Based on Lennon (1997), Castro et al. (2008), and Evans et al. (2015), we classify it as a Be2 V star.
The strong emission appears to be associated with the star, as opposed to being gaseous in origin. K6 is not near the known H II regions (see Figure 3). Moreover, the lack of typical nebular emission lines further suggests it is not from an H II region. Leščinskaitė et al. (2022) identified K6 as a strong Hα emitter based on the Subaru/Suprime-Cam narrowband imaging, postulating that it might be an emission star.
Due to the strong emission in the Balmer lines, the likely presence of a disk and the high χ 2 in the photometry, we adopt here the alternative fitting methods introduced in Section 4.3.
Using this adapted spectroscopic method, we determine that the star has T log eff ( ) = 4.43 ± 0.02 dex and a v i sin = 220 ± 10 km s −1 . This is in line with an early-type Be star (Rivinius et al. 2013). The high v i sin is a key indicator that the star has to be a rotating MS star as opposed to a BHeB star. An evolved star would break apart at such high v i sin . As previously mentioned, the star is among the bluest in the UV CMD but among the reddest in the optical CMD.
Using our composite photometry fitting (see Section 4.3), we determine that T log 4.48 eff star 0.03 0.04 ( ) =  dex, which is in line with our spectroscopic measurement. The disk T log eff disk ( ) = 3.73 0.03  dex. This serves primarily as an estimate, since we use the model of a cool star to mimic the disk spectra. Figures 10 and 11 show the median photometric model along with samples from the posterior. In Figure 10 we can clearly see that the UV photometry is fitted well by a rotating MS B-type star. In the NIR, the fit is not as good. However, since we are not actually fitting a disk model, this could explain the discrepancy. Overall, the fit yields a good χ 2 value (see Table 6).
We considered other explanations for the unusual photometric SED. The first possibility could be dust giving rise to NIR excess. However, given how blue the star is in the UV CMD, we believe dust is unlikely. A second possibility is an unresolved red stellar companion. We tried fitting two stellar models simultaneously. However, one of the stars always yields nonphysical parameters (i.e., a temperature and radius combination that is not expected for any star). Given that there is no evidence for a red star in the spectra, and we are not able to fit a physically reasonable star to the photometric data, this scenario is unlikely.
We note here that this does not mean that the star is necessarily a single star, simply that there is no evidence for a secondary polluting our current data.
Based on the present evidence, we suggest this is a Be star. This provides the first spectroscopic confirmation for Be stars in Leo A.
K4: K4 is a Be1 V type star. Due to the strong emission in the Balmer lines and the high χ 2 (>300) in the BEAST fitting, we postulate it is a star surrounded by a disk. Another possibility is an unresolved red stellar companion. For this particular star, we cannot completely exclude this possibility, since we do not have the NIR photometry to determine the nature of the excess. However, given that for other Be stars in our sample we find the same characteristics as for K4, we postulate that extinction occurring in the photometry is due to a disk.
K5: K5 is a Be2 V type star. Similar to K4, we assume that the star is a Be star whose photometry is affected by extinction/reddening due to the disk, as opposed to extensive dust or a red companion.
B5: B5 is a Be3 V type star. The 2D spectra show extended emission, which suggest it is due to nebular emission from the nearby H II region. However, the Hα/ Hβ emission ratio suggests that the Hα emission cannot be purely gaseous. While our data cannot exclude a binary influencing the star (e.g., we need time series spectra), the spectroscopic features and SED fit are consistent with a circumstellar disk.
K8: K8 is a Be2 V type star. Similarly to K6, the star is a Be star whose photometry is affected by extinction/reddening due to the disk as opposed to extensive dust or a red companion.
B8: B8 is a shell Be3 V type star. Hα and Hβ are in emission and both show a double-peaked line profile. The distinct Balmer profiles suggest that the star is most likely a shell star (i.e., a Be star viewed edge-on). Despite the low v i sin , the strong emission and spectral shape support the classification of the star as a shell star.

The Utility of Panchromatic Photometry for Massive Star Studies
We have shown the first extensive and quantitative comparison between stellar parameters for massive stars derived from optical spectroscopy and UV-optical-NIR broadband photometry at low metallicity. Figure 12 shows the comparison between photometric and spectroscopic T log eff ( ), g log( ), and  L L log( ) for stars that have χ 2 < 100 in both methods. We find that a median absolute deviation (MAD) between the two T log eff ( ) of <0.01 dex and a median difference of −0.03 dex (1600 K). The 16th and 84th percentiles of the difference span 0.05 dex (i.e., ±0.025 dex around the median), indicating that the observed differences are slightly smaller than scatter in the data. For g log( ), the MAD is 0.18 dex and the median difference is 0.06 ± 0.20 dex. Lastly, for  L L log( ), we find an MAD of 0.05 dex and a median difference of −0.03 ± 0.06 dex. The clear outlier in each panel is B12. B12 is among the lowest S/N stars in our spectroscopic sample (S/N∼ 12) and has modest emission, though not as strong as other stars, which could explain the larger differences between the photometric and spectroscopic parameters. In general, these findings show that for "normal" O/B type MS stars in our sample (i.e., those with no or little emission), photometric and spectroscopic parameters are consistent with each other.
This comparison clearly shows that panchromatic photometric studies can yield reliable stellar parameters at least for "normal" single stars, which is encouraging especially when spectra are not available. For example, most of the LUVIT sample consist of galaxies too far away for acquiring optical spectra with the current generation of telescopes. Poor photometric fits with the BEAST are often due to stellar activity (e.g., binarity, disks) that are not captured by single star models. We find that with full UV through NIR broadband SEDs it is possible to identify such candidates by their poor χ 2 values, providing an invaluable way to identify unusual objects for future follow-up (e.g., with JWST, ELT, GMT, and TMT). We note that the BEAST is in development of including more models.

Stars in H II Regions
Leo A hosts several known H II regions (see Figure 1 of van Zee et al. 2006), two of which are included in the HST and the spectroscopic data set. Previous analysis of these H II region ionized gas abundance patterns have placed indirect constraints on the properties of the ionization fields, and in turn, the underlying types of ionizing sources (e.g., the spectral type of the OB stars; Skillman et al. 1989;van Zee et al. 2006;Ruiz-Escobedo et al. 2018). There have been a few direct studies of Figure 10. Models for K6 from adapted composite SED fitting (see Section 4.3). The data (black dots) are matched by scaling and combining the stellar models. The star (blue line) and the disk (red line) are combined for the composite SED (black line). The stellar and the disk models are both varied to find the median SED. We note that the error bars of the data are smaller than the dots representing the data point. the massive stars that power H II regions at such low metallicities (e.g., Garcia et al. 2014;Camacho et al. 2016;Evans et al. 2019;Telford et al. 2021b). In most cases, these studies have found that a singular star is responsible for powering most of the H II region.
In this section, we discuss the results of our stellar spectroscopic and SED fitting in the context of literature nebular line analysis of the H II regions in Leo A. Figure 13 shows the two H II regions that are covered by our UV and optical HST data along with their corresponding optical CMDs. Stars located approximately within each H II region are color-coded by their age from our SED fitting, while the CMD of the entire Leo A field population is shown for reference in gray. For clarity of discussion, we label the two H II regions: Region A and Region B. We have six band coverage for Region B, and four band (two UV and two optical) coverage for Region A.    (2) spectroscopic temperature of the star; (3) photometric temperature of the star; (4) photometric radius of the star; (5) rotational velocity; (6) reddening value from the SED fitting; (7) photometric temperature of the disk; (8) photometric radius of the disk; (9) present-day mass picked for the combined SED fitting; (10) heliocentric corrected radial velocity from full spectral fitting; (11) heliocentric corrected radial velocity measurements using Gaussian fitting; and (12) reduced χ 2 value from the composite photometric fitting. Columns (2), (5), (10), and (11) are the stellar parameters derived from the adapted spectroscopic analysis. Columns (3), (4), (6), (7), (8), and (12) are the stellar parameters derived from the composite photometric analysis.
The brightest star in Region A is K2. From our spectra and photometry, we find that K2 is an O9.7V type star with a T log eff ( ) = 4.50 ± 0.01 dex, which is very similar to what van Zee et al. (2006) suggested. Furthermore, K2 has an age of t log 7.0 0.07 ( ) =  (yr), M ini = 16.81 ± 0.21M e , and A V = 0.06 ± 0.02. For K2, we are able to calculate the intrinsic ionizing photon production rate using the BEAST (Choi et al. 2020). We derive a median rate of 10 47.6±0.1 photons s −1 . Using the Hα value reported in van Zee et al. (2006), we compute the photoionization rate required to produce the Hα luminosity (Hummer & Storey 1987;Storey & Hummer 1995;Kennicutt 1998;Choi et al. 2020). We derive a rate of 10 47±0.05 photons s −1 . This is consistent with the value derived for K2 from the BEAST. This suggests that K2 could be responsible for a majority of the ionization in the H II region. Furthermore, this finding supports the, thus far, low extinction from dust along the line of sight, since we would otherwise expect our derived stellar value to be larger than the value reported in literature.

Region B
Region B is the right region in Figure 13. It is described as the east H II region in Ruiz-Escobedo et al. (2018), who also noted that it is faint and appears to have low extinction (c (H β ) = 0.08 ± 0.05). Ruiz-Escobedo et al. (2018) detected only the hydrogen Balmer series in emission as well as [O II] λ3727Å, [N II] λ6583Å, and [S II] λ6717, 6731Å. Since they were only able to obtain an upper limit for [O III]λ5007Å, they determined that a star with T log eff ( ) 4.45 is necessary and an ionizing photon production rate of 10 48.0 photons s 1 is required to produce the observed H II region size of 1.5 pc, assuming an electron density of 100 cm −3 .
Star K1 is located in Region B. From our spectra, we find K1 is an O9V type star. We find T log eff ( )= 4.49 ± 0.01 from our spectroscopic analysis and T log eff ( )= 4.52 ± 0.02 from SED fitting. Both are in reasonable agreement with the estimate of T log eff ( )4.45 from Ruiz-Escobedo et al. (2018). Our SED fitting shows that K1 has t log 6.9 0.08 ( ) =  (yr), M ini = 20.1 ± 1.29M e , and A V = 0.05 ± 0.01. The age and mass of this star are comparable to what is expected for an H II region, while the extinction is similar to what is derived from the nebular emission analysis. We find that B5 is in the H II region as well. However, B5 is a Be star, which makes the age and mass determination more difficult. For K1, we are able to calculate the intrinsic ionizing photon rate using the BEAST. We derive a median rate of 10 48.1±0.1 photons s −1 , which is consistent with the derivation from Ruiz-Escobedo et al. (2018). This supports that the majority of the ionization is likely coming from K1. Again, our results suggest that ). The bottom panel compares g log( ). The comparison between the two methods shows that panchromatic photometry can be used to recover reliable stellar parameters at least for "normal" single stars. dust along the line of sight in Leo A has to be low to negligible.
Overall, we find consistency between previous studies of the H II regions in Leo A and our study of massive stars in those regions. van Zee et al. (2006) considered whether the lack of hotter O-type stars can be explained by either a truncated initial mass function (IMF) or that the H II regions are evolved. To explore this, we consider the ages derived by the BEAST of the stars (K1, K2, and B5) in those regions ( Figure  13). Given that those stars are high-mass stars, they are already at the end of their lifetime on the MS. In Figure 13, we can clearly see that Region B is a more evolved H II region than Region A, since the average age of the massive stars in the region is older. However, we also find very few massive stars in both H II regions. Lastly, the results suggest that in both H II regions a majority of the ionization is coming from a singular O-star.

Field OB-stars
It is well established that most stars form in clustered environments (e.g., Lada & Lada 2003;Krumholz et al. 2019). Accordingly, OB-stars are expected to be found in or near H II regions or star clusters, as their short lifetimes imply they do not have time to wander to significant distances before ending their stellar lives.
However, a subset (20%-30%) of massive stars are found in the "field" with no clear association to clusters or H II regions (e.g., Dorigo Jones et al. 2020). The existence of those stars poses a dilemma for stellar kinematics, star formation theory, IMF sampling, etc. (e.g., Blaauw 1961;Clarke & Pringle 1992;Hoogerwerf et al. 2000;Oey et al. 2004;de Wit et al. 2005;Renzo et al. 2019). This has led to major differing theories to the formation of those stars: in situ star formation (Oey et al. 2013) and ejection Dorigo Jones et al. 2020). Recent studies (e.g., Oey et al. 2018;Dorigo Jones et al. 2020;Vargas-Salazar et al. 2020) favor the latter scenario in most cases, which is further broken down into two different type of ejections: a dynamical ejection scenario (DES) and a binary supernova scenario (BSS). While both require massive stellar binaries, in the DES, a close encounter with a binary system ejects the star (or binary) at high velocity out of the cluster or H II region. In the BSS scenario, a corecollapse supernova of the more evolved star ejects its OB companion or sometimes the entire system (Dorigo Jones et al. 2020).
The majority (∼70%, 13 stars) of the stars in our spectroscopic sample are not within d ∼ 40 pc of any H II region. To obtain a rough understanding of which mechanism is responsible for the position of the stars, we derive estimated required velocities from the age estimate yielded by the BEAST and distances from the H II region map (Figure 3). We summarize the derived values in Table 7. In this discussion, we take advantage of the coincidence that 1 pc Myr −1 is 0.98 km s −1 .
K3, K7, B4, and B10 can be considered runaway candidates, since their lower limit for velocity suggest that the stars could have been ejected. All four stars are at least 200 pc from the nearest H II region or cluster. They are among the youngest field stars (with the following ages in megayears: K3: 15.8  , and B10: 12.6 2.8 3.6  ). Although recent studies (e.g., Renzo et al. 2019;Schoettler et al. 2020;Vargas-Salazar et al. 2020) suggest a cutoff for runaways around ∼30 km s −1 , we still can consider those four stars runaway candidates, since we only report lower limits. Due to the lack of precise velocity measurements, it is currently not possible to distinguish whether the stars were ejected through DES or BSS. B11, B12, and B14 are more likely walkaway candidates (de Mink et al. 2014;Renzo et al. 2019;Schoettler et al. 2020;Vargas-Salazar et al. 2020  ). Their estimated velocities are of the order a couple of parsecs per megayear. This is in line with Note. The columns are: (1) star name; (2) initial mass (photometry); (3) initial mass (spectroscopy); (4) spectral type; (5) stellar age in megayears; (6) radial velocity; (7) approximate distances to the closest H II region or cluster for each star; and (8) approximate velocity derived from taking the value in column (7) and dividing it by the values in column (5). When available we list the initial mass; otherwise, we list the evolutionary mass. 1-The star has a really high χ 2 -value making the age unreliable and therefore the velocity unreliable as well.
intermediate-mass runaway velocities (Renzo et al. 2019). However, these types of walkaways remain rare (Renzo et al. 2019), and therefore are impossible to confirm without precise velocity measurements, much less to associate with different mechanisms. Five of our Be stars (K4, K5, K6, K8, and B8) and B7 are also in the field. Since we do not have age measurements, we can only present rough velocity ranges in Table 7. All five stars show velocity limits consistent with runaway stars. In the literature, Be stars in the field have been associated with BSS (Dorigo Jones et al. 2020). However, we cannot exclude that DES could be responsible for the ejection and velocity of these stars.
Further studies of OB-stars in Leo-A would likely provide a better understanding on the overall distribution of the mechanism responsible for runaway stars at extremely low metallicity.

Purity of the Blue Core Helium Burning Sequence
BHeB stars provide important constraints on stellar evolution and serve as excellent tracers of recent star formation (Langer & Maeder 1995;Dohm-Palmer et al. 2002;Dolphin et al. 2003;Larsen et al. 2011;Leščinskaitė et al. 2022). In particular, the extent of the blue loop and the ratio of blue to red core helium burning stars are sensitive to several poorly constrained aspects of massive stars physics (e.g., rotation, convective overshoot; Alongi et al. 1991;Ritossa 1996;McQuinn et al. 2011;Jie et al. 2015;Georgy et al. 2021;Eldridge & Stanway 2022;Farrell et al. 2022). The luminosities and spatial positions of BHeBs in nearby galaxies have been used to construct spatially resolved maps of star formation over the past few hundred megayears (Dohm-Palmer et al. 2002;McQuinn et al. 2012;Leščinskaitė et al. 2022).
However, some of the data and findings presented in this paper suggest that not all stars that reside in the BHeB region of an optical CMD are bona fide BHeBs. This has been previously suggested by  and Leščinskaitė et al. (2022). In our combined photometric and spectroscopic sample, we find that indeed, six stars in our sample out of seven that appear to be in the BHeB region on the optical CMD (K4, K5, K6, K8, B5, and B8) are, in fact, spectroscopically Be-type stars. This provides the first spectroscopic confirmation for Be stars in Leo A.
We compare the UV and optical CMDs to see if any other stars that appear to be in the MS region in the UV CMD are in the BHeB region in the optical. We compare stars that are brighter than 23.5 mag in the F475W filter. We further apply the following cuts: F475W-F814W >0.2, F475W-F814W < −0.15, and F275W-F475W <0.0. Lastly, we apply cuts by eye on stars that are on the MS region. These cuts are visualized in Figure 14. We find that ∼5%-10% of the stars in the BHeB region above 23.5 mag in F475W exhibit this behavior. This finding suggests that when taking the BHeB/RHeB-ratio from photometric surveys, we have to take into account that the BHeB region could be contaminated by Be stars. However, since we can disentangle the Be star candidates in the UV, corrections can be applied when multiwavelength coverage is available. Lastly, the findings propose a new avenue to potentially search for and detect Be star candidates through photometry.
In Figure 15, we compare our spectroscopic sample stellar types to other samples of massive stars with spectra at SMC and sub-SMC metallicities with published temperature and spectral type. We compare the four O-stars to the best fit of the Massey et al. (2009) sample and the Dufton et al. (2019) T log eff ( ) to spectral type relation. We further compare the stars to two stars in Sextans A: S3 Telford et al. (2021b) and OB521 Camacho et al. (2016). We note that two of our stars (K2 and B13: T log eff ( ) = 4.50 dex and SpT is O9.7V) fall on top of each other in the T log eff ( ) versus spectral type plot, which is why we only see three O-stars in the plot.
We find that our O-stars are consistent with the sub-SMC stars S3 and OB521. When comparing to the SMC trend lines, we see that our parameters are consistent with Massey et al.  The pink colored circles are stars that appear in the BHeB region both in the UV and the optical CMD. The blue colored circles are stars that are on the MS in the UV CMD, but appear on or near the BHeB region in the optical CMD. The cyan squares are the stars identified as Be stars in our spectroscopic sample. This shows candidate Be stars could possibly pollute the BHeB in the optical CMD, and that we can filter for candidate Be stars in panchromatic photometry studies. expectation (see, e.g., Mokiem et al. 2004). One possibility is the difference in the underlying models. While Telford et al. (2021b) and our work use the TLUSTY models, other groups use other models (e.g., CMFGEN, Hillier & Lanz 2001;FASTWIND, Santolaya-Rey et al. 1997), introducing at least some uncertainty.
For the photometrically determined values of T log eff ( ), we find that for two of the stars (K2 and B13), we derived hotter T log eff ( ), while for the other two (K1 and K7), the photometry yields the same T log eff ( ). We compare our B-stars to the SMC relation for B-stars from Dufton et al. (2019) for the dwarfs and to Dunstall et al. (2011) for the Be stars. For the B0-B1 stars, we find that we are cooler than the Dufton et al. (2019) derived relation. For stars after B1, the stars seem to be closer to the derived relation and hotter. We find that all of our Be stars are hotter than the derived relationship from the SMC for Be stars, both spectroscopically and photometrically. We did not find any sub-SMC B V stars to compare our sample to. While, for the O-stars, the slope between our points appears similar to the literature, we find that we have a flatter fall off for our B-stars. Since the B stars in our sample represent the first MS stars we have both spectral type and photometry for, they provide the first insight into the scale at sub-SMC. Overall, we find that at later spectral types T log eff ( ) decreases in our sample, as expected. However, there is a large spread, which is likely due to decreasing data quality. Increasing our sample of sub-SMC stars is crucial to obtain a full understanding of the relation between spectral type and T log eff ( ).

Current Observational Knowledge of Metal-poor Massive Stars
In Figure 16, we plot HR diagrams for the entire sub-SMC metallicity star sample in the literature with temperatures and luminosities that have been derived from spectroscopy. We note that recently Lorenzo et al. (2022) collected ∼150 spectra of metal-poor massive stars; however, currently there are no stellar parameters for this sample, and we cannot include them on this plot. We overplot the MIST stellar tracks  . The HR diagrams of massive sub-SMC metallicity stars in the Local Group with derived stellar parameters. We overplot MIST tracks. Top: all of the stars with temperatures and luminosities in the metal-poor galaxies: WLM (∼15% Z e ; aqua triangle), IC 1613 (∼14% Z e ; blue circles), and NGC 3109 (∼11% Z e ; purple-blue squares). The WLM stars are from Bresolin et al. (2006), Tramper et al. (2011, 2014, and Telford et al. (2021b); the IC 1613 stars are from Bresolin et al. (2007), and the NGC 3109 stars are from Evans et al. (2007), Tramper et al. (2014), andHosek et al. (2014). Bottom: all stars with temperatures and luminosities in the extremely metal-poor galaxies Sextans A (∼6% Z e ; purple triangle) and Leo P (∼3% Z e ; rose squares). As well as the stars in our spectroscopic sample (∼5% Z e ; pink circles). The Sextans A stars are from Kaufer et al. (2004), Camacho et al. (2016), and Telford et al. (2021b); the Leo P star is from Telford et al. (2021b). Our sample represents the first extended sample of MS stars at sub-SMC metallicities with derived stellar parameters. Furthermore it is the largest sample at sub-10% Z e with derived stellar parameters. at select masses for reference. We indicate the gas-phase metallicities for each galaxy after placing them on a common scale using Grevesse et al. (2010) for the solar abundance (12 log O H 8.69 ). The top panel shows metal-poor galaxies WLM (Z ∼ 15% Z e ; Bresolin et al. 2006;Tramper et al. 2011, 2014and Telford et al. 2021b, IC 1613 (Z ∼ 14% Z e ; Bresolin et al. 2007) and NGC 3109 (Z ∼ 13% Z e ; Evans et al. 2007;Tramper et al. 2014 andHosek et al. (2014)). The bottom panel shows extremely metal-poor galaxies Sextans A (Z ∼ 10% Z e ; Kaufer et al. 2004;Camacho et al. 2016;Telford et al. 2021b), and Leo P (Z ∼ 3% Z e ; Telford et al. 2021b), along with stars in our Leo A spectroscopic sample.
We clearly see that the literature thus far has populated mostly the very upper area of the HR diagram. This is not surprising, since the most luminous stars are more accessible to observe. In comparison, our sample extends to fainter luminosities, later types, and less evolved stars.
We double the number of OB-stars thus far with spectroscopically derived stellar parameters at sub-10% solar metallicity. It also contains the first sample of spectroscopically confirmed sub-SMC Be stars and the largest sample of sub-SMC B-type dwarf stars.

Initial and Present-day Masses across Analysis Methods
The mass of a star is among its most fundamental parameters. Here, we examine the masses we have measured from photometric SED fitting and from spectroscopy. SED fitting yields an initial mass of M ini and a present-day mass of M act . Spectroscopic fitting provides an initial mass of M 0 and present-day masses of M evo and M spec . As a reminder, M evo is determined by using T log eff ( ),  L L log( ) and evolutionary tracks, whereas M spec is determined using Equations (1)-(3) as explained in Section 2.2. The exact derivation of each of these masses is given in Sections 5.1 and 4. The different masses are listed in Tables 4 and 5 and comparison plots are presented in Figure 17. We note that the comparison includes stars that have χ 2 < 100 in both methods.
For the initial masses, M ini and M 0 , we find an MAD between the two masses of 0.59M e and a median difference of 1.1 ± 1.1M e . Overall, the M ini values are slightly larger than M 0 , but they are all roughly within the uncertainty except for B12 (which was pointed out as an outlier in Section 5.4). B12 shows a difference ∼3 times larger than the uncertainty and shows a significantly higher initial mass (M 0 ∼ 11 M e ) for spectroscopy. For the present-day masses of the stars in the sample, we find overall that M act and M evo are consistent with each other. We find an MAD between the two masses of 0.69M e and a median difference of 0.99 ± 1.42M e . Similar to the initial masses, the M act values are slightly larger than M evo , but within the uncertainty. We find B12 remains an outlier.
In Figure 17, we also compare M evo and M spec . We do not include M act , since the values agree well with M evo . In principle, M evo and M spec should be the same, and differences between them are used as a consistency check for evolutionary models. For example, as first noted by Herrero et al. (1992), there have been multiple studies (e.g., Herrero et al. 2002;Massey et al. 2005;Weidner & Vink 2010;Mahy et al. 2015;Markova & Puls 2015;McEvoy et al. 2015;Ramírez-Agudelo et al. 2017;Sabín-Sanjulián et al. 2017;Markova et al. 2018;Schneider et al. 2018;Castro et al. 2021) that discuss the tension between M evo and M spec for O-type stars. The tension rises from M evo being larger than M spec , although they should be the same. Most of these studies suggest that the discrepancy is the most significant at M > 40M e . We do not have any stars this massive, but we can examine differences in the ∼10 − 20M e range. For all of the stars in our sample, we find an MAD between the two masses of 3.0M e and a median difference of <0.01 ± 5.3M e . To get a sense of typical differences for stars that agree, we consider stars that agree within the uncertainty measurement, and we find an MAD between the two masses of 1.5M e and a median difference of 0.0 ± 2.2M e . We find three outliers in our comparison, K2 (∼1.5 times larger), K3 (∼1.5 times larger), and B14 (∼3 times larger). In Section 5.3, we noted that these stars show a slightly higherthan-expected g log( ) value, which is likely driving up the mass, in comparison to the stars' photometry. We note that these stars also show M spec > M 0 , which is not physical, and is likely due to a high g log( ) value or  L L log( ) value. It is interesting that K2 and K3 are outliers, since the K2 observations suggest it could be a binary candidate, while K3 has a peculiar spectra, which can be explained by a binary scenario. Overall, we find that for "normal" single stars, the evolutionary model derived masses and spectra derived masses are in agreement at the level of 1-2 M e . . We show all stars with χ 2 < 100 in both photometry and spectroscopy.

Finding Be Stars in Low-metallicity Environments
We provide the first confirmation of Be stars in Leo A. We also spectroscopically observed the lowest-metallicity Be stars to date. We confirm six Be stars (K4, K5, K6, K8, B5, and B8) and label two B stars as potential Be stars (B7 and B10). Excluding the last two stars, we find that for our spectroscopic sample, ∼45% of B-type stars are Be stars. This is high compared to the Be fraction in the Milky Way (≈17%; Zorec & Briot 1997), LMC (≈20%; Martayan et al. 2006), and SMC (≈30%; Martayan et al. 2007aMartayan et al. , 2007b, and appears in line with the expected trend of an increased Be star fraction at low metallicity due to weaker winds and faster rotation rates. However, our sample is also biased toward early spectral types, where the Be occurrence rate is higher; further study is needed to confirm this trend.
Furthermore, we show that Be stars in Leo A can be photometrically identified by a strong IR excess relative to stars near the main sequence in UV colors, even without spectroscopic follow-up. We show that the panchromatic SED allows us to put some constraints on both the stellar and disk parameters. While we explore the simplest parameters for the Be stars, the combination of our obtained spectra and photometry presents a unique opportunity to gain first insights on the disk physics at low metallicity. Lastly, the identification through an NIR excess is especially promising with JWST. HST UV and Hα imaging surveys can now be combined with the NIR imaging to identify and characterize Be stars.

Conclusions
We have presented the first combined stellar spectroscopy and HST photometry analysis of the extremely metal-poor OBstars in Leo A, a low-luminosity, dwarf galaxy at a distance modulus of μ = 24.40 mag. Our findings from the LUVIT HST, LRIS (KECK), and Binospec (MMT) observations are as follows: 1. We analyze the brightest OB-stars (F475W < 22 and F475W − F814W < 0.5) in Leo A using the LUVIT HST photometry. We spectroscopically observe and analyze a total of 18 stars within the photometric sample. We find four O-type dwarfs, eight B-type dwarfs, and six Be stars. 2. We derive stellar parameters of the photometric sample using the BEAST. We find that for good fits (χ 2 < 80) T log eff ( ) and g log( ) increase for bluer F475W−F814W colors as expected from stellar models. For stars on the MS region of the CMD, M ini trends as expected; the brightest stars (e.g., F275W <18 ) have the highest M ini with masses up to M ini ∼ 20M e . Similarly, A V matches findings from previous studies, which suggest little to no dust in Leo A (A V 0.08). However, for stars on the BHeB region of the CMD, the trends do not appear quite as clear, which we expand upon below. 3. We use The Payne and two different publicly available TLUSTY spectroscopic grids for O-stars and B-stars to measure spectroscopic parameters for 18 stars. We measure T log eff ( ), g log( ), V R , and v i sin for each star. We use MCMC sampling to obtain the uncertainties on all parameters. We find our spectroscopic sample has a T log eff ( ) range of 4.39 to 4.50, a g log( ) range of 3.69 to 4.64, and a v i sin range of <80 to 370. We find V R variation for K2 with Keck (−13 ± 2 km s −1 ) and Binospec (16 ± 4 km s −1 ), making it a viable binary system candidate. 4. The photometrically and spectroscopically derived parameters ( T log eff ( ), g log( ), and  L L log( )) are consistent to within 1σ for stars with χ 2 < 100 in both methods. We find an MAD of <0.01 dex for T log eff ( ), 0.18 dex for g log( ), and 0.05 dex for  L L log( ). When we find disagreement, it is usually due to stars with signs of emission, peculiar SED colors (i.e., IR excess), or high χ 2 values. This is a promising avenue to study more distant low-metallicity galaxies for which resolved stellar spectra cannot be obtained. 5. We compute and compare both initial masses and present-day masses with photometry and spectroscopy. The initial stellar masses (M ini (photometry) and M 0 (spectroscopy)) determined by the two techniques are consistent within 1σ. M ini ranges from 6.26M e to 20.1M e , and M 0 ranges from 8.0M e to 19 M e . Similarly, the present-day masses derived from the stellar parameters in combination with evolutionary tracks (M act (photometry) and M evo (spectroscopy)) are consistent within 1σ. M act ranges from 6.26M e to 19.99M e , and M evo ranges from 8.0M e to 19M e . We further compare M evo to the spectroscopic present-day mass (M spec ), and find that the stars are consistent to 1σ when excluding stars with M spec > 20M e (K2, K3 and B14). These stars show disagreements larger than 5M e . M spec ranges from 8.0M e to 27M e 6. We find that the gas, RGB, and stellar velocities are consistent within 1σ. We further are able to show that full spectral fitting is reliable at recovering RV measurements to higher precision. 7. Overall, we find that our study of massive stars in two H II regions is consistent with previous gas-phase studies of these H II regions in Leo A. The O-type stars (K1 and K2) in the H II regions provide the expected ionization power for each (10 48.1±0.1 photons s 1 and 10 47.6±0.1 photons s 1 , respectively), and suggest little dust along the line of sight. 8. We find a total of 13 isolated stars in our spectroscopic sample, i.e., they are not within d ∼ 40 pc of any H II region. For seven stars, we are able to use the derived ages and approximate distances to the nearest H II region to estimate velocities. Four stars yield velocities consistent with runaway stars (i.e., V  30 km s −1 ), while three are consistent with walkaway stars (i.e., V < 30 km s −1 ). 9. We confirm that six of our spectroscopically observed stars are Be stars; these are the first at sub-10% solar metallicity. We further find that these stars are photometrically contaminating the BHeB region in the optical CMD. We unexpectedly find that ∼5%-10% of stars in the BHeB region above 23.5 mag (F475W) in the optical CMD appear to be on the MS region in the UV CMD. 10. We find the observed sub-SMC metallicity late-type O dwarfs and our O dwarfs to be consistent in T log eff ( ) within 1σ for a given spectral type. We increase the sample of O dwarfs in the sub-10% solar metallicity regime with derived stellar parameters from five to nine. We provide the largest sub-SMC metallicity B dwarfs sample. We also spectroscopically observe the first sub-10% solar metallicity Be stars. We present the largest spectroscopic sample of MS stars at sub-SMC metallicity with derived stellar parameters.
M.G. thanks Miriam Garcia, Chris Evans, Sally Oey, Nate Bastian, Morgan Fouseneau, Dietrich Baade, and the referee for helpful discussions, feedback, and/or comments. M.G. acknowledges support of the "Schweizerische Studienstiftung" and UC Berkeley Cranor Fellowship. Support for this work was provided by NASA through grants GO-15275, GO-15921, GO-16162, GO-16717, AR-15056, AR-16120, and HST-HF2-51457.001-A from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. 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 most fortunate to have the opportunity to conduct observations from this mountain.
Part of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. Further, this data was made accessible by the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration.
Part of the data presented here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution.
This work is based on photometric observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This work made extensive use of NASA's Astrophysics Data System Bibliographic Services.

Appendix A The Effects of Stellar Model Choice on Photometric Fitting Results
We compare the stellar parameters derived from the BEAST using the different stellar models models it has available. They are: (i) MIST models with no rotation (MIST NR ; Choi et al. 2016); (ii) MIST models with rotation (MIST R ; Choi et al. 2017); and (iii) the PARSEC models that do not have rotation (PARSEC; Bressan et al. 2012). We limited this comparison to massive stars that fall within our spectroscopy color and magnitude limits (F475W < 22 and F475W − F814W < 0.5). This selects MS stars with M  5M e .
We compare parameters derived for the same stars from two different models. We divide the photometric SED fits into three groups based on the reduced χ 2 values. Fits with χ 2 < 20 for models being compared are considered good fits. Values of χ 2 80 in either of the evolutionary models are considered poor fits. All other fits are fits of decent quality. These categories are only meant to be used qualitatively to facilitate discussion. When comparing the agreement between different parameters, we only use stars deemed to have good or decent fits. Figure 18 shows a comparison between MIST R and MIST NR in the left column and the between MIST R and PARSEC in the right column. The y-axis is the difference between parameter values for the two models, and the x-axis is the MIST R parameter values.
For T log eff ( ) (top row), the MIST R versus MIST NR , we find an MAD of 0.01 dex and median difference of −0.01 dex. The 16th and 84th percentiles of the difference span 0.04 dex (i.e., ±0.02 dex around the median). This indicates that the differences are smaller than scatter in the data. In general, these findings show that for decent and good-fit stars, the rotational and nonrotational models are generally consistent with each other. We note that the scatter is larger for stars with T log eff ( )>4.0, where we find an MAD of 0.03 dex and scatter of 0.05 dex. As discussed in Choi et al. (2017), at this T log eff ( ), the MIST models change mass-loss prescriptions, which offers the possibility for increased scatter. For MIST R versus PAR-SEC, we find an MAD of 0.01 dex and a median difference of 0.00 ± 0.02 dex. The two models provide consistent values of T log eff ( ). The second row shows g log( ). For the MIST R versus MIST NR , we find an MAD of 0.05 dex and a median of −0.11 ± 0.09 dex. This implies that MIST NR has systematically larger g log( ) values than MIST R . It is possible that rotation may affect g log( ) (e.g., due to line broadening), which could contribute to the disagreement. For the MIST R versus PARSEC, we find an MAD of 0.02 and a median difference of 0.00 dex. The 16th and 84th percentiles of the difference span 0.08 dex (±0.04 dex), indicating that MIST R versus PARSEC are consistent with each other.
The third row shows M ini . For the MIST R versus MIST NR , we find an MAD of 0.46M e and a median of −0.78 ± 0.97M e . The two models are consistent, but MIST NR returns slightly larger values of M ini values for more massive stars. For the MIST R versus PARSEC, we find an MAD of 0.06M e and median difference of 0.05 ± 0.29M e . MIST R and PARSEC are consistent with each other. There is one clear outlier, which appears to be a good fit (i.e., χ 2 < 20) for both model fits.
Interestingly, for all other parameters, the fits are consistent with each other. We currently do not have spectra to further investigate this star.
The fourth row shows  L L log( ). For the MIST R versus MIST NR , we find an MAD of 0.04 dex and median of 0.15 ± 0.14 dex. In this case MIST R returns slightly higher luminosities than MIST NR . For the MIST R versus PARSEC, we find an MAD of 0.02 dex and median difference of −0.01 ± 0.09 dex. Showing that MIST R versus PARSEC are consistent.
The fifth row shows A log( ). For the MIST R versus MIST NR , we find an MAD of 0.10 dex and median of 0.20 ± 0.42 dex. We find the two models are consistent. For the MIST R versus PARSEC, we find an MAD of 0.00 dex and median difference of 0.00 ± 0.06 dex. The MIST R versus PARSEC models also yield values that are consistent.
The bottom row shows A V . For the MIST R versus MIST NR , we find an MAD of 0.04 dex and median of −0.01 ± 0.05 dex. Interestingly, most of the high A V values in either models are also associated with poor fits. The presence of larger A v values in a galaxy with minimal dust may be an indicator that we are not observing "normal" stars. Overall, the A V values from both models are consistent. For the MIST R versus PARSEC, we find an MAD of 0.02 dex and median difference Figure 18. The figure shows a direct comparison between the parameters derived from the BEAST using the different available evolutionary models. The left panels compare the MIST models with v/vcrit = 0.4 (MR) and v/vcrit = 0.0 (MNR). The right panels compare the MIST models with v/vcrit = 0.4 (MR) to the PARSEC models (P). We plot the values for the stars in both our spectroscopic and photometric sample. The top two panels compare T log eff ( ). The second two panels compare the g log( ). The third two panels compare M ini . The fourth panels compare  L L log( ). The fifth panels compare A log( ). The bottom panels compare the A V . Stars with χ 2 < 20 for both fits are pink (good fits), stars with χ 2 > 80 in either fit are cyan (poor fits), and all other stars are purple (decent fits). Detailed Discussion of Individual Stars

B.1. O-Type Stars
K2/B2: K2/B2 is the second brightest star in our sample (F475W VEGA = 20.24 and F336W VEGA = 18.57). Like K1, its locations on the UV and optical CMDs indicate it is likely an O star. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of a typical late-type O-star. We find emission line contamination in the Balmer series, which is not surprising given K2ʼs location in the middle of an H II region. The narrow emission lines are particularly clear for Hα and Hβ and weaker in the Hγ. The 2D spectrum shows that emission is spatially extended, consistent with a nebular origin. In most cases, we can only derive upper limits of EWs for the absorption features. Thus, we can only roughly estimate the spectral type. Following the EW classification process described for K1, we find that K2 is an O9.7 V type star in all of the schemes (Castro et al. 2008;Sota et al. 2011;Martins 2018).
When fitting the full spectrum using The Payne, we mask lines clearly affected by emission (e.g., the Hα feature and the center of the Hβ feature). The spectroscopic analysis yields T log eff ( )= 4.50 ± 0.01 dex and g log( ) = 4.10 ± 0.02 dex, which is consistent with a late O-type MS star. We find a v i sin = 95 ± 20 km s −1 , which is reasonable for a late O-type star. Using Equations (1)-(3), we derive M spec = 27 ± 6M e and M 0 = 17 ± 1M e , and  L L log 4.72 0.09 ( )=  dex. Overall, this agrees well with the EW measurements and the CMD location. We find the SED fitting yields T log eff ( )= 4.50 ± 0.01 dex and g log( ) = 3.85 ± 0.02 dex, which are also consistent with a late O-type star.
K4: K4 is among the brighter stars (F475W VEGA = 20.65 and F336W VEGA = 19.04) in our spectroscopic sample. On the CMD, it is located on the BHeB sequence in the optical, but appears on the MS in the UV-optical CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features. It is not located near a known H II region (Figure 3).
Several of the spectral features suggest a Be-type star: He I, Mg II, Si III, and Si IV are in absorption. The Balmer lines are strong in emission for Hβ and Hγ, and we observe emission features within absorption for Hδ and Hò. Based on the EWs, we classify K4 as Be1 V (Castro et al. 2008;Evans et al. 2015).
K4 only has measurements in four bands (two UV and two optical), since the stars falls outside of the NIR footprint. This makes the composite photometry method not viable, since we cannot constrain the disk. Using the adapted spectroscopy method, we determine that T log 4.45 0.01 eff ,star ( ) =  dex and v i sin = 339 ± 10 km s −1 . This is in line with an early-type Be star. The high v i sin is a key indicator that the star has to be a rotating MS star as opposed to a BHeB star. An evolved star would break apart at such high v i sin (e.g., Rivinius et al. 2013). Furthermore, the scenario is supported by Leščinskaitė et al. (2022), whose photometric analysis shows that K4 is a strong Hα emitter, similar to K6 and K8.
B4: B4 is the third brightest star in the sample (F475W = 20.40 and F336W = 18.75). Its locations in the optical and UV CMDs are typical of an MS star. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectra shows the characteristics of a late O-type/early B-type star: He I, He II, Si II, Si III, and Si IV in absorption and emission. B4 shows weak emission in the wings of Hβ, Hγ, Hδ, and Hò absorption lines. Hα is in absorption with weak narrow emission at the central wavelength. Based on the EW and the schemes by Lennon (1997), Evans et al. (2015), and Castro et al. (2008), we classify this as a B0V type star.
We model the full spectrum of B4, but mask out pixels clearly affected by emission lines, as they are not captured by the TLUSTY models. Our spectroscopic fitting yields T log eff ( )= 4.46 ± 0.01 dex and g log( ) = 3.75 ± 0.04 dex, which is consistent with an early B-type star. Using Equations (1)-(3), we derive M spec = 15 ± 2M e and M 0 = 15 ± 1M e , and  L L log 4.65 0.05 ( )=  dex. Overall, this aligns well with the EW measurements and the CMD location.
The results from SED fitting are consistent with the spectroscopic parameters. The SED fitting yields T log eff ( )= 4.50 ± 0.01 dex and g log( ) = 3.91 ± 0.01 dex. The photometric parameters are also consistent with a late-type O star. However, given the cooler T log eff ( ) from spectroscopy and the EW analysis from the spectra, we conclude the star is likely an early B-type star.
K5: K5 (F475W VEGA = 20.78 and F336W VEGA = 19.19), like the previous stars, is located on the BHeB sequence in the optical, but appears on the MS in the UV-optical CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectral features suggest an early Be-type star: He I is in absorption and emission. Mg II and Si III are in absorption. It is possible that He II is in absorption, but the features are nearly indistinguishable from the continuum. We do not detect any other absorption features for any other elements. This is likely because the fast rotation of the star has washed out these features. For the Balmer lines, we find Hβ in emission, and we note that the star is not near any H II regions. Hγ shows emission within the absorption profile. We further note that the Hδ feature appears to have emission nearly canceling the absorption feature. Based on the presence of Si IV and the He I lines, using the classification scheme of Castro et al. (2008) we determine the star is likely an early Be2 V star. The Be-star scenario is supported by Leščinskaitė et al. (2022), who found that K4 is a strong Hα emitter, similar to K6 and K8.
Strong emission in the Balmer lines and the high reduced χ 2 (>300) in the photometry lead us to adopt the adapted methods. We only have observations in four bands, making the composite photometric fitting not possible, since we cannot constrain the disk properties. Using the adapted spectroscopic method, we determine T log 4.46 0.02 eff ,star ( ) =  dex, and v i sin = 370 ± 90 km s −1 . This is consistent with an earlytype Be star.
B5: B5 (F475W VEGA = 21.5 and F336W VEGA = 20.18) is located on the BHeB sequence in the optical CMD, but appears on the MS in the UV-optical CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of a typical early-tomid-type Be star. He I, Si II, and Mg II are in absorption Hα is in emission, while Hβ shows a strong emission component within its absorption feature. The remaining Balmer lines are all in absorption with weak emission features. We estimate the stellar type of Be3 V using its EWs (Castro et al. 2008).
Balmer line emission and the large BEAST reduced χ 2 (>200) hint at the presence of a disk. Using the adapted spectroscopic method (Section 4.3), we determine T log eff ( ) = 4.40 ± 0.01 dex and v i sin = 165 ± 15 km s −1 . This is in line with an early-to-mid-type Be star. .46) appears to be an O-type star based on its CMD location. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
We observe both He I and He II in absorption. We further measure Mg II in absorption. Most Balmer lines appear in absorption; though the spectral region around Hα does not show a clear Gaussian absorption feature. However, due to the low resolution (R ∼ 1250), it is not clear if this is a data artifact or if there are other physical processes at work (e.g., stellar winds). The EW ratios suggest that K7 is an O9.7 (Sota et al. 2011) or O9.5 (Martins 2018) type star, with a luminosity class of V (Martins 2018). Lastly, using Castro et al. (2008), we determined the star to be an O9.5 V star. This fits well with the star's position in the CMD (see Figure 1).
The median parameters from full spectral fitting for K7 are T log eff ( )= 4.48 ± 0.01 dex and g log( ) = 4.04 ± 0.02 dex ( Table 4). Like the EW ratios, these values indicate the K7 is a late-type O star. The derived g log( ) is consistent with K7 being an MS star. We find a v i sin 97 20 =  km s −1 . Using Equations (1)-(3), we derive M spec = 13 ± 3M e and M 0 = 13 ± 1M e , and  L L log 4.38 0.08 ( )=  dex. These values are in agreement with the classification derived from the EW measurements and its location on the CMD. The spectroscopic and photometric parameters for K7 are in broad agreement, although we observe an elevated g log( ) value similarly to K1. The photometric parameters are T log eff ( ) = 4.52 ± 0.02 dex and g log( ) = 4.23 ± 0.05 dex; both analyses are consistent with this being an O-type star. We find that M 0 and M ini , as well as M spec and M act , are consistent with each other.
B7: B7 (F475W VEGA = 21.35 and F336W VEGA = 20.31) is among the bluest stars in the optical CMD, but in the UV, the star appears to be among the reddest stars. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The star shows emission resembling a P-Cygni profile in each of its Balmer lines. The clearest P-Cygni feature is the Hβ line, since Hα shows a strong emission feature. The central wavelength of Hα falls in the blue wing of the emission feature hinting at a P-Cygni profile. Hγ and Hδ on the other hand show narrow absorption features with very weak emission on the red side. We find He I, Mg II ,Si II, Si III, Si IV, and C I are in absorption (see Table 3). We also find He I in emission, as well as NO III, O II, and C II. Based on the EW ratios and following Castro et al. (2008), we determined the star to be an early B1type/Be1-type star.
The strong emission appears to be associated with the star, as opposed to being gaseous in origin. B7 is not near the known H II regions (see Figure 3). Due to the strong emission in the Balmer lines and the high reduced χ 2 (χ 2 > 300) in the photometry, we adopt here the alternative fitting methods introduced in Section 4.3. B7 only has measurements in four bands (two UV and two optical), since the stars falls outside of the NIR footprint. This makes the composite photometry method not viable, since we cannot constrain the NIR to determine the nature of the star's SED profile.
K8: K8 (F475W VEGA = 21.45 and F336W VEGA = 19.88), like the previous star, is located on the BHeB sequence in the optical CMD, but appears on the MS in the UV-optical CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The Keck spectrum of K8 is shown in Figure 8. The spectral features suggest an early-type fast rotating Be star. We find He I, Si III, and Si IV in absorption and He II in emission. The Balmer lines shows emission in Hβ. Hγ appears to show an emission feature within its absorption feature. Interestingly, the Hδ feature is a very shallow absorption feature, and lastly, we find Hò in absorption. All spectral features are quite shallow, indicating that there is significant rotation. Based on the detected upper limits for both He II and Si, we type the star as an early Be2 V. Similarly to K6, Leščinskaitė et al. (2022) identified K8 as a strong H-α emitter. Due to the stars similarity to K6, both in spectroscopy and photometry, we use the described alternative method of fitting the data (see Section 4.3).
Using our adapted spectroscopic method, we determined that T log eff ( ) = 4.44 ± 0.03 dex and v i sin = 290 ± 10 km s −1 . This is consistent with an early-type Be star. ( ) =  dex, which is in line with the spectroscopic method.
B8: B8 is a faint star (F475W VEGA = 21.5 and F336W VEGA = 20.18) and is the reddest of those we observed with MMT. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum contains several unusual features. We observe only the upper limit for He I in B8, and we find Mg II and Si II in absorption.
Hα and Hβ are in emission and both show a double-peaked line profile. For the other Balmer lines with combined absorption and emission features, we observe a widened wing and narrow absorption within the emission lines. It is likely that emission in the wings is making line profiles shallower. Using the EWs, we classify the star as a Be3 V (Castro et al. 2008). Lastly, the distinct Balmer profiles suggest the star is most likely a shell star (i.e., a Be star viewed edge-on).
The unusual pattern of emission and absorption in the Balmer lines and the high χ 2 (>200) in the BEAST fitting suggest B8 may have a circumstellar disk.
Using the adapted spectroscopic method (Section 4.3), we determine that the star has a T log eff ( ) = 4.42 ± 0.03 dex and v i sin = 135 ± 80 km s −1 . This is in line with an early-tomid-type Be star. Using the composite photometric method, we determine that  dex, which is in line with our previous findings. Interestingly, we would expect the v i sin to be higher given that we are likely looking at the star equator-on based on the doubled peaked profile in Hα. However, if the star is really viewed edge-on, the v i sin value is surprisingly low. However, the uncertainty on our measurement is quite large, suggesting that the star could have a higher v i sin value more in line with expectation (e.g., Iqbal & Keller 2013;Rivinius et al. 2013;Arcos et al. 2018). Leščinskaitė et al. (2022) identified B8 as a strong Hα emitter. Despite the low v i sin , the strong emission and spectral shape supports the classification of the star as a shell star.
B9: B9 is one of our fainter spectroscopic targets (F475W = 21.67, F336W = 20.25). Figure 1 shows that B9 is located between the MS and the BHeB sequences. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of an early B-type star. We find He I, Si II, Si III, and Mg II in absorption. The Hα appears in emission, while all other Balmer lines are in absorption. The Hα is marginally broader than what might be expected from purely gaseous emission, suggesting there may be a stellar contribution. The EW suggests that B9 is a B1.5 V in multiple classification schemes (Lennon 1997;Evans et al. 2015, andCastro et al. 2008).
We fit the full spectrum of B9 using The Payne and find T log eff ( ) = 4.45 ± 0.02 dex and g log( ) = 4.18 ± 0.08 dex, and v i sin = 115 ± 20 km s −1 . Using Equations (1)-(3), we derive M spec = 15 ± 5M e and M 0 = 12 ± 2M e , and  L L log( )= 4.26 0.04  , which are all in line with the rough classification described by the EW measurements and location in the CMD. The BEAST-based photometric parameters are T log eff ( ) = 4.49 ± 0.05 dex, g log( ) = 4.12 ± 0.68 dex, and M ini = 10.87 ± 6.13M e . The parameters are consistent with the spectroscopic measurement.
B10: B10 (F475W VEGA = 21.7 and F336W VEGA = 20.35) is the faintest star in our sample in the UV-optical CMD, but is brighter than several other stars in the optical-only CMD. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of an early B-type star with He I, Si III, and Si IV in absorption. Hα is in strong emission, while Hβ and Hδ show modest emission within their absorption feature. Hγ and Hò are both in absorption. B10 is not near any known H II region, suggesting the emission is not likely due to nebular contamination. Based on the observed EWs, we classify it as a B0.5 V (candidate Be) star.
The photometric parameters yield T log eff ( ) = 4.47 ± 0.03 dex, g log( ) = 4.02 ± 0.07 dex, and M 0 = 12.2 ± 1.21M e . B11: B11 (F475W VEGA = 21.90 and F336W VEGA = 20.25) is the second faintest star in the sample. Both in the opticalonly and the UV-optical sample, it lies on the MS. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of an early B-type star. We find He I, Si II, Si III, and Si IV in absorption. Most Balmer lines are in absorption with variable amounts of emission. For Hα, the strong emission nearly cancels out the absorption feature. Hβ and Hδ show nearly no emission in their absorption features. Although Hγ and Hò are predominantly in absorption, there is some emission in the blue wings. Taking the EWs, including upper limits at face value, we determine that the spectral type is a B2 type star (Castro et al. 2008;Evans et al. 2015).
Since B11 spectrum shows an S/N of <10, we cannot reliably recover the spectroscopic parameters. The photometric parameters are T log eff ( ) = 4.27 ± 0.04 dex, g log( ) = 3.81 ± 0.06 dex, and M ini = 6.26 ± 0.69M e . The photometry is consistent with the star's position in the CMDs as well as the spectral type derived from the EW.
Given that the S/N is ∼10, we refrain from making any further speculations as to the origin of the emission.
B12: B12 is among the faintest stars in our sample (F475W VEGA = 21.75 and F336W VEGA = 20.19). Its locations in the optical and UV CMDs are typical of an MS star. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of an early B-type star. He I, Si III, and Si IV are detected in absorption. The Balmer lines show absorption features with emission for Hα, Hβ, and Hγ. Overall, the emission within the star appears to be predominantly stellar, since the star is not near any know gas region. The emission feature overall appear quite weak in comparison to the other stars with emission in this sample. Lastly, based the He EW ratios, the presence of Si IV, and absence of Mg II, we use the classification scheme of Castro et al. (2008) to conclude that the star is likely a B1.5V star.
We fit the full spectrum of B12 using The Payne. The parameters for the spectral fitting are listed in Table 4. For B12, we find T log eff ( ) = 4.44 ± 0.05 dex and g log( ) = 4.00 ± 0.10 dex. Using Equations (1)-(3), we derive M spec = 8 ± 4M e and M 0 = 11 ± 2M e , and  L L log 4.07 0.05 ( )=  dex, which is all in line with the rough classification described by the EW measurements and location in the CMD. The photometric parameters yield T log eff ( ) = 4.34 ± 0.05 dex, g log( ) = 3.86 ± 0.08 dex, and M ini = 6.61 ± 1.21M e . B13: B13 (F336W VEGA = 20.25 and F475W VEGA = 21.9) is among the faintest stars in the optical-CMD, but is much brighter than several other stars in the UV-optical CMD. Its location on the UV and optical CMDs indicates it is likely a late-type O-type. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of a typical late-type O-star: He I and He II are in absorption. Most Balmer lines appear in absorption with light contamination, except Hα, which is fully in emission. The low S/N (S/N∼ 12) makes it challenging to definitely classify this star. It is likely that B13 is a late O-type MS star. However, it may also be consistent with an early B-type star; we cannot rule this out. Based on the scheme applied thus far for EW, we determined the star is an O9.7-B0 dwarf.
From the full spectral fitting of B13, we find T log eff ( ) = 4.50 ± 0.02 dex, and g log( ) = 4.64 ± 0.06 dex. Using Equations (1)-(3), we derive M spec = 19 ± 4M e and M 0 = 12 ± 1M e , and  L L log 4.07 0.08 ( )=  . M 0 and M ini are consistent, but M spec appears to be too large; once again, this is likely due to the high g log( ) value. The emission in the Balmer lines might affect the g log( ), since the parameter is the most sensitive to it. Nevertheless, these values are in line with the rough classification described by the EW measurements and location in the CMD.
SED fitting yields T log eff ( ) = 4.50 ± 0.01 dex and g log( )=4.37 ± 0.02 dex. Similarly to other O-stars, we find consistency in T log eff ( ) and a slight discrepancy in g log( ). M 0 and M ini are consistent with one another.
Given that the star shows narrow emission features, poor spectroscopic fitting, and extended emission in the 2D spectra, we postulate that the star must have some nebular contamination.
B14: B14 (F475W VEGA = 21.9 and F336W VEGA = 20.40) is the faintest star in our sample in the optical-only CMD, but is brighter than several stars in the UV-optical CMD. In both cases, it lies on the MS. Its spectrum is shown in Figure 8. We list its EW measurements in Table 3 and comment only on a subset of those features.
The spectrum shows the characteristics of an early-to-mid B-type star, with emission in that He I features showing both absorption and emission. We also find Mg II, Si II, and Si III in absorption and [O II] in emission. B14 shows that Hα is primarily present in absorption. Furthermore, we find that Hβ, Hγ, and Hδ all show some potential emission in their wings. However, due to the low S/N (S/N∼ 15) of the spectrum, the weak emission in the Balmer lines could be artifacts of the data. Based on the classification scheme previously employed for EWs, we determine the star is B3 star.
T log eff ( ) is consistent between the two methods, while g log( ) shows a discrepancy.