PyBlastAfterglow: GRB Afterglow Modeling

• PyBlastAfterglow is an open-source, modular code for simulating gamma-ray burst afterglows using a thin-shell model that self-consistently incorporates forward and reverse shocks.
• It numerically evolves the downstream electron distribution to compute observables such as light curves, spectra, and two-dimensional sky maps, bridging fast semi-analytic and full simulation approaches.
• The framework supports diverse jet structures, off-axis viewing, synchrotron self-absorption, and SSC effects, providing a versatile tool for benchmark studies and Bayesian inference applications.

Searching arXiv for PyBlastAfterglow and closely related afterglow modeling papers to ground the article in current literature. PyBlastAfterglow is an open-source, modular gamma-ray-burst afterglow code with a C++ core and Python interface for computing light curves, spectra, and sky maps of relativistic outflows, especially GRB jets. It is presented as occupying the middle ground between very fast but restrictive semi-analytic afterglow models and expensive hydrodynamic/MHD+radiation simulations. In its current formulation, it adopts a thin-shell dynamical model, includes forward and reverse shocks self-consistently, supports lateral structure and off-axis viewing, and implements synchrotron, synchrotron self-absorption, synchrotron self-Compton, pair production, and extragalactic background light attenuation, with outputs that include dynamical histories, observer-frame spectra, multi-band light curves, and imaging observables (Nedora et al., 2024).

1. Conceptual scope and development

The physical problem addressed by PyBlastAfterglow follows the standard decomposition of afterglow theory into two coupled modules: hydrodynamics of ejecta deceleration and radiative processes in the shocked plasma. In the standard external-shock picture, relativistic ejecta interact with the circumburst medium and produce a forward shock in the external matter, a reverse shock in the ejecta, and a contact discontinuity between the two shocked regions; the baseline afterglow model then identifies the observed emission primarily with synchrotron radiation from shock-accelerated electrons in the forward shock (Godet et al., 2012). PyBlastAfterglow encodes that backbone but is explicitly designed to extend beyond the simplest forward-shock-only implementation by including reverse-shock emission, structured jets, off-axis viewing, synchrotron self-absorption, synchrotron self-Compton, and sky-map generation (Nedora et al., 2024).

The code is also described, in an earlier application, as a semi-analytic afterglow model and as an extension of the version introduced in Nedora et al. 2021. In that use, it served as the central computational tool for kilonova afterglows with angle-dependent and velocity-dependent ejecta imported from numerical-relativity simulations, together with microphysical prescriptions and equal-arrival-time observer synthesis. That application showed that the framework is not limited to canonical long-GRB afterglows, but can also be adapted to merger ejecta, thermal-electron emission, and blast waves propagating through a medium modified by prior GRB ejecta (Nedora et al., 2022).

2. Dynamical formulation

PyBlastAfterglow models a jet as one or many angular blast-wave layers. Each angular segment is treated as a uniform thin shell containing forward-shock-shocked external medium, reverse-shock-shocked ejecta, and a contact discontinuity. The fluid stress-energy tensor is written as

$$T^{\mu\nu} = (\rho' c^2 + e' + p')u^{\mu}u^{\nu} + p'\eta^{\mu\nu},$$

with pressure related to internal energy through

$p' = (\hat{\gamma}-1)e',$

and an adiabatic index approximated as

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$

The shell total energy is then written as

$$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$

where $m$ is shell mass, $E_{\rm int}'$ is comoving internal energy, and $\Gamma_{\rm eff}$ is the effective Lorentz factor (Nedora et al., 2024).

The full forward-plus-reverse-shock system is treated through four standard regions: unshocked external medium, shocked external medium, shocked ejecta, and unshocked ejecta. Internal energies in the forward- and reverse-shock downstreams evolve through shock heating, adiabatic losses, and radiative losses,

$$dE_{\rm int;2}' = dE_{\rm sh;2}' + dE_{\rm ad;2}' + dE_{\rm rad;2}', \qquad dE_{\rm int;3}' = dE_{\rm sh;3}' + dE_{\rm ad;3}' + dE_{\rm rad;3}'.$$

For the forward shock,

$dE_{\rm sh;2}' = (\Gamma - 1)c^2 dm,$

where $dm$ is the swept-up external mass, while for the reverse shock

$p' = (\hat{\gamma}-1)e',$0

with $p' = (\hat{\gamma}-1)e',$1 the relative Lorentz factor between unshocked and shocked ejecta (Nedora et al., 2024).

Lateral spreading is imposed by prescription rather than by solving multidimensional lateral hydrodynamics. The comoving sound speed is

$p' = (\hat{\gamma}-1)e',$2

and a basic spreading law is

$p' = (\hat{\gamma}-1)e',$3

The default implementation adds a structure-aware correction intended to mimic two-dimensional jet behavior more closely than a simple conical expansion law (Nedora et al., 2024).

The explicit reverse-shock treatment is physically consequential because long-lived reverse shocks in stratified ejecta can drive plateaus, steep declines, bumps, re-brightenings, and a wide range of temporal decay indices, whereas forward-shock light curves are comparatively insensitive to ejecta stratification. That broader reverse-shock phenomenology was established in detailed mechanical-model calculations of long-lived RS blast waves, and provides a natural theoretical context for PyBlastAfterglow’s inclusion of self-consistent forward-plus-reverse-shock dynamics rather than an FS-only closure (1208.01558).

3. Shock microphysics and radiation modules

PyBlastAfterglow uses standard afterglow microphysical parameters $p' = (\hat{\gamma}-1)e',$4, $p' = (\hat{\gamma}-1)e',$5, and $p' = (\hat{\gamma}-1)e',$6, with separate parameter sets available for forward and reverse shocks. The downstream magnetic field is parameterized as

$p' = (\hat{\gamma}-1)e',$7

where $p' = (\hat{\gamma}-1)e',$8. Injected electrons satisfy

$p' = (\hat{\gamma}-1)e',$9

and the cooling timescale is

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$0

giving

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$1

The acceleration timescale is written as

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$2

which yields

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$3

The code assumes $\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$4 implicitly, so the accelerated-electron fraction is absorbed into $\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$5 (Nedora et al., 2024).

A major algorithmic feature is the numerical evolution of the downstream electron distribution. Rather than assuming a broken power law at every step, the code can solve

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$6

with

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$7

Adiabatic cooling is written as

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$8

and synchrotron cooling as

$\hat{\gamma} \approxeq \frac{4 + \Gamma^{-1}}{3}.$9

This kinetic equation is solved with a fully implicit, particle-number-preserving Chang–Cooper scheme on a logarithmic $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$0-grid from $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$1 to $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$2 (Nedora et al., 2024).

The synchrotron emissivity is evaluated numerically from the current $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$3. For $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$4, the code uses an analytic approximation to the angle-averaged synchrotron emissivity, while synchrotron self-absorption is computed from

$$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$5

The SSC module adopts a one-zone approximation, includes Klein–Nishina suppression, and treats only first-order inverse Compton scattering. Pair production is added through a $$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$6 source term in the electron kinetic equation and through attenuation of the escaping photon spectrum (Nedora et al., 2024).

The framework was subsequently extended, in kilonova-afterglow applications, to include both thermal and non-thermal electron populations. In that implementation, the total emissivity and absorption are

$$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$7

with the thermal part following a Maxwellian or Maxwell–Jüttner distribution. That extension showed that thermal electrons can dominate the early-time radio band and can generate double-peaked radio light curves at sufficiently high density, while the non-thermal power-law component dominates later (Nedora et al., 2022).

4. Jet structure, observables, and computational workflow

The code supports top-hat jets, Gaussian structured jets, and, according to the paper, arbitrary ejecta structures in principle. A typical workflow consists of specifying source and environment parameters, evolving blast-wave dynamics for each angular layer, computing microphysics and comoving spectra either analytically or with the numerical electron-evolution module, and then evaluating observables by equal-arrival-time-surface integration over the jet surface (Nedora et al., 2024).

Observer-frame radiation is obtained from the Lorentz-transformed emissivity and absorption coefficients,

$$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$8

with

$$E_{\rm tot} = \Gamma c^2 m + \Gamma_{\rm eff}E_{\rm int}',$$9

For a uniform slab, the emergent intensity is

$m$0

and the observer flux density is written as an angular integral over the jet surface. In the layered discretization, the code naturally handles on-axis and off-axis observers, Doppler boosting, equal-arrival-time surface effects, and counter-jet contributions (Nedora et al., 2024).

PyBlastAfterglow also computes two-dimensional sky maps. The visible shell is discretized into subrings and equal-solid-angle cells, projected onto the image plane, and post-processed in Python. The post-processing stage interpolates raw cell intensities by Delaunay triangulation, bins them on a regular $m$1-$m$2 grid, and derives image centroids and FWHM image sizes. This imaging pipeline was used in the paper’s GRB170817A application to compare structured-jet sky maps with VLBI centroid motion and apparent velocity constraints (Nedora et al., 2024).

The code’s output domain is correspondingly broad. It can generate blast-wave dynamical histories, downstream electron distributions, comoving synchrotron and SSC spectra, observer-frame spectra, multi-band light curves, sky maps, flux centroids, apparent proper motions, and image sizes. The physical emission range spans radio through optical and X-ray, up to gamma rays or TeV energies through SSC, and includes reverse-shock contributions at low frequencies and early times (Nedora et al., 2024).

5. Applications and use in inference

The paper presents several benchmark and science applications: top-hat forward-shock runs; top-hat forward-plus-reverse-shock runs; Gaussian structured forward-plus-reverse-shock jets; broadband modeling of GRB 190114C from soft X-ray to TeV with an SSC interpretation of the very-high-energy afterglow; and structured off-axis modeling of GRB170817A, including sky maps and centroid motion. For performance-oriented parameter studies, the paper reports a GRB 190114C grid of $m$3 runs using $m$4 CPU hours, while structured-jet full-physics GRB170817A-like runs with sky maps are stated to require about $m$5 CPU hour per run in the reported setup (Nedora et al., 2024).

The earlier kilonova-afterglow study broadened the astrophysical scope further. There, PyBlastAfterglow was used to propagate angle-dependent, velocity-stratified ejecta through both a uniform ISM and a GRB-modified circumburst medium with nonzero upstream velocity and pressure. That generalized blast-wave equation allowed the code to model a modified Mach number,

$m$6

and thus a shock/no-shock transition in a pre-cleared GRB wake. In that application, the code predicted early emission suppression at $m$7 days and only mild subsequent rebrightening when kilonova ejecta break through the GRB-processed environment (Nedora et al., 2022).

PyBlastAfterglow has also entered Bayesian inference through surrogate modeling. A 2025 study introduced fiesta, a Python package that trains machine-learning surrogates for GRB afterglow and kilonova models, including PyBlastAfterglow. In that work, a single PyBlastAfterglow calculation is described as potentially taking several minutes, which motivated surrogate emulation of full broadband spectral-flux surfaces $m$8 rather than fixed-filter light curves. For the Gaussian-jet PyBlastAfterglow surrogate, the training set size was $m$9, the test set size was $E_{\rm int}'$0, the preferred inference surrogate was a cVAE FluxModel, and the typical mismatch on held-out tests was usually $E_{\rm int}'$1 mag, although about $E_{\rm int}'$2 of test cases exceeded $E_{\rm int}'$3 mag somewhere along the light curve. That surrogate framework enabled the first Bayesian analyses of GRB afterglows with pyblastafterglow, including reanalyses of GRB170817A and GRB211211A (Koehn et al., 18 Jul 2025).

6. Validation, approximations, and relation to other frameworks

PyBlastAfterglow is validated in the 2024 framework paper against afterglowpy, jetsimpy, a simplified reference SSC code, and standard analytic expectations such as Blandford–McKee behavior, synchrotron spectra, and the synchrotron self-absorption frequency. The reported comparison is good before strong lateral spreading, with differences near spreading onset attributed to the fact that PyBlastAfterglow evolves independent one-dimensional blast waves with a spreading prescription, whereas jetsimpy uses a more self-consistent two-dimensional thin shell. The paper also introduces a faster comparison mode, denoted PyBlastAfterglow*, in which analytic synchrotron prescriptions are used without numerical electron evolution (Nedora et al., 2024).

The same paper is explicit about approximations. Its main physical simplifications are the thin-shell approximation, independent angular layers with prescription-based spreading rather than self-consistent lateral coupling, a simplified equation of state,

$E_{\rm int}'$4

a one-zone SSC treatment, a simplified reverse-shock upstream profile based on an exponential shell-density model and $E_{\rm int}'$5, no synchrotron-self-absorption heating of electrons, no diffusion or escape terms in the electron kinetic equation, and approximate late-time electron physics as $E_{\rm int}'$6. The paper quantifies energy conservation as typically within $E_{\rm int}'$7 during early evolution, worsening in the trans-relativistic stage to up to $E_{\rm int}'$8 for FS-only and up to $E_{\rm int}'$9 when the RS becomes relativistic (Nedora et al., 2024).

Within the broader afterglow-software landscape, PyBlastAfterglow occupies a distinct position. ScaleFit uses precomputed high-resolution two-dimensional relativistic hydrodynamic simulations, scaling relations, and lookup tables to generate light curves in about a millisecond, making it explicitly fitting-oriented rather than evolution-oriented (Ryan et al., 2013). DeepGlow is a neural emulator of BOXFIT that predicts single-frequency light curves on a fixed time grid, reaches $\Gamma_{\rm eff}$0 ms evaluation time, and delivers a $\Gamma_{\rm eff}$1 speedup over BOXFIT while retaining a few-percent average fidelity to the underlying RHD-based model (Boersma et al., 2022). jetsimpy adopts a reduced hydrodynamic model in which the blast wave is an infinitely thin two-dimensional surface that becomes effectively one-dimensional under axial symmetry, and it is explicitly designed for structured jets, flux centroid motion, and MCMC-oriented efficiency (Wang et al., 2024). VegasAfterglow, a later high-performance C++ framework with Python interfaces, extends the fast-afterglow paradigm toward arbitrary density profiles, arbitrary central-engine energy-injection histories, structured jets with arbitrary angular magnetization profiles, and analytic shock-jump solutions for arbitrary upstream magnetization $\Gamma_{\rm eff}$2 (Wang et al., 14 Jul 2025).

A complementary benchmark resource is the online library of afterglow light curves based on high-resolution two-dimensional RAM adaptive-mesh-refinement simulations with synchrotron radiation post-processing. That library provides synthetic light curves and broadband spectra over observer frequencies from radio to X-ray and observer times from hours to decades, and is explicitly intended for checking the accuracy of physical parameters derived from analytical model fits and for exploring observer-angle effects (Eerten et al., 2011). Against that landscape, the 2024 PyBlastAfterglow paper presents the code’s novelty not as a single dynamical prescription, but as the combination of forward-plus-reverse-shock thin-shell dynamics, numerical downstream electron evolution, synchrotron plus SSA plus SSC plus pair production plus EBL attenuation, structured-jet and off-axis equal-arrival-time imaging, and an open-source modular C++/Python implementation (Nedora et al., 2024).