GAMMA-RAY BURST AFTERGLOW BROADBAND FITTING BASED DIRECTLY ON HYDRODYNAMICS SIMULATIONS

Gamma-ray bursts (GRBs) are short intense flashes of gamma radiation produced by cataclysmic stellar events such as the collapse of the core of a massive star (Woosley 1993; Paczynski 1998; MacFadyen & Woosley 1999) or a neutron star–neutron star or neutron star–black hole merger (e.g., Eichler et al. 1989; Paczynski 1991). During these events a collimated relativistic outflow is produced that sweeps up the matter surrounding the GRB. Regardless of the original mass content or launching mechanism of this outflow (be it a fireball, Mészáros & Rees 1997, or a Poynting-flux-dominated jet, e.g., Drenkhahn 2002), the expanding blast wave will sweep up circumburst matter and will eventually start to decelerate. As the blast wave shocks the circumburst medium, broadband synchrotron radiation is produced by shock-accelerated electrons, giving rise to an afterglow signal that can be observed for up to days at X-ray and optical frequencies, and for up to years at radio frequencies.

Three kinds of parameters determine the shape of the observed afterglow light curves. First, the shock dynamics are set by the explosion energy and circumburst density. A second set of parameters captures the physics of synchrotron emission from shock-accelerated electrons. Finally, the observed flux depends on the parameters defined by a given observation: frequency, time, and observer angle.

Ever since the first afterglows were discovered (Costa et al. 1997; Groot et al. 1997), models based on synchrotron radiation from a decelerating relativistic blast wave have been successful in describing the broadband data (e.g., Wijers et al. 1997; Sari et al. 1998; Wijers & Galama 1999; Granot & Sari 2002). The synchrotron spectrum is typically described as consisting of a number of connected power-law regimes, with the critical frequencies connecting the regimes shifting over time and being determined by the basic spectral shape of synchrotron emission, synchrotron self-absorption, and cooling of the shock-accelerated electrons. In order to accurately model the late-time afterglow emission in the radio, afterglow models based on a non-relativistic rather than relativistic blast wave have been applied as well (e.g., Frail et al. 2000). Both the early-time ultrarelativistic and late-time non-relativistic stages of the evolution of the shock are self-similar, with solutions given by Blandford & McKee (1976), hereafter BM, and Sedov (1959), Von Neumann (1961), and Taylor (1950), hereafter ST, respectively. The intermediate stage can be approximated (e.g., Rhoads 1999; Huang et al. 1999), but this stage is complicated by the dynamics of jet decollimation. When jet decollimation is taken into account, a homogeneous jet surface that widens with the comoving speed of sound is often assumed (Rhoads 1999), while more recent jet spreading models (Granot & Piran 2012) do not take the radial structure of the outflow into account. The jet nature of the outflow ultimately reveals itself as a break in the light curve, the jet break, and a subsequent steepening of the light-curve slope. The physics of afterglow jets has been reviewed in, e.g., Piran (2005), Mészáros (2006), and Granot (2007).

Purely analytical models are severely limited in that they do not accurately capture many features of the jet dynamics (such as the aforementioned jet spreading and deceleration) and radiation. The simplifications inherent in a purely analytical approach lead to diverging predictions for a range of features such as the observed shape of the jet break, the size and shape of the counterjet (which was launched away from the observer and is only seen at late times, when relativistic beaming of the emitted radiation plays a lessened role), the nature and duration of the transition to the non-relativistic phase, and the effect of the orientation of the jet with respect to the observer. To gain a better understanding of these aspects, various authors have performed numerical relativistic hydrodynamics (RHD) simulations of afterglow jets, in one dimension (Kobayashi et al. 1999; Downes et al. 2002; Mimica et al. 2009; Van Eerten et al. 2010b), two dimensions (Granot et al. 2001; Zhang & MacFadyen 2006; Meliani et al. 2007; Ramirez-Ruiz & MacFadyen 2010; Van Eerten et al. 2011; Wygoda et al. 2011), and occasionally even three (Cannizzo et al. 2004). Over the past decade, these simulations have steadily increased in accuracy, mainly through the use of adaptive-mesh refinement (AMR) techniques, which locally increase resolution during simulation where needed, and are capable of resolving the six orders of magnitude difference between the initial width of the thin relativistic shell and the late-time outer radius of the decelerated jet. Recent high-resolution simulations have shown that relativistic jets spread sideways significantly slower than analytically expected and have a strongly inhomogeneous shock front (Zhang & MacFadyen 2009; Meliani & Keppens 2010; Van Eerten et al. 2011; Van Eerten & MacFadyen 2011a; Wygoda et al. 2011). Furthermore, they show that the transition from relativistic to non-relativistic expansion is a very slow process (Zhang & MacFadyen 2009; Van Eerten et al. 2010b) and that jet orientation strongly affects the jet break even for small observer angles (Van Eerten et al. 2010a), while for observers both on and off axes the observed jet-break time differs between different frequencies due to synchrotron self-absorption (Van Eerten et al. 2011).

An RHD jet simulation can be combined with a numerical synchrotron radiation calculation to yield a powerful tool to predict the evolution of the observable broadband afterglow spectrum in detail. The weaknesses of simplified analytical models are thus avoided, and local changes in fluid structure and arrival time effects are correctly accounted for. This calculation can be performed in different ways, for example, by summing over the emitted power of all fluid cells in the simulation in the case of an optically thin fluid (Downes et al. 2002; Nakar & Granot 2007; Zhang & MacFadyen 2009) or by fully solving the linear radiative transfer equations including synchrotron self-absorption (Van Eerten et al. 2010b; Van Eerten & MacFadyen 2011b).

An obvious drawback of simulation-based light curves compared to analytically calculated light curves is that calculating the former is a time-consuming process. A full jet simulation takes several thousand CPU-hours to complete. A purely analytical light curve, on the other hand, can be calculated almost instantaneously and can therefore be applied to iterative model fitting, where the procedure of minimizing χ² requires at least thousands of light-curve calculations with slightly differing explosion and radiation parameters. In this paper, we present a new method to use simulation results directly as a basis for iterative fitting of broadband data, which closes the gap between simulations and analytical models. This method should prove useful for further constraining the physics of GRB afterglows, thereby obtaining clues about the nature of the progenitor and the burst environment. This provides more accurate predictions for future surveys, including LOFAR (Rottgering 2006), SKA (Carilli & Rawlings 2004), or ALMA (Wootten 2003), and indirectly benefits gravitational wave predictions for LIGO (Abbott et al. 2009) and VIRGO (Acernese et al. 2008), where GRBs are potentially observable as electromagnetic counterparts (Nakar & Piran 2011; Van Eerten & MacFadyen 2011b). Finally, our method helps to establish a baseline for studies of the effect of more detailed models of the microphysics (see, e.g., Van Eerten et al. 2010b; Panaitescu et al. 2006; Filgas et al. 2011).

This paper is structured as follows. First, we briefly describe the numerical settings and code used for our RHD jet simulations of jets expanding into a homogeneous circumburst medium in Section 2. Our approach is made possible by two properties of decollimating and decelerating relativistic jets starting from a BM solution. First, the jet evolution is scale-invariant under rescaling of both explosion energy and circumburst density, which we discuss in Section 3. Second, for a given initial opening angle, the two-dimensional fluid profile of the blast wave evolves smoothly from a relativistic and purely radial outflow to a non-relativistic and spherical outflow. This implies that both the radial and lateral structure of the flow can be captured with sufficient accuracy by a low-resolution grid with specialized coordinates that can be determined a posteriori once the radial and angular extent of the jet at each moment in lab frame time are known from the high-resolution simulation. The dynamical evolution and condensed low-resolution description of jets with various opening angles are discussed in Section 4. In Section 5, we describe how the simulation results have been implemented in a broadband fitting code, which we apply to a case study, i.e., GRB 990510, in Section 6. We discuss our results in Section 7.

The source code of the broadband fit code will be made freely available on request or for download from http://cosmo.nyu.edu/afterglowlibrary. It can be run both on a single core and in parallel and allows the user either to quickly generate light curves and spectra for arbitrary explosion parameters or, when provided with a data set of observed fluxes in mJy, to perform a full broadband fit.

For this study a total of 19 jet simulations in two dimensions have been performed using the relativistic adaptive mesh (RAM) parallel RHD code (Zhang & MacFadyen 2006). The code employs the fifth-order weighted essentially non-oscillatory scheme (Jiang & Shu 1996) and uses the PARAMESH AMR tools (MacNeice et al. 2000) from FLASH 2.3 (Fryxell et al. 2000). For all jet simulations the BM solution for an adiabatic impulsive explosion is used in spherical coordinates to set the initial conditions. Instead of the full spherical solution, a conic section is used that is truncated at a different fixed opening angle for each simulation. The opening angles are listed in Table 1, along with the jet energy E_j for each simulation.

The jet energy E_j (the total for both jets) relates to the isotropic equivalent energy E_iso according to

$$\begin{equation} E_j = E_{{\rm iso}} (1 - \cos \theta _0) \approx E_{{\rm iso}} \theta _0^2 / 2. \end{equation} \tag{ 1 }$$

All jets expand into a homogeneous medium with number density n₀ = 1 cm⁻³ (mass density ρ₀ = 1 × m_p g cm⁻³, in terms of the proton mass m_p) and have an isotropic equivalent explosion energy E_iso = 6.25 × 10⁵¹ erg; note that due to the scale invariance of the simulations with respect to n₀ and E_iso, these can be scaled afterward to represent arbitrary values (see Section 3). All jets start at time t_b with fluid Lorentz factor γ_b = 25 directly behind the shock, ensuring that for all simulations γ_b > 1/θ₀. At this point the edges of the jet have not yet come into causal contact and lateral spreading has not yet set in, which allows us to use the spherically symmetric BM solution as the starting point.

The initial outer radii of the jets are given by R_b = 1.3102 × 10¹⁷ cm, determined from t_b and γ_b through R_b = ct_b(1 − 1/16γ²_b) (Equation (26) from BM, with c being the speed of light). The initial time t_b is determined from γ_b, E_iso, and n₀ using

$$\begin{equation} 17 E_{{\rm iso}} = 16 \pi m_p n_0 \gamma _b^2 c^5 t_b^3, \end{equation} \tag{ 2 }$$

which expresses conservation of the total energy (Equation (43) from BM).

For all simulations the stopping time is determined according to

$$\begin{equation} t_f = 10 \times t_{\rm NR} = 9700 \left(\frac{E_{{\rm iso}}}{10^{53} n_0} \right)^{1/3} {\rm days}. \end{equation} \tag{ 3 }$$

The time t_NR marks the time when the jet is analytically expected to transition from relativistic to non-relativistic flow and is determined from comparing the initial explosion energy to the total rest mass energy of the swept-up matter (Piran 2005). Numerical simulations have shown that in practice this transition takes more time to complete, which is why we have chosen to continue our simulations until two times the transition duration 5 × t_NR found numerically by Zhang & MacFadyen (2009). For E_iso = 6.25 × 10⁵¹ erg and γ_b = 25, it follows for all 19 simulations that t_b = 4.37 × 10⁶ s =50.6 days and t_f = 3.33 × 10⁸ s =3849 days; again, these values change under rescaling of E_iso or n₀.

All simulations use an equation of state with the adiabatic index as a function of comoving density and pressure changing smoothly from 4/3 for a relativistic fluid to 5/3 for a non-relativistic fluid (Zhang & MacFadyen 2009; Mignone et al. 2005).

2.1. Resolution and Refinement

What is straightforward, but perhaps obscured by the richness of features in the resulting light curves, and what has so far not been utilized in afterglow modeling, is that the full two-dimensional evolution of the jet is invariant under scaling of the initial explosion energy and circumburst medium density: independent of self-similarity, a more energetic jet (or a jet in a less dense environment) goes through exactly the same evolutionary stages as a jet with lower energy (or higher circumburst density), albeit that each stage occurs at a later time and larger radius.

The scalings can be understood as follows. The initial setup of the problem is determined completely by a limited number of parameters: E_iso, ρ₀ (≡ n₀ × m_p, with m_p the proton mass), θ₀, and c. The initial Lorentz factor γ_b is not included in the list because its precise value is arbitrary as long as γ_b > 1/θ₀. Only a limited number of independent dimensionless combinations of these parameters are possible, such as

$$\begin{equation} A = \frac{r}{ct}, \quad B = \frac{E_{{\rm iso}} t^2}{\rho _0 r^5}, \quad \theta, \quad \theta _0. \end{equation} \tag{ 4 }$$

Any dimensionless quantity that describes the local fluid conditions or global evolution of the jet, such as Lorentz factor γ, density ρ/ρ₀, or θ₉₅ (the angle with respect to the jet axis within which 95% of the jet energy is contained), can be expressed as a function of these parameters, i.e., γ(r, t, θ, E_iso, ρ₀, θ₀) → γ(A, B, θ, θ₀). Note that at both very early times (no lateral flow) and very late times (spherical flow) there is no dependency on θ or θ₀, which, together with the fact that A becomes a constant value both in the ultrarelativistic BM and non-relativistic ST limits (i.e., A → 1 and A → 0, respectively), accounts for the self-similarity of these limiting cases, with γ(r, t, θ) → γ(B), etc. The parameter B is identical to the ST self-similarity variable ξ⁻⁵ and is central to the self-similarity of the BM solution via Equation (2) above (where r ∼ ct), which fixes γ² and therefore the BM self-similarity variable χ(B, γ²(B)).

Even if we do not limit ourselves to either of the self-similar extremes and leave A and θ in place, we have a set of coordinates A, B, θ that is invariant under the following rescalings:

$$\begin{eqnarray} E^{\prime }_{{\rm iso}} & = & \kappa E_{{\rm iso}}, \nonumber \\ \rho _0^{\prime } & = & \lambda \rho _0, \nonumber \\ r^{\prime } & = & (\kappa / \lambda)^{1/3} r, \nonumber \\ t^{\prime } & = & (\kappa / \lambda)^{1/3} t. \end{eqnarray} \tag{ 5 }$$

The full scale-invariant hydrodynamics equations in terms of A, B, θ are provided for completeness in Appendix B. The direct implication of this is that we can determine the value of any dimensionless quantity for an explosion with E'_iso and ρ'₀ simply by probing a simulation with parameters E_iso, ρ at t, r, rather than at t', r'. Quantities that are not dimensionless, such as density ρ and internal energy density e, follow from ρ/ρ₀ = ρ'/ρ'₀, e/ρ₀c² = e'/ρ'₀c², etc.

Examples of scalings between jets with different E_iso and ρ₀ are given in Figures 1 and 2. The comoving fluid density and energy density profiles are drawn from the simulations presented in Van Eerten & MacFadyen (2011b), since the 19 simulations listed in Table 1 were set up to be unrelated via scaling. In the case of Figure 1, two simulations are compared for which λ = 1, whereas λ = 10³ in the case of Figure 2. In Figure 1, we set the outer radius of the left and right panels equal to their respective ct, resulting in both images being complete mirror images of each other, which confirms that both simulations are numerically indistinguishable. Although analytically equivalent, the two simulation runs in Figure 2 are no longer numerically identical, which leads to minor differences in (de-)refinement in the inner regions. However, those inner regions do not contribute to the observed flux and thus have no observable effect on the light curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have plotted a number of features of a subset of the 19 jet simulations in Figures 3 and 4. The jets follow the same evolution as those of earlier studies, such as Zhang & MacFadyen (2009), Van Eerten & MacFadyen (2011a), and Wygoda et al. (2011). The top plot of Figure 3 shows the evolution of the on-axis outer radius of the blast wave for jets with different θ₀ in the lab frame of the explosion. All jet simulations start out relativistically with identical γ_b, t_b, and blast wave radius R_b. Jets with a wide opening angle do not have a long spreading phase and undergo a single smooth transition from R ≈ ct to R ≈ 1.15(E_jt²/ρ₀)^1/5 (where the numerical factor 1.15 follows from the ST solution with adiabatic index 5/3). The figure shows that for a very wide jet with θ₀ = 0.5, the transition time is well approximated by the crossing point of the asymptotes for spherical outflow. However, this does not carry over to smaller angle jets, and the meeting point of the asymptotes for E_j(θ₀ = 0.05) severely underestimates the turnover point. The reason for this is that for narrower jets there is also an intermediate phase, where the jet decelerates due to lateral spreading. Although not as abrupt as originally predicted (Rhoads 1999, also discussed in Van Eerten & MacFadyen 2011a; Wygoda et al. 2011) and occurring throughout the transrelativistic phase of the jet, this leads to an extended period of jet deceleration in excess of the asymptotic jet deceleration of the ST phase and adds a second turnover point in the evolution of jet radius and jet velocity. As the bottom plot of Figure 3 indicates, this intermediate phase does not follow a simple power law. A full parameterization of the intermediate phase lies beyond the scope of this work and will not be required for model fitting based on the simulation results. In addition, due to inhomogeneity along the shock front, deceleration will be different for outflows along different angles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The early-time behavior visible in the bottom plot of Figure 3, where the peak Lorentz factor initially drops below its expected value in the BM regime but then moves back to the BM asymptote, is due to the resolution of the simulations. In the BM solution, the blast wave is extremely thin (ΔR_b ∼ R_b/2γ²_b) and we have not been able to achieve full convergence in two dimensions at early times (the issue is identical to that illustrated in Figure 5 of Zhang & MacFadyen 2009). By using the integrated values of the fluid quantities across a single fluid cell rather than the values at its central coordinate, we have ensured that all energy of the BM solution is accounted for during the initialization of the simulation. For this reason, the drop is only temporary. We emphasize that this is not due to lateral spreading of the jet, for if that were the case, the peak Lorentz factor would not have been able to recover.

Lateral spreading and the inhomogeneity of the shock front are illustrated in Figure 4. The plots show the time evolution of various characteristic angles θ₉₅, θ₇₅, θ₅₀ of the outflow, defined as the angles within which a fraction 0.95, 0.75, 0.50 of the volume-integrated rest-frame energy density τ is contained. The green dashed lines indicate the analytically expected onset of lateral spreading, when γ ∼ 1/θ₀. Because the plots show θ₉₅, ... rather than the outer edge angle θ_edge, these lines have been shifted according to $\theta (t) \rightarrow \sqrt{0.95} \theta _{{\rm edge}}$ (or $\sqrt{0.75} \theta _{{\rm edge}}$ or $\sqrt{0.50} \theta _{{\rm edge}}$). This means that if the simulations had followed the analytical estimate, the turnovers for all angles and fractions would have occurred along the green lines. Because the jets become strongly inhomogeneous along the shock front, the turnover is delayed for characteristic angles bounding smaller energy fractions: the jet energy stays concentrated near the jet axis even after the edges have begun to expand. Because narrower jets have smaller values of E_j, they will decelerate earlier, which leads to the behavior shown in the plots where θ₉₅, ... for small jets cross that of large jets with the same E_iso. Jets with the same E_j do not show this effect (the characteristic angle curves for wider jets would shift sufficiently far to the left in the plots if their energies were downscaled to match the E_j of the narrowest jet). For each figure, all curves at late times tend to the same fixed fractions of π/2 regardless of initial energy (i.e., the jet becomes spherical and homogeneous along the shock front). However, in our simulations these limits are only actually reached by the high-energy fractions. The early-time turnover behavior confirms the validity of our choice of starting γ_b > θ₀ (note that for the wide jets γ_b ≫ θ₀).

Each of the 19 simulations generates 3000 snapshots, varying in data size between ∼350 MB at early times and ∼40 MB at late times (when lower peak refinement levels are utilized). Therefore, it is currently not practically possible to load the complete output for a single simulation into computer memory at once, let alone the combined output of all 19 simulations. However, if we wish to use the simulations as a basis for iterative model fitting, we will need to be able to quickly access the fluid state at any requested point in time and space. We therefore need to summarize the simulation results in a way that adequately captures all aspects of the outflow but occupies only a relatively small amount of computer memory. In this study we have aimed for <1 GB in total, but the desired size in memory will be hardware-dependent.

The next three subsections describe the further technical aspects of producing light curves from the summarized "box," rather than the "grid" that originally contains the simulation data. In Section 5.4 we describe the details of the fit method.

5.1. "Box" Summary of Simulations

The phase space for each fluid variable is set by the variables E_iso, n₀, θ₀, r, θ, and t, and the local fluid state is fully specified by the fluid variables τ, ρ, v_r (radial fluid velocity), v_θ (angular fluid velocity), from which the other fluid quantities such as p (pressure), e, or γ can be derived. By storing e explicitly along with ρ, v_r, v_θ, while still including τ, even though the latter is not needed for the radiative transfer calculation, we ensure that we can explore the full fluid state at each point in time and space without having to use the equation of state. This implies that we will have to store five fluid variables along with the coordinates and sizes of the fluid cells they refer to (i.e., nine quantities in total). The scale invariance discussed in Section 3 implies that we only need to store a single value for E_iso and n₀. As mentioned already, we have used 19 possible values for θ₀. We calculate light curves for intermediate values of θ₀ by interpolating the appropriate results from the 19 simulations, as we demonstrate in Section 5.3. For the r, θ coordinates we use 100 entries each. For t we also use 100 entries rather than 3000, and as with the small subset of θ₀, we demonstrate in Section 5.3 that this number is sufficient and that the full light curve can be calculated using interpolation. All together, we now have the following number of entries in our summary of the fluid evolution: E_iso × n₀ × θ₀ × r × θ × t × variables = 1 × 1 × 19 × 100 × 100 × 100 × 9 = 171,000,000. We will refer to this summary as the "box," in order to distinguish it from grid-based fluid profiles from the simulations.

The reason that we only need as few as 100 cells in the r and θ directions has already been demonstrated partially by Figures 3 and 4 in the preceding section. These figures illustrate that, although high-resolution simulations were required in order to accurately determine the blast wave radius and lateral extent at each point in time, the evolution from BM to ST profile itself is smooth. This means that, once the large-scale properties of the outflow are known, it is possible to use this knowledge to select at each point in time new coordinates such that the key features of the outflow are resolved while ignoring parts of the outflow.

The box θ at each point in time will not cover the entire grid but only runs from 0 to θ_MAX ≡ θ₉₉/0.99 (i.e., slightly over the outer angle within which 99% of the volume integrated τ is contained). Anything outside these angles either is ISM or contributes a negligible amount of radiation. The lateral cells within this region are evenly spread. Alternatively we could have set the cell coordinates using θ₀₀, θ₀₁, etc., but there is no significant difference in the resulting light curves. The radial domain runs from 0 to the outer limit R_MAX of the blast wave at each box value of θ. Because the unshocked ISM is at rest, the outer boundary of the shock wave is readily determined using the criterion that γ > 1.000001. Even for extremely high resolution, this point will not exactly coincide with the peak of the blast wave. We therefore devote 10 cells to the region between R determined by the τ peak of the blast wave and the first shocked ISM cell. Of the remaining 90 radial box cells, 80 are devoted to resolving the blast wave. Since we know from the BM solution that the blast wave width is ΔR ≈ R/12γ² initially and from the ST solution that ΔR∝R eventually, we analytically determine the width of the blast wave by ΔR ≡ R/12γ², with γ² ≡ γ_b(t/t_b)^−3/2 + 1 (where the "1" has been added relative to BM Equation (24) to obtain the correct asymptotic behavior). Note that this is only approximately correct on-axis due to two-dimensional spreading and less accurate for the radial profile along the outer angles. This does not matter, however, as long as the approximation is sufficient to resolve the sharp feature of the blast wave. The final 10 radial cells are spaced between the origin and the back of the shock. All radial cells are equally spaced within their respective region. Figure 5 shows radial fluid profiles for the θ₀ = 0.05 simulation for different times and angles. Profiles for other values of θ₀ are similar.

Zoom In Zoom Out Reset image size

5.2. "Box" Interpolation

5.3. Light Curves

We calculate light curves and spectra from the box using the same method that we have used previously for grid-based light-curve calculations (Van Eerten et al. 2010a; Van Eerten & MacFadyen 2011b). The dominant radiation mechanism is assumed to be synchrotron radiation, and the broadband emission from each fluid cell is given by a series of connected power laws similar to those in Sari et al. (1998). The linear radiative transfer equations are solved simultaneously for a large number of rays. For each point in lab frame time t, the plane perpendicular to the direction of the observer and at a fixed distance R_EDS from the origin of the box (or grid), defined by R_EDS/c = t − t_obs, defines the area from which emission will arrive at exactly the same observer time t_obs. This plane is labeled the equidistant surface (EDS; see also Van Eerten et al. 2010b). In earlier work we have employed a procedure analogous to AMR for dynamically changing the number of rays through the EDSs that are followed simultaneously. For the blast waves of this study all EDS refinement would have occurred near the center of the EDS (defined as the intersection of the line from the grid origin to the observer and the EDS), in order to resolve the early-time BM profile, and the refinement level would have gradually decreased outward. Since this is essentially equivalent to base 2 logarithmic spacing, for the current study we use fixed logarithmic spacing between rays in the radial direction instead. The number of rays is evenly spaced in the angular direction. This approach is somewhat faster than the dynamical refining and requires less memory for bookkeeping. In the final step the rays are integrated over to yield the observed flux. Flux, observed frequency, and observer time are corrected for redshift z during the calculation.

In Appendix A, we give the exact expressions for the emission and absorption coefficients that are calculated while solving the radiative transfer equations through the evolving fluid. Their local values depend on a number of parameters that capture the microphysics behind the synchrotron radiative process as well as on the local fluid conditions. There are four such parameters: the power-law slope p of the shock-accelerated electrons, the fraction _B of magnetic energy relative to thermal energy, the fraction _e of downstream thermal energy density in the accelerated electrons, and the fraction ξ_N of the downstream particle number density that participates in the shock-acceleration process. By performing broadband fits on afterglow data using the box-based fit method from this paper, these parameters can be determined from the fits for the first time using the full blast wave evolution.

All simulations start at γ_b = 25. Before this time the outflow can be described by a conic section of the spherically symmetric BM solution. Because the observed flux at a given observer time contains emission from a wide range of emission times, initially including times for which γ > γ_b, the BM solution is used directly to determine the initial emission and absorption coefficients. Probing the BM solution between γ = 200 and γ_b using 1600 logarithmically spaced emission times has been found to be more than sufficient to capture the early-time emission.

We have performed various tests to check the resulting light curves from box-based calculations. First of all, we have compared box-based light curves to simulation-based light curves for the simulation underlying the box summary. An example is shown in Figure 6. At most, the difference between the two is on the order of a few percent (see bottom plots in figure). The exceptions are the early-time light curves for observers far off-axis, when the box approach smoothens out numerical noise more than the direct simulation calculation does.

Zoom In Zoom Out Reset image size

In Figures 7 and 8, we show a comparison to results of earlier studies. The close match between the box radio and optical light curves (solid and dashed lines for on- and off-axis observers, respectively) and the grid-based counterparts (thick dotted lines) in Figure 7 is remarkable given the differences in the methods by which they were obtained. The grid-based light curves are drawn from Van Eerten et al. (2010a) and are therefore based on a simulation with different refinement and resolution settings and strongly differing isotropic energy compared to the unscaled simulations of the current work. In addition to that, the earlier curves have been calculated with a completely different radiation algorithm, where a summation was done over the emitted power of each fluid cell (which therefore excludes the possibility of synchrotron self-absorption), rather than by employing the linear radiative transfer method used in the current work. On average (over the logarithmically spaced data points) the difference between the box light curves and the light curves from the earlier study shown in the figure is a factor of 1.15, with the biggest difference (a factor of 1.31) occurring for the optical light curves (both on and off axes) around 400 days. The fact that the off-axis light curves at some point cross the on-axis light curves and temporarily show higher flux levels is a result of relativistic beaming: at its most extreme it leads to the prediction of orphan afterglows, where the on-axis light curve remains effectively invisible relative to the off-axis light curve for observers at very high angles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 8, we show a comparison to light curves from Van Eerten & MacFadyen (2011b). These were obtained from simulations with lower resolution compared to the current simulations, which started from γ_b = 10 rather than γ_b = 25. This accounts for the early-time differences up to ∼20 days. Early-time differences aside, the average difference (over the logarithmically spaced data points) is a factor of 1.09, with the biggest difference (a factor of 1.23) occurring at late times for the two off-axis light curves.

Finally, in Figure 9 we show a comparison between optical and radio light curves based on a single opening angle simulation (0.2 rad) and light curves for the same opening angle reconstructed from interpolation between simulations with θ₀ = 0.175 rad and θ₀ = 0.225 rad. The figure shows that interpolated light curves generally match the light curves calculated from θ₀ = 0.2 rad. In this particular example the greatest discrepancy occurs between the low radio curves around 160 days, where the interpolated curve flux is briefly larger by a factor of 1.2. Note, however, that the figure represents an extreme scenario that does not occur in practice: by artificially removing θ₀ = 0.2 for the purpose of testing the interpolation, we have tested interpolation across Δθ₀ = 0.025 rad, whereas in practice the largest possible difference δθ₀ = 0.0125 rad. We conclude that the range of jet opening angles between θ₀ = 0.045 and θ₀ = 0.5 rad is adequately represented by our sample of simulations plus interpolation.

Zoom In Zoom Out Reset image size

The numerical errors between different simulations and box-based light curves given above should be compared to the difference between simulation/box-based light curves on the one hand and analytically calculated light curves on the other hand. Although the latter light curves have no numerical noise or resolution errors by definition, they have systematic errors due to the simplifications in the underlying assumptions for dynamics and radiating region that are far larger than the numerical noise in the former. In Van Eerten et al. (2010a) simulation results are compared to different analytical models and differences up to an order of magnitude in flux level and in time (for specific features such as the off-axis moment of peak flux) are seen, especially in the transrelativistic phase.

5.4. Fitting Methods

The "box"-based tool for GRB afterglow modeling presented in this paper can be applied to any broadband data set. Good spectral and temporal coverage is necessary to accurately determine all of the macro- and microphysical parameters. The broadband spectrum should span radio to X-ray frequencies to encompass all three characteristic frequencies of the synchrotron spectrum, and both early- and late-time data are needed to determine, for instance, the opening angle of the jet and the observer angle. We note that the tool allows for fitting of limited data sets by fixing some of the fit parameters.

To demonstrate the capabilities of our tool, we have selected the afterglow of GRB 990510. There are light curves available for this source in X-rays, in various optical bands, and at two radio frequencies. Historically, GRB 990510 was the first strong case for an achromatic jet break, and the broadband afterglow has been modeled by several authors, all using analytical expressions. Here, we show the results of our modeling with RHD simulations.

The modeling tool requires fluxes in mJy at all frequencies, which means that some conversions have to be applied to the optical and X-ray data. The radio data at 4.8 and 8.7 GHz were taken directly from Harrison et al. (1999). The X-ray count rates from Kuulkers et al. (2000) were converted to mJy using the conversion factors given in that paper. In our modeling we have used the V-, R-, and I-band data from Harrison et al. (1999), Israel et al. (1999), Stanek et al. (1999), Pietrzynski & Udalski (1999), Bloom et al. (1999), and Beuermann et al. (1999). We have corrected the optical magnitudes for Galactic extinction with E(B − V) = 0.20 (Schlegel et al. 1998) before converting the magnitudes into fluxes. Starling et al. (2007) have shown that the host galaxy extinction is negligible for GRB 990510 based on modeling of the X-ray to optical spectrum.

6.1. Fit Results

We have performed three different fits to the 205 data points of the GRB 990510 afterglow, for which we show the resulting fit parameters in Table 2 and the accompanying best-fit light curves in Figures 10–12. Figure 10 shows the fit results for a fit with fixed θ_obs = 0 and ξ_N = 1, which can be directly compared to the broadband fits performed by Panaitescu & Kumar (2002), who use analytical expressions. Figure 11 shows the best-fit light curves for a fit with ξ_N as a free parameter and fixed θ_obs = 0, and Figure 12 for a fit without any fixed parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Best-fit Model Parameters and 1σ Errors

Var.	PK02	On-axis, Fixed ξN	On-axis	Off-axis
θ0 (rad)	5.4+0.1− 0.6 × 10−2	7.5+0.2− 0.4 × 10−2	9.46−0.33+ 0.03 × 10−2	4.82+0.32− 0.04 × 10−2
Eiso (erg)	9.47+37− 2.27 × 1053	1.8+0.3− 0.1 × 1053	1.04+0.16− 0.02 × 1053	4.388+0.003− 0.605 × 1054
n0 (cm−3)	2.9+0.7− 0.9 × 10−1	3.0+0.4− 1.2 × 10−2	1.15+0.03− 0.19	1.115+0.258− 0.006 × 10−1
θobs(rad)	0 (fixed)	0 (fixed)	0 (fixed)	1.6+0.1− 0.1 × 10−2
p	1.83+0.11− 0.006	2.28+0.06− 0.01	2.053−0.006+ 0.007	2.089+0.013− 0.001
B	5.2+26− 2.9 × 10−3	4.6+0.8− 0.8 × 10−3	2.04+0.04− 0.30 × 10−3	1.36+0.19− 0.03 × 10−3
e	2.5+1.9− 0.4 × 10−2	3.73−0.68+ 0.07 × 10−1	6.8+0.6− 0.1 × 10−1	1.17+0.02− 0.12 × 10−2
ξN	1 (fixed)	1 (fixed)	5.4+0.6− 0.6 × 10−1	5.7+1.0− 1.7 × 10−2
χ2r	⋅⋅⋅	6.389	5.389	3.235
Ej	1.4+3.1− 0.3 × 1050	5.0+1.2− 0.8 × 1050	4.63+0.05− 0.10 × 1050	5.1+0.7− 0.8 × 1051

Notes. The column labeled PK02 lists the fit results from Panaitescu & Kumar (2002), translated from their units and 90% confidence intervals. For context we also include jet energy E_j results.

Download table as:  ASCIITypeset image

Table 2 shows a clear spread in best-fit parameters for the three fits we performed and also some differences with the results from Panaitescu & Kumar (2002). These differences can be mostly attributed to the fact that the synchrotron self-absorption frequency ν_a is not very well defined by this particular data set. The coverage at radio frequencies is fairly sparse compared to the optical and X-ray bands, and the flux uncertainties in the radio are also larger than at higher frequencies. For all three of our fits ν_a has a value around 10⁹ Hz at 1 day after the burst, while the peak frequency ν_m ∼ 10¹³ Hz at that time. The good coverage in the optical and X-ray bands enables an accurate determination of the cooling frequency ν_c, which is situated just above the optical bands, and the value of p. In contrast to the results from Panaitescu & Kumar (2002), we find p > 2 for all three of our fits rather than 1.8. Although our p > 2 lies within the error bar of their work, the converse is not true, which confirms p > 2, and that there is thus no need to include an additional high-energy cutoff on the relativistic particle distribution. The value of p has also been determined by several other authors based on optical and X-ray light-curve slopes (e.g., Harrison et al. 1999; Kuulkers et al. 2000; Panaitescu & Kumar 2001) or a combination of light curves and optical to X-ray spectra (Starling et al. 2008). In those studies the value for p falls in the range 2.1–2.2, consistent with the p value we have obtained in our fit without any fixed parameters.

It is clear from Table 2 that adding extra parameters introduces quite a strong variation in some of the parameters. In particular, the energy and circumburst density change from the on-axis to the off-axis fit with more than an order of magnitude, while in the off-axis fit the opening angle becomes a factor of two smaller by introducing a non-zero observer angle. This shows the importance of including the observer angle in broadband modeling of GRB afterglows, although this does require well-sampled light curves across the whole spectrum. More detailed studies of well-sampled broadband afterglows will be presented in a future paper.

In this paper, we present a method to directly fit light curves based on two-dimensional hydrodynamical jet simulations to broadband afterglow data. This provides a clear improvement over fits based on analytical models, which are not able to take complex features of the jet dynamics into account, such as the radial fluid profile, slow deceleration from ultrarelativistic to non-relativistic outflow, the sideways spreading, and resulting inhomogeneity along the shock front. The iterative fit procedure is possible because (1) we have shown that the jet evolution is scale invariant with respect to explosion energy and circumburst medium density, and (2) the results from high-resolution parallel RHD simulations can be summarized in a very compact form. The compressed version, a "box" summary, of the simulation "grid" data is possible once the blast wave lateral extent, radius, and radial width are known from the data at each emission time. The predicted flux is calculated for each data point for a given set of explosion and radiation parameters, and because the code takes into account both electron cooling and synchrotron self-absorption, it is applicable to the entire broadband afterglow spectrum.

In order to set up the boxes, a total of 19 RHD simulations were performed in two dimensions, with initial opening angles varying between 0.045 and 0.5 rad. Light curves for intermediate opening angles between different tabulated simulation results are obtained via interpolation at the fluid level. The simulations are run for a long time in order to ensure that late-time non-relativistic features such as the rise of the counterjet are covered as well. At very early times the outflow conforms to the self-similar Blandford–McKee solution, and this is used directly to calculate radiation from emission times before the starting point of the simulations. A comparison of the evolution of the 19 simulations reveals that, although hard to predict analytically, the evolution from strongly collimated BM jet to semi-spherical Sedov–Taylor blast wave is smooth. For small opening angles, the intermediate stage between the two asymptotes is more pronounced.

We present a number of tests for the resolution of box-based light curves by direct comparison to the simulations underlying the box summaries and to earlier work. The earlier light curves have been generated directly from simulations using different methods: by applying linear radiative transfer and by summation of the emitted power (in the optically thin case). The differences between box-based and simulation-based light curves are found to be a few percent at most, and the box summaries are found to correctly capture the shock profile at the different stages of blast wave evolution.

We have used GRB 990510 as a case study to demonstrate our method, because it has been observed over a wide range of frequencies. A simultaneous fit has been performed to data at two radio frequencies, three optical bands, and in X-rays. There are some substantial differences between our best-fit values and earlier-fit values found by Panaitescu & Kumar (2002), but also between the fits with various parameters fixed that we have performed. Most strikingly, including a non-zero observing angle in our modeling changes has a large impact on the obtained values for the blast wave energy and circumburst density. Furthermore, in contrast with Panaitescu & Kumar (2002), we find a value for the electron energy distribution index p > 2 and thus do not need to include a high-energy cutoff on the relativistic particle distribution.

The more accurate fits possible by the method of this paper help to further constrain the physics that shape the GRB afterglows. Explosion energy and circumburst density provide important clues to the nature of the progenitor star and burst environment. In addition, the types of fits to data sets covering a long time span that the method makes possible, where the complexities of the intermediate stage dynamics and the shape of the jet break are fully included, are necessary to establish a baseline for studies of the effect of more detailed models of the microphysics. The evolution of microphysical parameters, such as _B, is discussed in various papers (Van Eerten et al. 2010b; Panaitescu et al. 2006; Filgas et al. 2011).

Even without utilizing the possibility of fitting broadband data directly to simulation results, given the fact that the box summaries cover a large number of simulations not just directly but through scaling and interpolation, exploring the complete parameter space of impulsive jets in an ISM environment (with the inclusion of stellar wind environments being straightforward) should prove useful. It will allow us to test radiation mechanisms of physical interest other than pure synchrotron radiation in a realistic context. As long as there is no significant feedback on the dynamics of the outflow or dynamically relevant magnetic field involved, any generalization is possible even if it includes scattering. Examples of interest for different radiative processes are given by Giannios & Spitkovsky (2009) and Petropoulou & Mastichiadis (2009).

The source code of the broadband fit program described in this paper is publicly available at http://cosmo.nyu.edu/afterglowlibrary. The code has different settings: not only can it be used to fit box-based light curves to a data set and establish the uncertainties in the best-fit parameters, it can also be used to generate light curves and spectra directly for arbitrary frequencies, observer angles, observer times, and explosion and radiation parameters. This could be helpful for directly exploring how the different regions of the parameter space determine the shape of the afterglow, and for quickly creating light curves that can be expected or looked for in surveys (see, e.g., Nakar & Piran 2011; Roberts et al. 2011; Metzger & Berger 2012).

A number of practical improvements can be made to the code. In the case of very large data sets, even for the parallel version, where the data points to be calculated are evenly distributed on the cores, the total calculation time can become unwieldy. A remedy to this would be to no longer recalculate the flux independently for each data point but to calculate the flux at a fixed (large) number of values for observer time and frequency, chosen such to evenly cover the available data. The flux at the exact data point values can then be determined by interpolation between these values. Because the light curves are smooth, this will not impact the accuracy of the fit. For very long data sets (some GRBs have now been observed for several years, e.g., GRB 030329, Frail et al. 2000; van der Horst et al. 2008), more long-term simulation data can be added to the boxes. Boxes can also be generated for a stellar wind environment, since the same scale invariances apply. These applications and improvements will be presented in future work.

This research was supported in part by NASA through grant NNX10AF62G issued through the Astrophysics Theory Program and by the NSF through grant AST-1009863. The software used in this work was in part developed by the DOE-supported ASCI/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. We thank Erik Kuulkers for providing the X-ray fluxes. H.J.v.E. acknowledges hospitality from NASA/MSFC in Huntsville, where part of this research was conducted, and A.J.v.d.H. likewise acknowledges hospitality from New York University. Finally, H.J.v.E. and A.J.v.d.H. gratefully acknowledge Ralph Wijers, whose PhD supervision and research program at the University of Amsterdam have laid the foundation for their contributions to the current study. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

The expressions for the emission and absorption coefficient are drawn from Granot et al. (1999), based in turn on Sari et al. (1998), and for completeness we provide below the exact forms that we have implemented in our radiation code.³ The code solves the radiative transfer equation using

$$\begin{equation} \Delta I_\nu = (j_\nu - \alpha _\nu I_\nu) c \Delta t, \end{equation} \tag{ A1 }$$

where Δt is the lab frame time difference between two snapshots, j_ν is the emission coefficient, and α_ν is the absorption coefficient. The code subsequently calculates the flux by integrating over the area A of the emission plane, according to

$$\begin{equation} F_\nu = \frac{1+z}{d_L^2} \int d{A} I_\nu. \end{equation} \tag{ A2 }$$

Here, d_L is the luminosity distance and z is the redshift. Separating the coefficients into frequency-dependent and non-frequency-dependent components, j_ν = j_S × j_f and α_ν = α_S × α_f, we use for the non-frequency-dependent components

$$\begin{eqnarray} j_S &=& 9.6323 \frac{p-1}{3p-1}\frac{\sqrt{3} q_e^3}{8 \pi m_e c^2} \frac{\xi _N n^{\prime } B^{\prime }}{\gamma ^2 (1 - \beta \mu)^2}, \nonumber\\ \alpha _S &=& \frac{\sqrt{3} q_e^3 (p-1)(p+2)}{16 \pi m_e^2 c^2} \xi _N n^{\prime } B^{\prime } \gamma (1 - \beta \mu). \end{eqnarray} \tag{ A3 }$$

Here, q_e is the electron charge, m_e is the electron mass, n' is the comoving number density, B' is the comoving magnetic field strength, β is the fluid flow velocity as fraction of c, and μ is the angle between the outflow and the observer direction. The frequency-dependent part j_f is given by

$$\begin{equation} j_f = \left\lbrace \begin{array}{@{}ll} (\nu ^{\prime } / \nu ^{\prime }_m)^{1/3} & {\rm if }\ \nu ^{\prime } < \nu ^{\prime }_m < \nu ^{\prime }_c, \\ \ms (\nu ^{\prime } / \nu ^{\prime }_m)^{(1-p)/2} & {\rm if }\ \nu ^{\prime }_m < \nu ^{\prime } < \nu ^{\prime }_c, \\ \ms (\nu ^{\prime }_c / \nu ^{\prime }_m)^{(1-p)/2} (\nu ^{\prime } / \nu ^{\prime }_c)^{-p/2} & {\rm if }\ \nu ^{\prime }_m < \nu ^{\prime }_c < \nu ^{\prime }, \end{array} \right. \end{equation} \tag{ A4 }$$

when ν'_m < ν'_c and by

$$\begin{equation} j_f = \left\lbrace \begin{array}{@{}ll}(\nu ^{\prime } / \nu ^{\prime }_c)^{1/3} & {\rm if }\ \nu ^{\prime } < \nu ^{\prime }_c < \nu ^{\prime }_m, \\ \ms (\nu ^{\prime } / \nu ^{\prime }_c)^{-1/2} & {\rm if }\ \nu ^{\prime }_c < \nu ^{\prime } < \nu ^{\prime }_m, \\ \ms (\nu ^{\prime }_m / \nu ^{\prime }_c)^{-1/2}(\nu ^{\prime } / \nu ^{\prime }_m)^{-p/2} & {\rm if }\ \nu ^{\prime }_c < \nu ^{\prime }_m < \nu ^{\prime }, \end{array} \right. \end{equation} \tag{ A5 }$$

otherwise. Here, ν' denotes the comoving observer frequency. The synchrotron frequency ν'_m is given by

$$\begin{equation} \nu ^{\prime }_m = \frac{3 q_e}{4 \pi m_e c} (\gamma ^{\prime }_m)^2 B^{\prime }, \qquad \gamma ^{\prime }_m = \left(\frac{p-2}{p-1} \right) \left(\frac{\epsilon _e e^{\prime }}{\xi _N n^{\prime } m_e c^2} \right). \end{equation} \tag{ A6 }$$

The cooling break frequency ν'_c is estimated using a global cooling time, leading to

$$\begin{equation} \nu ^{\prime }_c = \frac{3 q_e}{4 \pi m_e c} (\gamma ^{\prime }_c)^2 B^{\prime }, \qquad \gamma ^{\prime }_c = \frac{6 \pi m_e c \gamma }{\sigma _T (B^{\prime })^2 t_c}, \qquad t_c \equiv t. \end{equation} \tag{ A7 }$$

Note that both γ'_m and γ'_c are comoving with the fluid, while γ_c in Sari et al. (1998) is in the rest frame.

For the self-absorption coefficient we ignore the effects of cooling, because in practice the self-absorption break frequency ν'_a ≪ ν'_c and because otherwise the limit of accuracy is set in practice by the use of a global cooling time rather than by any additional level of detail in the calculation of the absorption coefficient. We have

$$\begin{equation} \alpha _f = \frac{1}{\gamma ^{\prime }_m (\nu ^{\prime })^2} \left\lbrace \begin{array}{@{}ll} (\nu ^{\prime } / \nu ^{\prime }_m)^{1/3} & {\rm if }\ \nu ^{\prime } < \nu ^{\prime }_m, \\ \ms (\nu ^{\prime } / \nu ^{\prime }_m)^{-p/2} & {\rm if }\ \nu ^{\prime }_m < \nu ^{\prime }. \\ \end{array} \right. \end{equation} \tag{ A8 }$$

In our study we have solved the fluid equations using arbitrary values for explosion energy and density. For reference and in order to explicitly demonstrate the complete scale invariance of the simulations, we provide the fluid equations in terms of dimensionless parameters A, B, θ below.

The special relativistic fluid dynamics equations in spherical coordinates, assuming symmetry in angle ϕ in the x–y plane around the jet axis, are as follows:

$$\begin{eqnarray} &&\frac{\partial }{ct} \gamma \rho + \left(\frac{\partial }{\partial r} + \frac{2}{r}\right) \gamma \rho \beta _r + \frac{1}{r \sin \theta } \frac{\partial }{\partial \theta } \gamma \rho \beta _\theta \sin \theta = 0, \nonumber \\ &&\frac{\partial }{\partial ct} (h \gamma ^2 - p - \rho \gamma c^2) + \left(\frac{\partial }{\partial r} + \frac{2}{r}\right) (h \gamma ^2 \beta _r - \rho \gamma c^2 \beta _r)\nonumber\\ &&\qquad + \frac{1}{r \sin \theta } \frac{\partial }{\partial \theta } [(h \gamma ^2 \beta _\theta - \rho \gamma c^2 \beta _\theta) \sin \theta ] = 0, \nonumber \\ &&\frac{\partial }{ct} h \gamma ^2 \beta _r + \left(\frac{\partial }{\partial r} + \frac{2}{r}\right) \big(h \gamma ^2 \beta _r^2 + p\big) + \frac{1}{r \sin \theta } \frac{\partial }{\partial \theta } h \gamma ^2 \beta _r \beta _\theta = 0, \nonumber \\ &&\frac{\partial }{ct} h \gamma ^2 \beta _\theta + \left(\frac{\partial }{\partial r} + \frac{2}{r}\right) h \gamma ^2 \beta _\theta \beta _r + \frac{1}{r \sin \theta } \frac{\partial }{\partial \theta } \big(h \gamma ^2 \beta _\theta ^2 + p\big) = 0.\nonumber\\ \end{eqnarray} \tag{ B1 }$$

Here, β_r and β_θ are the fluid velocity components in the r and θ directions, respectively, in units of c and h the relativistic enthalpy density including rest-mass energy density. In terms of scale-free parameters A and B, the partial derivatives in r and ct can be expressed as

$$\begin{eqnarray} \frac{\partial }{\partial ct} & = & {-} \frac{A}{ct} \frac{\partial }{\partial A} + \frac{2B}{ct} \frac{\partial }{\partial B}, \nonumber \\ \frac{\partial }{\partial r} & = & \frac{1}{ct} \frac{\partial }{\partial A} - \frac{5B}{r} \frac{\partial }{\partial B}. \end{eqnarray} \tag{ B2 }$$

Combining these with the fluid equations above yields the scale-invariant forms

$$\begin{eqnarray} &&\left(-A \frac{\partial }{\partial A} + 2B \frac{\partial }{\partial B}\right) \gamma \frac{\rho }{\rho _0} + \left(\frac{\partial }{\partial A} - \frac{5 B}{A} \frac{\partial }{\partial B} + \frac{2}{A}\right) \gamma \frac{\rho }{\rho _0} \beta _r\nonumber\\ &&\qquad + \frac{1}{A \sin \theta } \frac{\partial }{\partial \theta } \gamma \frac{\rho }{\rho _0} \beta _\theta \sin \theta = 0, \nonumber \\ &&\left(-A \frac{\partial }{\partial A} + 2B \frac{\partial }{\partial B}\right) \frac{h \gamma ^2 - p - \rho \gamma c^2}{\rho _0 c^2}\nonumber\\ &&\qquad + \left(\frac{\partial }{\partial A} - \frac{5 B}{A} \frac{\partial }{\partial B} + \frac{2}{A}\right) \frac{h \gamma ^2 \beta _r - \rho \gamma c^2 \beta _r}{\rho _0 c^2}\nonumber\\ &&\qquad+ \frac{1}{A \sin \theta } \frac{\partial }{\partial \theta } \frac{(h \gamma ^2 \beta _\theta - \rho \gamma c^2 \beta _\theta) \sin \theta }{\rho _0 c^2} = 0, \nonumber \\ &&\left(-A \frac{\partial }{\partial A} + 2B \frac{\partial }{\partial B}\right) \frac{h \gamma ^2 \beta _r}{\rho _0 c^2} + \left(\frac{\partial }{\partial A} - \frac{5 B}{A} \frac{\partial }{\partial B} + \frac{2}{A}\right) \frac{h \gamma ^2 \beta _r^2 + p}{\rho _0 c^2}\nonumber\\ &&\qquad + \frac{1}{A \sin \theta } \frac{\partial }{\partial \theta } \frac{h \gamma ^2 \beta _r \beta _\theta }{\rho _0 c^2} = 0, \nonumber \\ &&\left(-A \frac{\partial }{\partial A} + 2B \frac{\partial }{\partial B}\right) \frac{h \gamma ^2 \beta _\theta }{\rho _0 c^2} + \left(\frac{\partial }{\partial A} - \frac{5 B}{A} \frac{\partial }{\partial B} + \frac{2}{A}\right) \frac{h \gamma ^2 \beta _\theta \beta _r}{\rho _0 c^2}\nonumber\\ &&\qquad + \frac{1}{A \sin \theta } \frac{\partial }{\partial \theta } \frac{h \gamma ^2 \beta _\theta ^2 + p}{\rho _0 c^2} = 0. \end{eqnarray} \tag{ B3 }$$

Quantities such as γρ/ρ₀, etc., are dimensionless and therefore (scale-invariant) functions of the scale-invariant dimensionless parameters A, B, and θ. For radial outflow and for limiting values of A, the fluid equations further reduce to self-similarity.