jetsimpy: Fast Hydrodynamic GRB Afterglow Simulator

• jetsimpy is a reduced hydrodynamic code that models relativistic blast waves as thin, axisymmetric surfaces to compute broadband synchrotron emission in GRB afterglows.
• The 2025 upgrade introduces radial stratification, enabling self-consistent reverse shock evolution and improved modeling of refreshed shocks and kilonova afterglows.
• By conserving energy through coupled forward and reverse shock dynamics, jetsimpy achieves thousands-fold speedups over full hydrodynamic simulations for inference-heavy analyses.

jetsimpy is a reduced hydrodynamic code and fast hydrodynamic afterglow simulator for gamma-ray burst afterglows that models a relativistic blast wave as a thin axisymmetric surface and computes broadband synchrotron emission for arbitrary jet structures. The original release supported arbitrary angular distributions of energy and Lorentz factor, reproduced full relativistic hydrodynamics for structured jets with lateral spreading and off-axis viewing, and was designed as a computationally efficient alternative to full relativistic hydrodynamic simulations in inference-heavy settings. The 2025 upgrade extends the framework to arbitrary axisymmetric two-dimensional ejecta by adding a stratified radial profile, thereby generating a self-consistent reverse shock and evolving the coupled forward–reverse shock system under global energy conservation rather than pressure balance at the contact discontinuity. In that form, jetsimpy is intended for early afterglows with reverse shocks, structured jets viewed on- and off-axis, refreshed shocks, kilonova afterglows, and profiles extracted from general relativistic magnetohydrodynamic simulations (Wang et al., 2024, Wang et al., 21 Jul 2025).

1. Origins and intended domain

jetsimpy was introduced to address a methodological gap between full relativistic hydrodynamics and semianalytic afterglow models. Full relativistic hydrodynamic simulations are accurate but computationally expensive, which makes broad parameter sweeps and Markov Chain Monte Carlo inference difficult. Semianalytic models are fast, but their fidelity degrades in the transitional spreading regime where self-similar solutions break down. jetsimpy was constructed as a reduced hydrodynamic model that approximates the blast wave as an infinitely thin two-dimensional surface and, under axial symmetry, reduces the dynamics to one dimension in polar angle.

The original code was formulated for structured jets with arbitrary angular energy and Lorentz-factor profiles, and it evolved the blast wave through ultrarelativistic, transitional, and Newtonian phases in ISM, wind, or smoothly varying media. It was calibrated to the Blandford–McKee and Sedov–Taylor limits, included an explicit coasting phase, and computed broadband synchrotron light curves, flux centroid motion, and image sizes. The code was described as achieving speedups by factors of thousands relative to 2D relativistic hydrodynamics, making millions of model evaluations feasible in MCMC workflows. Although developed with gamma-ray burst afterglows in mind, the formulation was presented as applicable to any relativistic, axisymmetric blast wave or jet for which thin-shell behavior is a reasonable approximation (Wang et al., 2024).

The later upgrade was motivated by cases in which the ejecta shell thickness and internal velocity stratification cannot be neglected. In those problems, the ejecta no longer behave as a purely angularly structured thin shell; instead, radial stratification generates a reverse shock and couples the forward-shock and reverse-shock regions dynamically. This extension enlarges the target problem class to include optical flashes, GW170817-like off-axis afterglows, refreshed shocks from long-lived energy injection, and kilonova afterglows (Wang et al., 21 Jul 2025).

2. Thin-shell hydrodynamic formulation

The original jetsimpy formulation represents the shocked region as an infinitely thin surface at radius $R(t,\theta)$. The forward shock sweeps the external medium, the blast is described by radially integrated conservation laws, and axial symmetry reduces the dynamics to a one-dimensional finite-volume Eulerian problem in polar angle. Lateral transport is not imposed through a phenomenological sound-speed prescription; rather, it arises from the polar momentum equation and the angular flux terms of the integrated Euler system. This construction is intended to capture jet-core deceleration, bow-shock formation at sharp edges, and the gradual evolution toward spherical symmetry.

In the upgraded model, the thin blast surface is bounded by the forward shock and reverse shock and separated internally by a contact discontinuity. The governing evolution is written as

$$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$

with

$$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$

Here $E_b$, $H_b$, and $P_b$ are the blast energy, enthalpy, and pressure integrated over the thin shell, while $M_{\rm fs}$ and $M_{\rm rs}$ are the shocked masses in the forward- and reverse-shock regions. The closure relations are written in terms of the shocked-fluid Lorentz factor $\Gamma$, the relative Lorentz factor $\Gamma_{34}$ between shocked and unshocked ejecta, and the corresponding masses. The code advances the shock and discontinuity radii through jump-condition-defined velocities for the forward shock, contact discontinuity, and reverse shock.

A notable feature of the 2025 formulation is its treatment of lateral expansion. Although the full integrated pressure includes both forward- and reverse-shock contributions, the production runs stabilize lateral dynamics by using forward-shock pressure only,

$$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$0

while retaining the reverse-shock contribution in the blast energy $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$1. The paper states that this yields robust jet spreading consistent with multidimensional simulations once the blast is mildly relativistic; an option using the full pressure closure is also provided. This places energy conservation, rather than pressure balance at the contact discontinuity, at the center of the dynamical scheme (Wang et al., 2024, Wang et al., 21 Jul 2025).

3. Angular structure and radial stratification

In its original form, jetsimpy accepted arbitrary, tabulated angular distributions of energy per solid angle and initial Lorentz factor. A common parameterization used in the upgraded paper is the power-law core-plus-wing structure

$$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$2

$$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$3

where $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$4 is the core energy, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$5 the core Lorentz factor, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$6 the core half-opening angle, and $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$7 the wing slope. More generally, arbitrary angular profiles can be tabulated and ingested directly.

The 2025 extension introduces radial stratification through two axisymmetric functions defined at each polar angle: $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$8, the isotropic-equivalent power of material launched at engine time $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$9, and $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$0, the isotropic-equivalent energy distribution over ejecta four-velocity $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$1. Their coupling is specified by

$$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$2

which guarantees that later-launched layers have lower $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$3 and therefore avoid internal collisions, and by the lab–launch time conversion

$$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$4

These relations determine the velocity structure of the ejecta layer reaching the reverse shock at a given $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$5.

The ejecta comoving density at the reverse shock is expressed as

$$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$6

with $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$7 obtained from $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$8, $$U = (E_b,\beta_\theta H_b,M_{\rm fs},M_{\rm rs})^T, \qquad F = \frac{c}{R_{\rm fs}}(\beta_\theta H_b,\beta_\theta^2 H_b + P_b,\beta_\theta M_{\rm fs},\beta_\theta M_{\rm rs})^T.$$9, and the same radial stratification functions. In practical terms, these tabulated functions can represent top-hat shells, thin shells that spread to $E_b$0 with homogeneous velocity stratification, and broad distributions for refreshed shocks, kilonova afterglows, and jet cocoons. Because $E_b$1 and $E_b$2 are tabulated per $E_b$3, the code can also ingest axisymmetric two-dimensional profiles extracted from GRMHD simulations once they are converted into $E_b$4, $E_b$5, $E_b$6, and $E_b$7 in physical units (Wang et al., 21 Jul 2025).

4. Radiation, geometry, and observables

jetsimpy computes emission by integrating Doppler-boosted comoving synchrotron emissivity over equal-arrival-time surfaces. In the reverse-shock region the specific intensity is written as

$E_b$8

and the observed flux density as

$E_b$9

with the Doppler factor

$H_b$0

Forward-shock emission is computed analogously using forward-shock thermodynamics and thickness. In the original code, the equal-arrival-time integration also supported image-domain observables such as flux centroid motion and image sizes, which were used in the GRB 170817A application.

The radiation model is the standard synchrotron prescription with microphysical parameters $H_b$1, $H_b$2, and $H_b$3. The 2025 paper states explicitly that the ejecta magnetization $H_b$4 is not included in the present work, corresponding to cold unmagnetized ejecta with $H_b$5, although mild reverse-shock magnetization can be mimicked phenomenologically through large $H_b$6 in the shocked ejecta. Characteristic frequencies $H_b$7, peak flux, spectral segments, and synchrotron self-absorption are implemented as in the original jetsimpy framework. A version-dependent distinction is important: the original code reports support for a deep-Newtonian correction, whereas the kilonova-afterglow applications in the 2025 reverse-shock paper presently use relativistic synchrotron formulas and do not include deep-Newtonian corrections when $H_b$8. This suggests that late-time mildly relativistic and Newtonian interpretations should be read with that approximation in mind (Wang et al., 2024, Wang et al., 21 Jul 2025).

5. Reverse shocks, refreshed shocks, and scientific results

A central result of the 2025 upgrade is the reformulation of reverse-shock dynamics under global energy conservation. For a power-law ambient density $H_b$9, the standard shell-thickness parameter is

$P_b$0

with $P_b$1. Thick shells satisfy $P_b$2, while thin shells satisfy $P_b$3. jetsimpy does not impose an analytic reverse-shock crossing time; the crossing occurs naturally when the reverse shock reaches the ejecta tail in the simulation. After that point, the code introduces an imagined zero-density tail with finite velocity to let $P_b$4 decrease smoothly while conserving energy and mass, avoiding numerical pathologies that would otherwise arise from abrupt removal of reverse-shock terms.

The paper contrasts this energy-conserving treatment with analytic pressure-balance models, which impose $P_b$5 at the contact discontinuity and use

$P_b$6

According to the paper, that relation tends to overestimate $P_b$7 and the reverse-shock energy density, especially in thin shells. The quantitative differences are substantial: in thick shells, the reverse-shock peak flux is approximately $P_b$8 dex below analytic predictions, and in thin shells, jetsimpy predicts reverse-shock peaks at least $P_b$9 dex fainter than in analytical thin-shell treatments. The paper attributes the thin-shell discrepancy to three compounded effects: neglect of the inherent velocity gradient in a homogeneously expanding thin shell, overestimation of $M_{\rm fs}$0 by pressure balance, and neglect of deceleration during crossing. It further states that these results naturally explain the rarity of bright optical flashes from reverse shocks.

The same radial formalism generates hydrodynamic energy injection when slower ejecta catch up with the decelerating blast. For kilonova afterglows the paper adopts

$M_{\rm fs}$1

and reports that the resulting forward-shock-dominated afterglow reproduces the analytic temporal slopes

$M_{\rm fs}$2

for $M_{\rm fs}$3 in the relevant spectral regime. In the example with $M_{\rm fs}$4, $M_{\rm fs}$5, and ISM density $M_{\rm fs}$6, the reverse-shock radio emission can exceed the forward shock for high $M_{\rm fs}$7, for example $M_{\rm fs}$8, while the forward-shock dynamics converge to the forward-shock-only case once all energy is deposited.

The upgraded code was also applied to the optical flash of GRB 990123. Using a thick-shell reverse shock and fixed values $M_{\rm fs}$9, $M_{\rm rs}$0, $M_{\rm rs}$1, and $M_{\rm rs}$2, the paper reports best-fit values $M_{\rm rs}$3, $M_{\rm rs}$4, $M_{\rm rs}$5, $M_{\rm rs}$6, $M_{\rm rs}$7, $M_{\rm rs}$8, $M_{\rm rs}$9, $\Gamma$0, and $\Gamma$1. The fit is described as showing good agreement with the optical and X-ray data, while minor discrepancies in the earliest X-ray points are attributed to global $\Gamma$2 approximations. The inferred $\Gamma$3 close to unity is noted in the paper as suggestive of at least mild magnetization in the shocked ejecta (Wang et al., 21 Jul 2025).

6. Benchmarks, implementation, and limitations

The original paper benchmarked jetsimpy against the moving-mesh 2D relativistic hydrodynamics code JET, against afterglowpy for Gaussian jets, and against BoxFit for top-hat jets in both ISM and wind environments. In ISM top-hat tests with $\Gamma$4, the code reproduced core evolution and spreading, although it developed a more extended low-energy tail at large angles than the hydrodynamic simulation; that tail was described as radiatively subdominant. In wind tests with $\Gamma$5, agreement remained good but somewhat less accurate than in ISM. For BoxFit comparisons in ISM with $\Gamma$6, on- and off-axis radio and optical light curves typically agreed within a factor of approximately $\Gamma$7 across phases, with early-time off-axis differences attributed to bow shocks present in hydrodynamics but absent in BoxFit startup profiles. The same paper then applied jetsimpy to GRB 170817A. A light-curves-only fit yielded $\Gamma$8 deg and $\Gamma$9 deg, whereas including centroid offsets favored $\Gamma_{34}$0 deg and $\Gamma_{34}$1 deg, consistent with the interpretation that superluminal motion strongly constrains jet opening angle and viewing geometry.

From the software perspective, jetsimpy is written in C++ with a Python interface and is publicly available at https://github.com/haowang-astro/jetsimpy; the 2024 paper also provides a Zenodo DOI, 10.5281/zenodo.11078596. The upgraded version solves the thin-shell forward-plus-reverse-shock blast as a two-dimensional surface in $\Gamma_{34}$2, using characteristic information and a stable lateral expansion prescription. The paper states that the additional operations introduced by radial stratification—solving for $\Gamma_{34}$3, root-finding for the launch-time relation, and density conversion—are subdominant in cost relative to the pre-existing angular hydrodynamics. Inputs include tabulated $\Gamma_{34}$4, $\Gamma_{34}$5, $\Gamma_{34}$6, $\Gamma_{34}$7, an external density profile, forward- and reverse-shock microphysics, and the reverse-shock cooling exponent $\Gamma_{34}$8. Outputs include multi-band light curves $\Gamma_{34}$9 and intermediate hydrodynamic variables such as $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$00, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$01, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$02, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$03, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$04, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$05, and $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$06.

The limitations are explicit. The framework assumes axisymmetry, cold ejecta, negligible initial ejecta pressure, and in the 2025 reverse-shock study $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$07. Rayleigh–Taylor instabilities at the contact discontinuity and internal turbulence are not resolved. The thin-shell approximation is excellent in the relativistic regime but becomes approximate in deep Newtonian flow, even though the original calibration to the Sedov–Taylor limit preserves the correct global scalings. Strongly magnetized ejecta, time- or space-dependent microphysics, and detailed radial structure are not modeled explicitly. The 2025 paper also emphasizes parameter degeneracies among $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$08, ambient density, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$09, $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$10, and $$\frac{\partial U}{\partial t} + \frac{1}{\sin\theta}\frac{\partial(F\sin\theta)}{\partial\theta} + S_{\rm fs} + S_{\rm rs} = 0,$$11, particularly when spectral breaks are not well constrained, and recommends multi-band, multi-epoch fitting, informative priors, and joint forward- plus reverse-shock modeling when interpreting early afterglow data. Taken together, these caveats define jetsimpy as a physically constrained but deliberately reduced framework whose principal value lies in fast, energy-conserving modeling of structured relativistic outflows over a broad astrophysical parameter space (Wang et al., 2024, Wang et al., 21 Jul 2025).