# Scylla. IV. Intrinsic Stellar Properties and Line-of-sight Dust Extinction Measurements toward 1.5 Million Stars in the SMC and LMC

Christina W Lindberg, Claire E Murray, Petia Yanchulova Merica-Jones, Caroline Bot, Clare Burhenne, Yumi Choi, Christopher J R Clark, Roger E Cohen, Karoline M Gilbert, Steven R Goldman, et al.

► To cite this version:

Christina W Lindberg, Claire E Murray, Petia Yanchulova Merica-Jones, Caroline Bot, Clare Burhenne, et al. Scylla. IV. Intrinsic Stellar Properties and Line-of-sight Dust Extinction Measurements toward 1.5 Million Stars in the SMC and LMC. The Astrophysical Journal, 2025, 982 (1), pp.33. 10.3847/1538-4357/adb4e8. insu-04994428

HAL Id: insu-04994428

https://insu.hal.science/insu-04994428v1

Submitted on 17 Mar 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Scylla. IV. Intrinsic Stellar Properties and Line-of-sight Dust Extinction Measurements toward 1.5 Million Stars in the SMC and LMC

Christina W. Lindberg<sup>1,2</sup> ID, Claire E. Murray<sup>1,2</sup> ID, Petia Yanchulova Merica-Jones<sup>2</sup> ID, Caroline Bot<sup>3</sup> ID, Clare Burhenne<sup>4</sup> ID, Yumi Choi<sup>5</sup> ID, Christopher J. R. Clark<sup>6</sup> ID, Roger E. Cohen<sup>4</sup> ID, Karoline M. Gilbert<sup>1,2</sup> ID, Steven R. Goldman<sup>2</sup> ID, Karl D. Gordon<sup>2</sup> ID, Alec S. Hirschauer<sup>2,7</sup> ID, Kristen B. W. McQuinn<sup>2,4</sup> ID, Julia C. Roman-Duval<sup>2</sup> ID, Karin M. Sandstrom<sup>8</sup> ID, Elizabeth Tarantino<sup>2</sup> ID, and Benjamin F. Williams<sup>9</sup> ID

Received 2024 September 17; revised 2025 February 4; accepted 2025 February 6; published 2025 March 13

## Abstract

By analyzing the spectral energy distributions (SEDs) of resolved stars in nearby galaxies, we can constrain their stellar properties and line-of-sight dust extinction. From the Scylla survey, we obtain ultraviolet to near-infrared photometry from Wide Field Camera 3 on board the Hubble Space Telescope for more than 1.5 million stars in the SMC and LMC. We use the Bayesian Extinction and Stellar Tool (BEAST) to analyze the multiband SEDs of these sources and characterize their initial masses, ages, metallicities, distances, and line-of-sight extinction properties (e.g., $A_V$, $R_V$). We apply quality cuts and perform validation simulations to construct a catalog of over 550,000 stars with high-reliability SED fits, which we use to analyze the stellar content and extinction properties of the SMC and LMC. We detect stars with masses as low as $0.6 M_\odot$. BEAST stellar age distributions show a jump in observed stars around 6 Gyr ago, which agrees with star formation histories. Extinctions ($A_V$) in both galaxies follow a log-normal distribution. We compare $A_V$ with ancillary gas and dust tracers like H I, $H\alpha$, and far-infrared (FIR) dust emission and find positive correlations on a field-by-field basis. We convert observed $A_V$ to predicted dust surface densities using the Draine et al. model and find $A_V$-based dust surface densities are a factor of ~2.5 lower than observed FIR-based dust surface densities, a correction factor similar to other studies.

Unified Astronomy Hirschaurer concepts: Magellanic Clouds (990); Hubble Space Telescope (761); Stellar populations (1622); Interstellar medium (847)

Materials only available in the online version of record: machine-readable table

## 1. Introduction

Studying resolved stellar populations in nearby galaxies provides a powerful window into the processes that drive galaxy formation and evolution. By observing individual stars over a range of wavelengths from ultraviolet (UV) to infrared (IR), we can infer their physical properties such as mass, age, metallicity, and the effects of dust extinction. These measurements enable a diverse range of science, including reconstructing galaxy assembly histories from stellar age-metallicity distributions (e.g., K. B. W. McQuinn et al. 2010); spatiotemporal mapping of star formation patterns (e.g., A. R. Lewis et al. 2015; B. F. Williams et al. 2017); direct constraints on stellar feedback processes (e.g., Y. Choi et al. 2020); mapping of dust properties via individual extinction measurements (e.g., C. W. Lindberg et al. 2024); and detailed studies of star cluster formation, dynamics, and mass functions by resolving individual members (e.g., L. C. Johnson et al. 2017).

Deriving the physical properties of resolved stars has often been done with detailed spectral line information. A high signal-to-noise stellar spectrum can reveal its precise chemical composition, temperature, surface gravity (e.g., Y.-S. Ting et al. 2019), rotational velocity (e.g., W. Huang et al. 2010), and magnetic field (e.g., F. Leone et al. 2017). However, these observations are expensive to obtain, and the required detailed modeling and analysis can be computationally complex. An alternative approach is to model the broadband spectral energy distribution (SED) of resolved stars (e.g., C. Conroy 2013). Although this approach cannot constrain the same level of detail on the atmospheric properties of a star as spectral modeling (e.g., O. G. Telford et al. 2024), broadband SED fitting is consistent with spectral-based fitting (M. Gull et al. 2022) and is observationally more efficient, making it ideally suited for analyzing large stellar surveys (e.g., $10^3-10^6$ sources) of nearby galaxies.

Although numerous SED modeling techniques exist, they generally focus on fitting Galactic sources that suffer minimal crowding compared to resolved star surveys in external galaxies (C. A. L. Bailer-Jones 2011; G. M. Green et al. 2014; R. Schönrich & M. Bergemann 2014; M. Ness et al. 2015).

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ The William H. Miller III Department of Physics & Astronomy, Bloomberg Center for Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA; christina.lindberg@live.com</footnote>

<footnote n="2">² Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="3">³ Observatoire Astronomique de Strasbourg, Université de Strasbourg, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854, USA</footnote>

<footnote n="5">⁵ NSF National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="6">⁶ AURA for the European Space Agency, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="7">⁷ Department of Physics & Engineering Physics, Morgan State University, 1700 East Cold Spring Lane, Baltimore, MD 21251, USA</footnote>

<footnote n="8">⁸ Department of Astronomy & Astrophysics, University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA</footnote>

<footnote n="9">⁹ Department of Astronomy, University of Washington, Box 351580, U. W., Seattle, WA 98195-1580, USA</footnote>
Overcoming these limitations is essential to take full advantage of the rich, but complex data sets provided by current and upcoming resolved stellar population surveys targeting nearby galaxies; e.g., M31: PHAT (J. J. Dalcanton et al. 2012), PHAST (B. F. Williams et al. 2021a), M33: PHATTER (B. F. Williams et al. 2021b), SMC: SMIDGE (P. Yanchulova Merica-Jones et al. 2017), LMC: HTTP (E. Sabbi et al. 2016), METAL (J. Roman-Duval et al. 2019), LUVIT (K. M. Gilbert et al. 2025). This motivated the development of the probabilistic Bayesian Extinction and Stellar Tool (BEAST; K. D. Gordon et al. 2016) to robustly measure stellar and dust properties of resolved stellar populations in nearby galaxies while precisely characterizing parameter uncertainties across all observed signal-to-noise regimes.

The LMC and SMC (collectively as the MCs) are nearby dwarf galaxies at 50 and 60 kpc, respectively (V. Scowcroft et al. 2016; G. Pietrzyński et al. 2019), whose close proximity provides an opportune chance to study resolved stars in external galaxies. These dwarf galaxies are likely on their first infall into the Milky Way (G. Besla et al. 2007; K. Bekki & S. Stanimirović 2009), but have been interacting with each other for several gigayears (E. Patel et al. 2020). Based on stellar orbital modeling, it is thought that these galaxies likely suffered a direct collision as recently as 140–160 Myr ago (Y. Choi et al. 2022; L. R. Cullinane et al. 2022), the effects of which are still kinematically, morphologically (e.g., D. L. Nidever et al. 2008; C. Mastropietro et al. 2009; K. Tsuge et al. 2019), and chemically (e.g., A. J. Fox et al. 2018; J. Roman-Duval et al. 2021) observable today, making them a valuable case study of dwarf-dwarf interactions and enhanced star formation.

The SMC and LMC have experienced relatively less chemical enrichment over their lifetimes compared to the Milky Way, with metallicities of $0.2 Z_⊙$ and $0.5 Z_⊙$, respectively (S. C. Russell & M. A. Dopita 1992). This makes them excellent laboratories for studying the effects of low metallicity on small-scale physics: stellar evolution, e.g., mass-loss rates in winds and ionizing radiation from massive stars (P. A. Crowther 2019; C. Evans et al. 2023; O. G. Telford et al. 2023, 2024; C. Hawcroft et al. 2024; V. Ramachandran et al. 2024; E. Tarantino et al. 2024); asymptotic giant branch (AGB) dust production (M. L. Boyer et al. 2012; S. R. Goldman et al. 2022); and the interstellar medium (ISM), e.g., dust extinction (K. D. Gordon et al. 2003, 2024), composition and dust-to-gas ratio (K. D. Gordon et al. 2014; J. Roman-Duval et al. 2014, 2017, 2021, 2022; C. J. R. Clark et al. 2023), and molecular gas formation (A. D. Bolatto et al. 2011; K. E. Jameson et al. 2016).

With new Hubble Space Telescope (HST) imaging from the Scylla survey (C. Murray et al. 2024), we have multiband photometry for nearly 100 fields across the SMC and LMC, sampling a diverse range of environments, metallicities, gas column densities, and radiation fields. Deep, high-resolution photometry means we can resolve stars down to subsolar masses (E. Sabbi et al. 2013), giving us unprecedented insight into the stellar content of these fields (R. E. Cohen et al. 2024a, 2024b), especially compared to crowding-limited ground-based surveys like the VISTA survey of the MC system (M. Cignoni et al. 2013) and APOGEE DR17 (S. R. Majewski et al. 2017; G. Zasowski et al. 2017; Abdurro'uf et al. 2022), as well as the flux-limited space-based Gaia survey (Gaia Collaboration et al. 2018). By characterizing resolved stars and their line-of-sight extinction in the MCs, we can gain information about the stellar and dust content in these galaxies.

In this work, we present new stellar fits and dust measurements for over 1.5 million stars with multiband SEDs in the MCs. In Section 2, we briefly describe the photometric catalogs used. In Section 3, we describe the BEAST analysis in detail, providing an overview of the stellar, dust, and observation models implemented in the fitting process, describe quality cuts for selecting high-reliability sources.

In Section 4, we use simulations to validate the accuracy of the BEAST as a function of filter coverage. In Section 5, we present stellar and dust measurements for all sources in the SMC and LMC. In Section 6, we discuss how these results compare with existing literature on stellar populations, star formation histories (SFHs), and alternative ISM tracers of the MCs. In Section 7, we summarize our findings.

## 2. Data

The Scylla survey is a pure-parallel program that obtained multiband HST Wide Field Camera 3 (WFC3) imaging of 96 fields in the SMC and LMC during Cycles 27-31 (C. Murray et al. 2024). All fields were observed in at least two optical bands (F475W and F814W) and, if scheduling allowed, additional filters were included in the following priority order: F336W, F275W, F110W, F160W, and F225W. UV filters were prioritized above IR because effects from dust extinction are strongest in the UV and measuring both sides of the Balmer jump (e.g., F336W combined with F475W) helps constrain stellar temperatures. To accurately measure both the stellar type and line-of-sight dust, we therefore only fit sources in fields with coverage in both the optical filters (F475W and F814W) and at least one UV filter, usually F336W, leaving us with 33 fields in the SMC and 24 fields in the LMC, the locations of which are shown in Figure 1 (blue).

In addition to the 57 Scylla fields with optical+UV coverage, we also include 32 fields from the Metal Evolution Transport and Abundance in the LMC (METAL) survey (J. Roman-Duval et al. 2019). METAL is a Cycle 24 HST (2016–2017) survey that obtained spectra toward 33 massive stars in the LMC while simultaneously obtaining parallel imaging using WFC3, similar to Scylla. Out of the 32 parallel WFC3 fields, 25 fields have seven-filter coverage with the same filters as in Scylla described above, one field has six-filter coverage, and six fields have four-filter coverage. The inclusion of these fields more than doubles the number of sources we can characterize in the LMC, the locations of which are shown in Figure 1 (red). Scylla and METAL fields analyzed in this work are summarized in Table 1.

In total, the data sample consists of 33 fields in the SMC and 56 fields in the LMC. Each field covers $160'' \times 160''$ (7.11 arcmin$^2$), except for fields with additional WFC3 IR coverage, in which case a field spans $123'' \times 137''$ (4.65 arcmin$^2$). At the distances of the LMC and SMC (50–60 kpc), the angular coverage of a single field translates to roughly 40–45 pc in physical scales. Calibrated images for both METAL and Scylla were retrieved from the Mikulski Archive for Space Telescopes (MAST) archive and processed using CalWF3 versions 3.3 or 3.4.1 for METAL, and versions 3.5.0, 3.5.1, 3.5.2, 3.6.1, or 3.6.2 for Scylla, based on the date of observations (J. Roman-Duval et al. 2019; C. Murray et al. 2024). For UVIS exposures, both surveys use images that have been corrected for charge transfer efficiency effects. The photometry for both METAL and Scylla was performed using the DOLPHOT package (A. E. Dolphin 2002; A. Dolphin
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Scylla and METAL field locations: maps of the peak brightness temperature of 21 cm emission in the SMC (top; N. M. Pingel et al. 2022) and the LMC (bottom; S. Kim et al. 1999), overlaid with the footprints of the 57 Scylla ($N \ge 3$ filters; blue) and 32 METAL fields (red). In each panel, regions where fields overlap are identified by dashed red bounding boxes and expanded at right. Each field is tagged by its “short” name number where the METAL fields start with “M” (see Table 1). At the distances of the LMC and SMC (50–60 kpc), the angular coverage of a single field (160° × 160°) translates to roughly 40–45 pc in physical scales.</figcaption></figure>

2016), using a similar photometry pipeline as the PHAT survey (B. F. Williams et al. 2014; C. Murray et al. 2024).

For all fields, we only analyze sources that pass the vgst (“very good star”) cuts (C. Murray et al. 2024), which require the following: (1) a source must be measured in all bands of a

field (i.e., if a field was observed in three filters, then a source should have measurements in all three filters, regardless of their quality); (2) a source must have a signal-to-noise ratio (SNR) >4 in at least one observed band; (3) to remove diffraction spikes and background galaxies, the squared sharpness,
<table><caption>**Table 1**
Field Summary</caption><thead><tr><th>Field Name<br>(short)<br>(1)</th><th>Field Name<br>(long)<br>(2)</th><th>Field Center R.A.<br>(°)<br>(3)</th><th>Field Center Decl.<br>(°)<br>(4)</th><th>N<sub>filters</sub><br>(5)</th><th>Guard Flag<br>(6)</th></tr></thead><tbody><tr><td>SMC_5</td><td>15891_SMC-3514se-8584</td><td>16.4523</td><td>-72.835</td><td>3</td><td>0</td></tr><tr><td>SMC_6</td><td>15891_SMC-3956ne-9632</td><td>15.9954</td><td>-72.1305</td><td>3</td><td>0</td></tr><tr><td>SMC_8</td><td>15891_SMC-3955ne-9818</td><td>16.0039</td><td>-72.1341</td><td>3</td><td>0</td></tr><tr><td>SMC_10</td><td>15891_SMC-3149ne-12269</td><td>15.0793</td><td>-72.1515</td><td>7</td><td>1</td></tr><tr><td>SMC_11</td><td>15891_SMC-8743se-11371</td><td>21.4361</td><td>-73.1212</td><td>3</td><td>0</td></tr><tr><td>SMC_13</td><td>15891_SMC-2983ne-12972</td><td>14.9313</td><td>-72.1748</td><td>5</td><td>0</td></tr><tr><td>SMC_14</td><td>15891_SMC-3669ne-13972</td><td>14.9314</td><td>-71.9391</td><td>3</td><td>1</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form in the online article.)

roundness, and crowding of a source must be equal to or less than 0.15, 0.6, and 0.2, respectively, in optical bands F475W and F814W; and (4) a source cannot be flagged as anything other than FLAG = 0 or FLAG = 2 in DOLPHOT, where 0 indicates a “good” PSF extraction, and 2 indicates too many bad or saturated pixels. We include sources with FLAG = 2 because a significant fraction of sources in certain fields are associated with bright stars of interest. We apply these same cuts to the METAL data for consistency. In total, we fit 623,763 sources in the SMC and 949,831 sources in the LMC, roughly 30% and 25% of the number of sources in the original catalogs, respectively.

## 2.1. Example Fields

Sources cover a wide range of colors and magnitudes. In Figure 2, we plot the optical color–magnitude diagram (CMD) of a seven-filter METAL field in the LMC, (14675_LMC-5536ne-9786 or M27¹⁰, top) and a three-filter Scylla field in the SMC (15891_SMC-3669ne-13972 or SMC_14, bottom). The locations of these fields were selected to represent the average Scylla and METAL fields, which will be preferentially located in the vicinity of star-forming regions (e.g., 30 Doradus for M27). Even though additional filters are available, both CMDs are plotted using F475W and F814W, since these (and F336W) are the default filters available in every Scylla and METAL field. We gray out all sources below the 50% completeness limits in either filter (F475W = 27.8 mag, F814W = 25.5 mag) to visualize the completeness limits calculated in C. Murray et al. (2024).

From the CMDs, we see there are brighter stars (F814W < 17 mag) in the LMC field compared to the SMC field. The lack of bright sources is the result of stars saturating during long exposures. For the majority of Scylla fields, we were only able to obtain standard long exposures due to scheduling restrictions as a parallel program. This means we were unable to obtain short “guard” exposures like METAL did, which would have allowed us to extract the photometry of bright sources. Consequently, we expect to observe few bright stars in the 56% of Scylla fields missing guard exposures.

Due to the large field of view of the MCs, we anticipate a small fraction of the sources in our fields to be contamination

from low-mass stars in the Milky Way. For both fields, we plot the CMD location of potential foreground sources from TRILEGAL as faint red stars (M. A. T. Groenewegen et al. 2002; L. Girardi et al. 2005; L. Girardi & P. Marigo 2007; E. Vanhollebeke et al. 2009). Most foreground sources are both redder and brighter than the main sequence in each field; however, there is some overlap with red clump stars. Foreground sources will be poorly modeled at a distance of 50–60 kpc, so we will remove these sources with quality cuts based on modeling agreement (Section 3.3).

In both fields, we observe a variety of CMD features such as a prominent main sequence, a track of red giant branch (RGB) stars, red clump (RC) stars, and pre-main-sequence stars. In the LMC METAL field, we also observe even brighter RGB stars overlapping with a potential AGB and red supergiant (RSG) stars. This track is likely absent in the SMC field due to the brighter saturation limit.

We overlay PARSEC isochrones to help us contextualize the range of masses and ages present in the fields. The isochrones span stellar ages of 10 Myr (dotted) to 10 Gyr (solid) (A. Bressan et al. 2012; Y. Chen et al. 2014, 2015; J. Tang et al. 2014; P. Marigo et al. 2017), assuming a distance of 60 kpc and 50 kpc and metallicities of 20% and 50% Z$_{⊙}$ for the SMC and LMC, respectively. In the SMC field, a 20% Z$_{⊙}$ isochrone is well matched with the observed upper main sequence, but is a poor fit to the RGB and red clump. Since RGB stars are several gigayears older than upper-main-sequence stars, we expect RGB stars to be more metal-poor (e.g., 10% Z$_{⊙}$). Previous Scylla papers reported similarly low metallicities for RGB stars (e.g., Figure 2 in R. E. Cohen et al. 2024b). RGB stars indicate the presence of older stellar populations; however, we also observe several sources in both fields that are likely younger than 100 Myr, implying both galaxies have ongoing star formation, which is consistent with recent SFH results (A. Mazzi et al. 2021; O. C. Jones et al. 2023; O. Nayak et al. 2024). For the 10 Myr isochrone (dotted), we note the CMD location of stars with initial masses ranging from 0.5–8 M$_{⊙}$. In the LMC field, there are a handful of sources that could potentially be more massive than 8 M$_{⊙}$. Based on the completeness limits, we should be able to detect stars with initial masses lower than 1.0 M$_{⊙}$.

In the LMC field, we apply 0.5 mag of extinction to the isochrones to improve alignment with the observations. In general, we expect the LMC to contain more dust than the SMC, as evidenced by the width of the RGB and red clump. Even with this adjustment, the 10 Gyr isochrone is still bluer than a significant fraction of sources in the RGB (~30%), a potential indication that these stars have more dust along

<footnote n="10">¹⁰ Long field names are formatted as “[PID] [galaxy] – [distance from galaxy center in arcsec] [cardinal direction] – [telescope rotation angle in deg (V3PA)],” where PID is the HST program ID under which the observation was taken (e.g., 14675_LMC-5945ne-34377). PID = 14675 denotes METAL fields, whereas all higher PIDs denote Scylla fields.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example CMDs: optical CMDs for one three-filter METAL field in the SMC (top; 15891_SMC-3669ne-139720 or SMC_14) and one seven-filter Scylla field in the LMC (bottom; 14675_LMC-5536ne-9786 or M27). All sources below the 50% completeness limits in either filter are shown in gray (F475W = 27.8 mag, F814W = 25.5 mag; C. Murray et al. 2024). Brighter sources (F814W < 17 mag) are absent in the SMC field due to a lack of guard exposures that capture brighter sources before they saturate, independent of the number of filters observed. PARSEC isochrones are overlaid, ranging from 10 Myr to 10 Gyr. For the SMC field, we apply no extinction, assume a distance of 60 kpc, and $Z = 0.2 Z_⊙$, while for the LMC field, we apply 0.5 mag of extinction, assume a distance of 50 kpc, and $Z = 0.5 Z_⊙$. In both fields, we show the direction of a 1 mag reddening vector (black arrow). The location of potential foreground sources from TRILEGAL are plotted as faint red stars. Example SEDS: inset plots (top right) show SEDs for five random sources in each of three distinct regions of the CMD: the upper main sequence (blue), the red clump (magenta), and the lower main sequence (yellow). Our ability to accurately model stars and their extinction is limited by the number of bands observed (i.e., more bands provide greater leverage).</figcaption></figure>

their line of sight, as initially shown in Figure 15 of C. Murray et al. (2024). Compared to the total internal extinction, we anticipate limited amounts of foreground extinction for both

fields with an estimated $A_V = 0.1–0.2$ mag for the SMC and $A_V = 0.2–0.5$ mag in the LMC (P. B. W. Schwering 1988; L. Staveley-Smith et al. 1997; B. Q. Chen et al. 2022).
<table><caption>**Table 2**
Beast Model Parameters</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Description</th><th>Min</th><th>Max</th><th>Resolution</th><th>Prior</th></tr></thead><tbody><tr><td>log(A)</td><td>years</td><td>stellar age</td><td>6.0</td><td>10.13</td><td>0.1</td><td>flat SFR</td></tr><tr><td>log(M<sub>ini</sub>)</td><td>M<sub>☉</sub></td><td>initial stellar mass</td><td>-0.8</td><td>2.0</td><td>variable<sup>a</sup></td><td>Kroupa IMF</td></tr><tr><td>log(Z)</td><td>Z<sub>☉</sub></td><td>stellar metallicity</td><td>-2.1</td><td>-0.3</td><td>0.3</td><td>flat</td></tr><tr><td>A<sub>V</sub></td><td>mag</td><td>dust column</td><td>0.01</td><td>10.0</td><td>0.05</td><td>flat</td></tr><tr><td>R<sub>V</sub></td><td>...</td><td>dust average grain size</td><td>2.0</td><td>6.0</td><td>0.5</td><td>peaked at ~2.7<sup>b</sup></td></tr><tr><td>f<sub>A</sub></td><td>...</td><td>dust mixture coefficient</td><td>0.0</td><td>1.0</td><td>0.1</td><td>peaked at 1<sup>b</sup></td></tr><tr><td>d [LMC]</td><td>kpc</td><td>distance</td><td>40</td><td>60</td><td>2.5</td><td>flat</td></tr><tr><td>d [SMC]</td><td>kpc</td><td>distance</td><td>47</td><td>77</td><td>2.5</td><td>flat</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Initial stellar mass resolution is variable and depends on sampling of stellar evolutionary models.

<sup>b</sup> R<sub>V</sub> and f<sub>A</sub> priors are set by 2D prior, described in Figure 7 of K. D. Gordon et al. (2016). Weighted average of the 2D prior is R<sub>V</sub> = 3.55 and f<sub>A</sub> = 0.73.

Although both fields have similar optical CMDs, our ability to characterize individual SEDs is highly dependent on the number and range of filters observed. In the inset plots (Figure 2, top right), we show the SEDs for five random sources in three separate regions: (1) the upper main sequence (blue); the RGB/RC (magenta); and the lower main sequence (yellow). The CMD locations of these sources are also indicated in the main panel. For the seven-filter LMC field, we see that the upper-main-sequence stars (blue) peak in the UV, while the RGB/RC stars (magenta) are also bright but peak in the optical bands and drop off sharply in the UV. Meanwhile, lower-main-sequence stars (yellow) are several orders of magnitude dimmer in all bands and have much larger photometric uncertainties (shaded). For the three-filter SEDs in the SMC field, we observe similar fluxes; however, the overall shapes of the SEDs are much less apparent. Therefore, we anticipate that SED parameter fits for sources in fields with fewer filters (e.g., three bands) will be less constrained than fields with more filters (e.g., four or more bands). We investigate the effects of filter number on BEAST constraints in Section 4.

Going forward, we will refer to the Scylla and METAL fields analyzed in this paper collectively as Scylla fields for simplicity.

**3. Analysis**

We use the BEAST (K. D. Gordon et al. 2016) to analyze the Scylla photometric catalogs. The BEAST uses a Bayesian probabilistic approach to model the intrinsic properties of individual stars and dust along each sight line based on their multiwavelength SEDs. While the BEAST is adept at analyzing individual stars, it is optimized for large-scale photometric surveys of resolved stellar populations, where factors like crowding and photon noise significantly influence flux measurements and photometric biases across filters. The BEAST is also designed to incorporate priors to accommodate previous knowledge about the expected stellar and dust properties from relevant studies (e.g., mass distribution/initial mass function, IMF, age distribution/SFH, distance ranges). The BEAST fitting procedure is explained in detail in K. D. Gordon et al. (2016), where the BEAST was used to interpret observations from the Panchromatic Hubble Andromeda Treasury (J. J. Dalcanton et al. 2012). Since then, the BEAST has been applied to other galaxies in the Local Volume and is established as a reliable SED fitting tool (Y. Choi et al.

2020; D. Van De Putte et al. 2020; P. Yanchulova Merica-Jones et al. 2025).

In Section 3.1, we provide a brief summary of the BEAST grid construction and the stellar, dust, and observational models employed for version 1.0 of the Scylla fitting. Details regarding the range, resolution, and priors for each parameter in the grid are summarized in Table 2. In Section 3.2, we describe the summary statistics the BEAST provides from its fitting procedure and present an example fit toward two stars: one with noticeable extinction and one without (Figure 3). In Sections 3.3, we present quality cuts to isolate high-reliability sources.

**3.1. Predicted SEDs**

To characterize individual stars, the BEAST compares the observed SED of each star with predicted SEDs. These predicted SEDs are constructed using a range of inputs from the stellar and dust models with corresponding priors, producing a 7D grid (one for each parameter; see Table 2) of extinguished SEDs to which we can compare the data. To construct this grid, first, the BEAST generates a grid of stellar spectra across a range of initial masses, ages, and metallicities (details in Section 3.1.1). Then, varying amounts of extinction are applied to each stellar spectrum based on the dust extinction models (defined by three dust parameters detailed in Section 3.1.2), producing an even larger grid of extinguished spectra. The grid is then duplicated at a range of distances, creating a 7D parameter grid of extinguished stellar spectra. We then integrate these spectra over HST transmission curves for each observed filter to generate model photometric measurements (i.e., simulated SEDs) for every combination of input stellar and dust parameters.

Lastly, to account for systemic observational effects in the actual photometry, we use artificial star tests to quantify the average photometric bias as a function of flux in each filter (details in Section 3.1.3). These biases are added to the predicted SEDs so we can accurately compare with the observed SEDs.

We note that the number of data points does not restrict the number of model parameters that can be fit since the BEAST is not fitting SEDs with a set of basis functions that can describe any shape (e.g., polynomials). Rather, we are fitting with a limited set of shapes (i.e., dust extinguished stellar SEDs), which means that even with seven stellar and dust parameters, there are SED shapes that cannot be fit well (e.g., unresolved binaries, stars with complex atmosphere models like AGB
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Example SED fits from the BEAST: SED fits of two stars in an LMC field (16786_LMC-5127ne-3118)—a massive star with moderate amounts of extinction (blue) and a red giant branch (RGB) star with almost no extinction (red). The observed fluxes are shown with black dots. The best-fit unextinguished stellar SED (dotted), extinguished SED (dashed), and bias-corrected SEDs (solid) are shown for comparison (dashed and solid lines are almost entirely overlapping). Model SEDs are constructed using the median (±1σ) values of each marginalized stellar and dust parameter. In the panels below, we show the distributions of these 1D marginalized parameters, noting the median ±1σ values. For both stars, there are several combinations of stellar and dust parameter values that can produce the observed flux. Primary parameters such as initial mass ($M_{ini}$), age, and extinction ($A_V$) are more tightly constrained compared to secondary parameters like metallicity ($Z$), distance, total-to-selective extinction ($R_V$), and extinction curve shape ($f_A$), which have a weaker effect on the observed flux.</figcaption></figure>

stars, nonstellar emission sources). Additionally, certain parameters (e.g., mass, age, extinction) will have a much larger influence on the resulting SED shapes than others (e.g., stellar metallicity, dust composition).

### 3.1.1. Stellar Model

The stellar model SED grid is constructed by integrating stellar evolutionary and atmospheric models. The primary stellar parameters consist of the initial stellar mass ($M_{ini}$), age ($A$), and metallicity ($Z$), from which we can derive stellar parameters such as the stellar effective temperature ($T_{eff}$), surface gravity ($g$), and luminosity ($L$). For the stellar evolutionary models, we use PARSEC stellar evolutionary tracks sourced from the CMD 3.7 web interface<fnref n="11" /> of A. Bressan et al. (2012), assuming a P. Kroupa (2001) IMF,

and encompassing a metallicity range approximately spanning from 0.008–0.5 $Z_\odot$ (J.-K. Lee et al. 2005). We use stellar atmosphere grids developed by F. Castelli & R. L. Kurucz (2003) and T. Lanz & I. Hubeny (2003, 2007) for local thermodynamic equilibrium (LTE) models and non-LTE models, respectively. Where the two grids overlapped, we preferred the T. Lanz & I. Hubeny (2003, 2007) models (see K. D. Gordon et al. 2016).

It is important to acknowledge that the BEAST fitting process has recognized limitations as it presently lacks models for certain stellar types including binaries, AGB stars, young stellar objects, and other rare stellar species. These stellar types require complex modeling and are therefore not usually included in common stellar evolutionary models. These omissions will impact the accuracy of fits for any unrepresented stellar types. For example, AGB stars are often characterized as highly extinguished RSGs since these two stellar types are relatively cop spatial in a CMD (C. W. Lindberg et al. 2024). In

<footnote n="11">¹¹ http://stev.oapd.inaf.it/cgi-bin/cmd</footnote>
addition, unresolved binaries will often appear redder and/or more luminous than a singular primary star due to the presence of a lower-mass secondary. However, without binary models, the BEAST will attempt to model these stars by applying increasing amounts of dust, resulting in an overestimation of extinction (e.g., Section 6.1.2). As a check, we can place upper limits on the expected extinction by observing the maximum extinction of brighter and older stellar populations, such as RGB stars (J. J. Dalcanton et al. 2023), which generally have a scale height larger than star-forming thin disks of galaxies (N. N. Patra 2019; X. Meng & O. Y. Gnedin 2021). For most of these unrepresented stellar types, there are often few SED models that can reproduce the observed magnitudes and colors, meaning that uncertainties on parameter fits will be unusually small relative to other common stellar types. This quirk gives us an opportunity to remove misfit stellar types by analyzing the uncertainties of the BEAST results, rather than just relying on the photometry.

### 3.1.2. Dust Model

One of the main science objectives of the BEAST is to quantify not only the amount of dust, but also how the nature of dust grains varies within galaxies. For most sight lines in the Milky Way (MW), variations in the shape of the extinction curve can be parameterized as a function of $R_V$, the ratio of total extinction, $A_V$, to selective extinction, $E(B-V)$. $R_V$ reflects the grain size distribution (J. A. Cardelli et al. 1989; E. L. Fitzpatrick 1999; K. D. Gordon et al. 2023). However, the majority of extinction curves measured toward the MCs do not display this dependence on $R_V$ (K. D. Gordon et al. 2003, 2024). Instead, notable distinctions are observed in the measured SMC extinction curves, particularly in the UV part of the spectrum, where the “SMC Average” curve exhibits the absence of a 2175 Å bump and a steeper rise in the far-UV. These variations are often attributed to the lack of metals in the MCs, but the large extinction curve variations within the MCs indicate other effects are important and may even be dominant (K. A. Misselt et al. 1999; K. D. Gordon et al. 2003, 2024). Regardless, extinction curves within the Local Group exhibit a wide range of properties falling between those of the MW and the SMC extinction curves (K. D. Gordon et al. 2003, 2024; G. C. Clayton et al. 2015; E. L. Fitzpatrick et al. 2019).

To address this variability, K. D. Gordon et al. (2016) proposed a mixture model for dust extinction where the coefficient $f_A$ represents the proportion of MW-type dust extinction and $1-f_A$ represents the proportion of SMC-type dust extinction ($A_\lambda$ signifies the extinction at wavelength $\lambda$):

$$ \frac{A_\lambda}{A_V} = f_A \left[ \frac{A_\lambda}{A_V} \right]_{\text{MW}} + (1 - f_A) \left[ \frac{A_\lambda}{A_V} \right]_{\text{SMC}} \quad (1) $$

where $[A_\lambda/A_V]_{\text{MW}}$ and $[A_\lambda/A_V]_{\text{SMC}}$ are the average MW and SMC curves. The mixture model smoothly interpolates between SMC-type dust ($R_V=2.74$; $f_A \sim 0$; K. D. Gordon et al. 2003) and MW-type diffuse ISM dust ($R_V=3.1$; $f_A \sim 1$; updated to E. L. Fitzpatrick et al. 2019). Future work will use the updated SMC Average curve (K. D. Gordon et al. 2024).

We analyze the Scylla data with this flexible dust model using three parameters: $A_V$, the total extinction in the V band; $R_V$, the total-to-selective extinction; and $f_A$, a mixture coefficient indicating the proportion of MW- and SMC Bar-type dust (e.g., range of 2175 Å bump strengths; K. D. Gordon

et al. 2016). BEAST extinction measurements will consist of both internal galactic extinction from the SMC/LMC and foreground extinction from the Milky Way because the BEAST uses individual stars as sight lines. Due to this foreground extinction, we expect the minimum potential extinction for extragalactic sources to be greater than zero. Other dust measuring techniques, such as quantifying the reddening vector of RGB stars on a CMD, only measure the internal extinction of a galaxy (J. J. Dalcanton et al. 2015; P. Yanchulova Merica-Jones et al. 2017).

Most dust column densities tend to exhibit power-law or log-normal distributions (J. J. Dalcanton et al. 2015; M. Lombardi et al. 2015; N. Schneider et al. 2015), with a tail of higher densities. Given the vast range of star-forming environments observed across the Scylla fields and the fact that we are probing the dust column densities with discrete sight lines, we want to ensure that we do not unintentionally bias BEAST fits against higher extinction. As such, we choose to adopt a flat prior for $A_V$. Meanwhile, for $R_V$ and $f_A$, we adopt the same 2D prior as described in K. D. Gordon et al. (2016), which has 1D projections that peak at ~3 and 1, respectively.

### 3.1.3. Observational Model

One of the main issues when comparing the observed photometry of stars to modeled stellar photometry is that we do not know how observing conditions might influence observations. Factors like extended background emission and nearby sources (e.g., binaries, background galaxies, other stars) can all bias what fluxes we observe for an individual source, or whether we even detect a source at all. In particular, stellar crowding due to neighboring sources can result in contaminated flux measurements, causing independent flux offsets in each filter. These effects can be significant in fields with high source densities, such as the center of M31, which are considered crowding-limited in HST data when stellar densities exceed 10 sources/arcsec² (see Figure 33 in B. F. Williams et al. 2014). We anticipate the impact of stellar crowding in Scylla fields will be minimal due to the relative proximity of the MCs—with most Scylla fields having stellar densities <1 sources/arcsec² (C. Murray et al. 2024). However, quantifying the observational effects in each Scylla field, it is still necessary if we want to compare observations with model SEDs.

To address these issues, we develop an observational model to simulate observational effects in the data. Specifically, the observational model quantifies the average bias between model SEDs and their simulated observations and average observational uncertainty across a range of flux bins for each filter. To simulate observations, we employ artificial star tests (ASTs) using the following system: (1) photometry is simulated for a large number of stars across a range of fluxes in each observed filter (usually several thousand); (2) for each star, its intrinsic flux is converted to electron counts that are added to the original observations of a field at a random location within the images; and (3) observations are re-reduced using the same photometry pipeline as the original observations. Using this process, we can determine whether a star that was injected into the image was recovered or not, helping us quantify what fraction of stars we might be missing at each magnitude (i.e., how “complete” our photometry is at increasing magnitudes). If a star is observed, we can then measure the photometric bias, i.e., how much its observed flux changed compared to its
injected flux, and the photometric uncertainty, i.e., the standard
deviation of the measured fluxes. Using this information, we
calculate the average bias and uncertainty of sources in each
filter as a function of the observed flux. A summary of the
observational model can be found in C. Murray et al. (2024).

We require ASTs to cover a wide range of fluxes to ensure
that we are able to quantify observational effects for any
realistic stellar model (realistic in this case means having at
least one band within the range of observed fluxes). To do this,
we estimated the total flux range in each observed band and
split that range into 40 flux bins. Model SEDs are then
randomly selected from the BEAST-predicted SEDs until there
are at least 50 model SEDs in each bin, resulting in a few
thousand SEDs per field. By selecting stars in this manner, we
assume bands are independent. Realistically though, we do
expect correlation between observed bands; however, quantify-
ing these correlations would require far more ASTs. In other
words, fewer samples are needed to cover *N* 1D spaces, where
*N* is the number of observed bands, as opposed to one *N*-
dimensional space. By assuming the bands are independent, we
thus minimize the number of ASTs needed to construct the
observational model.

3.2. Fitting SEDs

For each star, we quantify the likelihood of each model by
directly comparing the synthetic photometry with the observed
photometry and use the photometric uncertainties from the
observational model as the uncertainties in the calculations,
producing a grid of likelihoods. This likelihood grid is
multiplied by a similarly sized grid of priors, providing 7D
posterior probability distribution functions (pPDFs) for all
relevant models. We marginalize this 7D pPDF to obtain 1D
pPDFs for each input parameter.

The BEAST reports a sparse sampling of the full 7D pPDFs,
the 1D marginalized pPDFs themselves, and 16th, 50th, and
84th percentile values from the 1D pPDFs (i.e., median ±1σ)
for each parameter. The BEAST also reports the expectation
value (Exp)—the weighted average of the posterior, and the
parameter values for the single best model fit from the whole
grid (Best) along with the corresponding χ² of said Best fit.

In Figure 3, we show example SED fits for two stars in the
catalog, one with significant extinction and one without. In the
main upper panels, the black dots show the observed flux of
each star across the six observed filters, and overlaid in shaded
colors are the simulated SEDs based on the median ±1σ fits of
the marginalized parameters. In the panels below, we plot the
marginalized 1D pPDFs for each BEAST parameter, and note
the 16th, 50th, and 84th percentile fits based on the 1D pPDFs.
Going forward, all parameter measurements (X) will be based
on the median values of their 1D pPDFs (Xp50), and parameter
uncertainties will be based on the 1σ deviations defined as:
(Xp84 - Xp16)/2.

From the marginalized 1D pPDFs, we see there are several
parameters that exhibit double-peaked probability distributions,
specifically M<sub>ini</sub>, age, and A<sub>V</sub>. In these scenarios, there are
multiple combinations of parameter values that could produce
SEDs similar to the observations. Bimodalities in the parameter
fits arise when there are observational degeneracies between
luminosity, distance, and extinction. However, for most 1D
pPDFs, there is generally only one major mode that is often
aligned with the median, and even if there are multiple modes

of similar probability (e.g., A<sub>V</sub> for the red star), these modes
tend to be clustered in the same region of the parameter range.

Not all Scylla fields were observed with the same number of
filters. In Figure 4, we show how BEAST fits for a single star
varies when observed with seven near-UV (NUV)–near-IR
(NIR) filters (blue) versus three optical filters (red). While both
extinguished stellar model SEDs agreed with the observed
photometry (black dots), in the parameter subfigures, we see
notable differences in the marginalized 1D pPDFs. Primary
parameters, like initial mass (M<sub>ini</sub>), age, and extinction (A<sub>V</sub>),
were all recovered in a similar range, although the seven-filter
BEAST fits have smaller uncertainties. Meanwhile, secondary
parameters like metallicity (Z), distance, R<sub>V</sub>, and f<sub>A</sub> all span a
wide range of values and are generally in poor agreement. In
Section 4, we investigate this effect of filter coverage in more
detail.

3.2.1. Fitting Procedure

To perform the fitting, we made use of high-performance
computing (HPC) resources of the Advanced Cyberinfrastruc-
ture Coordination Ecosystem: Services & Support super-
computer program supported by the National Science
Foundation. The resources were allocated on the Bridges-2
HPC platform administered by the Pittsburgh Supercomputer
Center (PSC). The total resources used amounted to 4.5 million
service units (SUs) on Bridges-2 Regular Memory (RM) nodes,
and 220,000 SUs on the Bridges-2 Extreme Memory (EM)
nodes. EM nodes were used for tasks that required a large
amount of memory (e.g., generating BEAST model grids and
observational models, and trimming the BEAST model grids),
while the RM nodes were used for general-purpose computing
such as the fitting of each individual source. These resources
were also used for overhead tasks, such as code benchmarking
and occasional re-fitting.

Overall, we fit 949,831 stars in the LMC and 623,763 stars
in the SMC, located across 56 and 33 individual fields in the
LMC and SMC, respectively. To produce the final Scylla
catalogs presented in this paper, we trim the BEAST parameter
files consisting of fits for each individual star by post-
processing them with the quality cuts described in the next
Section.

3.3. High-reliability Catalog Quality Cuts

Quality cuts in the Scylla photometric catalog were less
restrictive compared to other photometric surveys (e.g.,
J. J. Dalcanton et al. 2012; B. F. Williams et al. 2014), with
only single-band cuts on SNR and no restrictions on crowding
(for more details, see Section 2). This was done to ensure that
as broad a range of potential sources were characterized by the
BEAST and that no rare or interesting objects were excluded.
With BEAST fits for all of these sources, we now define
additional quality cuts to establish a high-reliability catalog
(HRC) flag to distinguish between well-fit sources and poorly
fit sources/potential contaminants.

BEAST likelihoods were calculated using all available
wavelength information in each field. However, some Scylla
fields fit with the BEAST only have three-filter coverage, so we
want to ensure that any quality cuts can be applied across all
fields equally. Between the optical and the UV filters shared in
all fields (F814W, F475W, and F336W), the UV is the best
band for determining whether extinguished models are actually
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** High-reliability quality cuts. Left panel: observed F336W flux vs. the BEAST-modeled F336W flux of all sources in the SMC and LMC, colored by their signal-to-bias ratio (Δ<sub>F336W</sub>). Right panel: same as left panel, but sources are colored by their photometric signal-to-noise ratio (SNR) in F336W. For visual clarity, we limit the color bars to 3.0 for both plots. We exclude any sources with Δ<sub>F336W</sub> < 3 or SNR<sub>F336W</sub> < 3 from the high-reliability catalogs.</figcaption></figure>

dozen) per field with abnormally high $A_V$ estimates and low $A_V$ uncertainties spatially coinciding with saturated sources in the images. These sources have $A_V$ SNRs (SNR($A_V$) = 2 × $A_{Vp50}/(A_{Vp84} - A_{Vp16})$) much higher than all other sources in the catalogs: diffraction spike sources have $A_V$ SNRs from 10–20, while all other catalog sources have $A_V$ SNRs 1–5. Removing the top 1% of high $A_V$ SNRs sources removes all potential contamination from diffraction spikes, except for in one field (16786_SMC-4451ne-16362) where we opt to place the $A_V$ SNR manually at 10. We remove the top 1% of high $A_V$ SNR sources for each field after placing the previous two quality cuts. By isolating potential diffraction spikes via an SNR cut rather than a cut on $A_V$ directly, we are able to retain the majority of high $A_V$ sources in each field.

Cumulatively across all fields in the SMC and LMC, these quality cuts remove roughly 59% and 68% of sources from the BEAST catalogs. The final high-reliability catalogs contain 254,468 and 306,290 in the SMC and LMC, respectively. In the complete BEAST catalogs, we include a column called [HRC] to identify whether a source passes [1] or fails [0] these quality cuts.

In Figures 6 and 7, we show CMDs in the UV (left), optical (middle), and IR (right) of sources in the SMC and LMC, respectively, that pass (top row) and fail (lower rows) each of the HRC cuts. Due to the magnitude overlap between these populations, we plot each population separately. In the upper row, we show the density distribution of sources that pass the HRC cuts compared to all of the sources that failed (gray). Inversely, in the lower rows, we show the density distribution of sources that failed each specific HRC cut (F336W bias, F336W SNR, $A_V$ SNR), compared with sources that passed. We apply each cut independently to the original catalog, meaning the combined fractions of excluded sources should not add up to 100%.

The F336W bias and F336W SNR cuts both remove a similar fraction of sources (52.6% and 52.8%, respectively) that are generally dimmer and span a broader range of colors. These cuts remove poorly modeled young stellar objects to the left of the lower main sequence (D. A. Gouliermis et al. 2012), and dim red background galaxies. In general, these two cuts remove sources below the average 50% completeness limit of $m_{F275W}$ = 26.2, $m_{F336W}$ = 26.9, $m_{F475W}$ = 27.8, $m_{F814W}$ = 25.5, $m_{F110W}$ = 24.7, and $m_{F160W}$ = 24.0 (dotted lines). Conversely, the $A_V$ SNR cut

only removes 1.8% of sources, however, these sources are generally much brighter and are likely either foreground stars blueward of the optical main sequence (as shown in Figure 2) or sources associated with bright diffraction spikes.

## 4. Validation

Although a significant fraction of fields in this study were observed with six to seven filters from the NUV-NIR, there are still many fields with only three to four filters of coverage in the UV and optical. To test the accuracy of the BEAST fitting method as a function of filter coverage, we construct simple simulations. First, we generate a grid of BEAST models using the same parameters as the BEAST model grid for the LMC (Table 2), with coarser resolution in all primary variables [log(A) = 1.0 dex, log(Z) = 0.6 Z<sub>☉</sub>, $A_V$ = 0.2 mag, $R_V$ = 1.0, $f_A$ = 0.25, d = 10 kpc] and sample 10<sup>4</sup> stars from this grid. Next, using the ASTs from field 14675_LMC-15562nw-33974 (M<sub>7</sub>), which was observed in all seven filters, we generate an observation model to estimate photometric biases and uncertainties, and fit the simulated stars using the full-resolution BEAST model grid for the LMC. We repeat this procedure for nearly all combinations of Scylla filter coverage, summarized in Table 3. After fitting, we apply the same high-reliability catalog cuts as described in Section 3.3 to the simulated catalogs, leaving us with several thousand high-reliability sources for analysis in each case.

Although the parameter range of the simulated stars is limited, these simulations give us an idea of what BEAST parameters we can reliably recover within the range of available Scylla and METAL filters. In Table 4, we quantify the accuracy of each BEAST parameter for each filter combination by calculating the rms error for all high-reliability sources. We compute the rms error for each BEAST parameter (x) as follows:

$$ \text{rms}(x) = \sqrt{\sum_{i=1}^{n_{\text{HRC}}} \frac{(x_{i, \text{out}} - x_{i, \text{in}})^2}{n_{\text{HRC}}} } $$

where $n_{\text{HRC}}$ is the number of sources that pass the high-reliability cuts out of the original 10,000 stars simulated (Table 4; rightmost column). In instances where parameter values span a large
SMC

<figure><img src="image_9.png" /><figcaption>**Figure 6.** High-reliability catalog quality cuts: UV (left panel), optical (middle panel), and infrared (right panel) CMDs of all SMC sources that pass the HRC quality cuts (top rows), and sources that failed (lower rows). The inverted catalogs are shown in gray for reference. The majority of poorly fit BEAST sources are dim and generally fall below the 50% completeness limits calculated from artificial star tests (dotted lines), except for the $A_V$ SNR cut, which targets brighter diffraction spikes.</figcaption></figure>

dynamic range (e.g., mass, age), we compute the rms errors using the logarithmic parameter values. From these values, we see that measurements of certain parameters only improve marginally with additional filters (e.g., age, distance, $f_A$), whereas others improve significantly (e.g., mass, $A_V$, $R_V$).

To visualize these trends, in Figure 8, we plot the differences between the input BEAST parameter values used to generate each simulated star and the output BEAST median values for

five of the most common filter combinations, ranging from three to seven filters.<fnref n="12" /> We show both the stellar and dust parameters, and the derived stellar parameters ($T, g, L$).

With three-filter coverage, we can reliably measure most stellar masses and older ages ($\log(A) > 8$ yr). Despite large

<footnote n="12">¹² Exact filters are detailed in Table 3. We omit analyzing filter combinations that only occur in singular instances.</footnote>
LMC

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Figure 6 but for the LMC: UV (left panel), optical (middle panel), and infrared (right panel) CMDs of all LMC sources that pass the HRC quality cuts (first rows), and sources that failed (lower rows). The inverted catalogs are shown in gray for reference.</figcaption></figure>

uncertainties, we also reliably measure dust extinction ($A_V$), although, there is systematic bias toward higher magnitudes, the potential source of which we will address later on in Section 6.1.2. Most other BEAST parameters like metallicity ($Z$), $R_V$, $f_A$, and distance, are poorly recovered with three-filter coverage.

With additional filter coverage, we observe significant improvement in the measurement of certain parameters, such as metallicity ($N_{\text{filters}} \ge 6$), extinction ($N_{\text{filters}} \ge 4$), and $R_V$ ($N_{\text{filters}} \ge 6$). These filter requirements align with previous

BEAST papers, which required a minimum four-band detection for stars in NGC 4214 (Y. Choi et al. 2020). Both $f_A$ and distance remain poorly recovered, even with seven filters of observations.

Figure 8 reveals systematic errors present in the BEAST output, notably overestimation at low input values (e.g., age, $R_V$, $f_A$, and distance) and underestimation at high input values (e.g., $M_{\text{ini}}$, $Z$, $R_V$, $f_A$, and distance). For a parameter like initial mass with an IMF-based prior, the default expectation is that most observed sources will be lower mass. Therefore, if a
<table><caption>**Table 3**
Filter Combinations</caption><thead><tr><th>Name</th><th>No. Filters</th><th>F225W</th><th>F275W</th><th>F336W</th><th>F475W</th><th>F814W</th><th>F110W</th><th>F160W</th><th>N<sub>LMC</sub></th><th>N<sub>SMC</sub></th></tr></thead><tbody><tr><td><b>A</b></td><td><b>7</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>28</td><td>4</td></tr><tr><td><b>B</b></td><td><b>6</b></td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>2</td><td>6</td></tr><tr><td>C</td><td>6</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td><b>X</b></td><td>2</td><td>2</td></tr><tr><td>D</td><td>5</td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td><b>X</b></td><td>3</td><td>2</td></tr><tr><td>S5</td><td>5</td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td>...</td><td>0</td><td>0</td></tr><tr><td><b>E</b></td><td><b>5</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td>...</td><td>2</td><td>4</td></tr><tr><td>...</td><td>4</td><td>...</td><td><b>X</b></td><td>...</td><td><b>X</b></td><td><b>X</b></td><td>...</td><td><b>X</b></td><td>1</td><td>0</td></tr><tr><td>F</td><td>4</td><td>...</td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td><b>X</b></td><td>2</td><td>3</td></tr><tr><td><b>G</b></td><td><b>4</b></td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td>...</td><td>7</td><td>4</td></tr><tr><td><b>H</b></td><td><b>3</b></td><td>...</td><td>...</td><td><b>X</b></td><td><b>X</b></td><td><b>X</b></td><td>...</td><td>...</td><td>9</td><td>8</td></tr></tbody></table>

**Note.** This table summarizes the unique filter combinations for the Scylla and METAL fields used in this work. The filter combination names are carried over from C. Murray et al. (2024). The number of fields with each combination in the LMC ($N_{LMC}$) and SMC ($N_{SMC}$) is indicated in the right-hand columns. The bolded filter combinations are plotted in Figure 8.

<table><caption>**Table 4**
BEAST Simulation Validation</caption><thead><tr><th>Name</th><th>No. Filters</th><th>Δ log(M<sub>ini</sub>)<br>(M<sub>☉</sub>)</th><th>Δ log(A)<br>(yr)</th><th>ΔZ<br>(Z<sub>☉</sub>)</th><th>ΔDistance<br>(kpc)</th><th>ΔA<sub>V</sub><br>(mag)</th><th>ΔR<sub>V</sub></th><th>Δf<sub>A</sub></th><th>N<sub>stars</sub></th></tr></thead><tbody><tr><td><b>A</b></td><td><b>7</b></td><td>0.1821</td><td>1.3014</td><td>0.1326</td><td>15.6115</td><td>0.2285</td><td>0.4728</td><td>0.1714</td><td>8174</td></tr><tr><td><b>B</b></td><td><b>6</b></td><td>0.1892</td><td>1.3099</td><td>0.1333</td><td>15.7091</td><td>0.2384</td><td>0.4808</td><td>0.1827</td><td>8168</td></tr><tr><td>C</td><td>6</td><td>0.2013</td><td>1.3304</td><td>0.1461</td><td>15.7633</td><td>0.3366</td><td>0.4756</td><td>0.1760</td><td>8198</td></tr><tr><td>D</td><td>5</td><td>0.2076</td><td>1.3388</td><td>0.1471</td><td>15.8588</td><td>0.3451</td><td>0.4828</td><td>0.1870</td><td>8187</td></tr><tr><td>S5</td><td>5</td><td>0.2127</td><td>1.3477</td><td>0.1489</td><td>15.9908</td><td>0.3196</td><td>0.5870</td><td>0.1984</td><td>8211</td></tr><tr><td><b>E</b></td><td><b>5</b></td><td>0.2195</td><td>1.3582</td><td>0.1550</td><td>16.0989</td><td>0.3892</td><td>0.7586</td><td>0.2038</td><td>8231</td></tr><tr><td>F</td><td>4</td><td>0.2349</td><td>1.3848</td><td>0.1742</td><td>16.1338</td><td>0.4982</td><td>0.5069</td><td>0.1968</td><td>8201</td></tr><tr><td><b>G</b></td><td><b>4</b></td><td>0.2289</td><td>1.3749</td><td>0.1586</td><td>16.1850</td><td>0.4067</td><td>0.7917</td><td>0.2136</td><td>8232</td></tr><tr><td><b>H</b></td><td><b>3</b></td><td>0.2864</td><td>1.4516</td><td>0.1830</td><td>16.3361</td><td>0.6458</td><td>0.9012</td><td>0.2245</td><td>8229</td></tr></tbody></table>

**Note.** rms errors for multiband BEAST fits compared to inputs for 10⁴ simulated stars. Table 3 outlines filters used in each filter combination.

more-massive star is poorly constrained, its mass will be underestimated. Cumulatively, this results in a systemic underestimation, which we see for BEAST initial mass values greater than one solar-mass. The same principle applies to BEAST age estimates, which assumes a flat SFH prior, resulting in an overestimation of ages for young stars. At high metallicities, underestimation is likely due to edge effects in the model gridding. For example, stars simulated with $Z = 0.5 Z_{☉}$ will intrinsically be fit with some level of scatter, both greater than and less than $0.5 Z_{☉}$. However, in the absence of any models with $Z > 0.5 Z_{☉}$, only the lower scatter will be captured in the distributions, resulting in a systematic underestimation at grid edges. This edge effect is also present in the $R_V$, $f_A$, and distance results, where lower values are systematically overestimated, and high values are systematically underestimated.

In all filter scenarios, distances remain poorly constrained, with distances systematically trending toward 50 kpc regardless of input values, despite providing a flat distance prior. At this time, it is unclear what the source of this systematic trend is. These fits could potentially be improved by providing a more finely spaced distance grid, or by iteratively updating distance fits by collectively modeling the distance distribution of all sources in a field (which will be possible with the MegaBEAST.<fnref n="13" />

Due to the variable quality of parameters based on the number of filter coverage, some discussions will focus on parameter-optimized subsets of the observed fields. For

example, Z and $R_V$ measurements are only well constrained in the six-to-seven-filter cases, so when discussing Z and $R_V$ results in detail, we will only use measurements from fields with six or seven filters. Similarly, discussions of extinction results will emphasize data from fields observed with at least four filters to avoid systematic bias from three-filter fits.

## 5. Results

In this Section, we present the BEAST fits in detail, describe any notable systematics, and discuss the differences between fits in the two galaxies. The BEAST catalogs are made publicly available as High-Level Science Products on MAST.<fnref n="14" /> For all parameter values, we use the reported median values (p50), as opposed to the singular best-fit model values (Best) or the expectation values (Exp). We opt to use the p50 values since these measurements are calculated using the entire posterior grid, meaning we can also estimate parameter uncertainties from the width of the marginalized 1D pPDF.

In Figure 9, we plot the 2D distributions of the median p50 values for all primary stellar (initial mass, age, metallicity, distance) and dust ($A_V$, $R_V$, $f_A$) parameters for all high-reliability sources in the LMC (blue) and SMC (red). Correlations reflect trends in the stellar population. Along the diagonal, we plot the 1D histograms of all of the parameters individually, along with the corresponding median and 1σ values for the two galaxies.

<footnote n="13">¹³ https://github.com/BEAST-Fitting/megabeast</footnote>

<footnote n="14">¹⁴ https://archive.stsci.edu/hlsp/scylla</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Validation: the median ±1σ difference between the BEAST-modeled (“out”) and simulated (“in”) parameters vs. the simulated parameters for 10⁴ stars simulated in three-, four-, five-, six-, and seven-filter cases, with high-reliability cuts applied. For initial masses ($M_{ini}$), we show both the full input parameter range, and a zoom-in on sources with $M_{ini} \le 2 M_{\odot}$ (dashed). Additional filters improve BEAST SED fits for nearly all parameters, especially Z, $A_V$, and $R_V$.</figcaption></figure>

**Initial Mass and Age.** Similar stellar initial mass ($M_{ini}$) and age (log A) distributions are observed in the stellar populations of the two galaxies, with median and standard deviation initial mass estimates of 0.94$^{+0.31}_{-0.16}$ (0.92$^{+0.30}_{-0.11}$) $M_{\odot}$ and ages of 4.38$^{+1.83}_{-2.70}$ (5.67$^{+1.89}_{-3.39}$) Gyr in the SMC (LMC), respectively. There is a strong correlation between initial mass and age that is expected, since higher-mass stars have shorter life spans than low-mass stars. We anticipate few BEAST fits of massive stars since these stars are often saturated in the photometry (see Figure 2 in Section 2 for further details).

**Metallicity.** Although the metallicity grid was the same for both galaxies, stars in the LMC are on average fit to be more metal-rich than the SMC, with median metallicity fits of 14$^{+26}_{-9}\%$ $Z_{\odot}$ ([M/H] = −0.71$^{+0.44}_{-0.90}$, assuming $\alpha$/Fe = 0.2) in the LMC versus 6$^{+16}_{-3}\%$ $Z_{\odot}$ ([M/H] = −1.07$^{+0.65}_{-1.37}$) in the SMC. The 1D metallicity histogram for the LMC features several spikes, particularly at 20% and 50% $Z_{\odot}$. The metallicity grid was originally constructed to span from ~1% solar to ~50% $Z_{\odot}$, based on the prior knowledge of the metallicities in the MCs. However, the pile-up of LMC sources at ~50% $Z_{\odot}$ indicates that fits would likely have improved had higher metallicities been included. Additionally, simulations from Section 4 show that metallicity is systematically underestimated when fitting with only optical and UV filters (i.e., three to five filters), meaning many of the sources in Figure 9 would likely have been fit as higher metallicity had additional IR filters been available. Future survey fittings and releases will account for this issue and expand the metallicity range.

**Distance.** The distance grid was intentionally different for the two galaxies, with the LMC grid ranging from 40–60 kpc, and the SMC grid ranging from 47–77 kpc. This was done to account for the prior distance estimates for the galaxies, with the LMC at 49.89 kpc and the SMC at 61.94 kpc, respectively (R. de Grijs et al. 2014; R. de Grijs & G. Bono 2015). We expect a broad distance distribution for both galaxies, especially for the SMC, which is quite extended along the line of sight (reviewed in C. E. Murray et al. 2024).

For the LMC, BEAST distance estimates span the full grid with a median distance fit of 49.2$^{+3.9}_{-2.7}$ kpc; however, the distribution of distances is slightly asymmetric, skewed toward larger distances. This distance distribution agrees with depth measurements from S. Subramanian & A. Subramaniam (2009), who found an average depth of 4.0 ± 1.4 kpc in the bar region and 3.4 ± 1.2 kpc in the disk. We believe the BEAST-derived distance distribution might be slightly smaller, since we measure both young and old stars, whereas S. Subramanian & A. Subramaniam (2009) only measured the distances to red clump stars, which are generally older and would therefore have a larger scale height.

In the SMC, distance estimates also span the full grid, with a median distance fit of 61.2$^{+6.8}_{-6.7}$ kpc; however, the distance distribution is not symmetric, and there is evidence of bimodality in the distribution. These results corroborate other studies that show the SMC comprises two star-forming systems separated by several kiloparsecs along the line of sight (S. Subramanian et al. 2017; C. E. Murray et al. 2024). This multicomponent structure could potentially be the result of tidal
<figure><img src="image_12.png" /><figcaption>**Figure 9.** BEAST parameter results for all HRC sources in the LMC and SMC: 2D distributions of the median p50 values for all BEAST primary stellar parameters—initial mass ($M_{ini}$), age (A), metallicity (Z), distance, and dust parameters—$A_V$, $R_V$, $f_A$, for all high-reliability sources in the LMC (blue) and SMC (red). One-dimensional histograms of all parameters are shown along the diagonal along with the corresponding median and $1\sigma$ values for the SMC and LMC, respectively.</figcaption></figure>

stripping due to a recent near impact (~5 kpc) with the LMC roughly 150 Myr ago (Y. Choi et al. 2022).

While median BEAST distance fits do agree with existing distance measurements, validation tests (Section 4) showed that

individual distance fits are often poorly recovered, spanning ±5–10 kpc in uncertainty. One potential explanation for this paradox is that individual stars might not contain enough information to recover distances; however, as an ensemble with
thousands of sources, we can recover signal, even if individual sources have low SNR. As such, we caution anyone wanting to use distance measurements for individual stars to be careful in how they choose to interpret these measurements.

*Extinction.* Both galaxies exhibit low amounts of extinction ($A_V$) along the line of sight toward stars, with median extinction values of ~0.31$^{+0.26}_{-0.17}$ (0.42$^{+0.45}_{-0.24}$) mag in the SMC (LMC), respectively. These extinction measurements include contributions for Milky Way foreground dust. The extinction distributions in both galaxies follow log-normal distributions, with tails leading toward higher extinctions. The log-normal distribution in the LMC is broader than in the SMC.

*Dust Properties.* The effects of extinction are more influential on the observed SEDs than other dust parameters like $R_V$ and $f_A$. As a result, the influence of these secondary dust parameters is far harder to constrain without high-quality observations farther into the UV (e.g., F275W and F225W). When observations lack enough information to inform the likelihood, the posterior distribution will default to the prior distribution. The BEAST model grid adopted the same 2D prior for the secondary dust parameters $R_V$ and $f_A$ as described in Figure 7 in K. D. Gordon et al. (2016), with a peak in the marginalized $R_V$ distribution around ~2.7 and $f_A$ at 1.0. Given the resolution of the parameter gridding, the weighted average of this 2D prior falls near [3.55, 0.73], which is near what most of the observations indicate in both galaxies with median $R_V$ fits of 3.38$^{+0.18}_{-0.25}$ (3.35$^{+0.28}_{-0.24}$) and median $f_A$ fits of 0.78$^{+0.04}_{-0.04}$ (0.78$^{+0.05}_{-0.05}$) in the SMC (LMC). Lack of variation across both galaxies and the similarity to the 2D prior indicates most $R_V$ and $f_A$ measurements are poorly constrained in most fields. Individual fields with coverage in multiple UV filters (F275W, and F336W) could potentially constrain variations in $R_V$, something future Scylla science publications will investigate.

## 5.1. Filter Variation

There are considerable variations in the BEAST results and uncertainties when fitting sources with three versus four versus seven bands of observations. In this Section, we investigate these variations on an individual Scylla field and comment on their application to the broader results of the paper. We use 14675_LMC-15562nw-33974 (M6) to investigate the effects of filter variation on BEAST results. This field is located in the upper-west corner of the LMC and was originally observed in seven filters. Although this field spatially coincided with two other Scylla fields, LMC_8 and M5, with three- and six-filter coverage, respectively, we find that minor differences in the photometry between each field make it difficult to isolate the effects of filter coverage. For a true apples-to-apples comparison, we instead chose to refit the M6 field using data from only three- and four-filter combinations (H and G in Table 3). We apply HRC cuts to all BEAST catalogs, retaining 1646, 1670, and 1614 sources in the three-, four-, and seven-filter catalogs, respectively, out of the total 5342 sources with BEAST fits.

In Figures 10 and 11, we plot the 2D distribution of BEAST parameters for all HRC sources in the three-filter catalog (green), four-filter catalog (red), and seven-filter catalog (blue). By overlaying the 1D and 2D distributions, we can see that the distributions of certain BEAST parameters vary considerably depending on the number of filters. Most notably, stellar metallicity is measured to be much lower when fit with only three filters, compared to seven filters. Meanwhile, the median extinction tends to be double when measured with three filters

compared to seven. From these trends, we surmise that any three-filter fields with low-metallicity estimates are likely incorrect and that median galactic estimates of extinction in Figure 9 are likely overpredicted due to contributions from three-filter fields. From Figure 11, we see that the inclusion of even one additional filter (F275W) increases metallicity measurements for nearly all sources. In addition, median extinction measurements are lower; however, there is still a slight overestimation compared to seven-filter measurements.

# 6. Discussion

With high-reliability stellar and dust fits for over 550,000 sources, we can now analyze what BEAST results tell us about the stellar and dust content of the clouds and how it compares to existing literature. In this Section, we first discuss the stellar content of these Scylla fields, see how BEAST-fit mass estimates compare with isochrone-predicted masses, and discuss the potential influence of binary populations (Section 6.1). Second, we compare the BEAST age distributions with existing SFH measurements (Section 6.2). Third, we quantify our sensitivity to $A_V$ as a function of stellar luminosity, characterize the log-normal extinction distributions in each galaxy, and compare the average extinction values with other ISM tracers (Section 6.3).

## 6.1. Stellar Content

A galaxy consists of many types of stellar populations that can be used to address broader science questions about the galaxy itself: e.g., IMF (D. R. Weisz et al. 2013), SFH (A. R. Lewis et al. 2015), disk structure (B. F. Williams et al. 2012), metallicity distribution (D. Gregersen et al. 2015), total column extinction (J. J. Dalcanton et al. 2015)—as well as more targeted questions about physical processes in galaxy evolution: e.g., AGB stars as a source of dust production (S. R. Goldman et al. 2022) or formation of stellar clusters (L. C. Johnson et al. 2017).

From CMDs alone, we can identify several stellar populations based on their locations within color–magnitude space. In Figures 12 and 13, we plot optical CMDs for all high-reliability sources in the LMC and SMC, respectively. The sources in the CMDs are split up by the number of filters to help compare how variable filter coverage influences the measured BEAST parameters. The filter numbers were chosen such that similar numbers of sources were visualized in each CMD.

In the first row of Figures 12 and 13 (a,f,k), we show the density distribution of all source n fields with three-, four-, and six-to-seven-filter coverage in each galaxy, both as colored 2D histograms and outlined in logarithmic levels using kernel-density estimates (KDEs). These KDE outlines are propagated to all other rows to provide a reference for where the majority of stars are located. Both galaxies feature a prominent main sequence, running from the top left (F475W–F814W = −0.5, F814W = 18 mag) down to the lower right (F475W–F814W = 1.5/1.8, F814W = 25 mag). Perpendicular to this structure is the RGB, which extends upward from F814W = 23 mag to 17 mag. Along this branch is the red clump (RC) at ~19 mag, which can be seen as an overdensity along the RGB both in the SMC and LMC.

While the morphology of the RGB and RC can be intrinsically complex (L. Girardi 2016), their overall distribution on a CMD contains valuable information about a galaxy,
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Filter impact on BEAST parameter results (three vs. seven): 2D distribution of BEAST parameters for an example field in the LMC (14675_LMC-15562nw-33974; M6). The field was originally observed and fit with seven filters (blue), and BEAST fits were repeated using three filters (F336W, F475W, and F814W; green) for comparison. Some BEAST parameters, like initial mass ($M_{ini}$) and age, are measured to have similar marginalized distributions with both three- and seven-filter coverage. Other parameters like metallicity ($Z$), however, differ significantly depending on the number of filters. In particular, the median extinction ($A_V$) measured with three-filter coverage is nearly double the median measured with seven-filter coverage.</figcaption></figure>

especially extinction, which can cause considerable broadening of these CMD features. This is especially notable in the LMC, where the RC structure appears to be elongated down and to the right. Based on extinction coefficients for these optical bands, we can predict the reddening effect of some magnitude of extinction, also referred to as a reddening vector. This reddening vector for $A_V = 1$ mag is shown in the

bottom-right column. The width of the LMC RC, compared to the orientation of the reddening vector shows that the elongation of the LMC RC is almost certainly the result of dust extinction.

In all subsequent rows of Figures 12 and 13 (subfigures b-e, g-j, and l-o), we plot the spatially binned medians for primary BEAST parameters.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Filter impact on BEAST parameter results (four vs. seven): 2D distribution of BEAST parameters for an example field in the LMC (14675_LMC-15562nw-33974; M6). The field was originally observed and fit with seven filters (green), and BEAST fits were repeated using four filters (F275W, F336W, F475W, and F814W; red) for comparison. The inclusion of F275W increases BEAST metallicity measurements significantly compared to three-filter fits (Figure 10).</figcaption></figure>

*Initial Mass.* In Figures 12 and 13 (b, g, l), we plot the average initial mass estimates. The smallest/largest initial mass estimates for sources in the SMC and LMC are $0.58_{-0.04}^{+0.05}/24.39_{-6.74}^{+1.60}M_{\odot}$ and $0.55_{-0.02}^{+0.03}/51.27_{-2.99}^{+14.29}M_{\odot}$, respectively; however, for visual purposes, we limit the color range to 0.6–6 $M_{\odot}$. The main sequence spans a wide range of initial masses, whereas the RGB and RC are primarily composed of sources with initial masses ~1 $M_{\odot}$, consistent with

expectations (e.g., L. Girardi 2016; D. M. Nataf et al. 2021). Fields with six-to-seven-filter coverage contain more high-mass stars than fields with three-to-four-filter coverage. This is likely the result of multiple-variable-like source saturation without guard exposures and spatial location of six-to-seven-filter fields near star-forming regions.

*Stellar Age.* In Figures 12 and 13 (c, h, m), we plot the average ages of stars across the CMD. We find that age is spatially
<figure><img src="image_15.png" /><figcaption>**Figure 12.** CMD BEAST trends in the LMC: optical CMDs for all high-reliability sources in the LMC from fields with three-filter coverage (left; *n* = 80, 399), four-filter coverage (middle; *n* = 56, 446), and six-to-seven-filter coverage (right; *n* = 145, 170). Top columns (*a*, *f*, *k*): density distribution of sources as colored 2D histograms, and outlined in logarithmic levels using kernel-density estimates (KDEs). KDE density outlines (white) are propagated to all other rows. All other rows: spatially binned (100 × 100) medians for notable BEAST parameters: initial mass (*M*<sub>ini</sub>), age, metallicity (*Z*), and extinction (*A*<sub>V</sub>). In panel *o*, we show the direction of a 1 mag reddening vector with a black arrow. UV and IR versions of the six-to-seven-filter sources are provided in the Appendix.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** CMD BEAST trends in the SMC (the same as Figure 12): optical CMDs for all high-reliability sources in the SMC from fields with three-filter coverage (left; *n* = 74, 270), four-filter coverage (middle; *n* = 61, 031), and six-to-seven-filter coverage (right; *n* = 69, 225). Notable differences from Figure 12: extinction color range has been reduced to improve dynamic range. UV and IR versions of the six-to-seven-filter sources are provided in the Appendix.</figcaption></figure>

correlated with initial mass, i.e., higher-mass stars ($M > 4 M_\odot$) tend to be <100 Myr old, whereas lower-mass ($M < 1 M_\odot$) tend to be >10 Gyr old. For higher-mass stars in the SMC, there is a gradient in age, with redder stars being more evolved.

*Metallicity*. In Figures 12 and 13 (d, i, n), we plot the average metallicities and find notable differences between the distributions depending on the number of filters in both the SMC and the LMC. In the SMC, there is a gradient in metallicity from
3% to 20% $Z_\odot$, perpendicular to the main sequence, where redder sources are more metal-rich. This metallicity gradient could explain why the single-metallicity isochrones plotted in Figure 2 did not perfectly match the observed population distribution on the CMD. By comparison, the median metallicity for most sources in the RGB and RC is ~10% $Z_\odot$, similar to spectroscopic studies from APOGEE of RGB stars in the LMC (J. T. Povick et al. 2024). In the LMC, the median metallicity for main-sequence stars is more homogeneous (20%-30% $Z_\odot$), and both MS and RGB/RC metallicities are on average much greater than the SMC, 25% $Z_\odot$ compared to 10% $Z_\odot$. In both galaxies, high-mass stars ($M > 4 M_\odot$) show the highest amount of enrichment, at 50% $Z_\odot$, whereas low-mass stars ($M < 0.7 M_\odot$) are the least enriched ($<3\% Z_\odot$).

**Extinction.** In Figures 12 and 13 (e, j, o), we plot the average extinction and show the direction of a 1 mag reddening vector. The LMC contains far greater amounts of extinction (0.3–3.0 mag) than the SMC (0.2–1.0 mag). The elongated RC and RGB in both galaxies show a gradient in extinction, ranging from 0.1–0.7 mag in the SMC, which is in agreement with existing RGB $A_V$ estimates in the SMC of $\langle A_V \rangle = 0.41 \pm 0.09$ mag (P. Yanchulova Merica-Jones et al. 2021), and 0.1–4.0 mag in the LMC. Surprisingly, the three-filter SMC contains lower-main-sequence sources with extinctions greater than the maximum extinction experienced by the RGB: $A_V = 1.7$ mag for the lower main sequence versus $A_V = 0.7$ mag for the RGB. In Section 6.1.2, we discuss the credibility of these $A_V$ estimates and show how binary stars could influence these results.

### 6.1.1. Evolutionary Track Mass Estimates

Due to the deep high-resolution photometry from HST, we resolve subsolar mass sources. To validate mass estimates from the BEAST, we compare the CMD distribution of BEAST initial mass estimates in the SMC with PARSEC stellar evolutionary tracks (A. Bressan et al. 2012). We opt to only make this comparison in the SMC since the LMC has much more extinction and a larger range of metallicities, which complicates this comparison.

In Figure 14, we plot the spatially binned optical CMD of all high-reliability sources in the SMC, colored by the median initial mass estimates from the BEAST. This Figure is the same as Figure 13(b), (g), and (l), combined, except the color bar now shows discrete mass bins to aid in the comparison. PARSEC stellar evolutionary tracks are overlaid, with age ranging from 10 Myr (log(Age) = 7 yr) to 6.3 Gyr (log (Age) = 9.8 yr) for stars with three initial masses: 0.7, 1.0, and 2.0 $M_\odot$. These models assume a metallicity of 20% $Z_\odot$, a distance of 60 kpc, and apply no extinction.

The evolutionary tracks for stars with initial masses of 0.7 and 1.0 $M_\odot$ are in agreement with the BEAST-predicted initial mass estimates. The 2.0 $M_\odot$-track also overlaps with BEAST estimates; however, the track seems to overlap more with slightly higher-mass stars ($M = 2.5–3.0 M_\odot$), and there is less agreement given the limited spatial coverage of the CMD. There are two potential reasons for the lack of alignment between higher-mass stars and their evolutionary tracks: (1) it could be that higher-mass stars are located in dusty environments, which would make the observed stars dimmer and redder than the stellar evolutionary track; (2) it could be that the higher-mass stars have metallicities >20%, which would make

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Lower-mass limit: spatially binned optical CMD of all high-reliability sources in the SMC, colored by median initial mass. PARSEC stellar evolution tracks for stars with initial masses of 0.7, 1.0, and 2.0 $M_\odot$ are overlaid, assuming no extinction and $Z = 20\% Z_\odot$ (A. Bressan et al. 2021). We compare stellar evolutionary tracks with the SMC, because it has a smaller range of metallicities and less extinction than the LMC. BEAST initial mass estimates are in broad agreement with evolutionary tracks.</figcaption></figure>

them dimmer and redder for a given mass on the main sequence. From Figure 13, we know main-sequence stars with 2.0 $M_\odot$ in the SMC exhibit a gradient of metallicities, ranging from 3%-30% $Z_\odot$, while extinction varies from 0.1–0.5 mag, so both of the effects noted above are likely responsible for the offset between the observed stars and the stellar evolutionary tracks at higher masses. Overall, the BEAST initial mass estimates agree with stellar evolutionary tracks.

### 6.1.2. Binary Stars

As noted in earlier discussions of extinctions, parts of the lower main sequence in Figure 13(e) seemed to experience double as much extinction as the RGB: $A_{V,max} = 1.7$ mag for the lower main sequence versus $A_{V,max} = 0.7$ mag for the RGB. Given their bright luminosities and extended scale heights relative to most star-forming galaxy disks, RGB sources are often used to quantify the total column extinction of external galaxies (e.g., J. J. Dalcanton et al. 2015; P. Yanchulova Merica-Jones et al. 2021). The $A_V$ range of the RGB in the SMC should capture the average range of extinctions within the galaxy. It is therefore unexpected to observe dim lower-main-sequence stars with more than double the extinction, according to BEAST fits.

One potential explanation for this abnormally high extinction is the lack of binary models in the BEAST (adding binaries is in progress). The binary fraction of stars (and subsequent mass ratios) strongly depends on mass and metallicity and is still an active field of research. For most binaries, if the mass of the secondary is considerably less than the primary (assuming similar ages and metallicities), the difference between the combined flux versus the primary flux alone will be negligible since $L \propto M^3$. For example, a mass ratio of 10% will roughly translate to a luminosity ratio of 0.1%. However, for binaries with larger mass ratios (e.g., mass ratio ~50%), the combined increase in magnitude and color can be prominent on a CMD.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Effects of binaries: spatially binned optical CMD of all high-reliability sources in the SMC, colored by median BEAST extinction. In the lower left, we show the direction of a 1 mag reddening vector as a black arrow. PARSEC isochrones for log(Age) = 9.8 yr with Z = 20% Z⊙ and no extinction are overlaid in blue (A. Bressan et al. 2012). Applying $A_V = 1$ mag produces the isochrone shown in red. Assuming all stars are binaries with an equal-mass companion (binary mass ratio = 1) produces the isochrone shown in orange. The BEAST does not currently model include binary models; however, both effects are capable of populating the redder lower main sequence.</figcaption></figure>

To illustrate this effect, in Figure 15, we overlay PARSEC isochrones (log(Age) = 9.8 yr, $Z = 20\% Z_⊙$, distance of 60 kpc) with $A_V = 0$ mag (blue), $A_V = 1$ mag (red), and $A_V = 0$ mag but assuming equal-mass binary companions (mass ratio = 1; orange). Both the extinct isochrone (red) and the binary isochrone (orange) overlap the high-extinction sources in the lower main sequence. However, along the RGB, the effects of binaries are degenerate with single-star isochrones. We note that this overlap will vary depending on the binary mass ratio and the filter combinations.

It is possible that the high-extinction sources in the lower main sequences of the SMC are actually the consequence of misclassified binary stars. Future BEAST modeling of binary stars will amend these misclassifications; however, due to the scarcity of these stars, they have relatively little impact on the average extinction of a field; the median extinction per field is on average $0.050 \pm 0.035$ mag and $0.058 \pm 0.053$ mag less dusty in the SMC and LMC, respectively, if we exclude these dimmer sources ($F814W > 22$ mag). As such, we do not believe these sources will have a significant impact on the overall BEAST statistics.

## 6.2. Star Formation History

The MCs have a long interaction history, which is thought to have generated bursts of star formation when close encounters occurred (D. R. Weisz et al. 2013; E. Patel et al. 2020; R. E. Cohen et al. 2024a, 2024b). Ongoing or recent star formation ($\le 100$ Myr) can be traced via observables such as $H\alpha$ or UV radiation from massive stars (D. R. Weisz et al.

2012; K. B. W. McQuinn et al. 2015; M. Cignoni et al. 2019). However, tracing the SFH of galaxies requires studying all stellar populations.

Stellar modeling shows how stars with different initial masses and metallicities evolve and change color/magnitude over the course of their lives. On a statistical level, these stellar models allow us to figure out what history of star formation would have produced the distribution of stars observed in a CMD in the present day. Critical to this analysis is the assumption of an IMF, the rate of enrichment (i.e., age-metallicity relation), and the stellar evolution models employed.

The SFHs of the SMC and LMC have been extensively studied. Most recently, (R. E. Cohen et al. 2024a, 2024b) used two-band optical imaging from the Scylla survey (C. Murray et al. 2024) and archival HST imaging to constrain the SFHs across different regions in the MCs. They calculated the SFHs for each field by modeling the observed CMDs using the software code MATCH (A. E. Dolphin 2002). One of the main metrics measured from these studies is $\tau_X$, the amount of time it takes to form X% of the total stellar mass (e.g., $\tau_{50}$ is the lookback time by which 50% of the cumulative stellar mass of a galaxy has been formed). These metrics are used to make it easier to compare SFHs between galaxies of different masses.

The BEAST is not designed to infer the SFH of galaxies, as it is only able to observe the age distribution of stars in the present day (since older high-mass stars will be missing). However, we can still compare the observed age distributions with SFH measurements to validate BEAST age estimates.

In Figure 16, we plot the cumulative distributions of the BEAST stellar age for each Scylla field in the SMC (left; red) and LMC (right; blue), and highlight the cumulative distribution for all sources combined. Markers note the 50th ($\odot$), 75th ($\times$), and 90th ($\Delta$) percentiles of the cumulative stellar age distributions based on all HRC sources compared to the median $\pm 1\sigma\tau_{50}$, $\tau_{75}$, and $\tau_{90}$ measurements for each Scylla and METAL field in the SMC and LMC (R. E. Cohen et al. 2024a, 2024b). We also plot the total cumulative SFH for each galaxy as a gray line and show the cumulative distribution for a constant star formation rate as a dashed line.

The BEAST stellar age distributions show the majority of observed stars are younger than 6–7 Gyr, with 77.5% (93.9%) of stars in the SMC and 57.1% (78.8%) of stars in the LMC being younger than 6.0 (7.0) Gyr, respectively. These findings are in strong agreement with existing measurements of the SFR, which found a burst of star formation around 5–6 Gyr ago. In the SMC, the median observed stellar age distribution is in strong agreement with the cumulative mass fraction for the first 6 Gyr. However, we observe a divergence beyond this age limit that is likely caused by either: (1) observational bias of late-stage intermediate-mass stars (i.e., white dwarf stars are likely too dim to be observable); or (2) observational bias of low-mass stars beyond the completeness limits. In the LMC, there is less agreement between the observed stellar age distribution and the cumulative mass fraction; however, the median 50th, 75th, and 90th age percentiles are still within $1\sigma-2\sigma$ of the corresponding $\tau$ measurements.

Comparing $\tau_{50}$ and $A_{50}$ comes with caveats. Importantly, BEAST age distributions are calculated using observed present-day stellar populations, meaning the stellar age distributions will be missing stars that have already died or stars that are too dim/bright to be observable. By comparison,
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Stellar age vs. star formation history (SFH): cumulative stellar age distributions of sources in each Scylla field in the SMC (left; red) and LMC (right; blue). The 50th (○), 75th (×), and 90th (△) percentiles of the cumulative stellar age distributions for all HRC sources in the LMC and SMC are marked. We overlay median ±1σ τ<sub>50</sub>, τ<sub>75</sub>, and τ<sub>90</sub> measurements for all Scylla and METAL field from R. E. Cohen et al. (2024a, 2024b; black) and the total cumulative SFH for the SMC and LMC from gray). For comparison, the distribution of a constant star formation rate is shown as a dashed line. Overall, the distributions of age estimates from the BEAST are in broad agreement with independent SFH estimates, especially for the SMC.</figcaption></figure>

SFH fits take observational incompleteness into account. Additionally, SFH fits usually include some consideration for binary stars, whereas the BEAST assumes a binary fraction of 0 since it only uses single-star models. Intrinsically, SFH fits like τ<sub>50</sub> are a measurement of the time it took to form half of the cumulative stellar mass ever formed, whereas A<sub>50</sub> only measures the median age of the stars observable in the present day. In general, A<sub>50</sub> should be systematically younger than τ<sub>50</sub>, especially for further lookback times. With all of this in mind, it is rather impressive that the stellar age distributions from the BEAST are in such strong agreement with robust SFH measurements, especially in the SMC. These results ultimately indicate that, for stellar populations, the BEAST fits are capable of accurately recovering the ages of stars.

## 6.3. Dust Properties

Measuring dust extinction through resolved stellar SED fitting gives us valuable new insight into the gas content of nearby galaxies. In this Section, we examine the BEAST $A_V$ measurement more closely and discuss upcoming mapping techniques, define detection limits for high extinction, quantify the average extinction distribution for each galaxy, and compare dust measurements with other ISM tracers available.

### 6.3.1. Mapping Dust Extinction

Mapping the distribution and composition of dust structures in galaxies is critical for tracing the multiphase structure of the ISM. Most dust maps are either constructed using IR dust emission, which relies on assuming a dust emissivity coefficient, or use stellar types like RGB stars to measure a reddening vector in a stellar population (J. J. Dalcanton et al. 2015; P. Yanchulova Merica-Jones et al. 2017), although this method limits the number of viable stellar probes. By comparison, BEAST extinction measurements provide high-resolution probes of the ISM using all of the stars in a field, not just certain stellar types. In upcoming Scylla publications, we will harness these discrete probes of extinction to infer the underlying distribution of dust within individual fields. For the

purposes of this paper, the main goal is to: (1) quantify how much extinction we can probe with the observed Scylla sources; and (2) compare extinction with other tracers of the ISM.

### 6.3.2. Extinction Detection Limit

The maximum extinction we can probe is dependent on both the unextinguished luminosity of stars and the wavelengths observed. Brighter stars can probe larger magnitudes of extinction before the sources fall below observational detection limits. Meanwhile, the effects of extinction are wavelength-dependent and are more prominent in the UV, affecting the F225W, F275W, and F336W bands the most. To quantify the $A_V$ detection limit in the Scylla fields, we use ASTs to test how stellar recovery varies as a function of simulated extinction and intrinsic stellar luminosity.

In Figure 17 (left), we plot the model extinctions and intrinsic luminosities for 24,039 ASTs simulated for seven-filter Scylla field 15891_SMC-1588ne-12105. Sources were generated using a grid of extinctions ranging from 0.0–10.0 mag with a grid spacing of 0.5 mag, whereas luminosities are calculated from the initial stellar mass (which is sampled from a Kroupa IMF) and log-ages ranging from 6.0–1.0 yr with a grid spacing of 1.0 dex. ASTs were processed in the same pipeline as described in Section 3.1.3. For each AST, we note whether the star was observed in F275W (cyan) or not (magenta). We find that sources dimmer than $1 L_\odot$ are generally not detected under any observing conditions, a finding corroborated by the lower-mass limit of ~0.7 $M_\odot$ found in Section 6.1.1. Conversely, stars with $\log(L) > 4 L_\odot$ are too bright and saturated with long Scylla exposure times. However, some of these bright stars become observable as larger amounts of extinction are applied.

ASTs are generated to span a wide range of observed fluxes and are not inherently representative of the majority of sources observed in the Scylla fields. To understand which stellar populations are observable with larger extinctions, we compare the photometric overlap between the observed sources in the Scylla field and the detected ASTs. In Figure 17 (right), we plot
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Extinction detection limit. Left panel: 24,039 artificial star tests (ASTs) from a six-filter Scylla field in the SMC (15891_SMC-1588ne-12105) plotted as a function of model extinction and extinguished luminosity. Sources are colored by whether they were detected in the bluest filter F275W (cyan) or not (magenta). With the given exposure times, most artificial stars are undetectable; however, some sources are too bright and saturate. Right panel: optical CMD of ASTs overlaid with the actual observed sources in the example Scylla field, with all sources in gray, and high-quality sources in black. Sources in the upper main sequence and RGB are easily observable with 3 mag of extinction.</figcaption></figure>

ASTs by their extinguished model photometric in an optical CMD, maintaining the coloring showing where ASTs were detected in F275W (cyan) and not (magenta). We overlay all photometric sources originally observed in the Scylla field (gray), as well as the sources that pass the HRC quality cuts (black).

In summary, most luminous stars in the upper main sequence and RGB (log(L) ≈ 2 L<sub>☉</sub>) are detectable with 3 mag of extinction, well below current extinction estimates for the MCs (A<sub>V</sub> = 0.46 and 0.55 mag in the SMC and LMC, respectively; D. Zaritsky et al. 2002, 2004). Massive stars (log(L) = 3–4 L<sub>☉</sub>) provide the greatest opportunity to probe extreme levels of extinction (A<sub>V</sub> = +3–8 mag). Low detection limits for sources dimmer than 1 L<sub>☉</sub> once again evidence that high-extinction low-mass sources discussed in Section 6.1.2 are likely just misclassified binary stars given that these sources should, in theory, not be detectable with A<sub>V</sub> = 2 mag.

### 6.3.3. Log-normal Extinction Distribution

The ISM is a multiscale structure with volume-filling atomic gas that cools and coalesces to dense molecular structures that eventually form stars. As a result, the distribution of gas and dust densities in a given area of a galaxy is often well described by a log-normal distribution (E. Vazquez-Semadeni 1994; E. C. Ostriker et al. 2001; J. J. Dalcanton et al. 2015; M. Lombardi et al. 2015; N. Schneider et al. 2015), where turbulence is being injected into the medium, and only a small fraction of the ISM exists at higher densities. Although we cannot measure the continuous density distributions of the ISM via extinction directly, we can still probe this structure with individual stellar sight lines.

We observe log-normal A<sub>V</sub> distributions in both the SMC and LMC, as initially illustrated in Figure 9. This is despite the fact that BEAST priors for A<sub>V</sub> were flat and covered a wide range of magnitudes (0–10 mag). In Figure 18, we characterize these distributions by calculating the best-fit log-normal functions for each galaxy. The SMC is best characterized by a log-normal distribution with μ = 0.33 and σ = 0.66 (median

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Log-normal extinction distribution: 1D histogram of median BEAST extinctions (A<sub>V</sub>) for all high-reliability sources in the SMC (red) and LMC (blue) in fields observed with at least four filters. Dashed (SMC) and dotted (LMC) lines show the best-fit log-normal functions to each distribution.</figcaption></figure>

A<sub>V</sub> = 0.29 mag), whereas the LMC is dustier with μ = 0.44 and σ = 0.76 (median A<sub>V</sub> = 0.40 mag). These findings are similar to previous extinction studies in the MCs, which found log-normal distributions with median values of 0.46 mag in the SMC (see Figure 19 in D. Zaritsky et al. 2002) and 0.55 mag in the LMC (see Figure 9 in D. Zaritsky et al. 2004).

### 6.3.4. Comparison with ISM Tracers

The various phases and composition of gas say a lot about the environments present in galaxies. Atomic hydrogen gas (H1) is the most prevalent phase of gas in the ISM, and is especially dominant in dwarf galaxies like the SMC (A. Leroy et al. 2007; K. E. Jameson et al. 2016). We observe neutral
atomic hydrogen via 21 cm emission, which allows us to
measure the abundance of atomic hydrogen and its radial
velocity. As gas cools and condenses, shielding from ionizing
radiation allows for the formation of molecules that can be used
to trace denser regions prerequisite to star formation. While
molecular hydrogen, H₂, is the most prevalent molecule, its
weak emission means it is often underdetected. As a
compromise, carbon monoxide (CO) is often used as an
alternative molecular gas tracer since its dipole moment
produces rotational transitions (e.g., *J* = 1–0 at λ = 2.6 mm)
observable at low temperatures (~5 K), although its usage
requires significant assumptions (e.g., S. C. Madden et al.
1997; I. A. Grenier et al. 2005; M. G. Wolfire et al. 2010;
K. E. Jameson et al. 2018; S. C. Madden et al. 2020). Within
star-forming regions, UV radiation from young stars photo-
ionizes nearby hydrogen, which recombines and de-excites,
producing Hα radiation in the process (λ = 656.28 nm).

Dust is well mixed with atomic and molecular gas (F. Boulanger et al. 1996; B. L. Rachford et al. 2009), meaning dust can be used to trace the combined atomic and molecular gas column densities along lines of sight. Dust can either be observed through extinction, as we do in this paper via the BEAST, or through emission in the far-infrared (FIR) as dust particles are heated by the interstellar radiation field/ nearby stars and re-emit energy in the FIR/submillimeter regime (K. D. Gordon et al. 2014).

In the following Section, we look at how three separate ISM
tracers correlate with dust extinction across the fields: atomic
gas via H I emission, ionized gas via Hα emission, and dust via
FIR emission. We chose to not make comparisons with
molecular gas tracers since, to date, there are no surveys that
observe the entire LMC or SMC. These tracers were all
observed at different resolutions, so we opt to convolve all
emission measurements to the resolution of a single HST Scylla
field: ~160" for fields with only optical and UV coverage with
WFC3 UVIS, or ~130" for fields with IR coverage from
WFC3 IR; and only compare with BEAST extinction values
averaged across each field. In Figure 19, we show how each of
these tracers correlates with the median BEAST-derived
extinction measurements.

*Atomic Hydrogen.* To trace atomic hydrogen (H I) across the SMC, we use 21 cm emission from the recent GASKAP-HI collaboration (N. M. Pingel et al. 2022). These H I maps were constructed using new observations from the ASKAP telescope (A. W. Hotan et al. 2021) combined with single-dish observations from the GASS survey (N. M. McClure-Griffiths et al. 2009; P. M. W. Kalberla et al. 2010; P. M. W. Kalberla & U. Haud 2015), producing H I emission data cubes with an angular resolution of 30", or ~10 pc assuming a distance of 50–60 kpc. For the LMC, we use combined H I emission from the Australia Telescope Compact Array and Parkes receiver (S. Kim et al. 2003) with an angular resolution of 60". We calculate the median H I emission across each Scylla field.

In Figure 19(a), we find a linear correlation between H I and
extinction, up until $A_V > 0.8$ mag, after which H I plateaus in
LMC fields. These fields with higher average extinction are all
located within the 30 Doradus star-forming region of the LMC
and are likely dominated by molecular gas. Other extinction
surveys find similar trends when comparing $A_V$ and H I-based
$A_V$ measurements at 40" resolution in the LMC (see Figure 12
in Y. Choi et al. 2018).

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Extinction vs. ISM tracers, median BEAST $A_V$ compared to other ISM tracers per Scylla field: (a) atomic hydrogen ($N(H	ext{i})$; S. Kim et al. 2003; N. M. Pingel et al. 2022), (b) ionized hydrogen (Hα; J. E. Gaustad et al. 2001), and (c) FIR dust emission ($\Sigma_D$; C. J. R. Clark et al. 2021). Uncertainties denote the standard deviation of extinction or tracer values found within each field. Panel (c) shows the theoretical dust surface density as a function of $A_V$ based on Equation (7) in B. T. Draine et al. (2014; black dashed line), and the same equation but with a 2.5 correction factor as noted in J. J. Dalcanton et al. (2015; gray dashed line). A 0.1 mag offset is included along the x-axis to account for foreground extinction.</figcaption></figure>

In the Milky Way, the transition from H I-to-H₂ was
observed to occur at around $A_V = 0.3$ mag (H. Liszt 2014),
meaning, for fields below this transition threshold, the
relationship between H I column density and $A_V$ should become linear. Beyond this threshold, the relationship should flatten as $A_V$ continues to trace the total column density of neutral gas with contributions from H$_2$, while H I remains constant. The transition threshold from H I-to-H$_2$ should be independent of metallicity, so this inconsistency between the Milky Way and the LMC is likely due to under-resolving the molecular regions within Scylla fields in this analysis. Further work at subfield resolutions is needed to resolve these inconsistencies.

*Ionized Hydrogen*. Hα broadly traces the amount of ionizing radiation from young stars, and should therefore trace regions of recent star formation. To trace warm ionized interstellar gas, we use Hα emission from the Southern H-Alpha Sky Survey Atlas (SHASSA; J. E. Gaustad et al. 2001), which observed the southern hemisphere sky using the Cerro Tololo Inter-American Observatory (CTIO) at 48" pixel-resolution, or ~13 pc assuming a distance of 50–60 kpc. Similar to H I, we calculate the median Hα emission for each Scylla field.

Figure 19(b) shows a logarithmic correlation between $A_V$ and Hα emission for fields in the LMC, whereas in the SMC, there is very little correlation. Fields with higher levels of extinction tend to be correlated with higher levels of Hα emission in the LMC. This is not surprising given that these high-extinction fields are all located in 30 Doradus, the most active star-forming region in the Local Group (E. Sabbi et al. 2013). Perhaps more surprising is the fact that fields with $A_V < 0.5$ mag span several orders of magnitude in Hα emission, an indication that either: (1) dustier star-forming structures constitute a relatively small area within any given field and are therefore not contributing significantly to the median extinction; (2) any dusty surrounding material has been destroyed by the ionizing radiation (e.g., J. M. D. Kruijssen et al. 2014). Once again, research on subfield resolutions will be necessary to constrain the two scenarios.

*FIR Dust Emission*. To compare dust extinction with FIR dust emission-derived dust surface densities, we use a new combination of FIR emission observations of the MCs (C. J. R. Clark et al. 2021). These maps were constructed using re-reduced high-resolution Herschel data from the Herschel Inventory of the Agents of Galaxy Evolution (M. Meixner et al. 2013) combined with lower-resolution data from all-sky surveys like Planck (Planck Collaboration et al. 2011), the Infrared Astronomical Satellite (G. Neugebauer et al. 1984), and the Cosmic Background Explorer (N. W. Boggess et al. 1992; R. F. Silverberg et al. 1993), which capture extended dust emission.

In Figure 19(c), we plot dust surface density ($\Sigma_D$) as a function of the median $A_V$ for each Scylla field. We compare this correlation with previous estimates of FIR-derived dust surface densities as described by Equation (7) in B. T. Draine et al. (2014):

$$A_V = 0.74 \left( \frac{\Sigma_D}{10^{-1} M_{\odot} \text{ pc}^{-2}} \right) \text{mag.}$$

We calculate the theoretical dust surface density based on the range of observed BEAST $A_V$ and overlay the results as a dashed black line, including a 0.1 mag offset along the x-axis to account for potential foreground extinction. Although the relation between $A_V$ and $\Sigma_D$ is similar, there is a notable offset. This is not the first instance of discrepancy found: J. J. Dalcanton et al. (2015) found a factor of ~2.5 difference between their extinction inferred from CMD fitting versus dust

surface densities from FIR dust emission modeling in M31. Similarly, Planck Collaboration et al. (2016) found a factor of ~2 factor difference between B. T. Draine et al. (2007) extinction estimates versus quasar $A_V$ measurements, and P. Yanchulova Merica-Jones et al. (2021) found a factor of 1.8 difference between CMD extinction measurements in the SMC versus dust surface densities derived measured from FIR emission modeling. When we apply a correction factor of 2.5 to $\Sigma_D$ (gray dashed line), the model is now within $1\sigma$ of the observed emission, indicating that BEAST $A_V$ estimates are in fact in agreement with existing FIR-derived dust surface density measurements.

For all ancillary ISM tracers, we find positive correlations with the median BEAST-derived extinction measurements in the MCs. However, in the SMC, these correlations are much less pronounced owing to the smaller amounts of dust in the SMC and resolution effects. For both galaxies, these correlations indicate that BEAST-derived extinction measurements are broadly capturing the average ISM content in each Scylla field. Future Scylla publications will expand this investigation by comparing the BEAST-derived extinction parameters ($A_V$, $R_V$, and the presence of the 2175 Å bump) at the original resolutions of the ancillary ISM tracers, and use extinction to quantify the fraction of CO-dark molecular gas as a function of metallicity between the LMC and SMC. On their own, extinction measurements can be used to construct high spatial resolution extinction maps.

## 7. Summary

In this work, we present stellar and line-of-sight dust fits to over 1.5 million sources in the MCs. From these initial fits, we define quality cuts to extract 550,000 high-reliability sources in the SMC and LMC (257,037 and 309,359 sources, respectively). Sources were observed in multiple HST filters spanning from the UV to the IR as part of the Scylla survey (C. Murray et al. 2024). Our data consists of 33 fields in the SMC and 56 fields in the LMC observed with HST WFC3 as part of the Scylla (C. Murray et al. 2024) and METAL (J. Roman-Duval et al. 2019) surveys. Scylla fields were observed with filters in the following priority order when scheduling allowed: F475W, F814W, F336W, F275W, F110W, F160W, and F225W. We only fit fields with at least three-filter coverage.

We use the BEAST (K. D. Gordon et al. 2016) to characterize the multiband SEDs (three- to seven-filter coverage), providing us with probabilistic estimates of initial mass, age, metallicity, distance, extinction ($A_V$), total-to-selective extinction ($R_V$), and extinction curve (Milky Way-like or SMC-like without 2175 Å) for each individual star. We compare these measurements to existing literature on the stellar content of the MCs, their SFH, and the ISM. Our main findings are as follows:

1. *Stellar Mass*. We observe stars with initial masses as low as $0.6 M_{\odot}$, and find that initial mass estimates are in agreement with PARSEC stellar evolutionary tracks (A. Bressan et al. 2012; see also Section 6.1.1; Figures 12, 13, and 14).

2. *Binary Stars*. Binary stars are a potential source of high-extinction misclassifications, especially for low-mass sources along the main sequence where binary models populate the same CMD space as 1 mag extinction.
However, even if these potential binary sources were misclassified as highly extinguished, they would only increase the average extinction in each galaxy by ~0.05 mag (Section 6.1.2; Figure 15).

3. *Star Formation History.* BEAST stellar age distributions are in strong agreement with SFH measurements for the past 6 Gyr (R. E. Cohen et al. 2024a, 2024b). The majority of observed stars are fit to be younger than 6–7 Gyr, with 77.5% (93.9%) of stars in the SMC and 57.1% (78.8%) of stars in the LMC being younger than 6.0 (7.0) Gyr, respectively (Section 6.2; Figure 16).

4. *Extinction Limits.* Most upper-main-sequence and RGB/RC stars are detectable with up to 3 mag of extinction. Subsolar mass stars are not detectable with more than ~0.5 mag of extinction (Section 6.3.2; Figure 17).

5. *Extinction Distribution.* Extinction distributions for all sources in both the SMC and LMC follow log-normal functions, with $\mu = 0.33$ and $\sigma = 0.66$ (median $A_V = 0.29$ mag) in the SMC, and $\mu = 0.44$ and $\sigma = 0.76$ (median $A_V = 0.40$ mag) in the LMC (Section 6.3.3; Figure 18).

6. *Gas Tracers.* Field-averaged BEAST $A_V$ measurements are positively correlated with other ISM tracers like atomic hydrogen (H I) and ionized hydrogen (H$\alpha$). $A_V$ is linearly correlated with H I until ~0.8 mag. Above this limit, observed fields are mainly located in the 30 Doradus region of the LMC, meaning the gas content is likely predominantly in the molecular phase. $A_V$ is logarithmically correlated with H$\alpha$ emission, but low-$A_V$ fields ($A_V < 0.5$ mag) span several orders of magnitude in H$\alpha$ emission ($10^2-10^5$ Rayleigh; Section 6.3.4; Figure 19).

7. *FIR-based Dust Surface Densities.* We convert observed BEAST $A_V$ measurements to predicted dust surface densities using the B. T. Draine et al. (2014) model and find that the $A_V$-based dust surface densities are a factor of ~2.5 lower than observed FIR-based dust surface densities. This correction factor is similar to other studies, both in the SMC (e.g., P. Yanchulova Merica-Jones et al. 2021), and other galaxies like M31 (e.g., J. J. Dalcanton et al. 2015; see also Section 6.3.4; Figure 19).

8. *Future Work.* The results presented in this paper are predominantly focused on galaxy-/field-scale statistics of BEAST fits. In upcoming Scylla ISM publications, we plan to focus on subfield scales by: (1) constructing parsec-scale extinction maps; (2) investigating the source

of subfield variations in $R_V$; and (3) constraining the fraction of CO-dark molecular gas in the SMC and LMC.

*Facilities:* HST (WFC3/IR), HST (WFC3/UVIS), XSEDE (J. Towns et al. 2014).

*Software:* DOLPHOT (A. E. Dolphin 2002), BEAST (K. D. Gordon et al. 2016), astropy (Astropy Collaboration et al. 2013, 2018), numpy (C. R. Harris et al. 2020), scipy (P. Virtanen et al. 2020), matplotlib (J. D. Hunter 2007).

## Acknowledgments

This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs 15891, 16235, and 16786. This research has made use of NASA's Astrophysics Data System. All of the data presented in this paper were obtained from MAST at the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909/8ads-wn75. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts.

This project made use of data from the Southern H-Alpha Sky Survey Atlas (SHASSA), which is supported by the National Science Foundation

This work used Bridges-2 at Pittsburgh Supercomputing Center (PSC) through allocation \#AST190023P from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grant Nos. 2138259, 2138286, 2138307, 2137603, and 2138296. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. ACI-1548562. Specifically, it used the Bridges-2 system, which is supported by NSF award No. ACI-1928147, at the PSC.

## Appendix
UV and IR BEAST CMDs

In Figure 20, we show the UV and IR versions of the optical CMDs in Figures 12 and 13.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** UV and IR versions of Figures 12 and 13: UV (F275W and F336W) and IR CMDs (F110W and F160W) for all high-reliability sources in the SMC (top; $n_{UV}$ = 139, 551 and $n_{IR}$ = 59, 231) and LMC (bottom; $n_{UV}$ = 188, 555 and $n_{IR}$ = 137, 023). Left columns (a, b): density distribution of sources as colored 2D histograms, and outlined in logarithmic levels using kernel-density estimates (KDEs). KDE density outlines (white) are propagated to all other columns. All other columns (c-j): spatially binned (100 × 100) medians for notable BEAST parameters: initial mass ($M_{ini}$), age, metallicity (Z), and extinction ($A_V$).</figcaption></figure>

## ORCID iDs

Christina W. Lindberg https://orcid.org/0000-0003-0588-7360

Claire E. Murray https://orcid.org/0000-0002-7743-8129

Petia Yanchulova Merica-Jones https://orcid.org/0000-0002-9912-6046

Caroline Bot https://orcid.org/0000-0001-6118-2985

Clare Burhenne https://orcid.org/0009-0005-0339-015X

Yumi Choi https://orcid.org/0000-0003-1680-1884

Christopher J. R. Clark https://orcid.org/0000-0001-7959-4902

Roger E. Cohen https://orcid.org/0000-0002-2970-7435

Karoline M. Gilbert https://orcid.org/0000-0003-0394-8377

Steven R. Goldman https://orcid.org/0000-0002-8937-3844

Karl D. Gordon https://orcid.org/0000-0001-5340-6774

Alec S. Hirschauer https://orcid.org/0000-0002-2954-8622

Kristen B. W. McQuinn https://orcid.org/0000-0001-5538-2614

Julia C. Roman-Duval https://orcid.org/0000-0001-6326-7069

Karin M. Sandstrom https://orcid.org/0000-0002-4378-8534

Elizabeth Tarantino https://orcid.org/0000-0003-1356-1096

Benjamin F. Williams https://orcid.org/0000-0002-7502-0597

## References

Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bailer-Jones, C. A. L. 2011, MNRAS, 411, 435

Bekki, K., & Stanimirović, S. 2009, MNRAS, 395, 342

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668, 949

Boggess, N. W., Mather, J. C., Weiss, R., et al. 1992, ApJ, 397, 420

Bolatto, A. D., Leroy, A. K., Jameson, K., et al. 2011, ApJ, 741, 12

Boulanger, F., Abergel, A., Bernard, J. P., et al. 1996, A&A, 312, 256

Boyer, M. L., Srinivasan, S., Riebel, D., et al. 2012, ApJ, 748, 40

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
Castelli, F., & Kurucz, R. L. 2003, in IAU Symp. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray (Cambridge: Cambridge Univ. Press), A20

Chen, B. Q., Guo, H. L., Gao, J., et al. 2022, MNRAS, 511, 1317

Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068

Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525

Choi, Y., Dalcanton, J. J., Williams, B. F., et al. 2020, ApJ, 902, 54

Choi, Y., Nidever, D. L., Olsen, K., et al. 2018, ApJ, 866, 90

Choi, Y., Olsen, K. A. G., Besla, G., et al. 2022, ApJ, 927, 153

Cignoni, M., Cole, A. A., Tosi, M., et al. 2013, ApJ, 775, 83

Cignoni, M., Sacchi, E., Tosi, M., et al. 2019, ApJ, 887, 112

Clark, C. J. R., Roman-Duval, J. C., Gordon, K. D., Bot, C., & Smith, M. W. L. 2021, ApJ, 921, 35

Clark, C. J. R., Roman-Duval, J. C., Gordon, K. D., et al. 2023, ApJ, 946, 42

Clayton, G. C., Gordon, K. D., Bianchi, L. C., et al. 2015, ApJ, 815, 14

Cohen, R. E., McQuinn, K. B. W., Murray, C. E., et al. 2024a, ApJ, 975, 42

Cohen, R. E., McQuinn, K. B. W., Murray, C. E., et al. 2024b, ApJ, 975, 43

Conroy, C. 2013, ARA&A, 51, 393

Crowther, P. A. 2019, Galax, 7, 88

Cullinane, L. R., Mackey, A. D., Da Costa, G. S., et al. 2022, MNRAS, 510, 445

Dalcanton, J. J., Bell, E. F., Choi, Y., et al. 2023, AJ, 166, 80

Dalcanton, J. J., Fouesneau, M., Hogg, D. W., et al. 2015, ApJ, 814, 3

Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18

de Grijs, R., & Bono, G. 2015, AJ, 149, 179

de Grijs, R., Wicker, J. E., & Bono, G. 2014, AJ, 147, 122

Dolphin, A., 2016 DOLPHOT: Stellar photometry, Astrophysics Source Code Library, ascl:1608.013

Dolphin, A. E. 2002, MNRAS, 332, 91

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Evans, C., Marcolini, W., Bouret, J.-C., & Garcia, M. 2023, ExA, 56, 537

Fitzpatrick, E. L. 1999, PASP, 111, 63

Fitzpatrick, E. L., Massa, D., Gordon, K. D., Bohlin, R., & Clayton, G. C. 2019, ApJ, 886, 108

Fox, A. J., Barger, K. A., Wakker, B. P., et al. 2018, ApJ, 854, 142

Gaia Collaboration, Helmi, A., van Leeuwen, F., et al. 2018, A&A, 616, A12

Gaustad, J. E., McCullough, P. R., Rosing, W., & Van Buren, D. 2001, PASP, 113, 1326

Gilbert, K. M., Choi, Y., Boyer, M. L., et al. 2025, ApJS, 276, 8

Girardi, L. 2016, ARA&A, 54, 95

Girardi, L., Groenewegen, M. A. T., Hatziminaoglou, E., & da Costa, L. 2005, A&A, 436, 895

Girardi, L., & Marigo, P. 2007, in ASP Conf. Ser. 378, Why Galaxies Care About AGB Stars: Their Importance as Actors and Probes, ed. F. Kerschbaum, C. Charbonnel, & R. F. Wing (San Francisco, CA: ASP), 20

Goldman, S. R., Boyer, M. L., Dalcanton, J., et al. 2022, ApJS, 259, 41

Gordon, K. D., Clayton, G. C., Decleir, M., et al. 2023, ApJ, 950, 86

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Gordon, K. D., Fitzpatrick, E. L., Massa, D., et al. 2024, ApJ, 970, 51

Gordon, K. D., Fouesneau, M., Arab, H., et al. 2016, ApJ, 826, 104

Gordon, K. D., Roman-Duval, J., Bot, C., et al. 2014, ApJ, 797, 85

Gouliermis, D. A., Schmeja, S., Dolphin, A. E., et al. 2012, ApJ, 748, 64

Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2014, ApJ, 783, 114

Gregersen, D., Seth, A. C., Williams, B. F., et al. 2015, AJ, 150, 189

Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005, Sci, 307, 1292

Groenewegen, M. A. T., Girardi, L., Hatziminaoglou, E., et al. 2002, A&A, 392, 741

Gull, M., Weisz, D. R., Senchyna, P., et al. 2022, ApJ, 941, 206

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hawcroft, C., Sana, H., Mahy, L., et al. 2024, A&A, 688, A105

Hotan, A. W., Bunton, J. D., Chippendale, A. P., et al. 2021, PASA, 38, e009

Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605

Hunter, J. D. 2007, CSE, 9, 90

Jameson, K. E., Bolatto, A. D., Leroy, A. K., et al. 2016, ApJ, 825, 12

Jameson, K. E., Bolatto, A. D., Wolfire, M., et al. 2018, ApJ, 853, 111

Johnson, L. C., Seth, A. C., Dalcanton, J. J., et al. 2017, ApJ, 839, 78

Jones, O. C., Nally, C., Habel, N., et al. 2023, NatAs, 7, 694

Kalberla, P. M. W., & Haud, U. 2015, A&A, 578, A78

Kalberla, P. M. W., McClure-Griffiths, N. M., Pisano, D. J., et al. 2010, A&A, 521, A17

Kim, S., Dopita, M. A., Staveley-Smith, L., & Bessell, M. S. 1999, AJ, 118, 2797

Kim, S., Staveley-Smith, L., Dopita, M. A., et al. 2003, ApJS, 148, 473

Kroupa, P. 2001, MNRAS, 322, 231

Kruijssen, J. M. D., Longmore, S. N., Elmegreen, B. G., et al. 2014, MNRAS, 440, 3370

Lanz, T., & Hubeny, I. 2003, ApJS, 146, 417

Lanz, T., & Hubeny, I. 2007, ApJS, 169, 83

Lee, J.-K., Rolleston, W. R. J., Dufton, P. L., & Ryans, R. S. I. 2005, A&A, 429, 1025

Leone, F., Scalia, C., Gangi, M., et al. 2017, ApJ, 848, 107

Leroy, A., Bolatto, A., Stanimirović, S., et al. 2007, ApJ, 658, 1027

Lewis, A. R., Dolphin, A. E., Dalcanton, J. J., et al. 2015, ApJ, 805, 183

Lindberg, C. W., Murray, C. E., Dalcanton, J. J., Peek, J. E. G., & Gordon, K. D. 2024, ApJ, 963, 58

Liszt, H. 2014, ApJ, 780, 10

Lombardi, M., Alves, J., & Lada, C. J. 2015, A&A, 576, L1

Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141

Madden, S. C., Poglitsch, A., Geis, N., Stacey, G. J., & Townes, C. H. 1997, ApJ, 483, 200

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77

Mastropietro, C., Burkert, A., & Moore, B. 2009, MNRAS, 399, 204

Mazzi, A., Girardi, L., Zaggia, S., et al. 2021, MNRAS, 508, 245

McClure-Griffiths, N. M., Pisano, D. J., Calabretta, M. R., et al. 2009, ApJS, 181, 398

McQuinn, K. B. W., Skillman, E. D., Cannon, J. M., et al. 2010, ApJ, 721, 297

McQuinn, K. B. W., Skillman, E. D., Dolphin, A. E., & Mitchell, N. P. 2015, ApJ, 808, 109

Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 62

Meng, X., & Gnedin, O. Y. 2021, MNRAS, 502, 1433

Misselt, K. A., Clayton, G. C., & Gordon, K. D. 1999, ApJ, 515, 128

Murray, C., Lindberg, C. W., Yanchulova Merica-Jones, P., et al. 2024, ApJS, 275, 5

Murray, C. E., Hasselquist, S., Peek, J. E. G., et al. 2024, ApJ, 962, 120

Nataf, D. M., Cassisi, S., Casagrande, L., Yuan, W., & Riess, A. G. 2021, ApJ, 910, 121

Nayak, O., Nally, C., Hirschauer, A. S., et al. 2024, ApJ, 975, 262

Ness, M., Hogg, D. W., Rix, H. W., Ho, A. Y. Q., & Zasowski, G. 2015, ApJ, 808, 16

Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJL, 278, L1

Nidever, D. L., Majewski, S. R., & Butler Burton, W. 2008, ApJ, 679, 432

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Patel, E., Kallivayalil, N., Garavito-Camargo, N., et al. 2020, ApJ, 893, 121

Patra, N. N. 2019, MNRAS, 484, 81

Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Natur, 567, 200

Pingel, N. M., Dempsey, J., McClure-Griffiths, N. M., et al. 2022, PASA, 39, e005

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, A1

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A132

Povick, J. T., Nidever, D. L., Massana, P., et al. 2024, MNRAS, 533, 3685

Rachford, B. L., Snow, T. P., Destree, J. D., et al. 2009, ApJS, 180, 125

Ramachandran, V., Sander, A. A. C., Pauli, D., et al. 2024, A&A, 692, A90

Roman-Duval, J., Bot, C., Chastenet, J., & Gordon, K. 2017, ApJ, 841, 72

Roman-Duval, J., Gordon, K. D., Meixner, M., et al. 2014, ApJ, 797, 86

Roman-Duval, J., Jenkins, E. B., Tchernyshyov, K., et al. 2021, ApJ, 910, 95

Roman-Duval, J., Jenkins, E. B., Tchernyshyov, K., et al. 2022, ApJ, 928, 90

Roman-Duval, J., Jenkins, E. B., Williams, B., et al. 2019, ApJ, 871, 151

Russell, S. C., & Dopita, M. A. 1992, ApJ, 384, 508

Sabbi, E., Anderson, J., Lennon, D. J., et al. 2013, AJ, 146, 53

Sabbi, E., Lennon, D. J., Anderson, J., et al. 2016, ApJS, 222, 11

Schneider, N., Bontemps, S., Girichidis, P., et al. 2015, MNRAS, 453, L41

Schönrich, R., & Bergemann, M. 2014, MNRAS, 443, 698

Schwering, P. B. W. 1988, PhD thesis, Netherlands Univ. Leiden

Scowcroft, V., Freedman, W. L., Madore, B. F., et al. 2016, ApJ, 816, 49

Silverberg, R. F., Hauser, M. G., Boggess, N. W., et al. 1993, Proc. SPIE, 2019, 180

Staveley-Smith, L., Sault, R. J., Hatzidimitriou, D., Kesteven, M. J., & McConnell, D. 1997, MNRAS, 289, 225

Subramanian, S., Rubele, S., Sun, N.-C., et al. 2017, MNRAS, 467, 2980

Subramanian, S., & Subramaniam, A. 2009, A&A, 496, 399

Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287

Tarantino, E., Bolatto, A. D., Indebetouw, R., et al. 2024, ApJ, 969, 101

Telford, O. G., Chisholm, J., Sander, A. A. C., et al. 2024, ApJ, 974, 85

Telford, O. G., McQuinn, K. B. W., Chisholm, J., & Berg, D. A. 2023, ApJ, 943, 65

Ting, Y.-S., Conroy, C., Rix, H.-W., & Cargile, P. 2019, ApJ, 879, 69

Towns, J., Cockerill, T., Dahan, M., et al. 2014, CSE, 16, 62
Tsuge, K., Sano, H., Tachihara, K., et al. 2019, ApJ, 871, 44

Van De Putte, D., Gordon, K. D., Roman-Duval, J., et al. 2020, ApJ, 888, 22

Vanhollebeke, E., Groenewegen, M. A. T., & Girardi, L. 2009, A&A, 498, 95

Vazquez-Semadeni, E. 1994, ApJ, 423, 681

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2013, MNRAS, 431, 364

Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, ApJ, 744, 44

Williams, B. F., Bell, E. F., Boyer, M. L., et al. 2021a, HST Proposal, 29, 16778

Williams, B. F., Dalcanton, J. J., Bell, E. F., et al. 2012, ApJ, 759, 46

Williams, B. F., Dolphin, A. E., Dalcanton, J. J., et al. 2017, ApJ, 846, 145

Williams, B. F., Durbin, M. J., Dalcanton, J. J., et al. 2021b, ApJS, 253, 53

Williams, B. F., Lang, D., Dalcanton, J. J., et al. 2014, ApJS, 215, 9

Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191

Yanchulova Merica-Jones, P., Gordon, K., Sandstrom, K., et al. 2025, ApJ, 978, 144

Yanchulova Merica-Jones, P., Sandstrom, K. M., Johnson, L. C., et al. 2017, ApJ, 847, 102

Yanchulova Merica-Jones, P., Sandstrom, K. M., Johnson, L. C., et al. 2021, ApJ, 907, 50

Zaritsky, D., Harris, J., Thompson, I. B., & Grebel, E. K. 2004, AJ, 128, 1606

Zaritsky, D., Harris, J., Thompson, I. B., Grebel, E. K., & Massey, P. 2002, AJ, 123, 855

Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198