PyStaff: Python Stellar Absorption Fitting

• PyStaff is a Python FSF code that fits integrated-light spectra using a single SSP model to derive stellar and gas kinematics, IMF slope, and elemental abundances.
• It employs Bayesian MCMC sampling with emcee and multiplicative Legendre polynomials to robustly match observed spectra while mitigating continuum mismatches.
• Comparative studies show PyStaff yields precise results for pure SSPs but may introduce biases when applied to galaxies with composite star formation histories.

PyStaff, short for Python Stellar Absorption Feature Fitting, is a Python full spectral fitting (FSF) code for integrated-light galaxy spectroscopy that, in the configuration examined by the comparative study of elliptical galaxies in the Fornax cluster, derives stellar kinematics, gas kinematics, and stellar population parameters under a single simple stellar population (SSP) assumption. In that study, PyStaff is evaluated alongside ALF, Starlight, and pPXF using optical plus near-infrared spectra of NGC 1399 and NGC 1404, with particular emphasis on the robustness of inferred low-mass IMF slope, age, metallicity, and 19 elemental abundances (Lonoce et al., 22 Aug 2025).

1. Functional scope and scientific role

PyStaff is designed to fit integrated-light galaxy spectra and recover several classes of parameters from a common forward model. The quantities explicitly fitted in the study are stellar radial velocity $V$, stellar velocity dispersion $\sigma$, gas velocity $V_{\rm gas}$, gas velocity dispersion $\sigma_{\rm gas}$, SSP age, total metallicity $[Z/H]$, 19 elemental abundances, and a low-mass IMF slope over $0.08$–$1.0\,M_\odot$ (Lonoce et al., 22 Aug 2025).

The elemental abundance vector comprises Fe, O, C, Ca, Na, N, Mg, Ti, Si, K, V, Cr, Mn, Co, Ni, Cu, Sr, Ba, and Eu. Within the study, PyStaff is used to measure radial trends in these quantities for the two brightest galaxies in the Fornax cluster. The spectra combine optical coverage from $4040$–$6640$ Å, with two slightly different optical range configurations tested, and near-infrared coverage in the Ca triplet region, $8143$–$\sigma$0 Å (Lonoce et al., 22 Aug 2025).

PyStaff is operated in a full spectral fitting mode in which each “good” pixel contributes to the likelihood. Continuum mismatches between observed and model spectra are not interpreted as stellar population information directly; instead, they are absorbed through multiplicative polynomials. This places PyStaff within the FSF class rather than line-index-based approaches, and the comparative study treats it as one of four widely used codes for IMF-sensitive stellar population analysis in elliptical galaxies (Lonoce et al., 22 Aug 2025).

The implementation referenced in the study is available on GitHub at https://github.com/samvaughan/PyStaff. A plausible implication is that PyStaff is intended not only for kinematic extraction but for joint chemo-dynamical inference in high-quality galaxy spectra, provided the SSP approximation is defensible.

2. Population model, abundance treatment, and IMF parameterization

In the setup tested, PyStaff is a parametric FSF code. It assumes that the light at a given radius is described by a single SSP characterized by one age, one total metallicity, one low-mass IMF slope, one abundance pattern for 19 elements, and one stellar kinematic state $\sigma$1 (Lonoce et al., 22 Aug 2025). It does not, in this use case, fit composite star formation histories or multiple age components. That distinction is central to its contrast with Starlight and pPXF, which fit superpositions of SSPs.

PyStaff uses the Conroy et al. (2018) SSP models, with wavelength coverage from $\sigma$2–$\sigma$3 Å, spectral resolution $\sigma$4 km s$\sigma$5, an age grid up to 13.5 Gyr, and a metallicity grid spanning $\sigma$6 to $\sigma$7 dex (Lonoce et al., 22 Aug 2025). Non-solar abundances are implemented through response functions from Conroy et al. (2018), computed over the same wavelength range and resolution for a Kroupa IMF. PyStaff applies a Taylor expansion around solar abundance values to interpolate or extrapolate these response functions; the comparison paper explicitly contrasts this with ALF’s linear interpolation.

The abundance fitting bounds are typically $\sigma$8 to $\sigma$9 dex, with two stated exceptions: $V_{\rm gas}$0 is allowed over $V_{\rm gas}$1 to $V_{\rm gas}$2, and $V_{\rm gas}$3 over $V_{\rm gas}$4 to $V_{\rm gas}$5. Some abundance parameters are expressed relative to Fe, for example $V_{\rm gas}$6 (Lonoce et al., 22 Aug 2025).

For the IMF, the adopted SSP models parameterize the mass function as a broken power law,

$V_{\rm gas}$7

with $V_{\rm gas}$8 applying over $V_{\rm gas}$9–$\sigma_{\rm gas}$0, $\sigma_{\rm gas}$1 over $\sigma_{\rm gas}$2–$\sigma_{\rm gas}$3, and a fixed Salpeter slope $\sigma_{\rm gas}$4 above $\sigma_{\rm gas}$5 (Lonoce et al., 22 Aug 2025). Because previous work showed strong degeneracies between $\sigma_{\rm gas}$6 and $\sigma_{\rm gas}$7, the comparison paper fixes $\sigma_{\rm gas}$8 and treats them as a single IMF slope parameter,

$\sigma_{\rm gas}$9

In PyStaff as used there,

$[Z/H]$0

with Salpeter retained above $[Z/H]$1 (Lonoce et al., 22 Aug 2025).

Within this parameterization, “top-heavy” denotes $[Z/H]$2, whereas “bottom-heavy” denotes $[Z/H]$3. PyStaff fits only one IMF slope per radial spectrum. It does not permit distinct IMF slopes for hypothetical subcomponents, so any genuinely composite stellar population is compressed into a single “effective” SSP description.

3. Inference procedure and technical configuration

PyStaff borrows its kinematic handling from pPXF methodology, but the actual parameter inference is performed through Bayesian posterior sampling with the MCMC ensemble sampler emcee (Lonoce et al., 22 Aug 2025). The likelihood is effectively a pixel-by-pixel $[Z/H]$4 comparison between observed and model fluxes using the per-pixel uncertainties:

$[Z/H]$5

The corresponding posterior is written as

$[Z/H]$6

where the priors $[Z/H]$7 are defined implicitly by the allowed parameter bounds (Lonoce et al., 22 Aug 2025).

To absorb continuum-shape differences from flux calibration or instrumental effects, the spectrum is divided into four wavelength sections:

Section	Wavelength range
Range 1	4040–4975 Å
Range 2	4975–5843 Å
Range 3	5843–6640 Å
Range 4	8143–8600 Å

An alternative slightly shifted “Range 2” configuration is also tested. In each section, PyStaff fits a multiplicative Legendre polynomial to the model spectrum, with order

$[Z/H]$8

The absolute continuum is therefore not used directly as a stellar population constraint (Lonoce et al., 22 Aug 2025).

For computational efficiency and reduced per-pixel noise, the spectra are resampled to 1.25 Å per pixel. The comparison paper notes this resampling as a contributor to PyStaff’s very narrow posterior distributions, which appear as very precise outputs (Lonoce et al., 22 Aug 2025).

The code also models gas emission as a constrained component. Balmer lines H$[Z/H]$9, H$0.08$0, H$0.08$1, and H$0.08$2 are tied by fixed ratios,

• H$0.08$3 H$0.08$4
• H$0.08$5 H$0.08$6
• H$0.08$7 H$0.08$8

Forbidden lines such as [O III], [N I], and [N II] are included with fixed doublet ratios, with one amplitude fitted per species. For example, [O III] 4960 Å is fixed to $0.08$9 [O III] 5008 Å (Lonoce et al., 22 Aug 2025). Gas kinematics are constrained relative to the stellar solution: $1.0\,M_\odot$0 lies within $1.0\,M_\odot$1 km s$1.0\,M_\odot$2 of $1.0\,M_\odot$3, and $1.0\,M_\odot$4 within $1.0\,M_\odot$5 km s$1.0\,M_\odot$6. The stellar bounds listed for the fits are $1.0\,M_\odot$7 km s$1.0\,M_\odot$8 and $1.0\,M_\odot$9 km s$4040$0 (Lonoce et al., 22 Aug 2025).

The MCMC configuration used in the study comprises 100 walkers, 10,000 steps, and a burn-in of 5,000 steps. Reported parameter values are the means of the converged posterior distributions, with uncertainties combining statistical errors from the chains and systematic variations from different wavelength-range fits (Lonoce et al., 22 Aug 2025).

4. Simulation results: accuracy, precision, and failure modes

The central methodological conclusion of the comparison study is explicit: codes such as ALF and PyStaff, both operating under an SSP assumption, return greater precision and accuracy only when the underlying population is a pure SSP; when the star formation history is more complex, they return erroneous results (Lonoce et al., 22 Aug 2025).

In the SSP simulations, where the mock spectra are generated from old, metal-rich populations with a single IMF slope, all four codes recover age and $4040$1 well. For the IMF slope, PyStaff and ALF recover the correct $4040$2 with high accuracy and high precision, whereas Starlight and pPXF yield broader and often biased IMF posteriors. The paper further reports that pPXF shows double-peaked IMF distributions in top-heavy cases, while Starlight tends to shift $4040$3 by approximately $4040$4 in such cases. For the abundance vector, PyStaff and ALF retrieve nearly all 19 elements accurately in the SSP tests, with only small deviations noted, including Na in ALF (Lonoce et al., 22 Aug 2025).

The behavior changes sharply in the composite SFH simulations. Three classes are described: CSP, with 80% mass in an old, metal-rich component and 20% in a younger, metal-poor component sharing the same IMF; CSP2, with the same age–metallicity structure but different IMF slopes in the two components; and CSP3, a 50%–50% mixture of two old components with very different metallicities and IMF slopes (Lonoce et al., 22 Aug 2025).

For CSP models with the same IMF but different ages and metallicities, PyStaff tends to recover an age between the true component ages, such as around 10 Gyr rather than 13.5 and 7 Gyr, and near-solar metallicity rather than the true high- and low-metallicity values. In bottom-heavy cases the IMF can still be recovered reasonably, but in top-heavy cases the recovered IMF slope is biased by approximately $4040$5–$4040$6. The error bars remain small, producing high precision but misleading accuracy (Lonoce et al., 22 Aug 2025).

For CSP2 and CSP3, where the components have different IMF slopes, the biases become more severe. Age, metallicity, and IMF are all reported to lie between the true component values, and even the main-component age and metallicity are badly biased. The elemental abundances are generally recovered well only under the explicitly unrealistic assumption that both components share the same abundances (Lonoce et al., 22 Aug 2025).

By contrast, Starlight and pPXF, although noisier and often more biased in SSP-only tests, are better at exposing multi-component structure. Starlight is reported to recover the ages, metallicities, and IMF slopes of both components accurately, albeit with broad error distributions, while pPXF detects two components but often returns biased metallicities and IMF distributions that can become double-peaked and extreme (Lonoce et al., 22 Aug 2025).

A recurrent theme is therefore the distinction between precision and model adequacy. In PyStaff, narrow posteriors do not by themselves demonstrate that the SSP model is physically correct. In composite populations, they can instead reflect a tightly constrained but non-physical compromise solution.

5. Radial inferences for NGC 1399 and NGC 1404

Applied to radial spectra extending to approximately $4040$7, PyStaff and ALF yield consistent SSP-based descriptions for both NGC 1399 and NGC 1404 (Lonoce et al., 22 Aug 2025). The most prominent IMF result is a radial profile that is approximately flat, with $4040$8 over $4040$9–$6640$0. The study describes this as a “flat super-Salpeter IMF”:

$6640$1

Under the SSP interpretation, both galaxies thus appear consistently bottom-heavy in the low-mass regime, with no significant radial IMF gradient (Lonoce et al., 22 Aug 2025).

The inferred age structure is similarly simple in the PyStaff solution. NGC 1399 is assigned an age of approximately 13.5 Gyr across most radii, and NGC 1404 is described as old with a roughly flat age profile, aside from some spread in outer bins. Metallicity shows a negative gradient in both systems, from approximately $6640$2 dex in the center to approximately $6640$3 dex near $6640$4 (Lonoce et al., 22 Aug 2025).

PyStaff also returns radial trends for 19 elemental abundances. The study reports mild gradients overall, with several notable patterns: Na shows a strong negative gradient, Mg a positive gradient in both galaxies, and Mn a negative gradient. Other elements, including O and N, show code-dependent differences. The paper emphasizes that differences between PyStaff and ALF in Na and some other elements are attributed to extrapolation behavior and degeneracies rather than to a fundamental failure of either code (Lonoce et al., 22 Aug 2025).

The same galaxies admit a different interpretation when fitted with Starlight and pPXF. Those non-SSP fits identify two distinct old components in both NGC 1399 and NGC 1404. The main component is old ($6640$5 Gyr), has super-solar metallicity around $6640$6 dex, and has an IMF that is top-heavy in the center and becomes more bottom-heavy with radius. The secondary component is also old, somewhat younger in some fits, has sub-solar metallicity $6640$7 dex, and has a roughly constant bottom-heavy IMF at the upper end of the allowed range, $6640$8–$6640$9. Its fraction increases with radius, reaching approximately 50% around $8143$0 (Lonoce et al., 22 Aug 2025).

In the authors’ interpretation, PyStaff’s flat super-Salpeter IMF and negative metallicity gradient can be understood as a weighted compromise imposed by the SSP assumption. The metallicity gradient, for example, is interpreted not as a gradient within a single component but as the changing mixture of a metal-rich and a metal-poor population whose individual metallicities are themselves approximately flat in the non-parametric fits (Lonoce et al., 22 Aug 2025).

6. Comparative placement and methodological use

The comparison paper positions PyStaff relative to ALF, Starlight, and pPXF primarily through their SFH assumptions and the resulting bias–scatter tradeoff (Lonoce et al., 22 Aug 2025).

Code class in the study	SFH assumption	Reported behavior
PyStaff and ALF	Fixed SSP	Precise; accurate for SSPs; biased for composite SFHs
Starlight and pPXF	Free SFH as SSP combinations	Noisier; better at revealing secondary components

Within this framework, PyStaff and ALF are described as strong tools for IMF and abundance work when a galaxy is known, or strongly suspected, to be close to an SSP. Under those conditions they provide precise IMF slopes and detailed abundance patterns. The study’s simulations support that use case directly for age, metallicity, IMF slope, and nearly all elemental abundances (Lonoce et al., 22 Aug 2025).

For systems with potentially complex star formation histories, the paper recommends a different workflow. A non-parametric FSF code, particularly Starlight and in some contexts pPXF with caution, should first be used to check for secondary components. If the inferred SFH is close to a pure SSP, PyStaff can then be applied to refine IMF and abundance measurements. If multiple components are present, PyStaff’s outputs should be interpreted as effective SSP parameters rather than as the physical properties of any single component (Lonoce et al., 22 Aug 2025).

This comparative result addresses a common misconception in full spectral fitting: apparently small formal uncertainties are not, by themselves, evidence of physical fidelity. In the PyStaff use case studied here, narrow posteriors can coexist with strong systematic bias whenever the single-SSP assumption is violated. The final astrophysical interpretation adopted by the authors for NGC 1399 and NGC 1404 is therefore based primarily on Starlight, which they identify as showing the least bias in their composite-SFH simulations (Lonoce et al., 22 Aug 2025).