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

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 ~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.


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 starforming 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 ∼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 ob-tained 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 ∼a few Mpc), measurements of stellar metallicities have been achieved using spectroscopic observations of blue and red supergiants (e.g., Davies et al. 2010Davies et al. , 2015Davies et al. , 2017;;Kudritzki et al. 2012Kudritzki et al. , 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. 2011Colucci et al. , 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(Larsen et al. , 2018;;Hernandez et al. 2017Hernandez et al. , 2019Hernandez et al. , 2021Hernandez et al. , 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 ∼30,000) optical observations, but has since been successfully applied to spectra at lower resolutions R ∼ 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σ uncertainties induced changes in their inferred metallicities of order ∼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σ 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 use1 .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.

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 precomputed 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 Ta-  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 bles 3 and 4), making it unnecessary to interpolate to a finer grid.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. Gordon2 (Gordon et al. 2022) and extinction by K. Barbary3 (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 ⊙ 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 αelement abundance [α/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 α 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 where y n and σ 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 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.
Note-Summary of the main characteristics of some frequently-used SSP model suites.This is not a complete list of what SESAMME can accept as the input model suite, but showcases the flexibility of SESAMME to fit the user's preferred set of model assumptions and parameter spaces.
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.

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 Star-burst99.We use the wavelength-dependent emission coefficients given in Ferland (1980) and Aller (1984), which together account for free-free, bound-free, and twophoton emission from H and He + .These values assume typical Case B conditions (densities and temperatures of order ∼100 cm −3 and 10 4 K, respectively, as well as an ionizing photon escape fraction of ∼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 associ-ated with a given stellar population.For BPASS model spectra, for example, we use the tables of hydrogenionizing 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.

Inputs for SESAMME
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 ∼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α absorption profile prior to modeling with SESAMME; for very low-redshift clusters, this could include a Milky Way foreground component in addition to Lyα absorption intrinsic to the target.If the Lyα absorp-tion profile is sufficiently broad, the wings can suppress the continuum around and the apparent strength of N v λλ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 ∼ 3000 optical spectra from Keck/DEIMOS vs. R ∼ 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 ⊙ Å−1 .For example, if the user wishes to use BPASS v2.3 models with a given [α/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 [α/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).

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 fourdimensional walkers at each of the S steps in the MCMC chain generated by emcee -in other words, an array of dimension 4 × N × 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 ∼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 inde-pendent, 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× 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 ∼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 ∼5 iterations per second with the nebular continuum generation turned on (or ∼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.

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 [α/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 α 1 = -1.30,high-mass slope α 2 = -2.35,and a maximum mass of 300 M ⊙ ("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 ±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 ± 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 ∼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 (∆ Å ∼5 -20) at the expected locations of strong ISM emission and absorption features (e.g., OI λ1302, CII λ1334, SiIV λλ1393,1402, CIV λλ1548,1550, and geocoronal Lyα).
Columns 6 -9 of Table 3 show the 16, 50, and 84 th 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 ∼10 -15% uncertainties on the cluster age and mass, stellar metallicity, and average reddening.
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 ∼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.

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 ∼ 30%, and the inferred cluster mass is also overestimated by 0.2-0.3dex 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 χ 2 = 0.97, while the optimal model with the Calzetti curve yields χ 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 ∼ 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.

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 com- pare 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 metalpoor 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 ∼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 stel-lar 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-thanexpected 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.

RE-EXAMINING YMCS IN M83 AND NGC 1313
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 LMCmetallicity 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 lowresolution version of each COS spectrum.The highresolution 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 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).
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 weightcombining 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 ∼5 -7 per resel (0.06 Å for G130M and 0.07 Å for G160M) at 1310 Å. 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 ∼6.8 at the native 0.06 Å/resel resolution of COS (S/N ∼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.
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α 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α absorption profile.This procedure is similar to that used in James et al. (2014) and Hernandez et al. (2020Hernandez et al. ( , 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 λ1247, Si iii λ1294, C ii λ1324, O iv λ1342, and S v λ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 σ s = (σ 2 model − σ 2 COS ).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 ± 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× the literature value.The reddening and amplitude parameters are left essentially unconstrained.
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 [α/Fe] and a broken power law IMF with low-mass slope α 1 = -1.30,high-mass slope α 2 = -2.35,and a maximum mass of 300 M ⊙ .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 ∼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σ 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 singleburst models of a fixed mass of 10 6 M ⊙ , assuming the high mass-loss Geneva stellar tracks and a double-power law IMF with low-mass slope α 1 = -1.30,high-mass slope α 2 = -2.3, and a maximum mass of 100 M ⊙ .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; ∼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).Note-a Assumes a Calzetti et al. (2000) extinction curve with RV = 4.05.
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 ∼6 -8 Myr and reddening E(B − V ) ∼0.28, regardless of whether we use BPASS or Starburst99 SSP models.This stands in contrast to the light-weighted age of ∼30 Myr and slightly lower reddening of ∼0.2 mag determined by Sirressi et al. (2022), who used Star-burst99 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 (∼ 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 linearlyuniform age grid of our Starburst99 model suite (∆age = 1 Myr) compared to that of BPASS (∆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).

DISCUSSION AND SUMMARY
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 VAN-DELS 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 α enhancement at a fixed overall metallicity mass fraction Z (i.e., for a fixed Z, increases in [α/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 restframe 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, 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 λ1240 on the left and Civ λ1550 on the right).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λλ1238, 1242 wind feature, which is known to be present only in fairly young (≲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 ⊙ , 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 ⊙ ) and high (Z = 2Z ⊙ ) 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 ∼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 supersolar metallicity environments.(See also the median model for cluster M83-3 in Fig. A.1, which fails to re-produce 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 * ∼ 10 4 − 10 5 M ⊙ , 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 ⊙ .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 Comparison of N vλλ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 ∼5 Myr, while N v is seemingly absent in all but the youngest and most metal-poor BPASS models we show here.
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 noiseadded BPASS v2.3 stellar population models, show that SESAMME can successfully recover the known age and metallicity of the model, even at moderate signal-tonoise ratios ∼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.

Figure 1 .
Figure1.Visualization of the format of a SESAMME model cube, generated from BPASS v2.3 models with binary physics and a solar α/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⊙ Å−1 .For clarity, only the first ten age intervals and three of the 13 possible metallicities in BPASS are shown.
Figure2.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 ∼6.8 at the native 0.06 Å/resel resolution of COS (S/N ∼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 fromHernandez et al. (2019).The 2D and 1D histograms have been lightly smoothed for clarity.

Figure 4 .
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α and Oi λ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 λ1240 in the bottom left and Civ λ1550 in the bottom right.

Figure 5 .
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 Figure4, 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 λ1240 on the left and Civ λ1550 on the right).

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).
Figure 7.Comparison of N vλλ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 ∼5 Myr, while N v is seemingly absent in all but the youngest and most metal-poor BPASS models we show here.
Figure A.1.Continuation of Fig.4showing 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α 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.

Table 1 .
List of extinction curves available for SESAMME

Table 2 .
Example List of Stellar Population Models

Table 3 .
Cardelli et al. (1989)operties of synthetic validation spectra a Values are equivalent to AV = 0.31 and 3.1 for theCardelli et al. (1989)extinction law with RV = 3.1. b

Table 4 .
Inferred properties of YMCs from HST/COS FUV spectra a