#  jetsimpy: A Highly Efficient Hydrodynamic Code for Gamma-ray Burst
Afterglow

Hao Wang Department of Physics and Astronomy, Purdue University, 525
Northwestern Avenue, West Lafayette, IN 47907, USA  Ranadeep G. Dastidar
Department of Physics and Astronomy, Purdue University, 525 Northwestern
Avenue, West Lafayette, IN 47907, USA  Dimitrios Giannios Department of
Physics and Astronomy, Purdue University, 525 Northwestern Avenue, West
Lafayette, IN 47907, USA  Paul C. Duffell  Department of Physics and
Astronomy, Purdue University, 525 Northwestern Avenue, West Lafayette, IN
47907, USA

(Received xxx; Revised xxx; Accepted xxx)

###### Abstract

Gamma-ray Burst afterglows are emissions from ultra-relativistic blastwaves
produced by a narrow jet interacting with surrounding matter. Since the first
multi-messenger observation of a neutron star merger, hydrodynamic modeling of
GRB afterglows for structured jets with smoothly varying angular energy
distributions has gained increased interest. While the evolution of a jet is
well described by self-similar solutions in both ultra-relativistic and
Newtonian limits, modeling the transitional phase remains challenging. This is
due to the non-linear spreading of a narrow jet to a spherical configuration
and the breakdown of self-similar solutions. Analytical models are limited in
capturing these non-linear effects, while relativistic hydrodynamic
simulations are computationally expensive which restricts the exploration of
various initial conditions. In this work, we introduce a reduced hydrodynamic
model that approximates the blastwave as an infinitely thin two-dimensional
surface. Further assuming axial-symmetry, this model simplifies the simulation
to one dimension and drastically reduces the computational costs. We have
compared our modeling to relativistic hydrodynamic simulations, semi-analytic
methods, and applied it to fit the light curve and flux centroid motion of
GW170817. These comparisons demonstrate a high level of agreement and validate
our approach. We have developed this method into a numerical tool, jetsimpy,
which models the synchrotron GRB afterglow emission from a blastwave with
arbitrary angular energy and Lorentz factor distribution. Although the code is
built for GRB afterglow in mind, it is applicable to any relativistic jets.
This tool is particularly useful in Markov Chain Monte Carlo studies and is
provided to the community.

Gamma-ray bursts, Hydrodynamics, Jets

††journal: ApJS††software: astropy (Astropy Collaboration et al., 2013, 2018,
2022), Matplotlib (Hunter, 2007), Bilby (Ashton et al., 2019), pymultinest,
(Buchner et al., 2014)

xxx

##  1 Introduction

Gamma-ray Bursts (GRBs) are the most energetic catastrophic events in the
Universe. Their tremendous explosive power provides a platform for scientific
study of fundamental physical processes in extreme physical environments. One
of the intriguing aspects of GRBs is their long-term, multi-waveband afterglow
following the prompt emission. GRB afterglows are radiation from the ultra-
relativistic blastwaves produced by the interaction of a narrow jet with the
surrounding external medium. They are the promising sites to study the
hydrodynamical interactions and microphysical processes associated with some
of the most extreme macroscopic bulk motions in the Universe. Recently, it has
also been demonstrated that GRB afterglows could serve as probes to study the
Universe’s expansion rate (see Bulla et al. 2022 for a review), thereby
helping to unravel the nature of the Hubble tension (see Di Valentino et al.
2021 for a review).

To model the evolution of a narrow GRB jet, the angular energy profile is
often approximated by a so-called “top-hat” model, where the energy is
uniformly distributed within a narrow core, and sharply drops to zero at the
edge. This model has been proven highly successful in accounting for the GRB
phenomenology when the jet is observed from an on-axis direction, namely, when
the observer’s line of sight lies within the jet core (see Kumar & Zhang 2015
for a review). This is a natural consequence of the selection bias, since at
cosmological distances, it is only possible to observe a GRB within the jet
core, where the prompt emission is bright enough for detection (e.g.,
Beniamini & Nakar 2019). In this situation, the less energetic tail outside
the core becomes subdominant in observations, making the study of a smoothly
varying angular energy profile less compelling.

However, the revolutionary discovery of the first multimessenger observation
of a binary neutron star merger (Abbott et al. 2017a, b, c; see Margutti &
Chornock 2021 for a review), the famous event GW170817, has revealed the
importance of a more sophisticated jet model: the so-called “structured jet”.
A structured jet is characterized by the previously mentioned smoothly varying
angular energy profile, in which a less energetic tail extends far beyond the
core. This model is motivated by the argument that GW170817 was very likely
observed from an off-axis direction (Troja et al., 2017; Margutti et al.,
2017; Haggard et al., 2017; Alexander et al., 2018; Granot et al., 2018;
Lazzati et al., 2018; Margutti et al., 2018; Mooley et al., 2018; Troja et
al., 2018; Fong et al., 2019; Hajela et al., 2019; Lamb et al., 2019; Wu &
MacFadyen, 2019; Troja et al., 2020). In this case, the radiation from the
energetic core is subdominant at early times, as it is beamed away due to
ultra-relativistic motion. The less energetic tail close to the line of sight,
on the other hand, becomes important in the modeling of early emission. As the
jet decelerates, the line of sight eventually enters the beaming cone, and the
energetic core gradually becomes visible. This process leads to a slowly
increasing luminosity in the early phase of the afterglow, which effectively
explains the rising light curves observed in the GW170817. In contrast to a
structured jet, a top-hat jet model predicts a sharp increase in the early
phase for a misaligned jet, which is at odds with observations from this event
(e.g., Nakar & Piran 2018). Additionally, the apparent superluminal motion
observed in this event (Mooley et al., 2018; Ghirlanda et al., 2019; Mooley et
al., 2022) further supports an off-axis observation of a structured jet. It
has been suggested that future multimessenger observations of neutron star
mergers will predominantly be off-axis (e.g., Troja et al. 2017; Saleem et al.
2018; Kathirgamaraju et al. 2018; Duque et al. 2019; Beniamini et al. 2019;
Wang & Giannios 2021), indicating that the GW170817 scenario is likely to be
common.

Various works have suggested that the rising slope of off-axis GRB afterglows
helps to constrain the observing angle between the jets axis and the line of
sight (e.g., Finstad et al. 2018; Gill & Granot 2018; Wu & MacFadyen 2019;
Ryan et al. 2020). Such a constraint becomes particularly important when
considering the combined view with gravitational wave observations. The
observing angle, which is also assumed to be the inclination angle of the pre-
merger neutron star binary, is degenerate with the luminosity distance
inferred from gravitational wave data. A better constrained observing angle
helps break the degeneracy and provides a better distance measurement, making
it a powerful tool in the standard siren study of the Hubble constant (e.g.,
Guidorzi et al. 2017; Hotokezaka et al. 2019; Wang & Giannios 2021; Wang et
al. 2023; Gianfagna et al. 2023). This is particularly helpful in the current
context, where there is growing evidence of tension between different
measurements of the Hubble constant (see Di Valentino et al. 2021 for a
review). So far, all studies utilizing this method have been model-dependent,
as they assume a specific type of jet profile. The jet profile is either
approximated by an analytical expression, such as Gaussian jets (e.g., Gill &
Granot 2018; Ryan et al. 2020), power-law jets (e.g., Mooley et al. 2022), and
boosted fireball jets (e.g., Wu & MacFadyen 2018, 2019), or adopted by the
output of another numerical simulation (e.g., Wang & Giannios 2021). It is
unclear how the assumed jet profile affects the inferred observing angle.
Recently, Gianfagna et al. (2023) found a moderate difference in the observing
angle inferred from Gaussian and power-law jets, which motivates further
investigation into this effect. Therefore, it is necessary to develop a
reliable afterglow model that allows for the flexibility of various jet
structures.

Modeling the GRB afterglow with an arbitrary jet structure presents
challenges. Important features of off-axis afterglow can be well modeled by
analytic methods (Beniamini et al., 2020, 2022), but the accurate modeling of
late-time spreading effect (Rhoads, 1999; Sari et al., 1999) remains a
challenge. This complexity stems from the fact that spreading, driven by
internal pressure, is a non-linear hydrodynamic process. Consequently,
analytical methods struggle to accurately describe this phase. Several
attempts, including approximations based on sound speed expansion (e.g.,
Rhoads 1999; Ryan et al. 2020; Nedora et al. 2023b) or curvature effects
(e.g., Granot & Piran 2012), have been made to model the spreading of a top-
hat jet. Some of them are even calibrated to numerical simulations (e.g.,
Duffell & Laskar 2018). However, it remains uncertain to what extend these
semi-analytical methods are applicable to structured jet models. To date, the
most accurate and reliable approach for modeling in these scenarios is
relativistic numerical hydrodynamic simulations.

Fully numerical hydrodynamic approches present their own challenges. In the
ultra-relativistic limit the swept-up material are compressed in an extremely
thin shell, typically with a width of approximately $R/\gamma^{2}$ where $R$
represents the shock radius and $\gamma$ the Lorentz factor. Numerically
resolving this thin structure is challenging, especially when the jet’s
Lorentz factor is exceptionally high. This substantial computational cost
restricts the use of Markov Chain Monte Carlo (MCMC) for parameter estimation
in models with varying jet structures. Consequently, there is an pressing need
to develop a model that balances the computational efficiency of semi-analytic
methods with the accuracy of full numerical hydrodynamics.

In this work, we propose an alternative hydrodynamic approach that
significantly reduces the computational cost of a full hydrodynamic
simulation. Our method bears resemblance to that of Kumar & Granot (2003) and
Lu et al. (2020) where the shocked region is approximated by an infinitely
thin two dimensional surface. By additionally assuming the axial-symmetry of
the jet, the simulation is effectively simplified into a one-dimensional
representation. While our method maintains the accuracy of numerical
hydrodynamics, it attains a computational efficiency that is comparable to
semi-analytic methods. We have developed the methods described in this work
into a numerical tool jetsimpy. The code simulates the evolution of the
blastwave and generate synthetic afterglow observables for a jet with
tabulated angular energy and Lorentz factor profiles, traveling within either
an interstellar medium (ISM) environment or stellar wind environment.

This paper is organized as follows. In §2 we describe our hydrodynamic model
of a structured GRB jet. In §3 we present our numerical results and the
comparisons to full numerical relativistic hydrodynamics. In §4 we briefly
introduce our modeling of the afterglow observables. In §5 we verify and
compare our methods to existing afterglow modeling tools. In Section 6, we
apply our model to the case of GW170817. Finally, in §7, we discuss and
summarize our methods.

##  2 Hydrodynamic Model

###  2.1 Basic assumptions

In this work, we assume the jet is axial-symmetric and the external medium
follows a power-law density distribution: $\rho_{0}\propto r^{-k}$. When the
ultra-relativistic jet moves through a cold external medium, it gives rise to
a forward-reverse shock structure. The region between these two shocks is
commonly referred to as the “blast”. The forward shock sweeps up the external
material, increasing the mass in the blast, while the reverse shock propagates
through the ejecta in the opposite direction. In this study, we focus on the
forward shock and assume that the ejecta is fully shocked and has sufficiently
cooled down. That is, we consider the blast region to consist of a hot, swept-
up material layer and a cold ejecta layer. We further ignore the impact of
radiative loss to the hydrodynamics. Throughout this section, variables
associated with the ejecta are labeled as “ej,” and those related to the
swept-up material are labeled as “sw”.

To describe the thermodynamic properties of the fluid, we consider two
reference frames: the comoving frame, which moves along with the fluid
elements, and the burster frame, which remains static relative to the burst
origin. All physical quantities measured in the comoving frame are indicated
with a prime (’). Quantities in the burster frame, on the other hand, are
denoted without the prime. The velocity $\beta$ and Lorentz factor $\gamma$
are always measured in the burster frame.

The thermodynamic properties of the fluid elements just behind the forward
shock can be determined by the shock jump condition (Uhm, 2011):

$\displaystyle\rho^{\prime}_{\rm sw}$ $\displaystyle=4\gamma\rho_{0}$ (1)
$\displaystyle e^{\prime}_{\rm sw}$
$\displaystyle=4\gamma(\gamma-1)\rho_{0}c^{2}$ (2) $\displaystyle
p^{\prime}_{\rm sw}$
$\displaystyle=\frac{4}{3}\beta^{2}\gamma^{2}\rho_{0}c^{2}$ (3)

In the last expressions, $\rho^{\prime}_{\rm sw}$, $e^{\prime}_{\rm sw}$, and
$p^{\prime}_{\rm sw}$ represent the mass density, energy density (rest mass
excluded), and pressure of the shocked material in the comoving frame,
respectively. The $\rho_{0}$ denotes the mass density of the ambient matter.
These relations are derived from a semi-analytic trans-relativistic equation
of state, which smoothly interpolates the adiabatic index between the ultra-
relativistic regime and Newtonian regime (Uhm, 2011).

In our model, the blast shell is approximated as an infinitely thin, two-
dimensional surface. Within this approximation, the thermodynamic properties
of both the hot, swept-up material and the cold ejecta are effectively
represented by Dirac delta functions. We further assume that the radial
average of the swept-up material is proportional to the values immediately
behind the forward shock. This allows us to directly determine the scale of
the delta function based on the shock jump conditions. Our approximations
remain valid due to the following considerations:

i) The blast shell is effectively thin. We can estimate the typical width of
the blast shell, $\Delta R$, by noting that the total mass behind the forward
shock equals the total swept-up mass at radius $R$. For example, in a constant
external medium, it indicates that $4\pi R^{2}\Delta R\gamma\rho^{\prime}_{\rm
sw}\approx 4\pi R^{3}\rho_{0}/3$, yielding an estimate of $\Delta R\sim
R/12\gamma^{2}$ in the burster frame. This estimation indicates that a thin
shell is an excellent approximation in the ultra-relativistic limit, hence
motivates the delta function approximation. While this assumption might not
hold as well in the Newtonian regime, the use of a trans-relativistic equation
of state ensures that the radially integrated mass and energy scale correctly
with the well-known Sedov-Taylor solution. As a result, our model continues to
provide order-of-magnitude estimations in this limit.

ii) The radial profile is self-similar. It is well known that an isotropically
expanding blastwave follows the Blandford-McKee self-similar solution in the
ultra-relativistic limit (Blandford & McKee, 1976), and the Sedov-Taylor self-
similar solution in the Newtonian limit (e.g., Sedov 1959). The self-similar
nature of these radial profiles implies that the radial average of
thermodynamic quantities is proportional to the values immediately behind the
forward shock, with this ratio remaining approximately constant over time.
This fact motivates us to approximate the radial average of the blast by the
values of shock jump conditions with a calibration coefficient. The
calibration coefficient helps to ensure that the radial integrated value
aligns with the exact self-similar solutions.

Based on these considerations, we first define the following radial integrated
values: the mass of swept-up matter and ejecta per solid angle

$\displaystyle M_{\rm sw}$ $\displaystyle=\int\gamma\rho^{\prime}_{\rm
sw}r^{2}dr,$ (4) $\displaystyle M_{\rm ej}$
$\displaystyle=\int\gamma\rho^{\prime}_{\rm ej}r^{2}dr,$ (5)

where $\rho^{\prime}_{\rm ej}$ is the density of the ejecta. Within the
infinitely thin shell approximation, the densities are assumed to be delta
functions. In this limit, the above integral implies the following relations

$\displaystyle\rho^{\prime}_{\rm sw}$ $\displaystyle=M_{\rm
sw}R^{-2}\gamma^{-1}\delta(r-R)$ (6) $\displaystyle\rho^{\prime}_{\rm ej}$
$\displaystyle=M_{\rm ej}R^{-2}\gamma^{-1}\delta(r-R)$ (7)

where $\delta(x)$ is the Dirac delta function. Similarly, we can approximate
the energy density and pressure using delta functions. According to the shock
jump condition, the density just behind the shock is $\rho^{\prime}_{\rm
sw}=4\gamma\rho_{0}$. We can construct delta functions for energy and pressure
by substituting all terms of “$4\gamma\rho_{0}$” in the shock jump conditions
with the above delta function. Given that the radial profile of the blast is
self-similar, the scale is proportional to this value with a calibration
coefficient. These considerations lead to the following approximations:

$\displaystyle\bar{e}^{\prime}_{\rm sw}$ $\displaystyle=s(1-\gamma^{-1})M_{\rm
sw}R^{-2}c^{2}\delta(r-R),$ (8) $\displaystyle\bar{p}^{\prime}_{\rm sw}$
$\displaystyle=\frac{1}{3}s\beta^{2}M_{\rm sw}R^{-2}c^{2}\delta(r-R).$ (9)

Our calibration coefficient “$s$” will be determined by matching our solutions
to analytic and numerical results later. These relations will serve as the
building blocks in the derivation of evolution equations.

###  2.2 Evolution equations

To derive the evolution equations, we start from the relativistic Euler
equations,

$\displaystyle T^{\mu\nu}_{;\mu}=0,$ (10) $\displaystyle J^{\mu}_{;\mu}=0,$
(11)

where $T^{\mu\nu}$ is the energy momentum tensor, and $J^{\mu}$ is the mass
current. These relations represents the conservation of energy, momentum, and
mass. In this work we assume the blast is perfect fluid. In a spherical
coordinate with metric tensor $g_{\mu\nu}={\rm
diag}(-1,1,r^{2},r^{2}\sin^{2}\theta)$, the energy momentum tensor and mass
current take the following form

$\displaystyle T^{\mu\nu}=h^{\prime}_{\rm
b}u^{\mu}u^{\nu}+g^{\mu\nu}\bar{p}^{\prime}_{\rm sw},$ (12) $\displaystyle
J_{\rm sw}^{\mu}=\rho^{\prime}_{\rm sw}u^{\mu},$ (13) $\displaystyle J_{\rm
ej}^{\mu}=\rho^{\prime}_{\rm ej}u^{\mu},$ (14)

where $h^{\prime}_{\rm b}$ is the enthalpy of the blast defined by

$h^{\prime}_{\rm b}=\bar{e}^{\prime}_{\rm sw}+\bar{p}^{\prime}_{\rm
sw}+\rho^{\prime}_{\rm sw}c^{2}+\rho^{\prime}_{\rm ej}c^{2}.$ (15)

Assuming axial-symmetry, the 4-velocity $u^{\mu}$ is

$u^{\mu}=(\gamma,\beta_{\rm r}\gamma,\frac{\beta_{\theta}\gamma}{r},0)$ (16)

where $\beta_{\rm r}$ and $\beta_{\theta}$ are the velocity components along
radial and polar directions respectively. In the definitions above, the swept-
up mass and ejecta in the mass current are defined independently, while the
energy momentum tensor is defined by the sum. This is because energy
conservation applies to the entire blast, whereas the mass conservation holds
for each component individually. By evolving the total energy momentum tensor,
we ensure the conservation of total energy, while avoiding the details of the
interaction between ejecta and swetp-up material.

In the spherical coordinate, the Eulerian equations can be explicitly
expressed by the following vector form

$\frac{1}{c}\frac{\partial{\bf A_{\rm t}}}{\partial
t}+\frac{1}{r^{2}}\frac{\partial(r^{2}{\bf A_{\rm r}})}{\partial
r}+\frac{1}{r\sin\theta}\frac{\partial({\bf
A_{\theta}}\sin\theta)}{\partial\theta}+{\bf s}=0,$ (17)

where the vector terms are defined by

$\displaystyle{\bf A_{\rm t}}=\begin{pmatrix}h^{\prime}_{\rm
b}\gamma^{2}-p^{\prime}_{\rm sw}\\\ h^{\prime}_{\rm
b}\gamma^{2}\beta_{\theta}\\\ \rho^{\prime}_{\rm sw}\gamma\\\
\rho^{\prime}_{\rm ej}\gamma\end{pmatrix},\ {\bf A_{\rm
r}}=\begin{pmatrix}h^{\prime}_{\rm b}\gamma^{2}\beta_{\rm r}\\\
h^{\prime}_{\rm b}\gamma^{2}\beta_{\theta}\beta_{\rm r}\\\ \rho^{\prime}_{\rm
sw}\gamma\beta_{\rm r}\\\ \rho^{\prime}_{\rm ej}\gamma\beta_{\rm
r}\end{pmatrix},$ $\displaystyle{\bf
A_{\theta}}=\begin{pmatrix}h^{\prime}_{\rm b}\gamma^{2}\beta_{\theta}\\\
h^{\prime}_{\rm b}\gamma^{2}\beta_{\theta}^{2}+p^{\prime}_{\rm sw}\\\
\rho^{\prime}_{\rm sw}\gamma\beta_{\theta}\\\ \rho^{\prime}_{\rm
ej}\gamma\beta_{\theta}\end{pmatrix},\ {\bf s}=\frac{1}{r}\begin{pmatrix}0\\\
h^{\prime}_{\rm b}\gamma^{2}\beta_{\theta}\beta_{\rm
r}-\frac{\cos\theta}{\sin\theta}p^{\prime}_{\rm sw}\\\ 0\\\ 0\end{pmatrix}.$

In these equations, the r-component momentum equation is not included. The
explanation will follow after we derive the radial integrated equations.

Now, we can apply the delta function approximation of the thermodynamic
quantities to eq. 17. After this, we multiply this equation by $r^{2}$ and
integrate it over the blast shell radius. For the sake of convenience, we
define the following integrated variables

$\displaystyle E_{\rm b}$
$\displaystyle=\frac{1}{c^{2}}\int_{0}^{R}(h^{\prime}_{\rm
b}\gamma^{2}-\bar{p}^{\prime}_{\rm sw})r^{2}dr$
$\displaystyle=s(1+\frac{1}{3}\beta^{4})\gamma^{2}M_{\rm sw}+(1-s)\gamma
M_{\rm sw}+\gamma M_{\rm ej},$ (18) $\displaystyle P_{\rm sw}$
$\displaystyle=\frac{1}{c^{2}}\int_{0}^{R}\bar{p}^{\prime}_{\rm
sw}r^{2}dr=\frac{1}{3}s\beta^{2}M_{\rm sw},$ (19) $\displaystyle H_{\rm b}$
$\displaystyle=\frac{1}{c^{2}}\int_{0}^{R}h^{\prime}_{\rm
b}\gamma^{2}r^{2}dr=E_{\rm b}+P_{\rm sw}.$ (20)

These variables represent the radial integrated energy, pressure, and
enthalpy, each normalized to the unit of mass. Utilizing these variables, we
can express the radial integrated Eulerian equations as follows:

$\frac{\partial{\bf U}}{\partial t}+\frac{1}{\sin\theta}\frac{\partial({\bf
F}\sin\theta)}{\partial\theta}+{\bf S}=0$ (21)

where the conserved variables ${\bf U}$, flux vector ${\bf F}$, and source
terms $\bf S$ are defined by

$\displaystyle{\bf U}=\begin{pmatrix}E_{\rm b}\\\ \beta_{\theta}H_{\rm b}\\\
M_{\rm sw}\\\ M_{\rm ej}\end{pmatrix},$ (22) $\displaystyle{\bf
F}=\frac{c}{R}\begin{pmatrix}\beta_{\theta}H_{\rm b}\\\
\beta_{\theta}^{2}H_{\rm b}+P_{\rm sw}\\\ \beta_{\theta}M_{\rm sw}\\\
\beta_{\theta}M_{\rm ej}\end{pmatrix},$ (23) $\displaystyle{\bf
S}=\frac{c}{R}\begin{pmatrix}-\frac{1}{c}\frac{\partial R}{\partial
t}\rho_{0}R^{3}\\\ \beta_{\theta}\beta_{\rm r}H_{\rm
b}-\frac{\cos\theta}{\sin\theta}P_{\rm sw}\\\ -\frac{1}{c}\frac{\partial
R}{\partial t}\rho_{0}R^{3}\\\ 0\end{pmatrix}.$ (24)

In these equations, additional source terms emerge in both the energy and mass
equations. These terms result from the movement of the integral boundary,
specifically the forward shock. The forward motion of the outer boundary leads
to the accumulation of the external medium, leading to the source terms in
these equations.

The most significant difference between our method and similar approaches in
previous works (e.g., Kumar & Granot 2003) is that the r-component equation,
specifically the radial momentum equation is not included. In full numerical
hydrodynamics, all momentum equations are required to solve for the velocity
of a fluid element in internal steps. However, in our approach, the total
velocity can be directly determined by eq. 18, given a complete set of $\bf
U$. This is because the shock jump condition is directly incorporated into the
derivation of the evolution equations, providing additional information. To
uniquely determine the velocity, we have chosen to exclude the radial momentum
equation due to the following reasons:

i) By reducing one partial differential equation, we can enhance computational
performance and simplify the eigenstructure of the Jacobian matrix (see
Appendix A).

ii) The integrated radial momentum equation averages the velocity across the
blast, which does not directly correspond to the shock jump condition. Since
most radiation occurs near the shock, the Lorentz factor just behind the shock
is more important than the radial average. Our method directly evolves these
post-shock values, offering greater flexibility in modeling radiation from the
blastwave.

iii) Employing Eq. 18 allows us to calibrate our method with self-similar
solutions through a constant coefficient. This calibration ensures that the
Lorentz factor at any given radius is consistent with exact self-similar
solutions. Such a calibration would not be possible if we substituted Eq. 18
with the radial momentum equation.

To complete the set of equations, we also require the evolution of the blast
shell radius $R(t,\theta)$. This can be derived by considering the speed of a
fluid element, given by $dR/dt=\beta c$. It’s important to note that in the
infinitely thin shell approximation, the velocity ($\beta$) difference between
the blast and the forward shock is essentially disregarded. While this is a
reasonable approximation in the ultra-relativistic limit, it becomes less
valid in the Newtonian limit. To address this discrepancy, we substitute the
velocity of the blast with the velocity of the forward shock:

$\beta_{\rm f}=\frac{4\beta\gamma^{2}}{4\gamma^{2}-1}$ (25)

This can be better understood by defining the radius of the blast as the
position of the forward shock. Now we can derive the radius evolution equation
by changing from a total derivative to a local derivative:

$\frac{\partial R}{\partial t}=\beta_{\rm f}c-\frac{\partial
R}{\partial\theta}\frac{\beta_{\theta}c}{R}$ (26)

In the numerical tests, we find that replacing $\beta$ by $\beta_{\rm f}$ in
this equation improves the accuracy of the solutions.

Now, the evolution of the blastwave can be solved by coupling eq.21 with
eq.26, provided with a set of appropriate initial conditions. These initial
conditions could be any physically reasonable angular distribution of initial
energy $E_{\rm b,0}$ and Lorentz factor $\gamma_{0}$. The mass of the ejecta
can then be estimated by $M_{\rm ej,0}=E_{\rm b,0}/\gamma_{0}$.

Similar to numerical hydrodynamics, these equations are non-linear and require
a well-designed approach to solve and maintain numerical stability. The
details of the numerical scheme are explained in Appendix A.

We notice that our hydrodynamic model introduces minor inconsistencies.
Specifically, the direction predicted for a fluid element by our model does
not always align with the shock normal direction. This inconsistency
originates from the infinitely thin shell approximation, which is insufficient
to reproduce all predictions of a full hydrodynamic model. On the other hand,
assuming that fluid elements move along the shock normal direction in the
context of infinitely thin surface approximation implies that the spreading of
a blastwave is completely determined by the geometry of the shock front,
instead of being governed by hydrodynamic equations. By accepting this minor
inconsistency, our approach allows us to derive the evolution equations from
first principles (i.e., the Euler equations), which should provide better
accuracy than a purely geometric approach.

###  2.3 Calibration with self-similar solutions

The calibration coefficient $s$ appearing in the equations is determined based
on comparisons with exact self-similar solutions. It is well-known that a
blastwave follows the Blandford-McKee self-similar solution in the ultra-
relativistic limit, and follows the Sedov-Taylor solution in the Newtonian
limit. Calibration can be achieved by aligning eq. 18 with the radial integral
of these self-similar solutions. It is also possible to calibrate the
coefficient in both limits and smoothly vary it in the transitioning phase.

To calibrate in the relativistic limit, we first revisit the Blandford-McKee
solutions (Blandford & McKee, 1976). In a power-law external density medium
described by $\rho_{0}(r)=Ar^{-k}$, the pressure and Lorentz factor profiles
are as follows:

$\displaystyle p^{\prime}_{\rm BM}(r)=\frac{2}{3}\rho_{0}c^{2}\Gamma_{\rm
f}^{2}f(\chi),$ (27) $\displaystyle\gamma_{\rm
BM}(r)^{2}=\frac{1}{2}\Gamma_{\rm f}^{2}g(\chi),$ (28)

where the self-similar variable and functions are defined by

$\displaystyle f(\chi)=\chi^{-(17-4k)/(12-3k)},$ (29) $\displaystyle
g(\chi)=\chi^{-1},$ (30) $\displaystyle\chi=1+2(4-k)(1-\frac{r}{R})\Gamma_{\rm
f}^{2}.$ (31)

In these equations $\Gamma_{\rm f}=\sqrt{2}\gamma_{\rm BM}(R)$ represents the
Lorentz factor of the forward shock111$\Gamma_{\rm f}$ does not precisely
match the aforementioned $\beta_{\rm f}$, as the equation of state used in the
Blandford-McKee solution is not trans-relativistic.. Based on the provided
profile, the energy per solid angle inside the forward shock can be integrated
as follows:

$E_{\rm b,BM}=\int_{0}^{R}4p^{\prime}_{\rm BM}\gamma_{\rm BM}^{2}r^{2}dr.$
(32)

The mass per solid angle can be simply calculated by integrating the external
density profile within the radius $R$:

$M_{\rm BM}=\frac{A}{3-k}R^{3-k}.$ (33)

To compare these results to eq. 18, we first define the rest mass and ejecta
mass excluded energy

$E_{\rm sw}=E_{\rm b}-\gamma M_{\rm ej}-M_{\rm sw}.$ (34)

In the ultra-relativistic limit, this value approaches $E_{\rm
sw}=\frac{4}{3}s\gamma^{2}M_{\rm sw}$. In the same limit, eq. 32 approaches
$E_{\rm b,BM}\propto\gamma^{2}M_{\rm BM}$. To align the two values, we define
the calibration coefficient in this limit by

$s_{\rm BM}(k)=\lim_{\gamma\to\infty}\frac{3E_{\rm b,BM}}{4\gamma^{2}M_{\rm
BM}}.$ (35)

We find that this coefficient depends on $k$ only.

We also calibrate our model with the Sedov-Taylor solution in a similar
manner. Since the exact solution is non-trivial, we numerically calculate the
radial integrated energy $E_{\rm ST}$ using the publicly available tool
sedov.tbz (Kamm & Timmes,
2007)222https://cococubed.com/research_pages/sedov.shtml. The mass integral
remains the same as above, i.e., $M_{\rm ST}=M_{\rm BM}$. In the Newtonian
limit, our model approaches $E_{\rm sw}=\frac{1}{2}(1+s)\beta^{2}M_{\rm sw}$,
while the Sedov solution approaches $E_{\rm b,ST}\propto\beta^{2}M_{\rm ST}$.
As a result, the calibration coefficient in this limit can be defined by

$s_{\rm ST}(k)=\lim_{\beta\to 0}\frac{2E_{\rm b,ST}}{\beta^{2}M_{\rm ST}}-1.$
(36)

We find the coefficient also depends on $k$ only. Both calibration
coefficients are non-trivial functions of $k$, and they are presented in Fig.
1.

Figure 1: The coefficients to calibrate our method with Blandford-McKee and
Sedov-Taylor self-similar solutions.

It can be seen that the two coefficients are in general different. To smoothly
transition between the two limits, we propose the following approximation for
the total calibration coefficient:

$s(k)=\frac{s_{\rm ST}+2s_{\rm
BM}\beta^{2}\gamma^{2}}{1+2\beta^{2}\gamma^{2}}.$ (37)

The factor 2 is chosen to better agree with hydrodynamics results in
respective limits.

We can further extend this result to apply in a more sophisticated external
medium by defining a “local” power index for the radial profile:
$k=-d\log\rho_{0}/d\log r$. We assume that the evolution of the calibration
coefficient in an environment with varying $k$ is quasi-static. This approach
is particularly useful in scenarios where the external medium transitions
smoothly from a wind-like to an ISM-like environment.

We caution here that the calibration in the Newtonian limit is tentative,
since the thin shell approximation will eventually breakdown at very low
speed. In this regime, the extensive radial profile leads to hydrodynamic and
radiation processes that differ significantly from those in the relativistic
regime. Therefore, our modeling in the Newtonian limit should be regarded as
providing only an order of magnitude estimation.

###  2.4 jetsimpy

The hydrodynamic model introduced above, along with the modelling of the
synchrotron emission from the electrons accelerated at the forward shock to be
discussed in subsequent sections, have been constructed into a numerical tool,
jetsimpy, and is provided to the community333https://github.com/haowang-
astro/jetsimpy. The code is written in C++ and features a flexible and user-
friendly Python interface. It is designed for arbitrary initial energy and
Lorentz factor distribution as initial conditions, and can evolve the jet in
ISM-like, wind-like, or mixed density profiles. jetsimpy drastically reduces
the computational cost by a factor of thousands compared to numerical
simulation. This efficiency is particularly beneficial for computing millions
of cases with varying angular profiles, such as in a MCMC study. Furthermore,
as we will demonstrate in the following sections, our tool maintains a level
of accuracy comparable to that of full numerical simulations, even in
situations where non-linear hydrodynamic effects are significant.

##  3 Hydrodynamic Results

In this section, we employ our model to simulate the evolution of a jet and
compare the results with those obtained from full numerical hydrodynamics. We
test our solutions for a top-hat jet moving in two different mediums: a
constant interstellar medium (ISM), where $k=0$, and a wind medium, where
$k=2$. Compared to structured jets, a top-hat jet exhibits more significant
non-linear hydrodynamic effects due to its sharp edge, such as the emergence
of a bow shock around the jet. This feature presents a challenging scenario to
test our code under extreme conditions.

The full numerical hydrodynamics simulations are conducted using the moving-
mesh code JET (Duffell & MacFadyen, 2013). In this code the grid lines move
along the radial directions at local speeds of fluid, thereby capable of
capturing flow structures with high resolution and high Lorentz factors. The
simulations are two-dimensional and are performed in spherical coordinates,
assuming axial-symmetry.

###  3.1 The jet evolution

To begin with, we present the typical behavior of jet evolution as predicted
by our model. A jet generally experiences four stages during its evolution
(Kumar & Zhang, 2015). Initially, the total mass of the blast shell is
dominated by the ejecta carried with the jet, and the jet coasts at a constant
velocity. After the accumulated mass exceeds the mass of the ejecta, the jet
begins to decelerate, following the Blandford-McKee solution. When the jet
becomes sub-relativistic, the fluid elements along the surface of the blast
wave become causally connected, leading to the jet starting its lateral
spreading due to internal pressure. Eventually, the jet transitions to a
Newtonian phase, and the shape of the blastwave becomes mostly spherical. The
subsequent evolution is best described by the Sedov-Taylor solution.

In Fig. 2 we show the evolution of proper velocity for isotropic blastwaves
with varying initial Lorentz factors moving in an ISM environment. The time is
normalized by a scale $r_{0}/c$, where $r_{0}=(3E_{\rm
iso}/4\pi\rho_{0}c^{2})^{1/3}$ is the radius at which the jet becomes sub-
relativistic. For comparison, the exact self-similar solutions are presented
as well. Our results align with both the exact Blandford-McKee and Sedov-
Taylor solutions in their respective limits and exhibit a smooth transition in
between.

Figure 2: The evolution of proper velocity for isotropic explosion with
various initial Lorentz factors ($\gamma_{0}$=10, 100, 1000 for green, orange,
blue lines, respectively). The time is scaled by $r_{0}/c$ which is the time
when the jet becomes sub-relativistic. The dashed lines are the exact
Blandford-McKee and Sedov-Taylor self-similar solutions.

In Figure 3, we show the evolution of angular structure for a top-hat jet with
opening angle 0.1 radian moving through ISM. The four colors in the figure
correspond to the snapshots of four phases of the jet’s evolution, as
indicated in the labels. During the deceleration phase (green lines), a bow
shock emerges at the edge of the jet, forming a sub-relativistic tail. As the
core decelerates to mildly relativistic speeds (orange lines), the jet begins
lateral expansion due to its internal pressure. In the Newtonian phase (red
line), the jet gradually becomes spherical.

Figure 3: The evolution of a top-hat jet with half opening angle of 0.1 radian
traveling in a constant density interstellar medium. Different phases of the
jet are represented by the corresponding colors: blue lines the coasting
phase, green lines the deceleration phase, orange lines the spreading phase,
and red lines the Newtonian phase. Left panel: normalized energy. Middle
panel: proper velocity. Right panel: radius scaled with $R(\theta=0)$.

###  3.2 Comparison with full numerical hydrodynamics

We further compare our model to full numerical hydrodynamic code JET, in
particular the evolution of jet structure during the spreading phase.

In the first case, we consider a top-hat jet with a half-opening angle of
0.052 radians (3 degrees) traveling through ISM. The initial condition for the
JET simulation is a Blandford-McKee radial profile. In Fig. 4 we show the
evolution of normalized energy (rest mass excluded) at three stages: the
initial condition (blue lines), the spreading phase (green lines), and the
Newtonian phase (red lines). The comparison reveals good agreement with the
core of the jet, and shows a more extended low energy tail outside the core.

Figure 4: The comparison between jetsimpy and the numerical hydrodynamic code
JET. The numerical setup is a top-hat jet with a half opening angle of 0.052
radians (3 degrees) traveling in a uniform external medium. We present the
comparison of three stages during the evolution: the blue lines that
represents the initial condition, the green lines that represent the
intermediate stage of spreading, and the red lines that represent the
Newtonian phase.

In the second comparison we consider a top-hat jet with a half opening angle
of 0.1 radians (6 degrees) traveling in a wind medium. The result is shown in
Fig. 5. This comparison also reveals overall good agreement, but is somewhat
less accurate than in the ISM case. The slightly increased discrepancy in the
wind medium test may indicate that our model is less accurate if the density
of the external medium varies steeply as a function of radius. Consequently,
we consider our model more accurate in the ISM environment as compared to the
wind medium.

Figure 5: Same as last figure, but for a top-hat jet with a half opening angle
of 0.1 radians (6 degrees) traveling in a wind medium.

The major discrepancy appearing in the two tests is that our model generates a
more extended tail at large angles. This discrepancy arises from the
infinitely thin shell approximation, where the entire radial profile is
projected onto a single surface. As a result, the regions inside and outside
of the edge of the jet becomes adjacent to each other, although they are in
fact spatially distant along the radial direction. This effect leads to energy
exchange between the regions that doesn’t occur in practice. Nevertheless,
this effect only appears far away from the jet core where the energy is orders
of magnitudes lower than the core. Therefore, we find it doesn’t practically
impact the afterglow modeling.

###  3.3 Computational performance

In addition to the good agreement, our code performs significantly faster than
full numerical hydrodynamics. The performance boost can be attributed to two
aspects. First, with one less dimension one can significantly reduce the
required number of grid cells of the simulation. Second, the
Courant–Friedrichs–Lewy (CFL) condition (see Appendix A) is more relaxed than
a full simulation. The CFL condition limits the simulation time step to the
signal-crossing time of the grid cells in order to maintain numerical
stability. In our approach, the simulation is performed along the polar
direction, so the signal speed is effectively the angular speed. As the
blastwave expands, the angular speed decreases due to the increasing radius.
Consequently, the time step can grow progressively larger as the jet evolves,
leading to an improved computational efficiency.

##  4 Afterglow modeling

In this section, we describe the modeling of GRB afterglows by post-processing
the numerical data obtained from our hydrodynamic models. In the current stage
of the code, we consider the modeling of light curves, flux centroid proper
motion, and image size.

To start with, we consider an observer who performs the actual observation,
referred to as the lab observer. It is assumed that the lab observer’s line of
sight lies within the x-z plane of the jet’s coordinate system, and is offset
from the z-axis by a polar angle $\theta_{\rm v}$. From the lab observer’s
perspective, the jet’s image is projected onto a plane perpendicular to the
line of sight. To model the size of this image, we establish a two-dimensional
Cartesian coordinate system, $\tilde{x}-\tilde{y}$, within this plane,
referred to as the line of sight (LOS) coordinate system. This LOS coordinate
system is defined so that the jet’s motion is oriented along the $\tilde{x}$
direction, and the $\tilde{y}$ axis aligns with the y-axis of the jet’s
coordinate system.

###  4.1 Relativistic effects

When modeling the radiation from a relativistic expanding shell, two major
relativistic effects must be considered: the Doppler beaming and the time
difference in arrival between different latitudes. Due to the spreading and
non-zero tangent velocity, the calculation of these effects involves some
complexities, which we will describe below.

We begin by addressing the arrival time difference effect. When the jet’s
expansion speed is ultra-relativistic, the emission of fluid elements from
various latitudes that appear to arrive simultaneously at the lab observer do
not actually get emitted at the same time in the burster frame. The geometric
shape that characterizes the surface where emissions reach the lab observer at
the same time is commonly referred to as the “equal arrival time surface,” and
is defined by the following equation:

$t-R\mu_{\rm r}/c=T_{\rm obs}/(1+z).$ (38)

In this equation, $\mu_{\rm r}$ is the cosine angle between the radial
direction and the line of sight, $T_{\rm obs}$ is the observed time of
radiation since the initial burst, and $1+z$ is the time dilation factor due
to cosmological red shift. The cosine angle can be calculated by the inner
product of radial norm $\hat{r}$ and line of sight norm $\hat{n}$:

$\displaystyle\mu_{\rm r}$ $\displaystyle=\hat{r}\cdot\hat{n}$
$\displaystyle=\cos\theta\cos\theta_{\rm v}+\sin\theta\cos\phi\sin\theta_{\rm
v}$ (39)

In the case of an isotropic shell expanding at a constant speed, the equal
arrival time surface is an ellipsoid. However, for a structured jet, the shape
becomes non-trivial and must be solved numerically. One must convert from
$T_{\rm obs}$ to $t$ by solving this equation in order to derive the actual
emission time in the burster frame.

To efficiently solve the equation, we first use a binary search to identify
the range of the solution within the discrete numerical data. The data is
further linear interpolated in this range to approximate $R(t,\theta)$. This
routine has significantly improved the performance of the code.

When considering lateral expansion, potential peculiar solutions could arise
if the normal vector of a relativistic blast wave significantly deviates from
the radial direction. This situation might lead to the coordinate speed
$\partial R/\partial t$ exceeding the speed of light. While this doesn’t
reflect a problem for hydrodynamics because it doesn’t represent the actual
speed of a fluid element, it does introduce numerical challenges in solving
the equal arrival time surface. This is because a fluid element from one
direction can “catch up with” a light emitted from another direction, which
leads to multiple solutions of equal arrival time surface at a given pair of
$T_{\rm obs}$ and $\theta$.

In this work, we assume a single solution and doesn’t employ specific
treatment to address this issue, because it hasn’t been actually observed in
any of the numerical setups we have tested thus far. This might be due to the
fact that spreading happens only in the mildly relativistic regime. However,
we leave this caveat to acknowledge that it could potentially present a
numerical challenge in certain scenarios.

Another relativistic effect is the Doppler beaming, which is characterized by
the Doppler factor

$\delta=\frac{1}{\gamma(1-\beta\mu_{\beta})}.$ (40)

Here, $\mu_{\beta}$ represents the cosine angle between the velocity of the
fluid element and the direction of the line of sight. Due to the non-zero
polar velocity, we cannot substitute $\mu_{\beta}$ with $\mu_{\rm r}$ as is
done in many semi-analytical models. The exact value of $\mu_{\beta}$ can be
calculated by taking the inner product of velocity norm $\hat{\beta}$ and line
of sight norm $\hat{n}$ like above. This calculation yields

$\displaystyle\mu_{\beta}$ $\displaystyle=\hat{\beta}\cdot\hat{n}$ (41)
$\displaystyle=(\frac{\beta_{\rm
r}}{\beta}\sin\theta\cos\phi+\frac{\beta_{\theta}}{\beta}\cos\theta\cos\phi)\sin\theta_{\rm
v}$ (42) $\displaystyle+(\frac{\beta_{\rm
r}}{\beta}\cos\theta-\frac{\beta_{\theta}}{\beta}\sin\theta)\cos\theta_{\rm
v}.$ (43)

Numerical tests have revealed minor differences between using the exact
expression and approximating it with $\mu_{\rm r}$. This can be attributed to
the fact that spreading primarily occurs when the jet reaches mildly
relativistic speeds, and beaming effect is relatively unimportant during this
phase.

###  4.2 Light Curves

We now proceed to model the light curves from a relativistic jet. Light curve
is one of the most important observational outcomes of a GRB afterglow. For
simplicity, here we consider a minimal assumption of synchrotron radiation
following Sari et al. (1998). We briefly summarize it as following.

We assume that the electrons accelerated by the blastwave follow a non-thermal
power-law distribution with a power index of $-p$. The minimal Lorentz factor
$\gamma_{\rm m}$ of electrons can be determined by assuming that a fraction
$\epsilon_{\rm e}$ of the internal energy is converted into electron kinetic
energy. The result is expressed as:

$\gamma_{\rm m}=\frac{p-2}{p-1}\frac{\epsilon_{\rm e}m_{\rm p}}{m_{\rm
e}}(\gamma-1)$ (44)

The relativistic electrons are assumed to move in a randomly oriented magnetic
field. The average magnetic field can be calculated by assuming that a
fraction $\epsilon_{\rm B}$ of the internal energy is converted into magnetic
energy:

$B^{\prime}=\sqrt{8\pi\bar{e}^{\prime}_{\rm sw}\epsilon_{\rm B}}$ (45)

where $n_{0}=\rho_{0}/m_{\rm p}$ is the number density of external medium.

The synchrotron radiation spectrum can be approximated by a multi-band broken
power-law, which depends on several characteristic break frequencies. In this
work, we consider two such breaks: the frequencies emitted by electrons with
the minimum Lorentz factor $\gamma_{\rm m}$ and the cooling frequency
associated with the Lorentz factor $\gamma_{\rm c}$. For simplicity, we do not
consider synchrotron self-absorption. The cooling Lorentz factor is calculated
via global cooling approximation

$\gamma_{\rm c}=\frac{6\pi m_{\rm e}\gamma c}{\sigma_{\rm T}B^{\prime 2}t}$
(46)

The break frequencies can be calculated by

$\nu^{\prime}_{\rm i}=\frac{3eB^{\prime}\gamma_{\rm i}^{2}}{4\pi m_{\rm e}c}$
(47)

where ${\rm i}=\\{{\rm m},{\rm c}\\}$ and $e$ is the electron charge. The
synchrotron emissivity in the comoving frame is now expressed by

$\displaystyle\epsilon^{\prime}_{\nu^{\prime}}=\epsilon^{\prime}_{P}$ (48)
$\displaystyle\times\begin{cases}(\nu^{\prime}/\nu^{\prime}_{m})^{1/3}&\nu^{\prime}<\nu^{\prime}_{m}<\nu^{\prime}_{c}\\\
(\nu^{\prime}/\nu^{\prime}_{m})^{-(p-1)/2}&\nu^{\prime}_{m}<\nu^{\prime}<\nu^{\prime}_{c}\\\
(\nu^{\prime}_{c}/\nu^{\prime}_{m})^{-(p-1)/2}(\nu^{\prime}/\nu^{\prime}_{c})^{-p/2}&\nu^{\prime}_{m}<\nu^{\prime}_{c}<\nu^{\prime}\\\
(\nu^{\prime}/\nu^{\prime}_{m})^{1/3}&\nu^{\prime}<\nu^{\prime}_{c}<\nu^{\prime}_{m}\\\
(\nu^{\prime}/\nu^{\prime}_{c})^{-1/2}&\nu^{\prime}_{c}<\nu^{\prime}<\nu^{\prime}_{m}\\\
(\nu^{\prime}_{m}/\nu^{\prime}_{c})^{-1/2}(\nu^{\prime}/\nu^{\prime}_{m})^{-p/2}&\nu^{\prime}_{c}<\nu^{\prime}_{m}<\nu^{\prime}\end{cases}$
(49)

The peak emissivity can be calculate by

$\epsilon^{\prime}_{\rm P}=\frac{\sqrt{3}e^{3}B^{\prime}n^{\prime}}{m_{\rm
e}c^{2}}$ (50)

where $n^{\prime}=\rho^{\prime}_{\rm sw}/m_{\rm p}$ is the number density of
the blast.

The total luminosity can be determined by integrating the emissivity over the
radiating volume. In our model, since the radiating region is approximated as
an infinitely thin surface, we need to estimate a typical shell width to
perform the volume integral. This estimation can be accomplished by relating
the mass per solid angle ($M_{\rm sw}$) with the mass per volume
($\rho^{\prime}_{\rm sw}$) and realize that $dV=r^{2}drd\Omega$. As such, the
typical shell width can be approximated as

$\Delta R\approx\frac{M_{\rm sw}}{R^{2}\rho^{\prime}_{\rm sw}\gamma}$ (51)

where we recover $\rho^{\prime}_{\rm sw}$ from the delta function
approximation to the value determined by the shock jump condition.

With the estimation of typical shell width, the total luminosity can be
estimated by the integral

$L_{\nu}=\int\epsilon^{\prime}_{\nu^{\prime}}\delta^{3}\Delta
R^{\prime}R^{2}d\Omega$ (52)

where $\Delta R^{\prime}=\Delta R\gamma$ is the shell width in the comoving
frame. $\nu=\nu^{\prime}\delta/(1+z)$ is the observed frequency. The integral
is performed in the jet spherical coordinate at lab observer’s time $T_{\rm
obs}$. The time of radiation $t$ for each fluid element is solved by the equal
arrival time surface discussed above.

The integral is carried out using a 2-dimensional adaptive integration
algorithm. Due to the relativistic beaming, the integrand exhibits a very
narrow peak. To efficiently and accurately perform the integration, we first
employ a minimization algorithm to identify the peak. After that, we rotate
the spherical coordinates to align the peak with the north pole. This approach
takes advantage of the fact that polar lines converge at the poles, which
naturally results in an oversampling of integration points around the peak,
leading to improved integration efficiency. These implementations offer a
significant enhancement in performance compared to a direct integration in the
jet coordinate system.

Given the luminosity, we can calculate the flux density by

$F_{\nu}=\frac{L_{\nu}(1+z)}{4\pi D_{\rm L}^{2}}$ (53)

where $D_{\rm L}$ is the luminosity distance.

The above minimal synchrotron radiation model is set as default in jetsimpy,
and is sufficient for afterglow modeling in most cases. Nevertheless, the
flexibility of the code allows for high level of customization with more
sophisticated modeling.

###  4.3 Centroid proper motion and image size

Another important observable of GRB afterglow is the flux centroid proper
motion, which provides a directly constraint of the Lorentz factor. The proper
motion of the jet image in the LOS plane can be defined by the offset of a
“flux centroid” from the burst origin. Since the image of the jet has an
diffusive scale, the overall offset can be defined by a weighted average over
the intensity of the image (e.g., Zrake et al. 2018; Fernández et al. 2022;
Mooley et al. 2022; Nedora et al. 2023a). Therefore, the flux centroid
position in the LOS plane can be calculated by

$\tilde{x}_{\rm c}=\frac{\int R\lambda_{\rm
r}\frac{dL_{\nu}}{d\Omega}d\Omega}{\int\frac{dL_{\nu}}{d\Omega}d\Omega}$ (54)

where $dL/d\Omega$ is defined in eq.52, and $\lambda_{\rm r}$ is the sine
angle between the radial direction and the line of sight. It can be calculated
by

$\lambda_{\rm r}=-\sin\theta\cos\phi\cos\theta_{\rm
v}+\cos\theta\sin\theta_{\rm v}$ (55)

The angular offset from the burst origin can be calculated by

$\tilde{\theta}_{\rm c}=\frac{\tilde{x}_{\rm c}}{D_{\rm A}}$ (56)

where $D_{\rm A}=D_{\rm L}/(1+z)^{2}$ is the angular diameter distance.

We can also calculate the image size of the jet in a similar manner. Following
the definitions adopted in previous works (e.g., Ryan et al. 2023), the image
size is calculated by taking the moments of the intensity distribution with
respect to the flux centroid. The Gaussian-equivalent image sizes along the
$\tilde{x}$ and $\tilde{y}$ directions are defined as follows:

$\displaystyle\tilde{\sigma}^{2}_{\rm x}$
$\displaystyle=\frac{\int(R\lambda_{\rm r}-\tilde{x}_{\rm
c})^{2}\frac{dL_{\nu}}{d\Omega}d\Omega}{\int\frac{dL_{\nu}}{d\Omega}d\Omega}$
(57) $\displaystyle\tilde{\sigma}^{2}_{\rm y}$
$\displaystyle=\frac{\int(R\sin\theta\sin\phi-\tilde{y}_{\rm
c})^{2}\frac{dL_{\nu}}{d\Omega}d\Omega}{\int\frac{dL_{\nu}}{d\Omega}d\Omega}$
(58)

Given the axial symmetry of the jet, we have $\tilde{y}_{\rm c}=0$. The
angular sizes are then calculated as:

$\displaystyle\tilde{\theta}_{\rm x}$ $\displaystyle=\frac{\tilde{\sigma}_{\rm
x}}{D_{\rm A}}$ (59) $\displaystyle\tilde{\theta}_{\rm y}$
$\displaystyle=\frac{\tilde{\sigma}_{\rm y}}{D_{\rm A}}$ (60)

##  5 Tests of Afterglow Modeling

In this section we test our afterglow modeling described in the previous
section by comparing them to existing tools. For the light curve calculation,
we compare our code with the widely used tools afterglowpy (Ryan et al., 2020)
and BoxFit (van Eerten et al., 2012). These tools represent the state-the-art
modeling of GRB afterglow with semi-analytical methods or full numerical
hydrodynamics. For the centroid proper motion and image size, we directly
compare them to the simulated imaging of the GRB jet.

###  5.1 Comparison to afterglowpy

The package afterglowpy (Ryan et al., 2020) is a widely-used GRB afterglow
modeling tool with semi-analytical spreading prescription. In this model a
structured jet is approximated by the combination of various concentric top-
hat jets moving in a ISM environment (see also Lamb & Kobayashi 2017; Lamb et
al. 2018; Fernández et al. 2022; Nedora et al. 2023b with a similar approach).
Each jet segment laterally expands at local sound speed. The spreading is
initial frozen, until the jets become mildly relativistic. This semi-
analytical approach provides a reasonable approximation of the spreading phase
at a much lower computational cost compared to hydrodynamic simulations.
However, because the interactions between the fluid elements are neglected, it
can’t account for all non-linear hydrodynamic effects or achieve the same
level of accuracy as a full simulation.

Our model shares the same trans-relativistic equation of state and thin shell
approximation as in afterglowpy, as well as a very similar radiation
calculation. The major differences are our more flexible initial angular jet
profile, the inclusion of the coasting phase, a different spreading
prescription, the calibration with self-similar solutions, and the support for
the wind environment. In fact, if all these effects are disabled, the two
models agree exactly with each other.

To compare the light curves for a structured jet, we assume the isotropic
equivalent energy and initial Lorentz factor follow a Gaussian distribution
with a typical core width $\theta_{\rm c}$:

$\displaystyle E(\theta)$
$\displaystyle=E_{0}\exp{\left[-\frac{1}{2}\left(\frac{\theta}{\theta_{\rm
c}}\right)^{2}\right]}$ (61) $\displaystyle\gamma(\theta)$
$\displaystyle=(\gamma_{0}-1)\exp{\left[-\frac{1}{2}\left(\frac{\theta}{\theta_{\rm
c}}\right)^{2}\right]}+1.$ (62)

The comparisons of light curves are presented in Fig. 6, with the parameters
provided in the caption. In this comparison, we keep the viewing angle fixed
at $\theta_{\rm v}=0.3$ and vary the half-opening angle $\theta_{\rm c}$. To
account for radiation from all directions, we set the truncation angle in
afterglowpy to $\theta_{\rm w}=\pi/2$.

Figure 6: Comparison of light curves to the publicly available code
afterglowpy. The jet structures are assumed to be Gaussian profile with core
width $\theta_{\rm c}$ labeled in the figure. The rest of the parameters for
the light curves are: $n_{0}=1$ cm-3, $E_{0}=10^{51}$ erg, $\gamma_{0}=1000$,
$\theta_{\rm v}=0.3$, $\epsilon_{\rm e}=0.1$, $\epsilon_{\rm B}=0.01$,
$p=2.5$, $\nu=10^{18}$ Hz, $d_{\rm L}=474.33$ Mpc, and $z=0.1$. Top Panel:
light curves with spreading effect and calibration disabled
($\gamma_{0}=\infty$ and $s=1$). Bottom Panel: light curves with these effects
enabled.

In the upper panel, we enable only the spreading effect, while disabling the
coasting phase and calibration (i.e., $\gamma_{0}=\infty$ and $s=1$). This
comparison is intended to isolate and reveal the differences caused solely by
the modeling of the spreading effect. We observe overall good agreement in all
cases, suggesting that non-linear hydrodynamic effects are not significant for
structured jets. However, we do notice minor discrepancies especially for
narrow jets. Specifically, our model predicts a slightly faster spreading
effect than afterglowpy for very narrow jets (e.g., the blue line). Such rapid
spreading is expected in these cases, as the shock front in narrow jets is
more causally connected, leading to an exponential expansion similar to the
prescription in Rhoads 1999. This fast spreading effect for extremely narrow
jets has also been observed in a number of numerical hydrodynamic studies
(e.g., Granot & Piran 2012; Nedora et al. 2023b; Govreen-Segal & Nakar 2023).
The minor discrepancy observed at later times is likely due to the spreading
prescription of afterglowpy where the counter-hemisphere is not well modeled.
The slopes eventually align with each other after a sufficient long time.

In the lower panel, we maintain the same numerical setup, with the addition of
enabling both the coasting effect and calibration. The initial Lorentz factor
of the core is assumed to be 1000. This finite initial Lorentz factor leads to
an initial rising phase following a slope of $F_{\nu}\propto t^{3}$, resulting
in a significant difference for the early time. The effect of calibration is
approximately a temporal offset, since the coefficient can be equivalently
understood as a scaling to the external density. We also find the calibration
leads to a slightly faster decay in the spreading phase. This effect is due to
the transition of the calibration coefficients between the self-similar
solutions in each limits. In the relativistic limit, this coefficient is less
than 1 (see Fig. 1), resulting in a predicted velocity that is faster than
uncalibrated predictions. Conversely, in the Newtonian limit, the coefficient
is greater than 1, leading to a lower velocity prediction. A a result, the
combined prediction of velocity decreases more rapidly in the transitional
regime than in either of the limits, which occurs exactly during the spreading
phase.

###  5.2 Comparison to BoxFit

We further compare to the light curves that have been post-processed using
full numerical hydrodynamics, specifically, the widely used tool BoxFit (van
Eerten et al., 2012). This code employs top-hat angular structures and a
Blandford-McKee radial profile for its initial setup. The simulation data are
then interpolated for post-processing the light curve and spectrum modeling.
BoxFit calculates the observed flux density through a ray-tracing radiative
transfer module. Because the code is based on numerical simulations, the
synthetic light curves are considered highly reliable. However the code only
works for a top-hat jet structure. Moreover, BoxFit supports a wind external
medium, making it a unique tool for comparisons in such scenarios.

In Figure 7, we present the comparisons of radio and optical light curves for
both on-axis and off-axis observations in an ISM environment. Considering that
non-linear hydrodynamic effects are more significant in narrower jets, we
performed a test to evaluate the reliability of our code under these extreme
conditions. For this purpose, the comparison was carried out using the minimum
jet opening angle allowed by BoxFit, which is 0.045 radians. To align our
numerical setup with BoxFit, we disable the coasting phase while enabling the
spreading effect and calibration. Additionally, we add a counter-jet to assess
the robustness of our model under Newtonian limits. As a reference, the light
curves produced by afterglowpy is presented as well. Details of the jet
parameters and light curves are provided in the figure caption.

Figure 7: Comparison of light curves to the publicly available codes BoxFit
and afterglowpy. We assume a top-hat jet with isotropic equivalent energy of
$10^{52}$ erg, opening angle of 0.045 radian (3 degrees), and a symmetric
counter-jet. The initial Lorentz factor is infinite and the external medium is
ISM-like. The rest of the parameters are: $n_{0}=1$ cm-3, $E_{\rm
iso}=10^{52}$ erg, $\epsilon_{\rm e}=0.1$, $\epsilon_{\rm B}=0.01$, $p=2.5$,
$d_{\rm L}=10^{28}$ cm, and $z=0$. The viewing angle $\theta_{\rm v}$ and the
frequency $\nu$ are labeled in the figure. Top panel: Radio and optical light
curves. Bottom panel: Relative differences as compared to BoxFit.

Not surprisingly, our results show good agreement with BoxFit, typically
within a factor of 50%. In particular, the agreement in the spreading phase is
slightly better than afterglowpy. The major discrepancies observed in this
comparison can be explained as follows. The significant difference at early
times for off-axis cases (brown and olive lines) is attributed to the bow
shock (see §3.1) generated outside the jet’s core. This effect is particularly
significant in a top-hat jet due to its very sharp edge. In BoxFit
simulations, these time periods are earlier than the startup time, and the
calculations are based on extrapolations assuming a Blandford-McKee radial
profile and top-hat angular profile, where bow shocks are not present. The
rising times of counter-jets also show a temporal offset. This suggests that
the assumption of an infinitely thin shell may not hold in the Newtonian
limit, as we have anticipated. However, the temporal difference remains within
a factor of 30%. Additionally, the early spreading phase for off-axis
observation (the olive line) shows a minor deviation, likely due to the more
extended tail discussed in §3.2, but the difference remains within an order of
unity. This comparison demonstrates that our model can achieve an overall
accuracy comparable to that of full numerical hydrodynamics, even when non-
linear effects are significant.

In Figure 8, we present the comparison of our model with BoxFit in a wind
environment. For simplicity, we do not include counter-jets in this
comparison, as the BoxFit simulations in these cases are performed in a
Lorentz boosted frame, which does not resolve counter-jets effectively. This
comparison also demonstrates overall agreement, typically within a factor of
unity, but is somewhat less precise than in ISM environments, as we have
expected. The deviations at early time are still due to the bow shocks emerged
outside the jet core. It is worth noting that in the wind environment, this
level of agreement is only attainable when calibration is enabled. Without
calibration, the results will exhibits a significant temporal offset. This
suggests that the calibration coefficient has a bigger impact on solutions in
wind environments as compared to ISM environments.

Figure 8: Same as last figure, but for a wind-like environmental medium
profile: $n_{0}=1(r/10^{17}{\rm cm})^{-2}$ cm-3. The counter-jet is absent in
this comparison.

###  5.3 Tests of Flux Centroid Proper Motion and Image Size Modeling

To validate our calculation of the centroid proper motion and image size, we
directly compare it to simulated images generated from the numerical data. The
procedure for calculating afterglow imaging is summarized as follows: For each
pixel of the image we consider a ray parallel to the line of sight but offset
from the origin by a corresponding distance. This ray is then intersected with
the equal arrival time surface. At each intersection point, we approximate the
observed intensity using $I_{\nu}=\epsilon^{\prime}_{\nu^{\prime}}\Delta
R^{\prime}\delta^{3}/4\pi$. The total intensity for a given pixel is
calculated as the sum of intensities from all intersecting points. The image
reflects the observed intensity distribution within the LOS plane.

Our results are presented in Fig. 9 with parameters detailed in the caption.
The images, arranged from top to bottom, correspond to three typical stages of
the jet: the coasting phase, the deceleration phase, and the spreading phase.
The flux centroid locations are marked with “+” signs, and the image sizes are
presented by the $1-\sigma$ contours of the Gaussian-equivalent intensity
distribution. Our modeling of the flux centroid and image size closely aligns
with the simulated images, thereby validating the accuracy of our method.

Figure 9: The image of a jet with Gaussian profile traveling in an ISM
environment. From top to the bottom are the images at coasting phase,
deceleration phase, and spreading phase. The intensity distributions are
normalized by maximum intensity. The markers “$\star$” and “$+$” represents
the burst origin and the flux centroid respectively. The white contour is the
$1-\sigma$ size assuming Gaussian equivalent intensity distribution. The
parameters for the calculation are: $n_{0}=1$ cm-3, $E_{0}=10^{52}$ erg,
$\gamma_{0}=1000$, $\theta_{\rm c}=0.1$ rad, $\theta_{\rm v}=0.3$ rad,
$\epsilon_{\rm e}=10^{-2}$, $\epsilon_{\rm B}=10^{-3}$, $p=2.2$, $\nu=3\times
10^{9}$ Hz, $d_{\rm L}=44$ Mpc, and $z=0.01$.

##  6 Application to GW170817

In this section, our model is applied to real-world data from the
gravitational wave event GW170817 (Abbott et al., 2017a, b, c). This event
represents the first multimessenger observation of neutron star mergers. It is
also the first GRB afterglow whose relativistic jet is evidently observed from
an off-axis direction (see Margutti & Chornock 2021 for a review). The multi-
waveband observation spanning a significant duration, together with the well-
measured apparent superluminal motion (Mooley et al., 2018; Ghirlanda et al.,
2019; Mooley et al., 2022), makes it an ideal object to apply our model.

To better fit the data of this event, we have made a slight adjustment to the
synchrotron radiation model described in §4.2, incorporating the correction of
the “deep Newtonian phase” at late stages. This adjustment is necessary when
the minimal Lorentz factor of electrons, $\gamma_{\rm m}$, as calculated by
Eq. 44, drops below 1, which happens when the velocity of the blastwave has
sufficiently decreased. In this scenario, electrons follow a power-law
distribution in momentum space, rather than in energy space. Consequently,
only a subset of electrons are relativistic and contribute to the emission in
the GHz waveband. This adjustment might be necessary, since a recent study
(Ryan et al., 2023) has shown it to yield improved fits for late-time data.
Here we follow the approach described in Sironi & Giannios 2013. Specifically,
if $\gamma_{\rm m}$ falls below 1, we fix it at 1 and reduce the number of
radiating electrons by a fraction of $f$:

$f=\frac{p-2}{p-1}\frac{\epsilon_{\rm e}m_{\rm p}}{\gamma_{\rm m}m_{\rm
e}}(\gamma-1).$ (63)

This correction leads to a late time flattening of the afterglow light curve.

We emphasize that the objective of this section is not to conduct a
comprehensive study of GW170817 but rather to demonstrate the application of
our code using real-world data. As such, we do not explore the effects of
varying jet structures or the potential influence of additional phenomena,
such as kilonova afterglows (Nakar & Piran, 2011; Piran et al., 2013;
Alexander et al., 2017; Hotokezaka et al., 2018; Kathirgamaraju et al., 2019),
that may emerge at later stages. Furthermore, gravitational wave data is not
incorporated into our analysis, given its irrelevance to afterglow modeling. A
comprehensive study of GW170817 using our newly developed code is left for
future research.

###  6.1 data fitting

To fit the data with our model, we perform a Bayesian parameter estimation.
For this analysis, multi-waveband light curve data up to 1000 days are taken
from Makhathini et al. (2021), where the data have been re-analyzed to
minimize discrepancies arising from different data processing techniques. The
original data references are detailed within their work. Beyond 1000 days, we
incorporate radio wave observations from the Very Large Array (VLA)
(Balasubramanian et al., 2021) and X-ray measurements from Chandra
(Balasubramanian et al., 2021, 2022; O’Connor & Troja, 2022). Additionally, we
take into account the flux centroid positions measured by the Very Large
Baseline Interferometry (VLBI) (Mooley et al., 2018; Ghirlanda et al., 2019)
and the Hubble Space Telescope (HST) (Mooley et al., 2022), with the HST
measurement presumed to mark the origin of the burst. Since there is a
correlated systematic error between the VLBI and HST measurements due to
different coordinate systems, for simplicity, this error is integrated into
the HST data point. The total uncertainty for each offset measurement is
calculated as the root mean square of the statistical and systematic errors.
Since gravitational wave data is excluded from our analysis, we set the
luminosity distance at $D_{\rm L}$=43.9 Mpc and the redshift at $z=0.0098$,
assuming cosmological parameters adopted from Planck (Planck Collaboration et
al., 2020). To simplify the analysis, data points representing upper limits
are omitted from all mentioned datasets.

In this study, we adopt a Gaussian profile for the angular distribution of the
jet, as described in the previous section. The absence of an obvious coasting
phase leads us to presume an infinitely large initial Lorentz factor. The lack
of evident late time re-brightening at the current stage motives us to exclude
the counter-jet component (e.g. Li et al. 2024). To fit the apparent
superluminal motion, we assume that the trajectory of the flux centroid within
the LOS plane originates from the coordinates $(RA_{0},Dec_{0})$ and with an
orientation angled of $\alpha$. Similar to previous studies (Mooley et al.,
2018; Ghirlanda et al., 2019; Mooley et al., 2022; Gianfagna et al., 2023;
Ryan et al., 2023), we place the coordinate origin at the centroid’s position
at 75 days. Under these assumptions, we identify the following set of free
parameters for our model fitting: $\vec{\vartheta}$=$\log_{10}n_{0}$,
$\log_{10}E_{0}$, $\theta_{\rm c}$, $\theta_{\rm v}$, $\log_{10}\epsilon_{\rm
e}$, $\log_{10}\epsilon_{\rm B}$, $p$, $RA_{0}$, $Dec_{0}$, $\alpha$.

According to the Bayes’ theorem, the posterior probability distribution of the
fitting parameters can be expressed as:

$p(\vec{\vartheta}|\vec{d}_{\rm F},\vec{d}_{\rm
c})\propto\mathcal{L}(\vec{d}_{\rm F}|\vec{\vartheta})\mathcal{L}(\vec{d}_{\rm
c}|\vec{\vartheta})\pi(\vec{\vartheta}),$ (64)

where $\vec{d}_{\rm F}$ and $\vec{d}_{\rm c}$ represent the datasets for flux
density and flux centroid offset, respectively. The likelihood $\mathcal{L}$
for both datasets is modeled as $\exp(-\chi^{2}/2)$, with $\chi^{2}$ being the
reduced chi-squared value calculated from the discrepancies between observed
data and the model predictions. It has been observed by Ryan et al. 2023 that
the late-time X-ray observations for this event, due to their low photon
counts, would be better modeled using a Poisson likelihood. This adjustment
has not been implemented in our current analysis but is subject to future
exploration. The prior distribution $\pi(\vec{\vartheta})$ for the parameters
is detailed in Table 1.

Table 1: Prior distributions of fitting parameters. Parameters | Distributions | Bounds
---|---|---
$\log_{10}(n_{0}/{\rm cm}^{-3})$ | Uniform | [-5, 0]
$\log_{10}(E_{0}/{\rm erg})$ | Uniform | [49, 57]
$\theta_{\rm c}$ [rad] | Uniform | [0.01, $\pi/2$]
$\theta_{\rm v}$ [rad] | Sine | [0, $\pi$/2]
$\log_{10}\epsilon_{\rm e}$ | Uniform | [-6, 0]
$\log_{10}\epsilon_{\rm B}$ | Uniform | [-6, 0]
$p$ | Uniform | [2, 2.5]
$RA_{0}$ | Uniform | [-10, 10]
$Dec_{0}$ | Uniform | [-10, 10]
$\alpha$ | Uniform | [-$\pi$, $\pi$]

In this analysis, we conduct two separate Markov Chain Monte Carlo (MCMC)
studies: one that incorporates all aforementioned datasets, and another that
excludes the centroid offset data. To sample the posterior probability
distribution, we utilize the Bilby package, a public tool for Bayesian data
analysis (Ashton et al., 2019). Within Bilby, we employ the pymultinest
sampler (Buchner et al., 2014), which is based on a Nested Sampling algorithm.
Nested Sampling offers an improved computational efficiency for high-
dimensional problems.

In Fig. 10, we present the light curves generated using the sample with the
highest posterior probability, together with the 68% confidence intervals for
the flux density. The fitting results for the centroid offset are shown in
Figure 11. We summarize the parameter estimation results in Table 2. The full
corner plots of posterior probability distributions are available in Appendix
B.

Figure 10: Best fitting light curves for multi-waveband observation of GW170817 with flux centroid data included or not. The data points are normalized to 3 GHz, 5.06$\times 10^{14}$ Hz, and 1keV respectively. The shadow regions represent 68% confidence interval of the flux density.  Figure 11: The posterior distribution of flux centroid position with respect to the measurement data. The contour lines represent 68% and 95% confidence region. The measurement at 8 days is assumed to be the burst origin.  Table 2: Parameter estimation results for the two MCMC fittings. Parameters | LC | LC + Offset
---|---|---
$\log_{10}(n_{0}/{\rm cm}^{-3})$ | $-0.65^{+0.32}_{-0.36}$ | $-1.33^{+0.91}_{-1.03}$
$\log_{10}(E_{0}/{\rm erg})$ | $51.86^{+0.41}_{-0.48}$ | $54.53^{+0.92}_{-0.95}$
$\theta_{\rm c}$ [deg] | $7.55^{+0.62}_{-0.59}$ | $2.84^{+0.29}_{-0.23}$
$\theta_{\rm v}$ [deg] | $50.20^{+3.85}_{-3.74}$ | $18.16^{+1.86}_{-1.48}$
$\log_{10}\epsilon_{\rm e}$ | $-1.49^{+0.42}_{-0.36}$ | $-4.13^{+0.95}_{-0.93}$
$\log_{10}\epsilon_{\rm B}$ | $-3.27^{+0.41}_{-0.35}$ | $-3.86^{+1.00}_{-0.91}$
$p$ | $2.12^{+0.01}_{-0.01}$ | $2.12^{+0.01}_{-0.01}$
$RA_{0}$ [MAS] | - | $-1.95^{+0.30}_{-0.29}$
$Dec_{0}$ [MAS] | - | $0.24^{+0.35}_{-0.33}$
$\alpha$ [rad] | - | $-0.02^{+0.11}_{-0.11}$

###  6.2 fitting results and discussions

Our fitting results largely align with previous studies which also account for
apparent superluminal motion (e.g., Mooley et al. 2022; Gianfagna et al. 2023;
Ryan et al. 2023). We observe notable differences in the parameter estimations
between the two MCMC studies, particularly regarding the jet opening angle
($\theta_{\rm c}$) and the observing angle ($\theta_{\rm v}$). Specifically,
when the centroid motion is excluded, the light curves are well fitted, but
this leads to a considerably large jet opening angle (7.55 degrees) and
observing angle (50.2 degrees). Conversely, including all datasets results in
a good fit for the centroid motion data, albeit slight discrepancies in the
light curves at later times. This approach yields smaller estimates for both
the jet opening angle (2.84 degrees) and the observing angle (18.16 degrees).
This finding is consistent with the aforementioned works.

The discrepancy observed between the two sets of fittings arises from the
respective characteristics of the data sets used to deduce the jet opening
angle and observing angle. In both approaches, the ratio $\theta_{\rm
v}/\theta_{\rm c}$ is derived from the rising slope of the light curves, but
the data features for determining the respective values differ between the two
fittings. In the fitting that excludes the centroid motion, the jet opening
angle is primarily inferred from the decay phase after the peak. During this
phase, the jet transitions to a sub-relativistic speed, and the flux density
reflects the jet’s total energy. Conversely, prior to the peak, the flux
density is determined by the isotropic equivalent energy of a localized jet
region due to beaming effects. Therefore, the jet opening angle is deduced
from the scale ratio between these two phases. In the fitting that
incorporates centroid motion data, the observing angle is predominantly
derived from the apparent transverse velocity, and the likelihood function
associated with this measurement has a greater influence than that of the
late-time light curve data. Our fitting indicates that the two ways of
inferring the jet opening angle and observing angles yield different
estimation values.

This deviation may suggest a discrepancy between the late-time light curve and
the observed superluminal motion. The late-time light curve’s shallow decay
suggests a wide jet, whereas the superluminal motion implies a narrow jet.
However, using observation angles to infer the Hubble constant through joint
fitting with gravitational wave data, previous studies (e.g., Gianfagna et al.
2023) have shown that a narrow jet scenario aligns more closely with the
current constraints set by Planck (Planck Collaboration et al., 2020) and
SH0ES (Riess et al., 2022). This result implies that a narrow jet scenario is
more plausible than a wide jet. Thus, the discrepancy might indicate an excess
at late times, potentially caused by an additional component, such as a
kilonova afterglow.

While the late-time deviation is also noted in the aforementioned works, the
excess observed in our results is slightly more significant. This can be
attributed to our model which predicts a steeper decline in narrow jets during
the spreading phase, in contrast to semi-analytic methods (see §5.1 for a
detailed discussion). This finding suggests that the jet in this event may be
so narrow that its spreading phase is best represented by a hydrodynamic
model.

This discrepancy might also arise from other factors which does not require an
additional component. For instance, Ryan et al. 2023 have utilized a Poisson
likelihood function for the X-ray data, which significantly enhances the
fitting results compared to the chi-square likelihood function. However, this
approach is unable to account for the excess observed as well in the radio
waveband, and the radio dataset employed in this study is insufficient to span
the late-time light curve. Additionally, certain non-standard processes in GRB
afterglows may also explain the shallow decay observed at late times, such as
complex jet structures or changes in microphysical parameters during the
evolution. A more thorough investigation that incorporates these systematics
will be conducted in future research.

##  7 Summary and Conclusion

In this work, we have developed a reduced hydrodynamic model for relativistic
structured jets. The blastwave is approximated as an infinitely thin two
dimensional surface. By additionally assuming the axial-symmetry of the jet,
the simulation is effectively simplified into a one-dimensional
representation. The radial integrate values are calibrated with the exact
Blandford-McKee and Sedov-Taylor self-similar solutions in the respective
limits. This approximation is motivated by the effectively thin blast shell
and the nature of self-similar radial profile. Our method maintains the
accuracy of numerical hydrodynamics, while attains a computational efficiency
that is comparable to semi-analytic methods.

We have extensively compared the numerical results of our methods to existing
tools. The comparison reveals a high level of agreement, especially with those
from numerical hydrodynamics. We have also applied it to the case of GW170817,
fitting the light curve and flux centroid motion data. The parameter
estimations mostly align with previous results, but our findings indicate a
more evident late-time excess, suggesting the presence of an extra component
or non-standard processes. The possible existence of the additional component,
which is also suggested in previous works, motivates a further investigation
of GW170817.

The methods discussed in this work have been developed into a numerical tool
jetsimpy. The code is able to simulate a jet with arbitrary tabulated angular
energy and Lorentz factor profiles, and evolve the jet within ISM environment,
wind-like environment, or mixed environment. The simulation data is then used
to calculate the light curve or spectrum, provided with a customized
emissivity model. Additionally, the tool provides the calculation of the
apparent superluminal motion and the image size of the jet. jetsimpy is
optimized to be used in MCMC where millions of light curves with varying jet
structures need to be generated in a reasonable time. It is also worth noting
that although our model is built for GRB afterglow in mind, it’s not limited
to this scenario but any relativistic hydrodynamic jets. As such, it will
serve as a powerful tool for the community in the era of multimessenger
astronomy.

## Acknowledgments

We thank Paz Beniamini for the useful discussion and suggestions. H. W. and D.
G. acknowledge support from the NSF AST-22107802, AST-2107806 and AST-2308090
grants.

##  Appendix A The Numerical scheme

Here, we provide a concise overview of our numerical approach to solving
hydrodynamic equations. The governing equations closely resemble those of
regular numerical hydrodynamics, except that they are coupled with with a
Hamilton-Jacobi equation.

###  A.1 Eigenvalues of the Jacobian Matrix

Similar to numerical hydrodynamics, our approach requires the eigenvalues of
the Jacobian matrix to formulate the numerical scheme. In our model, due to
the direct involvement of approximate trans-relativistic shock jump
conditions, the eigenvalues are non-trivial and do not align with those found
in regular relativistic hydrodynamics.

To begin with, we first address the solution of the Jacobian matrix. This is
achieved by representing each component as a function of the conserved
variables:

$\displaystyle\left|\frac{\partial{\bf F}}{\partial{\bf
U}}\right|=\frac{c}{R}\begin{pmatrix}0&1&0&0\\\
(1-\beta_{\theta}^{2})(1+\frac{\partial P_{\rm sw}}{\partial E_{\rm
b}})-1&2\beta_{\theta}&(1-\beta_{\theta}^{2})\frac{\partial P_{\rm
sw}}{\partial M_{\rm sw}}&(1-\beta_{\theta}^{2})\frac{\partial P_{\rm
sw}}{\partial M_{\rm ej}}\\\ -\beta_{\theta}\frac{M_{\rm sw}}{H_{\rm
b}}(1+\frac{\partial P_{\rm sw}}{\partial E_{\rm b}})&\frac{M_{\rm sw}}{H_{\rm
b}}&\beta_{\theta}-\beta_{\theta}\frac{M_{\rm sw}}{H_{\rm b}}\frac{\partial
P_{\rm sw}}{\partial M_{\rm sw}}&-\beta_{\theta}\frac{M_{\rm sw}}{H_{\rm
b}}\frac{\partial P_{\rm sw}}{\partial M_{\rm ej}}\\\
-\beta_{\theta}\frac{M_{\rm ej}}{H_{\rm b}}(1+\frac{\partial P_{\rm
sw}}{\partial E_{\rm b}})&\frac{M_{\rm ej}}{H_{\rm
b}}&-\beta_{\theta}\frac{M_{\rm ej}}{H_{\rm b}}\frac{\partial P_{\rm
sw}}{\partial M_{\rm sw}}&\beta_{\theta}-\beta_{\theta}\frac{M_{\rm
ej}}{H_{\rm b}}\frac{\partial P_{\rm sw}}{\partial M_{\rm ej}}\end{pmatrix}.$
(A1)

The four eigenvalues of the this matrix can be analytically solved:

$\displaystyle\omega_{1,2}$ $\displaystyle=\beta_{\theta}c/R,$ (A2)
$\displaystyle\omega_{3,4}$
$\displaystyle=\frac{\beta_{\theta}c}{R}\left(1-\frac{1}{2}\frac{M_{\rm
ej}}{H_{\rm b}}\frac{\partial P_{\rm sw}}{\partial M_{\rm
ej}}-\frac{1}{2}\frac{M_{\rm sw}}{H_{\rm b}}\frac{\partial P_{\rm
sw}}{\partial M_{\rm sw}}\right)$
$\displaystyle\pm\frac{c}{R}\sqrt{(1-\beta_{\theta}^{2})\left(\frac{\partial
P_{\rm sw}}{\partial E_{\rm b}}+\frac{M_{\rm ej}}{H_{\rm b}}\frac{\partial
P_{\rm sw}}{\partial M_{\rm ej}}+\frac{M_{\rm sw}}{H_{\rm b}}\frac{\partial
P_{\rm sw}}{\partial M_{\rm
sw}}\right)+\frac{1}{4}\beta_{\theta}^{2}\left(\frac{M_{\rm ej}}{H_{\rm
b}}\frac{\partial P_{\rm sw}}{\partial M_{\rm ej}}+\frac{M_{\rm sw}}{H_{\rm
b}}\frac{\partial P_{\rm sw}}{\partial M_{\rm sw}}\right)}.$ (A3)

The partial derivatives can be solved by differentiating eq. 18 and eq. 19,
yielding the following results:

$\displaystyle\frac{\partial P_{\rm sw}}{\partial E_{\rm b}}$
$\displaystyle=\frac{\frac{2}{3}sM_{\rm sw}}{\frac{2}{3}sM_{\rm
sw}(4\gamma^{4}-1)+(1-s)\gamma^{3}M_{\rm sw}+\gamma^{3}M_{\rm ej}},$ (A4)
$\displaystyle\frac{\partial P_{\rm sw}}{\partial M_{\rm sw}}$
$\displaystyle=\frac{1}{3}s\beta^{2}-\frac{\frac{2}{3}sM_{\rm
sw}\left[s\gamma^{2}(1+\frac{1}{3}\beta^{4})+(1-s)\gamma\right]}{\frac{2}{3}sM_{\rm
sw}(4\gamma^{4}-1)+(1-s)\gamma^{3}M_{\rm sw}+\gamma^{3}M_{\rm ej}},$ (A5)
$\displaystyle\frac{\partial P_{\rm sw}}{\partial M_{\rm ej}}$
$\displaystyle=-\frac{\frac{2}{3}s\gamma M_{\rm sw}}{\frac{2}{3}sM_{\rm
sw}(4\gamma^{4}-1)+(1-s)\gamma^{3}M_{\rm sw}+\gamma^{3}M_{\rm ej}}.$ (A6)

The eigenvalues correspond to all the wave speeds present in this system.
These values approach zero as $\beta_{\theta}\to 0$ and $\gamma\to\infty$.
This effect suggests that, in the ultra-relativistic limit, the fluid elements
are initially causally disconnected, with interactions across the blast
surface initially frozen. In the numerical scheme we introduce in the
following section, we do not require all the eigenvalues, but only the maximum
absolute value $\omega_{\rm max}=\max_{\rm j}(|\omega_{j}|)$.

###  A.2 Solving the evolution equations

The evolution equations in our model are hyperbolic conservation laws coupled
with a Hamilton-Jacobi equation. Each of them are solved with their respective
numerical schemes.

To solve the hyperbolic equations, we employ a second-order Godunov-type
finite volume method. The mesh is formed by unequally spaced (but static)
cells arranged along the polar direction, as detailed in the following
section. Up to second-order accuracy, the averages within each cell align with
the values at the cell centers. This alignment is especially beneficial for
our numerical scheme, as the blastwave radius to be evolved are also evaluated
at the cell centers. The numerical scheme can be expressed by a semi-discrete
form:

$\frac{d{\bf U}_{\rm i}}{dt}=\frac{{\bf F}_{\rm i-1/2}\sin\theta_{\rm
i-1/2}-{\bf F}_{\rm i+1/2}\sin\theta_{\rm i+1/2}}{\Delta V_{\rm i}}+{{\bf
S}}_{\rm i}.$ (A7)

where $i$ denotes the center of a cell, $i-1/2$ and $i+1/2$ are the left and
right edges of the cell. The volume of each cell is $V_{\rm i}=\cos\theta_{\rm
i-1/2}-\cos\theta_{\rm i+1/2}$.

For the reconstruction of the primary variables within a cell, we perform a
piecewise linear reconstruction approach with a Minmod slope limiter. The
slopes for the cells adjacent to the poles are reconstructed using ghost
cells, extrapolated based on a reflective boundary condition. The numerical
flux at the cell edges is determined using the Rusanov Riemann solver
(Rusanov, 1962):

${\bf F}_{\rm i+1/2}=\frac{1}{2}\left[{\bf F}_{\rm L}+{\bf F}_{\rm
R}-\omega^{\rm stencil}_{\rm max}({\bf U}_{\rm R}-{\bf U}_{\rm L})\right]$
(A8)

where ${\rm L}$ and ${\rm R}$ denotes the values at the left and right sides
of the edges. The coefficient $\omega^{\rm stencil}_{\rm max}$ is the maximum
$\omega_{\rm max}$ over the four cells around the edge. The time integration
is performed by a second order strong stability preserving Runge-Kutta method
with a time step $\Delta t$.

To maintain numerical stability, the time step must satisfies the following
Courant–Friedrichs–Lewy (CFL) condition

$\Delta t=c_{\rm cfl}\max_{\rm i}\left(\frac{\theta_{\rm i+1/2}-\theta_{\rm
i-1/2}}{\omega^{\rm sig}_{\rm i}}\right)$ (A9)

where $0<c_{\rm cfl}<1$ is the CFL number, $\omega^{\rm sig}_{\rm i}$ is the
signal speed of each cell.

In principle the signal speed $\omega^{\rm sig}_{\rm i}$ should match
$\omega_{\rm max}$. However, due to the aforementioned causality
considerations, we find $\omega_{\rm max}\approx 0$ at very early times. This
leads to a practically infinite time step at the beginning, which blows up the
scheme. To ensure a finite time step, we apply the following correction:

$\omega^{\rm sig}_{\rm i}=\omega_{\rm max,i}+0.05\frac{\beta_{\rm i}c}{R_{\rm
i}}.$ (A10)

In the ultra-relativistic limit the correction term corresponds to 5% speed of
light, and it vanishes in the Newtonian limit. In numerical tests we find this
correction provides an appropriate step size such that in coasting phase the
step size is sufficiently large to maintain good computational efficiency,
while not so large that misses the onsite of lateral expansion.

The radius evolution equation can be regarded as a non-linear Hamilton-Jacobi
equation $\partial R/\partial t+H(\partial R/\partial\theta)=0$, where the
Hamiltonian $H$ is defined by the negative right hand side of eq. 26.
Following Shu 2007, we adopt a second order monotone scheme within the
framework of finite difference method

$\frac{dR}{dt}=\hat{H}(u^{-},u^{+}),$ (A11)

where $\hat{H}$ is the numerical Hamiltonian, $u^{-}$ and $u^{+}$ are the left
and right biased estimation of slope. The numerical Hamiltonian is calculated
using a Lax-Friedrichs flux (Shu & Osher, 1989):

$\hat{H}(u^{-},u^{+})=H(\frac{u^{-}+u^{+}}{2})-\frac{1}{2}\alpha_{\rm
max}(u^{+}-u^{-}),$ (A12)

where $\alpha=|\partial H/\partial u|=|\beta_{\theta}c/R|$. The maximum is
evaluated over a stencil including the cell itself and the two adjacent cells.
The time integration is the same as the hyperbolic equations.

###  A.3 Initial condition and resolution

The initial conditions are defined by an arbitrary initial energy profile
$E_{\rm b,0}$ and a Lorentz factor profile $\gamma_{0}$. Correspondingly, the
ejecta mass profile is determined as $M_{\rm ej,0}=E_{\rm b,0}/\gamma_{0}$. In
a realistic scenario, a collimated jet has negligible energy far from the
core. However, for numerical reasons, it is necessary to avoid 0 energy. This
is achieved by adding an isotropic tail with sufficiently low energy and
speed. To avoid zero radius as well, we initially “coast” the jet from the
origin to a sufficiently small radius at a constant initial speed, and start
the simulation since then. The simulated data is subsequently interpolated for
the calculation of the afterglow.

The resolution of the simulation depends on the initial condition. Cells are
unequally sampled to adapt the energy and Lorentz factor distribution. For
simplicity, we uniformly distribute the cells within the jet core, and
logarithmic outside the core.

##  Appendix B Posterior Probability Distribution of the fitting to GW170817

Figure 12: The posterior probability distribution of the parameters in the
MCMC study without the centroid offset data.  Figure 13: The posterior
probability distribution of the parameters in the MCMC study with the centroid
offset data included.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2017c) —. 2017c, ApJ, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Alexander et al. (2017) Alexander, K. D., Berger, E., Fong, W., et al. 2017, ApJ, 848, L21, doi: 10.3847/2041-8213/aa905d
  * Alexander et al. (2018) Alexander, K. D., Margutti, R., Blanchard, P. K., et al. 2018, ApJ, 863, L18, doi: 10.3847/2041-8213/aad637
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Balasubramanian et al. (2021) Balasubramanian, A., Corsi, A., Mooley, K. P., et al. 2021, ApJ, 914, L20, doi: 10.3847/2041-8213/abfd38
  * Balasubramanian et al. (2022) —. 2022, ApJ, 938, 12, doi: 10.3847/1538-4357/ac9133
  * Beniamini et al. (2022) Beniamini, P., Gill, R., & Granot, J. 2022, MNRAS, 515, 555, doi: 10.1093/mnras/stac1821
  * Beniamini et al. (2020) Beniamini, P., Granot, J., & Gill, R. 2020, MNRAS, 493, 3521, doi: 10.1093/mnras/staa538
  * Beniamini & Nakar (2019) Beniamini, P., & Nakar, E. 2019, MNRAS, 482, 5430, doi: 10.1093/mnras/sty3110
  * Beniamini et al. (2019) Beniamini, P., Petropoulou, M., Barniol Duran, R., & Giannios, D. 2019, MNRAS, 483, 840, doi: 10.1093/mnras/sty3093
  * Blandford & McKee (1976) Blandford, R. D., & McKee, C. F. 1976, Physics of Fluids, 19, 1130, doi: 10.1063/1.861619
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Bulla et al. (2022) Bulla, M., Coughlin, M. W., Dhawan, S., & Dietrich, T. 2022, Universe, 8, 289, doi: 10.3390/universe8050289
  * Di Valentino et al. (2021) Di Valentino, E., Mena, O., Pan, S., et al. 2021, Classical and Quantum Gravity, 38, 153001, doi: 10.1088/1361-6382/ac086d
  * Duffell & Laskar (2018) Duffell, P. C., & Laskar, T. 2018, ApJ, 865, 94, doi: 10.3847/1538-4357/aadb9c
  * Duffell & MacFadyen (2013) Duffell, P. C., & MacFadyen, A. I. 2013, ApJ, 775, 87, doi: 10.1088/0004-637X/775/2/87
  * Duque et al. (2019) Duque, R., Daigne, F., & Mochkovitch, R. 2019, A&A, 631, A39, doi: 10.1051/0004-6361/201935926
  * Fernández et al. (2022) Fernández, J. J., Kobayashi, S., & Lamb, G. P. 2022, MNRAS, 509, 395, doi: 10.1093/mnras/stab2879
  * Finstad et al. (2018) Finstad, D., De, S., Brown, D. A., Berger, E., & Biwer, C. M. 2018, ApJ, 860, L2, doi: 10.3847/2041-8213/aac6c1
  * Fong et al. (2019) Fong, W., Blanchard, P. K., Alexander, K. D., et al. 2019, ApJ, 883, L1, doi: 10.3847/2041-8213/ab3d9e
  * Ghirlanda et al. (2019) Ghirlanda, G., Salafia, O. S., Paragi, Z., et al. 2019, Science, 363, 968, doi: 10.1126/science.aau8815
  * Gianfagna et al. (2023) Gianfagna, G., Piro, L., Pannarale, F., et al. 2023, arXiv e-prints, arXiv:2309.17073, doi: 10.48550/arXiv.2309.17073
  * Gill & Granot (2018) Gill, R., & Granot, J. 2018, MNRAS, 478, 4128, doi: 10.1093/mnras/sty1214
  * Govreen-Segal & Nakar (2023) Govreen-Segal, T., & Nakar, E. 2023, arXiv e-prints, arXiv:2311.09297, doi: 10.48550/arXiv.2311.09297
  * Granot et al. (2018) Granot, J., De Colle, F., & Ramirez-Ruiz, E. 2018, MNRAS, 481, 2711, doi: 10.1093/mnras/sty2454
  * Granot & Piran (2012) Granot, J., & Piran, T. 2012, MNRAS, 421, 570, doi: 10.1111/j.1365-2966.2011.20335.x
  * Guidorzi et al. (2017) Guidorzi, C., Margutti, R., Brout, D., et al. 2017, ApJ, 851, L36, doi: 10.3847/2041-8213/aaa009
  * Haggard et al. (2017) Haggard, D., Nynka, M., Ruan, J. J., et al. 2017, ApJ, 848, L25, doi: 10.3847/2041-8213/aa8ede
  * Hajela et al. (2019) Hajela, A., Margutti, R., Alexander, K. D., et al. 2019, ApJ, 886, L17, doi: 10.3847/2041-8213/ab5226
  * Hotokezaka et al. (2018) Hotokezaka, K., Kiuchi, K., Shibata, M., Nakar, E., & Piran, T. 2018, ApJ, 867, 95, doi: 10.3847/1538-4357/aadf92
  * Hotokezaka et al. (2019) Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, Nature Astronomy, 3, 940, doi: 10.1038/s41550-019-0820-1
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kamm & Timmes (2007) Kamm, J., & Timmes, F. 2007, Technical Report LA-UR-07-2849
  * Kathirgamaraju et al. (2018) Kathirgamaraju, A., Barniol Duran, R., & Giannios, D. 2018, MNRAS, 473, L121, doi: 10.1093/mnrasl/slx175
  * Kathirgamaraju et al. (2019) Kathirgamaraju, A., Giannios, D., & Beniamini, P. 2019, MNRAS, 487, 3914, doi: 10.1093/mnras/stz1564
  * Kumar & Granot (2003) Kumar, P., & Granot, J. 2003, ApJ, 591, 1075, doi: 10.1086/375186
  * Kumar & Zhang (2015) Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1, doi: 10.1016/j.physrep.2014.09.008
  * Lamb & Kobayashi (2017) Lamb, G. P., & Kobayashi, S. 2017, MNRAS, 472, 4953, doi: 10.1093/mnras/stx2345
  * Lamb et al. (2018) Lamb, G. P., Mandel, I., & Resmi, L. 2018, MNRAS, 481, 2581, doi: 10.1093/mnras/sty2196
  * Lamb et al. (2019) Lamb, G. P., Lyman, J. D., Levan, A. J., et al. 2019, ApJ, 870, L15, doi: 10.3847/2041-8213/aaf96b
  * Lazzati et al. (2018) Lazzati, D., Perna, R., Morsony, B. J., et al. 2018, Phys. Rev. Lett., 120, 241103, doi: 10.1103/PhysRevLett.120.241103
  * Li et al. (2024) Li, J.-N., Wang, Y.-Y., Wang, Y., et al. 2024, arXiv e-prints, arXiv:2401.17978, doi: 10.48550/arXiv.2401.17978
  * Lu et al. (2020) Lu, W., Beniamini, P., & McDowell, A. 2020, arXiv e-prints, arXiv:2005.10313, doi: 10.48550/arXiv.2005.10313
  * Makhathini et al. (2021) Makhathini, S., Mooley, K. P., Brightman, M., et al. 2021, ApJ, 922, 154, doi: 10.3847/1538-4357/ac1ffc
  * Margutti & Chornock (2021) Margutti, R., & Chornock, R. 2021, ARA&A, 59, 155, doi: 10.1146/annurev-astro-112420-030742
  * Margutti et al. (2017) Margutti, R., Berger, E., Fong, W., et al. 2017, ApJ, 848, L20, doi: 10.3847/2041-8213/aa9057
  * Margutti et al. (2018) Margutti, R., Alexander, K. D., Xie, X., et al. 2018, ApJ, 856, L18, doi: 10.3847/2041-8213/aab2ad
  * Mooley et al. (2022) Mooley, K. P., Anderson, J., & Lu, W. 2022, Nature, 610, 273, doi: 10.1038/s41586-022-05145-7
  * Mooley et al. (2018) Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018, Nature, 561, 355, doi: 10.1038/s41586-018-0486-3
  * Nakar & Piran (2011) Nakar, E., & Piran, T. 2011, Nature, 478, 82, doi: 10.1038/nature10365
  * Nakar & Piran (2018) —. 2018, MNRAS, 478, 407, doi: 10.1093/mnras/sty952
  * Nedora et al. (2023a) Nedora, V., Dietrich, T., & Shibata, M. 2023a, MNRAS, 524, 5514, doi: 10.1093/mnras/stad2128
  * Nedora et al. (2023b) Nedora, V., Dietrich, T., Shibata, M., Pohl, M., & Crosato Menegazzi, L. 2023b, MNRAS, 520, 2727, doi: 10.1093/mnras/stad175
  * O’Connor & Troja (2022) O’Connor, B., & Troja, E. 2022, GRB Coordinates Network, 32065, 1
  * Piran et al. (2013) Piran, T., Nakar, E., & Rosswog, S. 2013, MNRAS, 430, 2121, doi: 10.1093/mnras/stt037
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Rhoads (1999) Rhoads, J. E. 1999, ApJ, 525, 737, doi: 10.1086/307907
  * Riess et al. (2022) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7, doi: 10.3847/2041-8213/ac5c5b
  * Rusanov (1962) Rusanov, V. 1962, USSR Computational Mathematics and Mathematical Physics, 1, 304, doi: https://doi.org/10.1016/0041-5553(62)90062-9
  * Ryan et al. (2020) Ryan, G., van Eerten, H., Piro, L., & Troja, E. 2020, ApJ, 896, 166, doi: 10.3847/1538-4357/ab93cf
  * Ryan et al. (2023) Ryan, G., van Eerten, H., Troja, E., et al. 2023, arXiv e-prints, arXiv:2310.02328, doi: 10.48550/arXiv.2310.02328
  * Saleem et al. (2018) Saleem, M., Resmi, L., Misra, K., Pai, A., & Arun, K. G. 2018, MNRAS, 474, 5340, doi: 10.1093/mnras/stx3104
  * Sari et al. (1999) Sari, R., Piran, T., & Halpern, J. P. 1999, ApJ, 519, L17, doi: 10.1086/312109
  * Sari et al. (1998) Sari, R., Piran, T., & Narayan, R. 1998, ApJ, 497, L17, doi: 10.1086/311269
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics
  * Shu (2007) Shu, C.-W. 2007, High order numerical methods for time dependent Hamilton-Jacobi equations (World Scientific), 47–91, doi: 10.1142/9789812709066_0002
  * Shu & Osher (1989) Shu, C.-W., & Osher, S. 1989, Journal of Computational Physics, 83, 32, doi: 10.1016/0021-9991(89)90222-2
  * Sironi & Giannios (2013) Sironi, L., & Giannios, D. 2013, ApJ, 778, 107, doi: 10.1088/0004-637X/778/2/107
  * Troja et al. (2017) Troja, E., Piro, L., van Eerten, H., et al. 2017, Nature, 551, 71, doi: 10.1038/nature24290
  * Troja et al. (2018) Troja, E., Piro, L., Ryan, G., et al. 2018, MNRAS, 478, L18, doi: 10.1093/mnrasl/sly061
  * Troja et al. (2020) Troja, E., van Eerten, H., Zhang, B., et al. 2020, MNRAS, 498, 5643, doi: 10.1093/mnras/staa2626
  * Uhm (2011) Uhm, Z. L. 2011, ApJ, 733, 86, doi: 10.1088/0004-637X/733/2/86
  * van Eerten et al. (2012) van Eerten, H., van der Horst, A., & MacFadyen, A. 2012, ApJ, 749, 44, doi: 10.1088/0004-637X/749/1/44
  * Wang & Giannios (2021) Wang, H., & Giannios, D. 2021, ApJ, 908, 200, doi: 10.3847/1538-4357/abd39c
  * Wang et al. (2023) Wang, Y.-Y., Tang, S.-P., Jin, Z.-P., & Fan, Y.-Z. 2023, ApJ, 943, 13, doi: 10.3847/1538-4357/aca96c
  * Wu & MacFadyen (2018) Wu, Y., & MacFadyen, A. 2018, ApJ, 869, 55, doi: 10.3847/1538-4357/aae9de
  * Wu & MacFadyen (2019) —. 2019, ApJ, 880, L23, doi: 10.3847/2041-8213/ab2fd4
  * Zrake et al. (2018) Zrake, J., Xie, X., & MacFadyen, A. 2018, ApJ, 865, L2, doi: 10.3847/2041-8213/aaddf8