PySSED: Automated Stellar SED Fitting

• PySSED is a Python package for automated construction and fitting of stellar spectral energy distributions from multi-survey catalogue data.
• It performs tasks like cross-matching, quality control, dereddening, and model selection to derive effective temperatures and luminosities for Gaia DR3 stars.
• Its scalable all-sky pipeline processes millions of objects, enabling advanced studies such as machine-learning classification and excess detection.

Searching arXiv for papers on PySSED and closely related uses. arxiv_search(query="PySSED OR GASPS OR \"spectral energy distribution fitting\" PySSED", max_results=10) PySSED is a Python-based software package for automated construction and fitting of stellar spectral energy distributions from heterogeneous survey catalogues. In the current literature it is presented most explicitly as the computational engine behind the Gaia All-Sky Stellar Parameters Service (GASPS), where it performs catalogue collation, cross-matching, quality control, dereddening, model selection, and SED fitting for Gaia DR3 sources with secure parallax-derived distances (McDonald et al., 7 Jan 2026). Its primary products are effective temperatures and luminosities, but its workflow also yields model-versus-observed residual information, ultraviolet and infrared excess diagnostics, and feature sets used in downstream machine-learning studies (Cody et al., 2024). The name should be distinguished from PYSED, a separate Python package for molecular-dynamics phonon analysis based on spectral energy density methods (Liang et al., 1 May 2025).

1. Identity, scope, and nomenclature

PySSED is not the same thing as GASPS. GASPS is the broader project or service whose aim is to collate the SEDs of Gaia DR3 objects and fit those sources with secure parallax-derived distances, whereas PySSED is the Python software package that executes the automated catalogue collation, cross-matching, quality control, dereddening, model selection, and SED fitting. The GASPS paper refers specifically to “PySSED version 1.2” and treats it as the fitting pipeline rather than as a public portal or service layer (McDonald et al., 7 Jan 2026).

The software is designed for a regime in which stellar parameters must be inferred for very large numbers of sources from catalogue photometry rather than from homogeneous spectroscopy. The central scientific motivation is that temperature and luminosity are not directly observed in survey data. Instead, given distance and extinction, SED fitting estimates $T_{\rm eff}$ from the location of the SED peak and angular size from the long-wavelength tail, then converts that to luminosity using distance and the Stefan–Boltzmann law. The GASPS implementation is explicitly all-sky, highly automated, and catalogue-driven.

In the surrounding literature, PySSED also appears as an upstream astrophysical inference engine rather than as an end-user classifier or visualization-only tool. In the stellar-classification study, it supplies physically motivated features to an XGBoost model; in the planetary-nebula study, it is adapted to fit hot central-star continua and line-of-sight extinction (Cody et al., 2024, Csukai et al., 12 Sep 2025).

2. Catalogue assembly and data harmonization

PySSED begins from source identities or source positions and gathers multi-band photometry from major survey catalogues. In the GASPS workflow, Gaia DR3 provides the master source list, while external catalogues are accessed through CDS with astroquery. The paper emphasizes that the practical difficulty is not only fitting but constructing a reliable cross-matched photometric record for each source. PySSED is described as handling one-to-many and many-to-one associations, variable beam sizes, proper-motion correction, bad-data identification, and outlier rejection (McDonald et al., 7 Jan 2026).

The photometric coverage spans ultraviolet through optical and infrared. The cited catalogue set includes GALEX; Hipparcos/Tycho, Pan-STARRS DR1, APASS DR9, SDSS DR16, SkyMapper DR4, CMC15, DES DR2, TASS, and VST surveys in the optical; 2MASS, UKIDSS, and VISTA-based surveys in the near-infrared; and Akari, Spitzer, and WISE products including AllWISE, catWISE, and unWISE in the mid-infrared. In each sky field, Gaia is cross-matched against 25 other catalogues, and the model precomputation is stated to use 26 queried catalogues (McDonald et al., 7 Jan 2026).

Distance assignment is hierarchical. PySSED attempts to use geometric distances from Bailer-Jones et al. (2021), then direct inversion of Gaia DR3 parallaxes, and then direct inversion of Hipparcos parallaxes. Extinction is incorporated through the Lallement et al. (2022) $6\,\mathrm{kpc}\times6\,\mathrm{kpc$ 3D extinction cube, which is the default dereddening source. The code precomputes reddening estimates $A_\lambda/A_V$ at $A_V=0.1, 3.1, 10,$ and $31$ mag, because finite filter widths make the effective band extinction depend on the reddened spectral shape across the bandpass rather than on a single monochromatic correction (McDonald et al., 7 Jan 2026).

The software also regularizes heterogeneous survey uncertainties with catalogue-dependent noise floors. The GASPS paper gives typical floors of 10% for optical data, 20% for near-infrared data, and 25–33% for mid-infrared data. When multiple measurements exist in the same filter, they are merged if consistent within errors; otherwise the more precise value is used. This architecture makes PySSED as much a catalogue-engineering system as a model-fitting system.

3. Physical model selection and fitting procedure

The default stellar-atmosphere basis in PySSED is the BT-Settl grid. The authors selected BT-Settl because it provides a nearly complete grid across a wide range of temperatures, gravities, and metallicities, and because broadband SED fitting is sensitive mainly to large-scale spectral energy redistribution rather than to detailed line lists. For each atmosphere model, PySSED precomputes synthetic photometry by convolving the model spectrum with the filter transmission functions of the queried catalogues (McDonald et al., 7 Jan 2026).

Before the main fit, PySSED constructs an initial parameter guess. Metallicities and $\alpha$-enhancement are taken either from Gaia BP/RP spectra in DR3 or, if unavailable, from the fallback parametrization

$${\rm [Fe/H]} = \tan^{-1} Z, \qquad {\rm [\alpha/Fe]} = \frac{\tan^{-1} Z}{5}.$$

An initial surface gravity is then generated by fitting a blackbody, deriving a radius and luminosity, and inferring $g$ after assuming a stellar mass. The adopted mass prior is piecewise: a zero-age main-sequence mass from PARSEC for main-sequence stars; $1\,M_\odot$ for giants with $T<5500\,\mathrm{K}$ below the RGB tip; and, for AGB stars above the RGB tip, an average of a core mass from Blöcker (1993) and an initial mass from Casewell et al. (2009). The GASPS description adds that final parameters are usually insensitive to the selected $6\,\mathrm{kpc}\times6\,\mathrm{kpc$0 and composition, with temperature changes generally below 100 K even for fairly large departures from truth (McDonald et al., 7 Jan 2026).

The fitting itself is based on reduced $6\,\mathrm{kpc}\times6\,\mathrm{kpc$1, written $6\,\mathrm{kpc}\times6\,\mathrm{kpc$2, evaluated in logarithmic flux “to equalise effects from fluxes that are too faint and too bright.” The best-fit temperature is found with a Nelder–Mead minimizer. Surface gravity is then updated and the fit is repeated. Outlier handling is explicit and iterative: points or sets of points are rejected if their removal improves $6\,\mathrm{kpc}\times6\,\mathrm{kpc$3 by an average factor of 1.5 per removed datum, and up to five outliers can be treated together before the fit is considered stable. If the optimization converges outside the atmosphere grid, PySSED falls back to a blackbody fit, which is stated to yield poorer fits and less accurate parameters than a proper stellar-atmosphere fit (McDonald et al., 7 Jan 2026).

The principal fitted quantities are effective temperature and luminosity. The code also estimates angular scale and, through the adopted distance, a physical radius. Surface gravity and composition function primarily as atmosphere-selection coordinates rather than as the dominant science outputs.

4. All-sky deployment in GASPS

PySSED is engineered for bulk processing. In the GASPS deployment, the sky is partitioned into Galactic-latitude strips of $6\,\mathrm{kpc}\times6\,\mathrm{kpc$4, then into longitude intervals of roughly $6\,\mathrm{kpc}\times6\,\mathrm{kpc$5. Crowded regions are subdivided more finely, using $6\,\mathrm{kpc}\times6\,\mathrm{kpc$6 fields for $6\,\mathrm{kpc}\times6\,\mathrm{kpc$7 and $6\,\mathrm{kpc}\times6\,\mathrm{kpc$8 fields in the densest region $6\,\mathrm{kpc}\times6\,\mathrm{kpc$9. Search radii are bespoke per catalogue and are conservatively chosen to minimize false matches (McDonald et al., 7 Jan 2026).

The code processes at about one second per star. Running it over the entirety of Gaia DR3 required 8.6 CPU years and generated 2.8 TB of data, reduced to 520 GB if only fitted sources are retained. The resulting catalogue covers about 240 million objects, corresponding to about 13.3% of Gaia DR3. The main completeness limit is the requirement of distance precision $A_\lambda/A_V$0, with a secondary requirement of at least four photometric points (McDonald et al., 7 Jan 2026).

Population-level validation is one of the strongest empirical arguments for the pipeline. In a 33 million star subset dominated by low-extinction Galactic-pole sources, the expected HR-diagram structures are recovered: the main sequence from O to L types, the main-sequence turnoff, giant branch, white-dwarf cooling track, red clump, horizontal branch, RGB bump, AGB bump, and RGB tip. The giant branch also splits into metal-poor and metal-rich components, and that morphology is stated to correlate with Gaia metallicity and height above the Galactic plane (McDonald et al., 7 Jan 2026).

PySSED’s outputs in this deployment include fitted photometry, final stellar parameters, and residual diagnostics such as ratios of observed to modelled flux. Those residual products are operationally important because they support excess analysis rather than merely parameter estimation.

5. Derived products and scientific applications

A major use of PySSED is wavelength-dependent excess detection. Because the code assembles broad-band SEDs from the ultraviolet to the mid-infrared and compares them against photospheric models, it can identify ultraviolet excesses associated with hot companions such as accreting stars or young white dwarfs, and infrared excesses associated with circumstellar dust, mass loss, or accretion (McDonald et al., 7 Jan 2026).

In the stellar-classification study, PySSED is a central upstream component rather than the classifier itself. It retrieves coordinates from SIMBAD, queries VizieR for multi-wavelength photometry and other parameters, corrects the photometry for interstellar reddening, fits a stellar-atmosphere model, estimates effective temperature and luminosity, and records the differences between observed and model fluxes. The default PySSED output used there contained 442 feature columns before preprocessing, including original photometry, reddening-corrected photometry, model-fit photometry, effective temperature, luminosity, sky position, and additional parameters (Cody et al., 2024).

That study transformed the model residuals into per-band significance features using

$A_\lambda/A_V$1

where $A_\lambda/A_V$2 is the observed flux, $A_\lambda/A_V$3 is the modelled flux, and $A_\lambda/A_V$4 is the quoted observational uncertainty. After preprocessing and feature engineering, the dataset was reduced from 442 to 96 retained features. The paper reports that the main classifier reached an accuracy of about 0.7 with a macro F1 score of 0.61, and it explicitly concludes that PySSED’s residual spectral features carry information not contained in $A_\lambda/A_V$5 alone (Cody et al., 2024).

A second major adaptation appears in the study of central stars of planetary nebulae. There PySSED is modified away from its ordinary BT-Settl usage and coupled instead to TMAP non-LTE atmospheres in order to fit reddened hot-star continua. For 162 planetary nebulae with identified central stars, the authors use PySSED-built SEDs and fit reddened model spectra with Nelder–Mead $A_\lambda/A_V$6 minimisation, evaluating uncertainties through MCMC; they report extinctions accurate down to $A_\lambda/A_V$7 of $A_\lambda/A_V$8 mag, identify four new close binaries from the SEDs, and infer a binary fraction between 23% and 36% for the full sample (Csukai et al., 12 Sep 2025). This application shows that PySSED is not confined to its default atmosphere grid, but can function as a general SED-reduction and fitting framework when the model family is replaced.

6. Limitations, failure modes, and adjacent software

The dominant limitation identified in the GASPS deployment is extinction correction. The Lallement cube is coarse, with 50 pc resolution, and can overestimate extinction for stars lying just in front of mapped dust structures. In HR diagrams this appears as diagonal smearing away from expected loci, particularly for bright blue stars that remain visible to large distances and therefore encounter larger cumulative extinctions. Poor fits also arise for intrinsically variable or dusty stars near the RGB tip and on the AGB, because static stellar-photosphere models do not capture variability or circumstellar dust emission (McDonald et al., 7 Jan 2026).

The pipeline also exhibits recognizable grid and optimization artifacts. The GASPS paper notes vertical banding when poor SED coverage or noisy data cause the optimizer to settle on atmosphere grid points, and a diagonal artifact near $A_\lambda/A_V$9 associated with Gaia-derived starting temperatures and a model-grid temperature node. The paper is equally explicit that formal posterior uncertainty estimation is not part of the baseline all-sky workflow; the operating model is based on catalogue filtering, adopted noise floors, and reduced-$A_V=0.1, 3.1, 10,$0 optimization rather than on MCMC or evidence calculations (McDonald et al., 7 Jan 2026).

In the machine-learning setting, missingness remains severe. After preprocessing, the PySSED-derived dataset had 74% total missing data; only 27 of 96 features were present for at least 50% of samples. The authors therefore relied on XGBoost’s native missing-value handling rather than on imputation, and they suggest that a future direction is physics-informed filling of photometric gaps using the stellar-atmosphere models already used by PySSED (Cody et al., 2024). A second stated future direction is broader application of an “identical reduction” to Galactic globular clusters and resolved external galaxies where distances and foreground reddening are externally constrained; a third is simultaneous fitting of temperature, extinction, and eventually the extinction-curve slope

$A_V=0.1, 3.1, 10,$1

with the aim of mineralogical mapping of Galactic dust (McDonald et al., 7 Jan 2026).

Methodologically related software helps define PySSED’s niche. SEDBYS is another Python-based SED tool, but it is oriented toward young-star workflows, emphasizing archive querying, zero-point-based flux conversion, interactive inspection, infrared spectral retrieval, and bibTeX-oriented provenance capture rather than all-sky stellar-parameter fitting at Gaia scale (Davies, 2020). At the opposite extreme, the similarly named PYSED belongs to condensed-matter simulation and computes kinetic-energy-weighted phonon dispersion and phonon lifetimes from molecular-dynamics trajectories; despite the orthographic similarity, it is unrelated in data model, scientific domain, and methodology (Liang et al., 1 May 2025).