# Simultaneous Estimates of Star-cluster Age, Metallicity, Mass, and
Extinction (SESAMME) I: Presenting an MCMC Approach to Spectral Stellar
Population Fitting

Logan H. Jones Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Svea Hernandez AURA for ESA, Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Linda J. Smith
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218,
USA  Bethan L. James AURA for ESA, Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Alessandra Aloisi Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Søren
Larsen Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500
GL Nijmegen, The Netherlands

(Received 5/16/23; Accepted 9/28/23)

###### Abstract

We present the first version release of SESAMME, a public, Python-based full
spectrum fitting tool for Simultaneous Estimates of Star-cluster Age,
Metallicity, Mass, and Extinction. SESAMME compares an input spectrum of a
star cluster to a grid of stellar population models with an added nebular
continuum component, using Markov chain Monte Carlo (MCMC) methods to sample
the posterior probability distribution in four dimensions: cluster age,
stellar metallicity $Z$, reddening $E(B-V)$, and a normalization parameter
equivalent to a cluster mass. SESAMME is highly flexible in the stellar
population models that it can use to model a spectrum; our testing and initial
science applications use both BPASS and Starburst99. We illustrate the ability
of SESAMME to recover accurate ages and metallicities even at a moderate
signal-to-noise ratio (S/N $\sim 3-5$ per wavelength bin) using synthetic,
noise-added model spectra of young star clusters. Finally, we test the
consistency of SESAMME with other age and metallicity estimates from the
literature using a sample of HST/COS far-UV spectra towards young, massive
clusters in M83 and NGC 1313. We find that, on the whole, SESAMME infers star
cluster properties that are consistent with the literature in both low- and
high-metallicity environments.

Star clusters (1567); Ultraviolet spectroscopy (2284); Astronomy software
(1855)

††journal: ApJ

##  1 Introduction

An accurate understanding of the current and past lives of galaxies requires a
complete census of their elemental abundances. Innumerable works have examined
the global and spatially-resolved metallicities of star-forming galaxies
(SFGs), as well as metallicity gradients, usually through the analysis of
emission lines from H ii regions (e.g., Maiolino & Mannucci, 2019, and
references therein). These studies have also led to the development of global
relations between a galaxy’s metal content and other physical characteristics
(e.g., the mass-metallicity relation or MZR; Zaritsky et al., 1994; Tremonti
et al., 2004). Taken together, these relations and metallicity gradient
estimates are frequently used to investigate a galaxy’s star-formation and
merger history, the impact of galactic winds, and the flow of pristine gas
from the intergalactic medium onto galaxies (e.g., Ho et al., 2015; Kudritzki
et al., 2015).

However, relying solely on the emission from ionized gas to characterize the
chemical composition of SFGs may result in an incomplete picture of their
metal content. The nebular emission from H ii regions, for example, can only
be sustained for the $\sim$10 Myr period in which the most massive and UV-
bright stars can photoionize their surrounding gas. As such, the metallicity
of this gas is associated more with the most recent episodes of star formation
and may be more enriched than the surrounding interstellar medium (ISM; Kunth
& Sargent, 1986; Lebouteiller et al., 2013). Alternative perspectives on the
metal content of SFGs can be obtained through analyses of individual stars or
star clusters, which retain some memory of the chemical composition of their
birth cloud. For very nearby galaxies ($D\sim$a few Mpc), measurements of
stellar metallicities have been achieved using spectroscopic observations of
blue and red supergiants (e.g., Davies et al., 2010, 2015, 2017; Kudritzki et
al., 2012, 2014; Gazak et al., 2014; Hosek et al., 2014; Bresolin et al.,
2016), while other works use integrated-light (IL) spectroscopy of star
clusters and stellar populations (e.g., Halliday et al., 2008; McWilliam &
Bernstein, 2008; Colucci et al., 2011, 2012; Larsen et al., 2012; Steidel et
al., 2016; Chisholm et al., 2019).

A spectral fitting technique developed by Larsen et al. (2012, hereafter L12)
has been widely used in the literature to model the integrated light (IL) of
star clusters (Larsen et al., 2012, 2018; Hernandez et al., 2017, 2019, 2021,
2022; Asa’d et al., 2022; Gvozdenko et al., 2022). In brief, the method
involves first deriving stellar parameters of individual stars via color-
magnitude diagrams or, more typically, theoretical isochrones of a given age
and metallicity, from which one can synthesize stellar atmosphere models and
mock spectra for each star. The sum of these synthetic spectra creates an IL
spectrum of the cluster, which is then compared to the data. This process is
iterated until an optimal fit is obtained for the individual elemental
abundances across predefined wavelength bins. The L12 technique was originally
developed to study the properties of extragalactic star clusters with very
high spectral resolution ($R\sim$30,000) optical observations, but has since
been successfully applied to spectra at lower resolutions $R\sim 8000$
(Hernandez et al., 2017) and at both infrared (Larsen et al., 2018) and UV
wavelengths (Hernandez et al., 2019).

However, Hernandez et al. (2019) report a caveat to the L12 process that
becomes especially prevalent when dealing with UV observations. When applying
the L12 method to HST/COS FUV spectra of a sample of young massive clusters
(YMCs) in M83, they find that varying the assumed age of the input isochrone
can cause significant changes in the inferred stellar metallicities. For
clusters with relatively uncertain ages from, e.g., spectral energy
distribution (SED) fitting, dialing the assumed age up or down by the
$1\sigma$ uncertainties induced changes in their inferred metallicities of
order $\sim$0.2 – 0.3 dex on average and as high as 0.4 dex in individual
clusters. Additionally, the scatter of metallicity measurements from different
wavelength bins was consistently higher when assuming the $-1\sigma$ age. They
attribute these effects to the rapid evolution in the UV properties of massive
star clusters in the first 10 Myr of their lives, reinforcing the need for
precise age estimates for very young star clusters.

In this work, we introduce a complement to the L12 method in the form of a
Python tool for Simultaneous Estimates of Star-cluster Age, Metallicity, Mass,
and Extinction (SESAMME), which we release for public
use111https://github.com/astrolojo/SESAMME. We describe the structure and
computational performance of the code in Section 2, then apply the code to a
suite of synthetic FUV observations of YMCs in Section 3. Finally, we further
validate the usefulness of the code to real data by analyzing a small sample
of four YMCs with published estimates for their ages and stellar
metallicities, and show that SESAMME can reliably reproduce past measurements
using independent methods. We discuss our validation tests and immediate
future uses for the code in Section 5.

##  2 The Algorithm

Given an input spectrum of a star cluster and a set of model spectra for
simple stellar populations (SSPs), SESAMME samples the posterior probability
distribution function (PDF) for a four-parameter model to return constraints
on the physical and observable properties of the cluster. We further describe
the necessary inputs and resulting outputs in Sections 2.2 and 2.3.

The first two model parameters – cluster age and stellar metallicity mass
fraction $Z$ – typically have discrete values determined by the suite of
stellar population models used in a SESAMME run. Boundaries on the age and
metallicity priors can span the entire range used in a model suite (e.g., ages
between 6.0 $<$ log(age) $<$ 11.0 in v2.2 and v2.3 of the Binary Population
and Spectral Synthesis (BPASS) models) or can be restricted if a reasonable
constraint is already known (e.g., age estimates derived from SED fitting).

We note that, in the current implementation of SESAMME, the sampler does not
perform any sort of interpolation over the user-specified SSP model grid when
generating a model to evaluate. Instead, when the sampler chooses a random
value for the cluster age and metallicity, the code will round to the nearest
pre-computed value in the model grid. This choice was motivated by the non-
linear evolution of some spectral features with respect to age and metallicity
in young stellar clusters, which can lead to inaccurate UV line profiles when
using a naive interpolation of a coarser model grid. Additionally, the
uncertainties on the inferred cluster age and metallicity are typically
comparable to or larger than the bin size for each of these parameters (see
Tables 3 and 4), making it unnecessary to interpolate to a finer grid.

Table 1:  List of extinction curves available for SESAMME Type | Source | Default $R_{V}$
---|---|---
MW | Fitzpatrick (1999) | 3.1
| Cardelli et al. (1989) | 3.1
| Fitzpatrick & Massa (2007) | 3.1
| Gordon et al. (2023) | 3.1
LMC | Gordon et al. (2003) | 3.41
SMC | Gordon et al. (2003) | 2.74
Starburst | Calzetti et al. (2000) | 4.05
Table 2: Example List of Stellar Population Models Property |  BPASS v2.3 |  BPASS v2.2 |  Starburst99
---|---|---|---
Metallicity |  [10-5, 10-4, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.01, 0.014, 0.02, 0.03, 0.04] |  [10-5, 10-4, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.01, 0.014, 0.02, 0.03, 0.04] |  [0.001, 0.004, 0.008, 0.02, 0.04]
Age range (yr) |  106 – 1011 |  106 – 1011 |  Variable
Age step |  0.1 dex |  0.1 dex |  Variable
IMF options |  Double power law |  Single and double power laws; power law with exponential cutoff (Chabrier, 2003) |  Double power law
Abundance information |  [$\alpha$/Fe] $\in$ [-0.2, 0.0, 0.2, 0.4, 0.6]. Asplund et al. (2009) Solar abundances. |  [$\alpha$/Fe] = 0.0. Grevesse & Noels (1993) Solar abundances. |  Varies with stellar evolution tracks.
Notes |  Includes models with binaries and with single stars only. |  Includes models with binaries and with single stars only. |  User defines the age range and intervals at which spectra are printed.

The third parameter is the degree of reddening $E(B-V)$, in magnitudes, which
is then translated to extinction $A_{V}$ assuming a user-specified extinction
law. Sampled values of $E(B-V)$ are continuous within a user-specified range.
We use the Python packages dust_extinction by K. Gordon222https://dust-
extinction.readthedocs.io/ (Gordon et al., 2022) and extinction by K.
Barbary333https://extinction.readthedocs.io/ (Barbary, 2021) to implement the
many extinction model options available to SESAMME (see Table 1 for those that
are currently implemented). We note that although some of the extinction
models in these packages (e.g., the MW-type model from Gordon et al. (2023))
allow for variable values of $R_{V}$, SESAMME assumes that $R_{V}$ is fixed
for a given choice of model. The source code for SESAMME can be easily
modified by users who wish to implement additional extinction models or assume
alternative values of $R_{V}$ for existing models.

The final parameter is a continuous rescaling parameter log($A$), which can
account for the difference between the assumed stellar mass of the models
(e.g., initial mass $10^{6}M_{\odot}$ in BPASS) and the actual stellar mass of
the cluster under inspection, but only if the flux array of the spectrum is in
physical units. That is, if the input spectrum has been relatively flux
calibrated, such that the overall continuum shape and line profiles are
correct but the flux array is given in arbitrary units, then log($A$) is a
pure unit conversion factor to facilitate comparisons between the data and the
model. However, if i) the flux array is fully calibrated and ii) a distance to
the target is known, such that the spectrum can be converted to luminosity
units before analyzing with SESAMME, then the value of $A$ accounts for the
difference between the assumed mass of the models and the actual mass of the
cluster. his means that the mass estimate is thus degenerate with the assumed
distance to the cluster, but such a degeneracy is not unique to SESAMME. We
note that SESAMME does not explicitly take distance uncertainties into
account, which adds an additional systematic source of uncertainty in the
log($A$) parameter that is not captured by the modeling process. This means
that the errors on log($A$) inferred by SESAMME should be considered lower
limits on the true uncertainty. Finally, we note that mass estimates derived
using SESAMME are also limited by the completeness of the observations. For
example, if a star cluster is more extended than the spectroscopic aperture,
the $A$ value inferred by SESAMME will represent a lower limit in mass.

SESAMME has been explicitly tested using the SSP models from BPASS v2.2
(“Tuatara”; Stanway & Eldridge, 2018), BPASS v2.3 (“Broc”; Byrne et al.,
2022), and Starburst99 v7.0 (Leitherer et al., 1999) on both synthetic (Sec.
3) and empirical (Sec. 4) ultraviolet spectra of star clusters. The
distributed version of SESAMME v1.0 uses BPASS v2.3 with binary stars and a
Solar $\alpha$-element abundance [$\alpha$/Fe] = 0 as the default stellar
population model set. However, we stress that any set of stellar population
models can be used with SESAMME, as long as the model spectra can be formatted
appropriately (see Section 2.2). This flexibility allows users to use the SSP
suite that is preferred for their science goals and/or cross-examine the
impact of different physical assumptions on the age and metallicity
constraints inferred by SESAMME. For example, Starburst99 does not take the
physics of binary stars into account, while BPASS can model populations that
include binary systems, which can significantly affect the time evolution of a
stellar population’s UV luminosity. Similarly, BPASS v2.3 and v2.2 differ in
the assumed stellar atmosphere model grids, the choice of Solar chemical
composition, and the available range of options for $\alpha$ element
enhancement and initial mass function (IMF). We summarize the major
differences between these three SSP model codes in Table 2.

Evaluating the posterior probability of a model in SESAMME is straightforward
and follows standard practices set out by the authors of emcee (Foreman-Mackey
et al., 2013). Priors are uniform within boundaries set by the user before the
modeling run. The logarithm of the model’s likelihood is computed as

$-\frac{1}{2}\sum_{n}\frac{(y_{n}-m_{n})^{2}}{\sigma^{2}_{n}}+ln(2\pi\sigma^{2}_{n}),$

where $y_{n}$ and $\sigma_{n}$ are the flux and flux uncertainty arrays of the
observed spectrum at wavelength bin $n$, and $m_{n}$ is the emission model
(consisting of the SSP spectrum and an optional nebular continuum component –
see below) generated by the sampler in the same bin. The function that
calculates the log-probability also requires a pixel mask (see Section 2.2)
that removes undesirable wavelength bins from the computation.

###  2.1 Nebular Continuum

When generating samples from the pre-generated SSP model cube (described in
Section 2.2), SESAMME can also add in a nebular continuum component. Nebular
continuum emission can account for a substantial fraction of emission along
the line of sight to a star cluster, particularly at young ages and low
metallicities. In Section 3.3, we explore how the exclusion of nebular
continuum emission in the models affects the inferred properties and their
uncertainties.

When nebular continuum emission is “turned on,” the continuum is generated
with each sampling, simultaneous with the selection of a stellar population
model. Our approach is similar to the subroutine continuum in Starburst99. We
use the wavelength-dependent emission coefficients given in Ferland (1980) and
Aller (1984), which together account for free-free, bound-free, and two-photon
emission from H and He+. These values assume typical Case B conditions
(densities and temperatures of order $\sim$100 cm-3 and $10^{4}$ K,
respectively, as well as an ionizing photon escape fraction of $\sim$0) and a
number ratio He/H = 0.1 for the emitting gas. We note that, in this
formulation, the intensity and shape of the nebular continuum is agnostic to
the gas-phase metallicity.

The strength of the nebular continuum emission is determined by the incident
ionizing radiation field associated with a given stellar population. For BPASS
model spectra, for example, we use the tables of hydrogen-ionizing photon
rates $Q$ that are provided in the v2.2 and v2.3 releases. SESAMME then
interpolates the computed nebular continua to the wavelength grid of the
stellar population models, adds the nebular component to the stellar model,
and finally rescales the total model spectrum by the parameter $A$. We note
that this approach may not be robust for very low mass star clusters, where
the high-mass end of the IMF may be sparsely or stochastically sampled
(Orozco-Duarte et al., 2022). Lastly, we note that some codes for generating
stellar population models already include the nebular continuum in the model
spectra (e.g., Starburst99; Leitherer et al., 1999). In these cases, the
SESAMME user can turn off the nebular continuum feature.

###  2.2 Inputs for SESAMME

Figure 1:  Visualization of the format of a SESAMME model cube, generated from
BPASS v2.3 models with binary physics and a solar $\alpha$/Fe, in the form of
Astropy Tables. Within each slice of the model cube, the first column gives
the wavelength grid (which, in this case, has been truncated to only cover the
rest-frame UV), while the remaining columns give the light output at a fixed
log(age) at each wavelength bin in units of $L_{\odot}$ Å-1. For clarity, only
the first ten age intervals and three of the 13 possible metallicities in
BPASS are shown.

SESAMME uses emcee (Foreman-Mackey et al., 2013) for its statistical machinery
and requires minimal tuning before it can be run. For example, at the start of
every SESAMME run, the user specifies the size of an ensemble of Markov chains
or “walkers,” each of which undergo a semi-random walk in the parameter space.
In general, we recommend a minimum of $\sim$64 walkers, as we find that a
larger ensemble of walkers with moderately long MCMC chains provides better
constraints on the PDF than a small ensemble with long chains.

We also note that, if any of the four parameters have existing measurements
from the literature, then the (flat) priors on the relevant parameters can be
restricted before running SESAMME. For example, when applying SESAMME to IL
spectra of YMCs in M83 in Section 4, we restrict our age and metallicity
priors to a moderately narrow range of values for these parameters, since
estimates of the cluster ages and metallicities are known from Hernandez et
al. (2019). Finally, the user must also specify which extinction curve to
assume throughout the modeling process by setting the parameter use_ext_law,
as well as whether to include a nebular continuum emission component in the
modeling process. We list the available extinction curves in Table 1 and
revisit the topic of nebular continuum emission in Section 3.3.

In addition to setting these pre-run parameters, SESAMME also expects three
files as input:

1) A spectrum in ascii format with three labelled columns for wavelength,
flux, and flux uncertainty. For ease of comparison to models, SESAMME assumes
that the input spectrum has been corrected for i) redshift/radial velocity and
ii) Galactic extinction. The code generally does not require that ISM emission
and absorption features be modeled and removed from the spectrum before
running SESAMME, as these can be ignored during the modeling process through
the use of pixel masks. However, for FUV observations, we recommend that the
user models and subtracts the Ly$\alpha$ absorption profile prior to modeling
with SESAMME; for very low-redshift clusters, this could include a Milky Way
foreground component in addition to Ly$\alpha$ absorption intrinsic to the
target. If the Ly$\alpha$ absorption profile is sufficiently broad, the wings
can suppress the continuum around and the apparent strength of N v
$\lambda\lambda 1238,1242$, which can have small but significant effects on
the resulting age and metallicity posterior distributions.

Additionally, the spectral resolution of the chosen stellar population models
may be either too fine or too coarse, depending on the data the user wishes to
model (e.g., $R\sim 3000$ optical spectra from Keck/DEIMOS vs. $R\sim 20,000$
FUV spectra from HST/COS). We recommend that the user smooths and resamples
either the empirical spectrum or the model set, whichever is finer, to the
coarser wavelength grid to ensure a proper comparison between models and data.
We show an example of this in Section 4, where we use SESAMME to model HST/COS
FUV spectra (native resolution 0.06 Å per resolution element) with BPASS model
spectra (with fluxes given in bins of 1 Å).

2) A single FITS file (hereafter the “model cube”) containing the set of
stellar population models to compare with the empirical spectrum. The model
cube is generated by the user with a supplemental piece of Python code
(included in the SESAMME release) prior to running the main algorithm for
SESAMME. The cube generator collates a set of stellar population models into a
single multi-extension FITS file, where each FITS extension contains all model
spectra for a given metallicity. Each extension is formatted as a table in
which every column is the model spectrum of a fixed age (except for the first
column, which gives the wavelength grid in Å). By default, model fluxes are
assumed to be given in units of $L_{\odot}$ Å${}^{-1}.$ For example, if the
user wishes to use BPASS v2.3 models with a given [$\alpha$/Fe] value, the
resulting FITS file will have 13 extensions (one per stellar metallicity in
BPASS). We illustrate the format of the model cube in Figure 1, where we
render the BPASS v2.3 models with binary physics and a solar [$\alpha$/Fe]
value as an Astropy Table.

The generation of a model cube, rather than using individual model spectra as
provided by their respective developers, serves two purposes. First, the cube
generator allows the user to truncate the model spectra to only the
approximate wavelength regime that is relevant for their analysis (instead of,
e.g., the full BPASS spectra which give fluxes from 1 to 100,000 Å in steps of
1 Å). Limiting the wavelength coverage of the model spectra greatly reduces
the file size of the model suite and speeds up the MCMC process. Second,
collating all model spectra for a given IMF, abundance ratios, or stellar
tracks into a single file is faster and simpler to load, manipulate, and
examine. This is especially true for users who wish to examine the effects of
varying the assumed IMF or abundance ratios when modeling stellar populations,
where the number of individual files would otherwise quickly become unwieldy.

3) A NumPy boolean array object to act as a pixel mask, which blocks out
problematic regions of the spectrum to be modeled (e.g., chip gaps, poorly-
subtracted sky lines, or emission and absorption features from Galactic and
extragalactic ISM along the line of sight). We recommend creating and saving
the pixel mask prior to running SESAMME and loading it as an input file
(assumed to be in ascii format).

###  2.3 Model Evaluation, Outputs, and Performance

The basic unit of output from SESAMME is an EnsembleSampler object, which is
the fundamental vehicle for information in emcee. Once the MCMC run is
complete, the sampler describes the state of all $N$ four-dimensional walkers
at each of the $S$ steps in the MCMC chain generated by emcee – in other
words, an array of dimension $4\times N\times S$. The sampler must then be
further processed (using standard methods and functions described in emcee)
before further analysis of the posterior PDFs. First, we recommend discarding
a number of “burn-in” steps in which the walkers explore the parameter space
to find regions of maximum probability density. The specific number to discard
(usually $\sim$a few hundred to 1000) should be determined by the user after
examining the completed chain. Further, because consecutive samples in an MCMC
chain are not fully independent, many steps may be needed before the ensemble
ceases to retain information about past states (that is, to become
independent). We recommend “thinning” the chain by a factor (approximately
$0.5-1\times$ the autocorrelation time, which is computed with the emcee
function get_autocorr_time) that must also be determined after completing the
MCMC run. Typical autocorrelation times for each variable in SESAMME are of
order $\sim$a few hundred steps. We find that using a thinning factor that is
roughly equal to the average of the four autocorrelation times is usually
sufficient to achieve a good sampling.

From the validation tests using either synthetic or empirical UV spectra that
we describe in Sections 3 and 4, we find that the sampler usually settles onto
its preferred solution(s) relatively quickly, with only small gains in
precision when using very long chains ( $>50,000$). Chains as short as 10,000
steps are often sufficient for the walker ensemble to stabilize around regions
of maximum density, even if formal convergence is not reached (as determined
by estimates of the autocorrelation time; see the next section).

SESAMME’s computational performance is adequate for use on personal computers.
For reference, most of the development and testing of SESAMME v1.0 was done on
a laptop with a 2.3 GHz Intel Core i9 processor and 32 GB of DDR4 memory.
Using this setup, emcee runs at $\sim$5 iterations per second with the nebular
continuum generation turned on (or $\sim$15 iterations per second with no
nebular continuum); a single run of 50,000 steps takes approximately 150
minutes (50 minutes with no nebular continuum) to complete without CPU
parallelization.

##  3 Validating the Algorithm

###  3.1 Testing on synthetic spectra

To illustrate the extent to which SESAMME can recover accurate stellar ages
and metallicities even from spectra with moderate signal-to-noise ratio S/N,
we created a suite of mock cluster IL spectra based on models from BPASS v2.3.
We used the BPASS v2.3 binary models with [$\alpha$/Fe] = 0, beginning with
the pure stellar population model of a given age, $Z$, and mass. Then we added
an appropriately scaled nebular continuum component and reddened the combined
stellar+nebular spectrum assuming a Cardelli et al. (1989) extinction law with
$R_{V}=3.1$. Because our initial science applications with SESAMME will
involve far-UV spectra of YMCs from HST/COS, the eight mock spectra that we
generate span 1 dex in age, metallicity, and reddening, but are generally
fairly young (at most 30 Myr). The true model parameters for the eight mock
spectra are given in columns 2 – 5 of Table 3. All spectra assume a broken
power law IMF with low-mass slope $\alpha_{1}$ = -1.30, high-mass slope
$\alpha_{2}$ = -2.35, and a maximum mass of 300 $M_{\odot}$ (“135_300” in
BPASS v2.2 and v2.3).

We then added a small noise component to the flux in each wavelength bin to
“jitter” the spectrum away from an ideal stellar population model. The noise
amplitude per wavelength bin (which is distinct from the flux uncertainty in
that bin) was drawn from a uniform distribution such that the flux is changed
by up to $\pm$15% from its noiseless value. We similarly created a flux
uncertainty array by drawing from a normal distribution, scaled to the
noiseless model flux in each bin, such that the typical flux uncertainty is 30
$\pm$ 10%. We then convolve this flat noise profile with the throughput of the
combined COS FUV G130M+G160M gratings, which creates a noise array with
similar shape to a typical observation with COS. The final noise array is
scaled such that the S/N ranges from $\sim$2 to 15 per Å across the 1150 –
1800 Å range we use to approximate the COS wavelength coverage. Finally, we
create a pixel mask that blocks out the approximate locations of COS chip
gaps, as well as narrow wavelength ranges ($\Delta$Å $\sim$5 – 20) at the
expected locations of strong ISM emission and absorption features (e.g., OI
$\lambda$1302, CII $\lambda$1334, SiIV $\lambda\lambda$1393,1402, CIV
$\lambda\lambda$1548,1550, and geocoronal Ly$\alpha$).

Columns 6 – 9 of Table 3 show the 16, 50, and 84th percentiles of the samples
in each of the four marginalized distributions. SESAMME is generally
successful in recovering accurate values for all 4 parameters across the
parameter space that we explore in this work, with $\sim$10 – 15%
uncertainties on the cluster age and mass, stellar metallicity, and average
reddening.

Table 3:  Input and recovered properties of synthetic validation spectraa | Input parameters | PDF medians from SESAMME
---|---|---
Name | log(Age) | $Z$ | $E(B-V)^{b}$ | Mass | log(Age) | $Z$ | $E(B-V)$ | Mass
| (yr) |  | (mag) | (10${}^{5}M_{\odot}$) | (yr) |  | (mag) | (10${}^{5}M_{\odot}$)
1 | 6.5 | 0.002 | 0.1 | 1.00 | 6.49${}^{+0.10}_{-0.10}$ | 0.002${}^{+0.001}_{-0.001}$ | 0.10${}^{+0.01}_{-0.01}$ | 0.94${}^{+0.34}_{-0.23}$
2 | 6.5 | 0.002 | 1.0 | 1.00 | 6.46${}^{+0.08}_{-0.08}$ | 0.002${}^{+0.001}_{-0.001}$ | 1.00${}^{+0.01}_{-0.01}$ | 0.85${}^{+0.19}_{-0.19}$
3 | 7.5 | 0.002 | 0.1 | 1.00 | 7.48${}^{+0.07}_{-0.09}$ | 0.002${}^{+0.001}_{-0.001}$ | 0.09${}^{+0.02}_{-0.02}$ | 0.97${}^{+0.15}_{-0.13}$
4 | 7.5 | 0.002 | 1.0 | 1.00 | 7.50${}^{+0.10}_{-0.11}$ | 0.002${}^{+0.001}_{-0.001}$ | 1.01${}^{+0.02}_{-0.02}$ | 1.08${}^{+0.19}_{-0.17}$
5 | 6.5 | 0.020 | 0.1 | 1.00 | 6.51${}^{+0.04}_{-0.04}$ | 0.020${}^{+0.004}_{-0.003}$ | 0.10${}^{+0.02}_{-0.01}$ | 1.00${}^{+0.15}_{-0.14}$
6 | 6.5 | 0.020 | 1.0 | 1.00 | 6.54${}^{+0.08}_{-0.06}$ | 0.020${}^{+0.003}_{-0.002}$ | 1.01${}^{+0.01}_{-0.01}$ | 1.16${}^{+0.31}_{-0.19}$
7 | 7.5 | 0.020 | 0.1 | 1.00 | 7.47${}^{+0.06}_{-0.08}$ | 0.020${}^{+0.004}_{-0.003}$ | 0.11${}^{+0.03}_{-0.02}$ | 1.01${}^{+0.11}_{-0.11}$
8 | 7.5 | 0.020 | 1.0 | 1.00 | 7.49${}^{+0.05}_{-0.08}$ | 0.022${}^{+0.004}_{-0.003}$ | 1.00${}^{+0.02}_{-0.01}$ | 0.92${}^{+0.11}_{-0.11}$

One of the key caveats of the L12 technique is that its accuracy in
determining stellar metallicities requires an age-accurate choice of
isochrone, especially for very young clusters whose properties evolve
significantly over the course of $\sim$a few Myr. Indeed, avoiding such age
assumptions was one of the driving forces behind the creation of SESAMME.
However, the code and user must make several assumptions which may affect the
posterior PDF for a given SESAMME run. In the following sections, we explore
the sensitivity of SESAMME’s output to changes in the assumed extinction curve
and in the assumed prevalence of nebular continuum emission.

###  3.2 Varying the Extinction Law

The eight model spectra we describe above were reddened assuming a Milky Way-
like extinction curve with $R_{V}$ = 3.1, and for simplicity our nominal runs
with SESAMME assume the same curve in the modeling process (i.e., with the
parameter use_ext_law set to the Cardelli et al. (1989) curve). To examine the
impact of an improper choice of extinction curve on the inferred stellar
properties, we also ran the eight synthetic spectra through SESAMME with
use_ext_law set to the Calzetti et al. (2000) starburst extinction curve with
$R_{V}$ = 4.05. All other parameters within SESAMME (e.g., priors, chain
lengths, etc.) were identical to the baseline runs.

When using the Calzetti et al. (2000) law to model spectra that, in reality,
are subject to Milky Way-like extinction, SESAMME infers reddening values that
are higher than the true values by as much as $\sim 30$%, and the inferred
cluster mass is also overestimated by 0.2-0.3 dex on average. The median ages
increase by 0.1 dex at most, and the median metallicities are essentially
unaffected. Further, upon visually comparing the models that correspond to the
median PDF values in both the MW and the Calzetti case, it becomes clear that
assuming a Calzetti-like extinction curve provides an inferior “best fit” to
our synthetic spectra. For example, for spectrum 5 (young, high-metallicity,
and $E(B-V)$ = 1), the optimal model with MW-like reddening yields a reduced
$\chi^{2}=0.97$, while the optimal model with the Calzetti curve yields
$\chi^{2}=4.9$. These changes are driven by a combination of i) the shallower
slope at very blue wavelengths in the Calzetti et al. (2000) curve compared to
the Cardelli et al. (1989) curve, and ii) the lower S/N at the bluest and
reddest wavelengths of our synthetic spectra (such that the model fitting is
heavily weighted towards the highest S/N regions of the spectrum around $\sim
1500$ Å). These factors lead to a situation where the optimal model is more
heavily reddened (to better approximate the overall continuum shape), but
correspondingly more massive to compensate for the dearth of flux at all
wavelengths that results from the increased reddening. This test highlights
the importance of accurately selecting an extinction curve that best describes
the specific environmental conditions of the targets.

###  3.3 Effect of the Nebular Continuum

We remind users of SESAMME v1.0 that SESAMME’s default behavior is to add
nebular continuum emission to the stellar component during the modeling
process. However, for flexibility, we also allow users to turn nebular
continuum emission off.

To illustrate how the (non)inclusion of nebular continuum emission affects the
inferred properties of a star cluster, we re-ran the eight synthetic spectra
through SESAMME with the nebular continuum component set to zero at all
wavelengths during the modeling process (hereafter the “no-neb” case), with
all other parameters identical to the baseline run. In Figure 2, we compare
the median ages, reddening values, and masses that are inferred from the
baseline and no-neb runs; the median metallicity does not change significantly
between the baseline and no-neb cases and thus is not shown in the figure.
From these tests, we can conclude that the non-inclusion of nebular continuum
emission in the models has the following effects on each of the four
variables:

i) Age – Cluster ages in the no-neb case are slightly but systematically older
than in the baseline case, with the most extreme differences for very young
and metal-poor clusters (clusters 5 and 6).

ii) Metallicities – Differences in the inferred metallicities are negligible
within the uncertainties. However, the metallicity uncertainties seem to
broaden slightly at low $Z$ values when nebular continuum is not taken into
account.

iii) Reddening – SESAMME prefers models that are systematically more reddened
in the no-neb case than in the baseline case, with differences as large as 50%
for the youngest and most metal-poor clusters.

iv) Mass – For clusters younger than $\sim$10 Myr, SESAMME infers stellar
masses that are systematically higher (by up to a factor of 2) in the no-neb
case.

These behaviors are mostly expected due to the nature of nebular continuum
emission and how SESAMME treats its inclusion in the modeling. SESAMME’s
default modeling procedure includes a physically motivated amount of nebular
continuum emission – that is, for a given age, mass, and metallicity, the
stellar ionizing photon output $Q$ and hence nebular continuum strength is
fixed. In the FUV, the nebular continuum rises rapidly from 912 Å, peaks
around 1500 Å, then gently decreases to redder wavelengths. For young stellar
populations, this means that the nebular contintuum contributes a greater
fraction of UV light (relative to the stellar continuum) at longer
wavelengths. The result is a spectrum whose UV slope is shallower than one
with purely stellar emission, which can mimic the combined effects of
extinction and aging. In situations where nebular continuum is actually
present in the system under observation but is omitted from the modeling
procedure, the sampler will converge to a model that is both more massive (to
compensate for the dearth of flux at the red end) and more reddened (to
account for the flatter-than-expected spectral slope). In short, we can
conclude that the inclusion of a nebular continuum emission component in the
models has a non-negligible impact on the inferred properties of young,
massive star clusters and should not be ignored.

Figure 2:  Fractional residual difference of the ages (top), reddening
(middle), and masses (bottom) inferred for the eight synthetic FUV spectra
with (“nominal”) and without (“no-nebular”) a nebular continuum component in
the modeling. All points are colored according to the input age for that mock
cluster spectrum (blue for 3 Myr, gold for 30 Myr). Squares denote low-
metallicity clusters (spectra 1 – 4) and circles denote solar-metallicity
clusters (spectra 5 – 8).  Figure 3:  Corner plot for the SESAMME fit to the
UV spectrum of M83-8, using either Starburst99 (purple) or BPASS v2.3 (red) as
the assumed model suite. This spectrum has a S/N at 1310 Å of $\sim$6.8 at the
native 0.06 Å/resel resolution of COS (S/N $\sim$26.8 after binning to at the
1 Å resolution of BPASS). Red and purple dashed lines mark the median of each
parameter’s PDF for the two model sets. Grey dotted lines mark the BPASS age
and metallicity grid points; grid points for our SB99 models are omitted for
clarity. Black solid lines mark the best-fit age and metallicity from
Hernandez et al. (2019). The 2D and 1D histograms have been lightly smoothed
for clarity.

##  4 Re-examining YMCs in M83 and NGC 1313

Figure 4:  An example of SESAMME’s performance in modeling the FUV spectrum of
a metal-rich YMC. Top: The black curve shows the Hi absorption-corrected,
combined G130M+G160M spectrum of YMC M83-8 (Hernandez et al., 2019) from
HST/COS, after smoothing and resampling to the 1 Å resolution of the BPASS
models. Intervals in grey shaded regions were masked during the modeling
process and include ISM emission and absorption features, chip gaps, and
geocoronal Ly$\alpha$ and Oi $\lambda 1302$. The blue curve shows the median
model (see text for details), while the teal dotted curves show 50 random
samples from the cleaned MCMC chain. Middle: Model residuals for the 50
samples. Bottom: Zoom-ins to regions that include age-sensitive features: Nv
$\lambda 1240$ in the bottom left and Civ $\lambda 1550$ in the bottom right.

Figure 5:  An example of SESAMME’s performance using SB99 models to
characterize the FUV spectrum of a YMC. Top: The black curve shows the Hi
absorption-corrected, combined G130M+G160M spectrum of YMC M83-8 (Hernandez et
al., 2019) from HST/COS, now smoothed and resampled to the 0.6 Å resolution of
our SB99 models. As in Figure 4, intervals in grey shaded regions were masked
during the modeling process, and teal dotted curves show either 50 random
samples from the MCMC chain (top) or the residuals of those samples from the
observed spectrum (middle). In the bottom panels we show zoom-ins to age-
sensitive features (Nv $\lambda 1240$ on the left and Civ $\lambda 1550$ on
the right).

As a first science use of SESAMME, we also applied it to a preliminary sample
of four FUV spectra of YMCs from HST/COS. The sample includes two sightlines
towards YMCs in the massive, metal-rich spiral M83 (e.g., Bresolin et al.,
2016) and two towards YMCs in the LMC-metallicity galaxy NGC 1313 (e.g., Walsh
& Roy, 1997). For all clusters, we use the combined G130M+G160M data from COS.
The spectra towards M83 were originally obtained as part of HST program 14681
(PI: Aloisi), while the spectra towards NGC 1313 were obtained as part of HST
program 15627 (PI: Adamo).

For each YMC, we created both a high- and low-resolution version of each COS
spectrum. The high-resolution spectrum was used only for determining the
recession velocity of the cluster, while the low-resolution one was used as
input to SESAMME. We downloaded individual spectra for each object from the
Mikulski Archive for Space Telescopes (MAST) and reduced them using the
calibration pipeline  CALCOS v3.4.0. We further processed the individual
exposures using an IDL code by the COS Guaranteed Time Observer Team (Danforth
et al., 2010), which interpolates each spectrum onto a common wavelength grid
before weight-combining them. Lastly, we rebinned the final spectra by one COS
resolution element (1 resel = 6 pixels) that corresponds to the nominal point-
spread function. This process also rescales and stitches together the spectra
from the different COS gratings into a single, coadded file. These high-
resolution spectra have a typical S/N of $\sim$5 – 7 per resel (0.06 Å for
G130M and 0.07 Å for G160M) at 1310 Å.

We correct for Galactic reddening along the line of sight using the extinction
maps by Schlafly & Finkbeiner (2011), assuming a Cardelli et al. (1989)
extinction curve with $R_{V}=3.1$. We also fit and divide out a Ly$\alpha$
absorption profile, which includes a foreground component from Milky Way H i,
using the Python software VoigtFit (Krogager, 2018). We convolve the COS line
spread function profiles with the FWHM of the source in the dispersion
direction (as measured from the acquisition images), which accounts for
instrument-induced broadening of the Ly$\alpha$ absorption profile. This
procedure is similar to that used in James et al. (2014) and Hernandez et al.
(2020, 2021).

We then use the high-resolution spectrum to measure the recession velocity of
the cluster from photospheric lines. After fitting the local stellar continuum
around the relevant lines, we fit a series of Gaussians to 4 – 6 photospheric
lines across the wavelength range of the spectrum; these typically included C
iii $\lambda 1247$, Si iii $\lambda 1294$, C ii $\lambda 1324$, O iv $\lambda
1342$, and S v $\lambda 1502$. We take the average redshift of these lines as
the recession velocity of the cluster and use it to shift the COS spectrum to
the rest frame.

Proper comparisons with stellar population models must account for the
differing spectral resolutions of the COS spectra and the models (e.g., 1 Å
for BPASS). To do this, we smooth the spectrum with a Gaussian filter, such
that $\sigma_{s}=\sqrt{(}\sigma_{model}^{2}-\sigma_{COS}^{2})$. We then
resample the smoothed spectrum onto the model wavelength grid (e.g., in steps
of 1 Å to match the BPASS wavelength grid). Finally, we use the dereddened, H
i-corrected, and smoothed spectra as input to SESAMME.

We placed loose constraints on the age priors according to literature age
estimates for each cluster. For the YMCs in M83, the ages and age
uncertainties are derived from SED fitting (Hernandez et al., 2019) and have a
typical uncertainty of 50%. For example, Hernandez et al. (2019) estimated the
cluster M83-3 to be $5\pm 3$ Myr old from HST/WFC3 photometry; therefore, when
modeling this cluster with SESAMME, we constrain the age to be at most 10 Myr
old. For the YMCs in NGC 1313, the ages and age uncertainties are taken from
the CLusters in the UV as EngineS survey (CLUES; Sirressi et al., 2022), which
uses a two-SSP spectral fit to constrain the properties of the clusters. We
use their light-weighted ages as the comparison point for our work. Priors on
the stellar metallicity were determined in a similar manner and generally
spanned a range of 0.5 – 2$\times$ the literature value. The reddening and
amplitude parameters are left essentially unconstrained.

Table 4:  Inferred properties of YMCs from HST/COS FUV spectraa | BPASS | Starburst99
---|---|---
Name | log(Age) | $Z$ | $E(B-V)^{a}$ | Mass | log(Age) | $Z$ | $E(B-V)^{a}$ | Mass
| (yr) |  | (mag) | (10${}^{4}M_{\odot}$) | (yr) |  | (mag) | (10${}^{4}M_{\odot}$)
M83-3 | 6.53${}^{+0.43}_{-0.06}$ | 0.023${}^{+0.007}_{-0.002}$ | 0.27${}^{+0.02}_{-0.00}$ | 71${}^{+417}_{-3}$ | 6.48${}^{+0.09}_{-0.06}$ | 0.023${}^{+0.005}_{-0.003}$ | 0.25${}^{+0.02}_{-0.00}$ | 51${}^{+47}_{-1}$
M83-8 | 6.46${}^{+0.59}_{-0.08}$ | 0.010${}^{+0.005}_{-0.002}$ | 0.14${}^{+0.03}_{-0.01}$ | 2${}^{+14}_{-1}$ | 6.68${}^{+0.04}_{-0.14}$ | 0.011${}^{+0.003}_{-0.001}$ | 0.15${}^{+0.06}_{-0.003}$ | 3${}^{+0}_{-1}$
NGC 1313-1 | 6.45${}^{+0.606}_{-0.07}$ | 0.005${}^{+0.002}_{-0.001}$ | 0.18${}^{+0.02}_{-0.01}$ | 5${}^{+56}_{-1}$ | 6.50${}^{+0.28}_{-0.08}$ | 0.003${}^{+0.002}_{-0.001}$ | 0.12${}^{+0.01}_{-0.00}$ | 3${}^{+3}_{-0}$
NGC 1313-2 | 6.80${}^{+0.38}_{-0.11}$ | 0.003${}^{+0.001}_{-0.000}$ | 0.28${}^{+0.05}_{-0.05}$ | 8${}^{+24}_{-3}$ | 6.73${}^{+0.28}_{-0.05}$ | 0.004${}^{+0.000}_{-0.000}$ | 0.28${}^{+0.01}_{-0.01}$ | 5${}^{+9}_{-0}$

We then allow SESAMME to explore the BPASS parameter space using an ensemble
of 128 walkers in chains of length 5000 steps, with the nebular continuum
component turned on. For each cluster, we use the BPASS v2.3 binary models
with solar [$\alpha$/Fe] and a broken power law IMF with low-mass slope
$\alpha_{1}$ = -1.30, high-mass slope $\alpha_{2}$ = -2.35, and a maximum mass
of 300 $M_{\odot}$. We further assume a Calzetti et al. (2000) extinction
curve, since test runs using other extinction curves generally produce a
poorer fit to the data, particularly in the heavily star-forming environments
we are interested in. To generate the posterior PDFs for our four parameters,
we use emcee’s helper function get_chain to discard the first 300 – 800 steps
(the “burn-in”) of the chain, then thin the chain by a factor of $\sim$100 –
300. The exact value of the thinning factor is determined at the end of each
run using the emcee function get_autocorr_time. This leaves a collection of
several thousand independent samples from which we can construct posterior
PDFs for each of the four variables in the model. For each parameter, we take
the PDF median as the “best-fit” value and the 16th and 84th percentile values
as the 1$\sigma$ uncertainties. We show the 2D and 1D posterior PDFs as a
corner plot in Figure 3.

We show an example of SESAMME’s performance on the YMC M83-8 in Figure 4. In
addition to the smoothed and rebinned COS FUV spectrum, we also show the
“median model” inferred by SESAMME. This is the total model (stellar plus
nebular continuum) for which the four parameters are set to the 50th
percentiles of their respective, one-dimensional posterior PDFs (see columns 2
– 5 of the second row of Table 4). Further, we show 50 random draws from the
ensemble sampler as teal curves, which illustrate the range of models that
SESAMME considers to be a plausible fit to the data within the uncertainties.
We show the median models for the other three clusters in the Appendix (Fig.
A.1).

To test the self-consistency of SESAMME across different model sets, we then
performed this entire procedure a second time using a different set of stellar
population models. Specifically, we used Starburst99 single-burst models of a
fixed mass of 10${}^{6}M_{\odot}$, assuming the high mass-loss Geneva stellar
tracks and a double-power law IMF with low-mass slope $\alpha_{1}$ = -1.30,
high-mass slope $\alpha_{2}$ = -2.3, and a maximum mass of 100 $M_{\odot}$. We
generated models with $Z$ = 0.001, 0.004, 0.008, 0.02, and 0.04, each of which
includes model spectra in steps of 1 Myr for ages between 1 and 100 Myr.
Because we are interested in modeling FUV spectra from COS, we used only the
high-resolution, theoretical UV spectra (ifaspec files; $\sim$0.6 Å/pix) from
Starburst99. After smoothing and rebinning the COS spectra to the 0.6 Å/pixel
wavelength grid of the models, we re-ran the analysis described above for each
YMC. The ifaspec spectra already include a nebular continuum component in
addition to the stellar-only light, so we turned off the nebular continuum
feature within SESAMME. All other aspects of each MCMC run (priors, chain
length, etc.) were identical to the runs with BPASS models. We show the median
SB99 model and 50 random draws from the ensemble sampler in Figure 5.

We compare the results from SESAMME for each of the two model sets in Figure
6. For each cluster, the median age and metallicity inferred by SESAMME is
generally consistent (within the uncertainties) across the two model sets. We
also compare the ages and metallicities from SESAMME against literature
estimates for the same clusters. For the YMCs in M83, the ages are derived
from HST/WFC3 photometry in Hernandez et al. (2019), while the metallicities
are determined via isochrone fitting (Hernandez et al., 2019). For the YMCs in
NGC 1313, the literature ages and metallicities are the light-weighted
averages of the two-population SSP modeling in Sirressi et al. (2022).

For both BPASS and SB99 models, SESAMME’s agreement with the literature
estimates of the cluster ages is generally good within the uncertainties. The
major exception to this is the cluster NGC 1313-2, where SESAMME prefers
solutions with a cluster age around $\sim$6 – 8 Myr and reddening
$E(B-V)\sim$0.28, regardless of whether we use BPASS or Starburst99 SSP
models. This stands in contrast to the light-weighted age of $\sim$30 Myr and
slightly lower reddening of $\sim$0.2 mag determined by Sirressi et al.
(2022), who used Starburst99 model SSPs for their two-population fits.
Comparing their best fit for NGC 1313-2 (the purple curve in their Figure 11)
and the median model inferred from our SESAMME run with BPASS (the bottom
panel of Fig. A.1), we see that each model provides a reasonably good by-eye
match to the observed stellar continuum. The older model preferred by Sirressi
et al. (2022) is also consistent with the lack of strong Si iv and C iv wind
features, but does not explain the moderate N v emission seen in the COS
spectrum. Our preferred model(s), on the other hand, are young enough to
reproduce the N v profile, but also predict P-Cygni-like profiles for the Si
iv and C iv lines, which is not observed. This may be an indication that this
cluster indeed consists of two populations – a dominant component that is
$>20$ Myr old, and a small, very young component ($\sim 1-3$ Myr) that drives
the weak N v emission.

We also note that the age uncertainties when using Starburst99 are generally
smaller than when using BPASS, which is due in part to the finer and linearly-
uniform age grid of our Starburst99 model suite ($\Delta$age = 1 Myr) compared
to that of BPASS ($\Delta$log(age) = 0.1). Finally, the inclusion of binary
star physics in BPASS can cause a population to remain UV-luminous at slightly
older ages compared to single-star models. This may explain why SESAMME
considers older populations to be plausible solutions when using BPASS but not
Starburst99 (i.e., a tail to older ages or a secondary maximum in the age
posterior PDF, which leads to the asymmetric errorbars on the red points in
Figure 6).

Figure 6:  Comparison of the 16th, 50th, and 84th percentiles on the ages
(left) and stellar metallicities (right) of the YMCs in M83 and NGC 1313 as
determined by SESAMME (circles) and from other methods in the literature (teal
triangles). For each cluster, we performed two SESAMME runs – one with BPASS
v2.3 models (red) and one with Starburst99 (purple). The literature values for
the M83 YMCs are from Hernandez et al. (2019), while the values for the NGC
1313 YMCs are from Sirressi et al. (2022).

##  5 Discussion and Summary

Figure 7:  Comparison of N v$\lambda\lambda 1238,1242$ line profiles from
BPASS v2.3 (top) and Starburst99 (bottom) models. Red and purple curves
respectively show populations of age 2 and 5 Myr. Populations with a stellar
metallicity $Z=0.008$ are shown as solid lines, while $Z=0.04$ models are
shown with dashed lines. Some models have been shifted in luminosity for
visualization purposes. The vertical grey lines mark the central wavelengths
of the N v doublet components. Starburst99 predicts a significant P-Cygni line
profile even at high metallicity and ages $\sim$5 Myr, while N v is seemingly
absent in all but the youngest and most metal-poor BPASS models we show here.

One of the driving principles behind the design of SESAMME is that it should
be highly flexible in terms of what kinds of models it can accept as input,
which in turn enables it to be useful in a greater variety of science cases.
In particular, our focus on the rest-frame FUV in this work (and in a
forthcoming work in which we use SESAMME to study the homogeneity of stellar
and ionized-gas metallicities in nearby SFGs) is especially relevant for
ongoing studies of the stellar populations of galaxies at $z>3$, where the UV
output of YMCs is shifted into optical and IR bandpasses (e.g., the VANDELS
survey of Pentericci et al. (2018); Cullen et al. (2019) or the explosion of
recent results from JWST spectroscopy). For example, unlike past versions,
stellar population models in BPASS v2.3 are allowed to vary in their level of
$\alpha$ enhancement at a fixed overall metallicity mass fraction $Z$ (i.e.,
for a fixed $Z$, increases in [$\alpha$/Fe] are offset by decreases in [Fe/H]
and vice versa). This allows SESAMME to be used to constrain the properties of
stellar populations in a greater variety of environments, including clusters
in starbursting galaxies and other systems that are likely analogs of higher-
redshift objects.

On the other hand, stellar population models in BPASS v2.3 are only computed
for a single assumed IMF, compared to the nine supplied as part of the v2.2
release. Using BPASS v2.3 may thus present issues when modeling the IL of old
or intermediate-age clusters, where assumptions about the low mass end of the
stellar IMF may become important in the interpretation of spectra. However,
because the IMFs used in BPASS v2.2 do not vary much at the high mass end, the
choice of IMF should not greatly affect the modeling of rest-frame UV spectra,
which are dominated by the most massive stars.

Because SESAMME performs a full spectral fit to constrain the age and
metallicity of a cluster, it has the advantage of not being dependent on the
predicted strength of any single stellar absorption or wind feature, which may
be model-dependent. We illustrate the importance of this in Figure 7, where we
show a series of models centered around the N v$\lambda\lambda 1238,1242$ wind
feature, which is known to be present only in fairly young ($\lesssim$5 Myr)
populations. The lower panel of Figure 7 shows the N v profiles that are
predicted by Starburst99 for a single burst of mass 10${}^{6}M_{\odot}$,
assuming the high mass-loss Geneva stellar tracks and a double-power law IMF
with the default values for the slopes and mass cutoffs. It is clear that, for
a 2 Myr old population, Starburst99 predicts a stark P-Cygni profile at both
low ($Z=0.5Z_{\odot}$) and high ($Z=2Z_{\odot}$) metallicity. Once a
population reaches 5 Myr old, the absorption component of N v has weakened
significantly, but there are still hints of a P-Cygni profile across a range
of metallicities. In the upper panel of the same figure, we show the BPASS
v2.3 models with binary physics and solar abundance ratios for the same range
of ages and metallicities. Unlike in Starburst99, P-Cygni-like emission for N
v is seemingly only predicted for populations that are very young and fairly
metal-poor. The lack of obvious N v emission from a high metallicity model
population, even at ages $\sim$2 Myr, is somewhat unusual and could lead to an
inaccurate characterization of a YMC if N v is used as the primary age
diagnostic in solar or super-solar metallicity environments. (See also the
median model for cluster M83-3 in Fig. A.1, which fails to reproduce the N v
profile despite adequately tracing the continuum and other features like the
Si iv wind line.)

We also note that all of our tests (with both mock and real FUV IL spectra)
focus on relatively massive star clusters with $M_{*}\sim
10^{4}-10^{5}M_{\odot}$, such that the IMF should be well-sampled across the
full mass range. However, Orozco-Duarte et al. (2022) recently showed that
stochastic sampling of the IMF in low-mass star clusters can drastically
affect the cluster’s UV light output. For example, they report that if the
high-mass end of the IMF is not well-sampled, then a dearth of OB or Wolf-
Rayet stars can lead to highly uncertain cluster age determinations from UBVI
photometry, even in clusters as massive as $10^{5}M_{\odot}$. Although their
conclusions were based on photometry of massive star clusters, we expect that
a similar effect may be seen while modeling the FUV spectra of low-mass YMCs.
However, since the most frequently-used stellar population codes assume that
the IMF is sampled deterministically, this variance in the number of very
massive and UV-bright stars is difficult to account for.

Finally, we emphasize again that SESAMME is designed to model the IL spectra
of individual clusters using single SSP fitting (i.e., every star cluster
assumed to be the product of a single burst of star formation). As such, it is
not meant to be a spectral fitting tool for entire galaxies, for which one
must account for a potentially complex star-formation history. Although
SESAMME can be used to characterize YMCs in relatively distant galaxies, we
caution that any such observations must have sufficient spatial resolution to
isolate the YMC from its environs (to mitigate significant light contribution
from, for example, an underlying old population).

In summary, we have presented the v1.0 public release of SESAMME, a flexible
Python tool for Simultaneous Estimates of Star-cluster Age, Metallicity, Mass,
and Extinction. In this work, we described the structure and functionality of
SESAMME as a full spectral fitting tool for IL spectroscopy of star clusters,
including validation tests on both synthetic (Sec. 3) and empirical (Sec. 4)
FUV spectra of young, massive clusters. The test runs with synthetic YMC
spectra, built from noise-added BPASS v2.3 stellar population models, show
that SESAMME can successfully recover the known age and metallicity of the
model, even at moderate signal-to-noise ratios $\sim$5 Å-1. Our tests with
empirical FUV spectra from HST/COS, which included two sightlines towards YMCs
in M83 and two in NGC 1313, were similarly successful. The ages and
metallicities inferred by SESAMME for these clusters are both i) consistent
with literature estimates that use independent methods, and ii) internally
consistent across different assumed stellar population model sets. With the
utility of SESAMME thus established in the FUV, we plan to apply it to a
larger sample of extragalactic star clusters with IL spectroscopy from HST/COS
(Jones et al., in prep.). SESAMME is modular and open-source, and we encourage
contributions at all levels from the community.

We thank Claus Leitherer for helpful discussions regarding Starburst99 and
Russell Ryan for helpful discussions on the use of emcee. We also thank the
anonymous referee for their careful review and suggestions, which improved the
contents of this manuscript. SESAMME makes use of Astropy (Astropy
Collaboration et al., 2013, 2018, 2022), emcee (Foreman-Mackey et al., 2013),
dust_extinction (Gordon et al., 2022), extinction (Barbary, 2021), and NumPy
(van der Walt et al., 2011) in its computations, as well as Matplotlib
(Hunter, 2007) and Jupyter (Kluyver et al., 2016) for data visualizations.
This work additionally makes use of SciPy (Virtanen et al., 2020) in our
analyses of empirical HST/COS spectra and of v2.3 of the Binary Population and
Spectral Synthesis (BPASS) models as last described in Byrne et al. (2022) and
Stanway & Eldridge (2018). S. H. and B. L. J. acknowledge support from the
European Space Agency (ESA). Support for this work was provided by NASA
through grant AR-16130 from the Space Telescope Science Institute, which is
operated by AURA, Inc., under NASA contract NAS5-26555. 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
https://doi.org/10.17909/mnp4-f880 (catalog DOI: 10.17909/mnp4-f880).

## References

  * Aller (1984) Aller, L. H. 1984, Physics of thermal gaseous nebulae, doi: 10.1007/978-94-010-9639-3
  * Asa’d et al. (2022) Asa’d, R., Hernandez, S., As’ad, A., et al. 2022, ApJ, 929, 174, doi: 10.3847/1538-4357/ac5f3e
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of Astronomy and Astrophysics, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barbary (2021) Barbary, K. 2021, extinction: Dust extinction laws, Astrophysics Source Code Library, record ascl:2102.026.  http://ascl.net/2102.026
  * Bresolin et al. (2016) Bresolin, F., Kudritzki, R.-P., Urbaneja, M. A., et al. 2016, ApJ, 830, 64, doi: 10.3847/0004-637X/830/2/64
  * Byrne et al. (2022) Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329, doi: 10.1093/mnras/stac807
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chisholm et al. (2019) Chisholm, J., Rigby, J. R., Bayliss, M., et al. 2019, ApJ, 882, 182, doi: 10.3847/1538-4357/ab3104
  * Colucci et al. (2011) Colucci, J. E., Bernstein, R. A., Cameron, S. A., & McWilliam, A. 2011, ApJ, 735, 55, doi: 10.1088/0004-637X/735/1/55
  * Colucci et al. (2012) —. 2012, ApJ, 746, 29, doi: 10.1088/0004-637X/746/1/29
  * Cullen et al. (2019) Cullen, F., McLure, R. J., Dunlop, J. S., et al. 2019, MNRAS, 487, 2038, doi: 10.1093/mnras/stz1402
  * Danforth et al. (2010) Danforth, C. W., Keeney, B. A., Stocke, J. T., Shull, J. M., & Yao, Y. 2010, ApJ, 720, 976, doi: 10.1088/0004-637X/720/1/976
  * Davies et al. (2010) Davies, B., Kudritzki, R.-P., & Figer, D. F. 2010, MNRAS, 407, 1203, doi: 10.1111/j.1365-2966.2010.16965.x
  * Davies et al. (2015) Davies, B., Kudritzki, R.-P., Gazak, Z., et al. 2015, ApJ, 806, 21, doi: 10.1088/0004-637X/806/1/21
  * Davies et al. (2017) Davies, B., Kudritzki, R.-P., Lardo, C., et al. 2017, ApJ, 847, 112, doi: 10.3847/1538-4357/aa89ed
  * Ferland (1980) Ferland, G. J. 1980, PASP, 92, 596, doi: 10.1086/130718
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293
  * Fitzpatrick & Massa (2007) Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320, doi: 10.1086/518158
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gazak et al. (2014) Gazak, J. Z., Davies, B., Bastian, N., et al. 2014, ApJ, 787, 142, doi: 10.1088/0004-637X/787/2/142
  * Gordon et al. (2022) Gordon, K., Larson, K., McBride, A., et al. 2022, karllark/dust_extinction: NIRSpectralExtinctionAdded, v1.1, Zenodo, doi: 10.5281/zenodo.6397654
  * Gordon et al. (2023) Gordon, K. D., Clayton, G. C., Decleir, M., et al. 2023, ApJ, 950, 86, doi: 10.3847/1538-4357/accb59
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Grevesse & Noels (1993) Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, ed. N. Prantzos, E. Vangioni-Flam, & M. Casse, 15–25
  * Gvozdenko et al. (2022) Gvozdenko, A., Larsen, S. S., Beasley, M. A., & Brodie, J. 2022, A&A, 666, A159, doi: 10.1051/0004-6361/202243415
  * Halliday et al. (2008) Halliday, C., Daddi, E., Cimatti, A., et al. 2008, A&A, 479, 417, doi: 10.1051/0004-6361:20078673
  * Hernandez et al. (2020) Hernandez, S., Aloisi, A., James, B. L., et al. 2020, ApJ, 892, 19, doi: 10.3847/1538-4357/ab77c6
  * Hernandez et al. (2017) Hernandez, S., Larsen, S., Trager, S., Groot, P., & Kaper, L. 2017, A&A, 603, A119, doi: 10.1051/0004-6361/201730550
  * Hernandez et al. (2022) Hernandez, S., Winch, A., Larsen, S., James, B. L., & Jones, L. 2022, AJ, 164, 89, doi: 10.3847/1538-3881/ac7ebe
  * Hernandez et al. (2019) Hernandez, S., Larsen, S., Aloisi, A., et al. 2019, ApJ, 872, 116, doi: 10.3847/1538-4357/ab017a
  * Hernandez et al. (2021) Hernandez, S., Aloisi, A., James, B. L., et al. 2021, ApJ, 908, 226, doi: 10.3847/1538-4357/abd6c4
  * Ho et al. (2015) Ho, I. T., Kudritzki, R.-P., Kewley, L. J., et al. 2015, MNRAS, 448, 2030, doi: 10.1093/mnras/stv067
  * Hosek et al. (2014) Hosek, Matthew W., J., Kudritzki, R.-P., Bresolin, F., et al. 2014, ApJ, 785, 151, doi: 10.1088/0004-637X/785/2/151
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * James et al. (2014) James, B. L., Aloisi, A., Heckman, T., Sohn, S. T., & Wolfe, M. A. 2014, ApJ, 795, 109, doi: 10.1088/0004-637X/795/2/109
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in IOS Press, 87–90, doi: 10.3233/978-1-61499-649-1-87
  * Krogager (2018) Krogager, J.-K. 2018, VoigtFit: Absorption line fitting for Voigt profiles, Astrophysics Source Code Library, record ascl:1811.016.  http://ascl.net/1811.016
  * Kudritzki et al. (2015) Kudritzki, R.-P., Ho, I. T., Schruba, A., et al. 2015, MNRAS, 450, 342, doi: 10.1093/mnras/stv522
  * Kudritzki et al. (2014) Kudritzki, R.-P., Urbaneja, M. A., Bresolin, F., Hosek, Matthew W., J., & Przybilla, N. 2014, ApJ, 788, 56, doi: 10.1088/0004-637X/788/1/56
  * Kudritzki et al. (2012) Kudritzki, R.-P., Urbaneja, M. A., Gazak, Z., et al. 2012, ApJ, 747, 15, doi: 10.1088/0004-637X/747/1/15
  * Kunth & Sargent (1986) Kunth, D., & Sargent, W. L. W. 1986, ApJ, 300, 496, doi: 10.1086/163828
  * Larsen et al. (2012) Larsen, S. S., Brodie, J. P., & Strader, J. 2012, A&A, 546, A53, doi: 10.1051/0004-6361/201219895
  * Larsen et al. (2018) Larsen, S. S., Pugliese, G., & Brodie, J. P. 2018, A&A, 617, A119, doi: 10.1051/0004-6361/201832767
  * Lebouteiller et al. (2013) Lebouteiller, V., Heap, S., Hubeny, I., & Kunth, D. 2013, A&A, 553, A16, doi: 10.1051/0004-6361/201220948
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Maiolino & Mannucci (2019) Maiolino, R., & Mannucci, F. 2019, A&A Rev., 27, 3, doi: 10.1007/s00159-018-0112-2
  * McWilliam & Bernstein (2008) McWilliam, A., & Bernstein, R. A. 2008, ApJ, 684, 326, doi: 10.1086/589957
  * Orozco-Duarte et al. (2022) Orozco-Duarte, R., Wofford, A., Vidal-García, A., et al. 2022, MNRAS, 509, 522, doi: 10.1093/mnras/stab2988
  * Pentericci et al. (2018) Pentericci, L., McLure, R. J., Garilli, B., et al. 2018, A&A, 616, A174, doi: 10.1051/0004-6361/201833047
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Sirressi et al. (2022) Sirressi, M., Adamo, A., Hayes, M., et al. 2022, MNRAS, 510, 4819, doi: 10.1093/mnras/stab3774
  * Stanway & Eldridge (2018) Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75, doi: 10.1093/mnras/sty1353
  * Steidel et al. (2016) Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159, doi: 10.3847/0004-637X/826/2/159
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walsh & Roy (1997) Walsh, J. R., & Roy, J. R. 1997, MNRAS, 288, 726, doi: 10.1093/mnras/288.3.726
  * Zaritsky et al. (1994) Zaritsky, D., Kennicutt, Robert C., J., & Huchra, J. P. 1994, ApJ, 420, 87, doi: 10.1086/173544

##  Appendix A SESAMME Fits to Remaining Clusters

Figure A.1: Continuation of Fig. 4 showing the performance of SESAMME in
fitting FUV spectra of YMCs. Top: The black curve shows the H i absorption-
corrected, combined G130M+G160M spectrum of YMC M83-3 (Hernandez et al., 2019)
from HST/COS, after smoothing and resampling to the 0.6 Å resolution of our
Starburst99 models. Intervals in grey shaded regions were masked during the
modeling process and include ISM emission and absorption features, chip gaps,
and geocoronal Ly$\alpha$ and O i. The blue curve shows the median model (see
text for details), while the teal dotted curves show 50 random samples from
the cleaned MCMC chain. Bottom: Model residuals for the 50 samples.

Figure A.2: Continuation of Fig. A.1 for the subsolar metallicity YMC NGC
1313-1, as modeled using BPASS (top) and Starburst99 (bottom) SSPs.

Figure A.3: Continuation of Fig. A.1 for subsolar metallicity YMC NGC 1313-2,
as modeled using BPASS (top) and Starburst99 (bottom) SSPs.