Pysersic: Bayesian Galaxy Fitting

• Pysersic is a Python package that uses Bayesian inference with JAX acceleration to determine galaxy structural properties from imaging data.
• It fits parametric galaxy light profiles directly while accounting for PSF convolution and uses advanced Fourier-space rendering for improved efficiency.
• It has been validated across multiple observational regimes, making it a powerful tool for survey-scale analysis and wavelength-dependent morphology studies.

Searching arXiv for papers on pysersic and recent uses. pysersic is a Python package for determining galaxy structural properties via Bayesian inference, accelerated with jax. It was introduced as a framework for fitting parametric galaxy light profiles directly to imaging data, with explicit treatment of the point-spread function and posterior uncertainty, and it is designed for both interactive use and integration into survey-scale analysis pipelines (Pasha et al., 2023). Across subsequent work, it has been used for single-source fitting, multi-source fitting, wavelength-dependent morphology, catalog production, and methodological development for faster Fourier-space rendering (Pasha et al., 2023, Miller et al., 27 Aug 2025).

1. Software definition and scope

At its core, pysersic addresses a standard task in extragalactic astronomy: inferring structural parameters such as total flux, effective radius, Sérsic index, ellipticity, position angle, and centroid from galaxy images. The package is written in pure Python and built using jax, which provides just-in-time compilation, auto-differentiation, and execution on CPUs, GPUs, or TPUs. Inference is performed with numpyro, enabling gradient-based posterior estimation through methods such as No-U-Turn Sampling, as well as optimization-based point estimation and stochastic variational inference (Pasha et al., 2023).

The software was presented as a Bayesian alternative within the ecosystem of profile-fitting tools that also includes GALFIT, IMFIT, ProFit, GaLight, PetroFit, and PyAutoGalaxy. Its distinguishing emphasis is the combination of probabilistic programming, hardware acceleration, and direct access to full posterior distributions rather than only point estimates (Pasha et al., 2023). This suggests that pysersic is intended not merely as a fitting executable but as an inference framework whose outputs can be propagated into downstream population-level analyses.

The package has since been used in distinct observational regimes: HST imaging of post-starburst galaxies at $1Zhang et al., 2024), optical and ultraviolet morphology from JWST at $4Miller et al., 2024), low-mass galaxy size measurements across environments in the SAGA survey (Asali et al., 29 Sep 2025), and large morphology catalogs in Abell 2744 (Zhang et al., 27 Feb 2026). Those applications establish pysersic as a general-purpose structural fitting framework rather than one tied to a single instrument or redshift range.

2. Parametric model and computational methods

The standard profile used by pysersic is the Sérsic law,

$$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$

where $I(R)$ is the surface brightness at radius $R$, $I_e$ is the surface brightness at the effective radius $R_e$, $R_e$ is the half-light radius, $n$ is the Sérsic index, and $b_n$ is chosen so that half the total luminosity lies within $4Zhang et al., 2024, Pasha et al., 2023). In practical fits, the free parameters commonly include flux, $4Zhang et al., 2024, Asali et al., 29 Sep 2025).

A central computational issue in Sérsic fitting is rendering the model accurately enough for PSF convolution and likelihood evaluation. pysersic was described with multiple rendering strategies: oversampled real-space rendering, Fourier-space rendering via a sum of Gaussians, and a hybrid real-Fourier algorithm intended to improve speed and reduce aliasing while handling steep inner profiles (Pasha et al., 2023). Because jax supplies automatic differentiation, these rendering pipelines are compatible with gradient-based optimization and sampling (Pasha et al., 2023).

A later development replaced expensive approximations to the radial Fourier transform with a symbolic-regression emulator implemented in pysersic. In injection-recovery tests on 500 synthetic galaxies, the emulator achieved residual scatter in $4$$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$0 scatter in $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$1. The median MCMC runtime per galaxy decreased from 64 seconds to 23.5 seconds, a $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$2 speed-up (Miller et al., 27 Aug 2025). A plausible implication is that rendering, rather than only inference, has become a primary optimization target for survey-scale morphology pipelines.

3. Bayesian inference workflow

The generic pysersic workflow consists of providing an image cutout and, optionally, a PSF; configuring priors and the model; running inference; and inspecting best-fit parameters and posterior samples (Pasha et al., 2023). The package supports automatic priors derived from the image as well as user-specified priors, and the outputs include best-fit parameters, posterior samples, and diagnostic products (Pasha et al., 2023).

Published applications show how this abstract workflow is operationalized. In HST/WFC3 F110W imaging of DESI post-starburst galaxies, source extraction and deblending were performed with SEP, segmentation maps were used for masking, empirical PSFs were constructed to account for spatial variation and drizzling artifacts, and all unmasked sources plus the sky were fitted simultaneously in each cutout (Zhang et al., 2024). In the SAGA size-mass study, pysersic fits were performed independently in each of the $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$3, $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$4, and $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$5 bands, with contaminating sources masked using sep and priors informed by image statistics from photutils (Asali et al., 29 Sep 2025). In the Abell 2744 catalog, sources were fit independently in each band using empirical PSFs, simultaneous modeling of bright neighbors when needed, and two inference regimes: SVI-MVN for all $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$6 sources and NUTS additionally for $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$7 objects (Zhang et al., 27 Feb 2026).

The package’s Bayesian formulation matters scientifically because it makes uncertainty and parameter covariance first-class outputs rather than afterthoughts. That property was explicitly used in applications requiring error propagation into scaling-relation fits, injection-recovery validation, or comparison between subtle environmental or evolutionary effects (Asali et al., 29 Sep 2025, Miller et al., 2024, Zhang et al., 27 Feb 2026).

4. Multiband fitting and survey-scale extensions

pysersic has also been extended beyond independent single-band fits. In JWST UNCOVER and MegaScience analyses of galaxies at $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$8, a multiband module was developed from v0.1.4 onward to fit six NIRCam wide filters simultaneously. In that implementation, central position and position angle were fixed across bands, each band’s total flux was fitted independently, and effective radius, Sérsic index, and ellipticity were allowed to vary smoothly with wavelength through a 2nd-order B-spline with 3 knots (Miller et al., 2024). The joint model increased the precision of measured radii by $$I(R) = I_e \exp\left\{ -b_n \left[ \left(\frac{R}{R_e}\right)^{1/n} - 1 \right] \right\},$$9 (Miller et al., 2024).

Large-catalog work has stressed a different form of extensibility: scale. In Abell 2744, pysersic was used to derive spectrally resolved structural measurements for 29,608 sources over observed $I(R)$0, mainly from 20 JWST NIRCam broad and medium bands. Through quality evaluation criteria, robust structural parameters were inferred for more than $I(R)$1 of selected $I(R)$2 sources, and injection-recovery tests established fractional error below $I(R)$3 above $I(R)$4 (Zhang et al., 27 Feb 2026). The same study used the fitted size-wavelength relations to infer rest-frame UV, optical, and near-infrared sizes (Zhang et al., 27 Feb 2026).

Independent-band fitting remains important where direct wavelength coupling is not desired. In JADES imaging used to study galaxies at $I(R)$5, pysersic was applied band by band with stochastic variational inference and posterior checks against MCMC, producing continuum sizes, Sérsic indices, axis ratios, and position angles that were later compared across rest-frame FUV, NUV, and optical bands (Danhaive et al., 7 Oct 2025). The coexistence of these strategies indicates that pysersic functions both as a single-image fitter and as a substrate for higher-level modeling of wavelength-dependent morphology.

5. Scientific applications

The package’s scientific role is best understood through representative studies.

Study	Dataset or regime	pysersic contribution
DESI post-starburst galaxies (Zhang et al., 2024)	171 massive PSBs at $I(R)$6 with HST/WFC3 F110W	Structural fitting showed PSBs are $I(R)$7 dex smaller than quiescent galaxies, have $I(R)$8, $I(R)$9, and $R$0
SAGA Survey VI (Asali et al., 29 Sep 2025)	Low-mass galaxies across environments	Uniform fits showed SAGA satellites are larger than isolated galaxies by $R$1 to $R$2 dex, with similar Sérsic index distributions
UNCOVER high-$R$3 sizes (Miller et al., 2024)	995 galaxies at $R$4 from JWST	Simultaneous multiband fitting yielded $R$5 better radius precision and supported growth from 400 pc at $R$6 to 830 pc at $R$7 at $R$8
Abell 2744 global catalog (Zhang et al., 27 Feb 2026)	29,608 sources over 20 NIRCam bands	Produced quality-flagged posteriors and rest-frame size estimates over $R$9
Cosmic-noon quiescent galaxies (Prasal et al., 4 May 2026)	14 massive quiescent galaxies at $I_e$0	Morphological analysis found median $I_e$1 and mean $I_e$2 kpc

In the DESI post-starburst study, pysersic underpinned the conclusion that massive PSBs are more compact than the general quiescent population while having similar central surface mass densities, which the authors reconciled through later ex-situ growth via minor mergers or slight progenitor bias (Zhang et al., 2024). In SAGA, a uniform pysersic pipeline was necessary because the environmental size offsets under investigation were as small as $I_e$3–$I_e$4 dex, making method-dependent systematics a major concern (Asali et al., 29 Sep 2025).

At very high redshift, the package enabled direct population-level size-mass analyses. The UNCOVER work used joint multiband modeling to show rapid growth in average size for low-mass galaxies in the first billion years, while the JADES-based study in $I_e$5 galaxies found no significant differences between stellar continuum sizes at FUV, NUV, and optical wavelengths, with $I_e$6 kpc at $I_e$7 (Miller et al., 2024, Danhaive et al., 7 Oct 2025). These results illustrate that pysersic has been used both to detect wavelength dependence and to establish its absence.

6. Validation, limitations, and interpretive cautions

The package has repeatedly been validated against external methods and simulated data. In the DESI post-starburst analysis, cross-checks with single-component GALFIT fits found no significant systematic differences, with median $I_e$8 dex (Zhang et al., 2024). In the SAGA analysis, comparisons with Tractor, NSA, and SGA showed only very small global offsets below $I_e$9 dex, and $R_e$0 of pysersic sizes were within $R_e$1 dex of Tractor values (Asali et al., 29 Sep 2025). In the Abell 2744 catalog, injection-recovery tests showed systematic bias below $R_e$2 for $R_e$3 (Zhang et al., 27 Feb 2026).

At the same time, the literature emphasizes important limits. Many studies fit all galaxies with a single-component Sérsic model; this is standard for statistical work, but it cannot capture complex morphologies, barred systems, or strongly clumpy structures in detail (Asali et al., 29 Sep 2025, Danhaive et al., 7 Oct 2025). Some applications therefore added residual corrections or adopted two-component models. In the DESI post-starburst study, a residual-correction method based on adding the masked residual profile to the analytic Sérsic profile altered sizes only slightly, with $R_e$4 arcsec (Zhang et al., 2024). In the cosmic-noon quiescent-galaxy analysis, pysersic was used with a Sérsic-plus-exponential two-component model to extract bulge-to-total ratios alongside Sérsic indices and axis ratios (Prasal et al., 4 May 2026).

The most explicit caution concerns extremely compact JWST “little red dots.” Simulations showed that pysersic and GALFIT have difficulty recovering input parameters at $R_e$5, that estimated PSF fraction may be a more robust physically motivated description of compactness, and that almost all modeled LRDs with $R_e$6 cannot be differentiated from a point source regardless of intrinsic extent (Whalen et al., 25 Sep 2025). That work reported that $R_e$7 of observed LRDs are PSF-dominated and $R_e$8 show lower estimated PSF fractions or fitted sizes $R_e$9 pc, implying morphological diversity rather than a homogeneous population (Whalen et al., 25 Sep 2025).

A further caution arises from uncertainty calibration. In the Abell 2744 catalog, the width of the posterior distribution underestimated the true error by a factor of $R_e$0 for flux and size and $R_e$1 for Sérsic index in injection-recovery tests, and users were advised to inflate the reported $R_e$2 uncertainties accordingly (Zhang et al., 27 Feb 2026). This suggests that posterior precision in practical survey pipelines must be interpreted together with empirical calibration, especially near the detection limit or under model mismatch.