A Numerical implementation

Gamma-ray burst afterglows from trans-relativistic blast wave simulations

Abstract

We present a study of the intermediate regime between ultra-relativistic and nonrelativistic flow for gamma-ray burst afterglows. The hydrodynamics of spherically symmetric blast waves is numerically calculated using the amrvac adaptive mesh refinement code. Spectra and light curves are calculated using a separate radiation code that, for the first time, links a parametrisation of the microphysics of shock acceleration, synchrotron self-absorption and electron cooling to a high-performance hydrodynamics simulation. For the dynamics we find that the transition to the nonrelativistic regime generally occurs later than expected, that the Sedov-Taylor solution overpredicts the late time blast wave radius and that the analytical formula for the blast wave velocity from Huang et al. (1999) overpredicts the late time velocity by a factor 4/3. Also we find that the lab frame density directly behind the shock front divided by the fluid Lorentz factor squared remains very close to four times the unshocked density, while the effective adiabatic index of the shock changes from relativistic to nonrelativistic. For the radiation we find that the flux may differ up to an order of magnitude depending on the equation of state that is used for the fluid and that the counterjet leads to a clear rebrightening at late times for hard-edged jets. Simulating GRB030329 using predictions for its physical parameters from the literature leads to spectra and light curves that may differ significantly from the actual data, emphasizing the need for very accurate modelling. Predicted light curves at low radio frequencies for a hard-edged jet model of GRB030329 with opening angle 22 degrees show typically two distinct peaks, due to the combined effect of jet break, non relativistic break and counterjet. Spatially resolved afterglow images show a ring-like structure.

1 Introduction

Gamma-ray burst (GRB) afterglows can be explained from the interaction between an initially relativistic shock wave of hot fluid and the medium surrounding the burster. On passage of the shock electrons get accelerated to relativistic velocities (even with respect to the already relativistic local fluid flow) and small scale magnetic fields are generated. Under influence of the magnetic field, the electrons will produce synchrotron radiation, which will be seen by the observer. This model has been very successful when applied to broadband afterglow data, but thus far model predictions have been made using simplifying assumptions for the blast wave structure (approximating the blast wave width by a homogeneous slab, e.g. Wijers et al. 1997; Meszaros & Rees 1997; Sari et al. 1998; Rhoads 1999) or from analytical solutions in either the ultrarelativistic or the nonrelativistic regime (e.g. Granot et al. 1999; Gruzinov & Waxman 1999; Wijers & Galama 1999; Frail et al. 2000).

Since the beginning of the decade, fluid simulations have been performed to study afterglow blast waves and their resulting spectra (see Granot et al. 2001; Downes et al. 2002). More recent simulations have been used to address the specific theoretical issue of the visible effect of the blast wave encountering a density perturbation (Nakar & Granot, 2007; van Eerten et al., 2009). Very recently Zhang & MacFadyen (2009) studied the transition to the transrelativistic regime and the spreading of a collimated outflow, using an adaptive mesh technique for the fluid simulation. They made some simplifying assumptions for the radiation mechanism, when compared to the early analytical efforts (e.g. Granot et al. 1999), such as approximating the cooling time by the lab frame time and ignoring synchrotron self-absorption.

The aim of this paper is to present a theoretical and qualitative study of the transition regime between relativistic and nonrelativistic blast waves and the effect on the light curves and spectra at various wavelengths, using adaptive mesh relativistic fluid simulations for blast waves from an explosion in a homogeneous medium, while including all details of the synchrotron radiation mechanism that have been used for earlier analytical estimates. Also we present resolved afterglow images. We study spherical blast waves and sharp edged jets obtained by taking conic sections from a spherically symmetric fluid flow.

Obviously these simulations do not yet fully address the complete GRB afterglow picture of a realistic, 2D-dynamical jet, which we address in future work. However, some GRB afterglows have power law decays that last for months without a jet break, and thus may be (nearly) spherical. These are of course already addressed in the present work. Also, by studying the conic sections from spherical flows, we already address some aspects of jet behaviour, which allows us to probe some outstanding issues, such as whether the receding jet may lead to visible features in the late light curve, and whether a dynamical jet break must be truly achromatic. Finally, any fluid flow behaviour typical to higher-dimensional simulations, like lateral spreading of the jet, is best understood from a direct comparison to one-dimensional simulations and its effects on the light curve will in practice be modeled as a deviation from the heuristic description based on analytical approximations and one-dimensional simulations (i.e. as an additional smooth jet break). A companion paper is in preparation that will discuss the practical consequences for broadband afterglow data fitting from the underlying model from this paper.

This paper is organized as follows. In section 2 we discuss our radiation code and how it expands upon an approach outlined earlier in Van Eerten & Wijers (2009), hereafter EW09. A proper treatment of synchrotron radiation and shock wave generation of accelerated particles and small scale magnetic fields requires us to trace some additional quantities along with the fluid quantities.

In section 3 we provide the details of our simulations that assume typical GRB parameters. We show how the blast wave starts out in the ultrarelativistic regime and smoothly approaches the nonrelativistic regime. We discuss the consequences of different equations of state for the fluid and how our simulations differ from analytical approximations for the nonrelativistic regime. We show how the fluid lab frame density divided by the fluid Lorentz factor squared right behind the shock remains always close to four times that in front of the shock, even though we have differing adiabatic indices in both the relativistic and nonrelativistic regimes. Three additional quantities needed to be traced and we present results for the behaviour of these three: the accelerated electron number density, the magnetic field energy density and the accelerated particle distribution upper cut-off Lorentz factor. We explain how calculation of the latter especially is numerically challenging and how it shapes the spectrum beyond the cooling break.

In section 4 we take our results from section 3 and calculate spectra and light curves. We calculate spectra at 1, 10, 100, 1,000 and 10,000 days in observer time. We separately discuss the different factors contributing to the shape of the light curves: the equation of state, the evolution of the magnetic field and the evolution of the accelerated particle distribution.

We then turn to the specific case of GRB 0303029 in section 5. We take the explosion parameters that have been established for this burst by previous authors to set up a simulation. We qualitatively compare the resulting light curves to radio data at different wavelengths, assuming both a spherical explosion and a hard edged jet with opening angle of degrees. We provide spatially resolved radio images and make a qualitative prediction for the expected signal at radio wavelengths that will be observable with the next generation of telescopes, like lofar.

We discuss our results in 6. In the appendices we provide additional technical details on the numerical implementation of our approach and a discussion on the theoretical limitations and assumptions of our approach.

2 The radiation code

In this paper we follow the approach first outlined in EW09, where we calculate spectra and light curves from the output of a relativistic hydrodynamics (RHD) code using a separate radiation code. For the RHD simulations we use amrvac, a high performance code that includes adaptive-mesh refinement (AMR) (see Keppens et al. 2003; Meliani et al. 2007). amrvac calculates the evolution of the following conserved variables:

(1)

with the Lorentz factor, the proper density, the relativistic (i.e. including rest mass) enthalpy density, the three velocity, the pressure and the speed of light. In the entire paper, all comoving quantities will be primed.

In the second stage we use a radiation code to obtain the received flux for a given observer frequency, time and distance, from the local values of conserved variables at any contributing point in the fluid (we also use two auxiliary quantities, and , that amrvac stores as well in order to facilitate its calculation of the time evolution of the conserved variables). The radiation mechanism that is considered is synchrotron radiation and a number of parameters have been introduced in EW09 that capture the underlying radiation and shock microphysics. There are four of these ‘ignorance’ parameters. The fraction of the thermal energy that resides in the tangled-up magnetic field that is generated by the passage of a shock usually has a value around . The fraction of electrons that is accelerated into a relativistic power law distribution in energy also by the passage of a shock is usually of order unity in the relativistic regime. The thermal energy fraction captured by these electrons and minus the slope of the electron distribution .

The flux calculated by the radiation code is given by

(2)

Here denotes redshift, denotes the observer luminosity distance, the received power per unit volume, frequency and solid angle, the equidistant surface element given by the intersection of the fluid grid with that surface from which radiation is poised to arrive exactly at , and the emission time. Note that in this terminology, flux is defined per unit frequency. The integral

is effectively an integral over the entire radiating volume. is the angle between the local fluid velocity and the observer position, the fluid velocity in units of and the factor is a retardation effect due to the moving of the radiating source. The detailed dependency of the received power on the ignorance parameters and local fluid conditions is explained in EW09. However, in that paper only ultra-relativistic flows were addressed and in order to include subrelativistic and nonrelativistic flows as well, a number of features were added to our radiation code. Also we have added synchrotron self-absorption and the possibility to resolve the signal from the fluid into an image on the sky. We now have a generic radiation code that is capable of calculating the spatially resolved synchrotron radiation profile from an arbitrary fluid flow. The additional physics that we have included is explained below, with some of the practical numerical issues discussed separately in appendix A.

2.1 Realistic equation of state

In EW09 we applied a fixed adiabatic index equation of state (EOS)

(3)

where is the thermal energy density. In practice was always set to . However, when following a fluid from the relativistic regime (with flow velocities and thermal energy density dominating the rest mass energy density) down to the classical regime, this fixed adiabatic index becomes too restrictive. We therefore apply a Synge-like EOS (Synge, 1957) that results in an effective adiabatic index varying smoothly from to its classical limit :

(4)

where denotes the comoving energy density including rest mass, . This EOS has already been applied in amrvac (see Meliani et al. 2008, 2004). Also, because the radiation code reads both the conserved variables as well as from disc directly, it does not invoke any EOS itself, and no change in the radiation code was needed. The resulting effective adiabatic index is given by

(5)

The effect of an advanced EOS on the behaviour of the fluid is profound and we discuss this in detail in section 4.3.

2.2 Electron cooling

The shape of the observed spectrum from a single fluid cell, if electron cooling does not play a role, follows directly from the dimensionless function , first introduced in EW09. It has the limiting behaviour for small and for large . The received power depends on this shape and on the local fluid quantities via

(6)

Here denotes the lab frame number density (of all electrons, both accelerated and thermal). denotes the local comoving magnetic field strength, calculated from the thermal energy density after the passage of a shock. and denote electron mass and charge, respectively. The frequency is the synchrotron peak frequency, and it is related to the lower cut-off Lorentz factor of the power law accelerated electrons via

(7)

If cooling plays no role, the evolution of is completely adiabatic, which has as a consequence that the total fraction of the local thermal energy density residing in the power-law accelerated particle distribution remains fixed. will be related to throughout the downstream fluid according to

(8)

When cooling does play a role, however, this picture is changed. It now becomes necessary to introduce an upper cut-off Lorentz factor as well. In a single fluid element, no accelerated electrons with energies above will be found, because these have cooled to energies at or below . The temporal evolution of any electron Lorentz factor , and therefore of and as well, is given by

(9)

where is the Thomson cross section, the electron mass and the comoving time. The final term in this equation reflects synchrotron radiation losses, and if it is omitted only the adiabatic cooling term is left and it can be shown that this will result in the aforementioned fixed relation between and . In light of the previous subsection on the EOS, it may be worth noting that equation 9 is derived for a relativistic electron distribution with adiabatic index and that this remains valid even if the bulk of the fluid becomes nonrelativistic. After all, the power-law accelerated electrons are relativistic by definition.

The above has the following consequences for the simulations and radiation code. Because for low values of , the radiation loss term can be neglected next to the adiabatic expansion term, we will not apply equation 9 to and we will continue to calculate locally using equation 8 in the radiation code. For this is not an option and we numerically solve equation 9 in amrvac, resetting to a high value wherever we detect the passage of a shock. This reset implements the shock-acceleration of particles. In appendix A we discuss the numerical issues of this approach in some more detail. Also, in appendix B.2 we show that the gyral radius even for the high energy electrons contributing to the observed spectrum (within the frequency range under consideration, Hz) is orders of magnitude smaller then the relevant fluid scales.

Finally, we summarize the consequences of electron cooling for the spectrum discussed earlier in EW09. The received power from a fluid element is now given by

(10)

with the relations between the upper and lower cut-off Lorentz factors and their corresponding critical frequencies given by equation 7. is a generalisation of and the flux at frequencies above drops exponentially. That the resulting spectrum from the entire fluid does not show an exponential drop is due to the fact that there will always be some fluid elements contributing for which is still sufficiently high. The effect of this ‘hot region’ close to the shock front (with a size that depends on the observer frequency) on the composite synchrotron spectrum from a shock will be a steepening of the slope by instead. The cooling break is found at that frequency for which the width of the hot region becomes comparable to the width of the blast wave.

2.3 Magnetic field energy evolution

The magnetic field directly behind a shock has been parametrised using

(11)

Furthermore we assumed the number of magnetic flux lines threading a surface comoving with a fluid element to remain invariant, resulting in

(12)

For relativistic fluids this implies that the fraction remains fixed downstream, because . For a changing adiabatic index it is no longer possible to calculate a posteriori from , since the relation between the two is now no longer fixed. It becomes necessary to numerically solve in amrvac the equation

(13)

Like , we reset whenever a shock is encountered. The practical implementation of the evolving magnetic field is again discussed in appendix A. We note here that the assumption of frozen field lines is not essential, and that we can in principle include different magnetic field behaviour either by adding a source term to equation 13 (parametrising for example, magnetic field decay through reconnection) or by implementing a different equation entirely.

2.4 Changing fraction of accelerated particles

Although , the fraction of electrons accelerated by the passage of a shock is often assumed to be of the order unity for highly relativistic blast waves, it has to be lower at late times because otherwise there would not be enough energy available per accelerated electron to create a relativistic distribution (in other words, to ensure that ). We have implemented this change in our code by replacing user parameter by , that is, the fraction of electrons that is accelerated in the nonrelativistic limit. The fraction at the relativistic limit we set to one. Because is the most direct measure of how relativistic the fluid flow locally is, we have parametrised the simplest possible smooth transition between both limiting cases by

(14)

Whenever the passage of a shock is detected, amrvac resets the number density of accelerated electrons according to , with determined using the equation above. As with the magnetic field energy density, we now need to follow explicitly. Because is a number density, its evolution is described by a continuity equation, following

(15)

and is therefore easily implemented in amrvac.

2.5 Synchrotron self-absorption

In previous work we have solved equation 2 by first integrating over for a given emission time (and thus for a single snapshot), followed by an integration over . If we switch the order of the integrations then the integral over represents the solution to a linear radiative transfer equation without absorption, with the intensity given by

(16)

The integral over then represents a summation over all rays. The full linear radiative transfer equation including synchrotron self-absorption has the form

(17)

with and . The synchrotron self-absorption coefficient is given by

(18)

Here denotes the emitted power per ensemble electron and the electron number density for relativistic electrons accelerated to . Integrating over possible electron Lorentz factors yields by definition. These quantities are defined and explained in detail in EW09 (see also appendix A.3).

In this treatment of the self-absorption coefficient we only take into account transitions between already occupied energy levels of electrons, leading to the integration limits of equation 18 being exactly and . In this way we ignore stimulated emission arising from a population inversion below . This results in values of the absorption coefficient that are larger by a factor of when compared to Granot & Sari (2002).

In our radiation code, we now calculate the linear radiative transfer equation for each individual ray by not integrating over the two-dimensional surface (to get a single flux value from the collection of rays) until after the final snapshot has been processed. In addition to allowing us to include the effect of self-absorption, we now also get a spatially resolved signal from the fluid, showing the expected ring structure (extending predictions from Granot & Loeb 2001 to the nonrelativistic regime). We use an adaptive-mesh type approach to in order to ensure an adequate spatial resolution, see appendix A.4.

3 Fluid dynamics

In this section we describe the setup of our relativistic fluid simulations and compare the results against the theoretically expected behaviour.

3.1 Expected early and late time behaviour

Both the early and late time behaviour of the fluid can be described by a self-similar solution that is determined completely from the explosion energy and the circumburst number density .

At early stages, the Blandford-McKee (BM) solution (Blandford & McKee, 1976) for relativistic blast waves predicts the following relation between the shock front fluid Lorentz factor and the explosion time (, which is the same as the emission time ):

(19)

The density is related to the number density through the proton mass: . The shock radius is then given by

(20)

To lowest order is just , while the shock front fluid velocity . Further analytical equations for the fluid profile (in terms of pressure , Lorentz factor , number density , etc.), behind the shock front can be found in Blandford & McKee 1976.

At late stages the evolution of the blast wave is described by the Sedov-Taylor (ST) solution (Sedov, 1959; Taylor, 1950). For a fixed adiabatic index , the shock radius is now given by

(21)

which follows directly from dimensional analysis (except for the numerical constant). In this classical approximation, the speed of light does not appear. The shock front Lorentz factor is approximately one, while can be found from . Again analytical formulae for the fluid profile exist in the literature (Sedov, 1959).

At some point in time the evolution of the blast wave will no longer be adequately described by the BM solution but will become more and more dictated by the ST solution. An estimate for the turning point (Piran, 2005) can be made by equating the explosion energy to the total rest mass energy that is swept up:

(22)

Solving for returns the approximate radius at which the original explosion energy in the blast wave is no longer dominant over its rest mass energy.

Analytical estimates for the bulk fluid flow velocity including the intermediate regime also exist. One such example is found in Huang et al. (1999) and we discuss it in more detail as well as compare their prediction for right behind the shock front directly against our simulation results in section 3.3.1.

3.2 Setup of simulations

We have performed a number of simulations using the typical values for a GRB exploding into a homogeneous medium. We set up our simulations starting from the BM solution. The isotropic explosion energy erg, the medium number density cm. We have set the initial shock Lorentz factor to 10 (and the fluid Lorentz factor therefore , differing by a factor ). Although both amrvac and the radiation code are able to deal with far higher Lorentz factors, the focus for this research is on the transition to the nonrelativistic regime and for that purpose this relatively low Lorentz factor is sufficient. We have continued the simulations until the fluid proper velocity in the lab frame .

We have used both the advanced equation of state and a fixed adiabatic index at and . In the advanced EOS simulation we have also calculated the other quantities mentioned in the previous section: and . The value at the shock front for was set to the standard 0.01 and the non-relativistic limit for was set at . Sufficiently high values for at the shock front are chosen, generally on the order of .

In amrvac it is the number of refinement levels that determines the accuracy of the simulation. We have used 17 levels of refinement and 120 cells at the lowest refinement level. The grid was initially taken to run from cm to cm. The effective spatial resolution due to adaptive mesh refinement was therefore cm. This should be compared against the width of the blast wave at the start of the simulation, when it is the smallest. This is approximately equal to cm, for a starting shock Lorentz factor of 10.

Convergence of our results has been checked by performing simulations at different refinement levels and by simulations running for a shorter time on a smaller grid (thereby increasing the resolution). For the light curves and spectra we have used simulations with a shorter running time of days. At this stage the fluid velocity directly behind the shock is still six percent of light speed, but we have full coverage up to days in observer time. The corresponding grid size is cm to cm, leading to an effective resolution of cm. On a standard desktop PC1 the RHD simulations typically took a few days to complete and the radiation calculation a few hours.

3.3 Results

Blast wave velocity

Figure 1: at the shock front for the advanced EOS simulation (solid line), along with its asymptotic behaviour both at early and late stages. For comparison we have also plotted a prediction from Huang et al. (1999) (see text).

The solid line in figure 1 shows at the shock front for the advanced EOS simulation. The expected scaling behaviour at the early stage is dictated by and at the late stage by . We have plotted this asymptotic behaviour as well, setting the early stage scaling coefficient from the initial value at and the late stage scaling coefficient at (this point lies far to the right outside the plot). The shock velocity is shown to smoothly evolve from the BM solution to the ST solution. The meeting point of the asymptotes at days lies at . At this point for the fluid , so the fluid is still moving at a significant fraction of the speed of light.

According to eq. 22, the predicted radius for the transition to occur is parsec for the initial explosion energy and circumburst density that we have used, corresponding to a lab frame time days. We therefore conclude that the transition point from the relativistic to the nonrelativistic regime is far later than predicted by .

Also plotted in fig. 1 is the predicted value for from Huang et al. (1999), which we have implemented as follows. The starting point is

(23)

the differential equation proposed by the authors to depict the expansion of GRB remnants, simplified to the adiabatic case. Here denotes the rest mass of the swept-up medium and the mass ejected from the GRB central engine. Our approach starting from the BM solution is a limiting case where . The term was included by Huang et al. (1999) to incorporate a coasting phase. When solving eq. 23 we will use a very high () initial bulk fluid Lorentz factor and by assuming we converge on the limiting scenario used in our simulations. Eq. 23 can be analytically solved to yield

(24)

which (numerically) leads to once we apply

(25)

and

(26)

Here is measured in the simulation lab frame (i.e. it does not refer to observer time).

The resulting curve for initially lies below the simulation result, but ends up above at times the simulation value. The initial and final slopes for the analytical curve are correct by construction. We conclude that the approach from Huang et al. (1999) initially underestimates the BM phase and significantly overestimates the late stage flow velocity. The transition point between the relativistic and nonrelativistic regime also lies at an earlier time for the analytical curve, closer to the analytically predicted .

Blast wave radius

Figure 2: The resulting blast wave radii as a function of lab frame time for different simulations. The steady slope line shows the radius as predicted by the ST solution. The different simulations end up in the asymptotic regime with different radii: the ends up above the ST solution, the advanced EOS below the ST solution between the others and the lowest. The bottom curve shows the effective adiabatic index for the advanced EOS, minus . It starts at approximately zero at the left of the plot and proceeds to its asymptotic limit in the nonrelativistic case.

In figure 2 we plot the blast wave radius as a function of lab frame time for three different simulations: fixed adiabatic index at and and using the advanced EOS. Also we plot the radius as predicted from eq. 21 and, for the advanced EOS simulation, the difference between the effective adiabatic index and its relativistic limit . The latter illustrates how relativistic the fluid still is in terms of temperature (as opposed to flow velocity). At the intersecting point for the asymptotes at 1290 days the effective adiabatic index is already quite close () to its nonrelativistic limiting value. After 3800 days, when the time evolution of all the radii has become practically indistinguishable from from the ST solution we still see a difference between the different radii. At this time the ST radius is 1.358 parsec, the radius is 1.197 parsec, the radius is 1.388 parsec and the advanced EOS radius 1.313 parsec. Taking the advanced EOS radius as a standard, this implies that ST overpredicts the radius at this stage by 3.4 percent, overpredicts the radius by 5.7 percent and underpredicts the radius by 8.9 percent. Because all radii follow close to the same temporal evolution at this stage, these errors will only very gradually become smaller throughout the further evolution of the blast wave. Therefore, a derivation of the quantity from the radius using the Sedov-Taylor equation 21, is likely to overpredict its value by approximately 18 percent within any time interval of practical interest.

Density and energy profiles

Figure 3: Comoving number density profile. The profiles were taken at times corresponding to emission arrival times for the closest part of the shock front (i.e. with velocity directly towards the observer) at 10, 100, 1,000 and 10,000 days. Listed in increasing number and including the initial profile, these times correspond to lab frame times of 137, 387, 761, 2,227 and 12,583 days. Later times correspond to curves peaking further to the right in the plot.
Figure 4: Lab frame number density divided by . This effectively scales the shock profile to 4 at the shock front. The same lab frame times as in fig. 3 have been used.

In figure 3 we have plotted the comoving fluid number density profile (of the protons or the electrons, not both) at five different moments in time. Because later on we discuss spectra and light curves up to an observer time of 10,000 days, we have chosen emission times corresponding to arrival times of the shock front (using ) up to 10,000 days as well. The earliest fluid profile shows the initial conditions calculated from the BM solution when the shock Lorentz factor is 10. After some time, the number density at the shock front can be seen to tend to the value predicted from the shock-jump conditions for a strong classical shock, which is for the classical value of the adiabatic index .

What is shown in figure 4 is an interesting feature of the blast wave, which is that the lab frame density directly behind the shock divided by the squared fluid Lorentz factor directly behind the shock, throughout the entire simulation. This can be seen analytically to hold from the Rankine-Hugoniot relations in both the ultrarelativistic and nonrelativistic case, even though the adiabatic indices have different fixed values, from

(27)

which can be viewed as the relativistic generalization of the classical compression ratio and holds for arbitrary . When we use an advanced EOS, where we let smoothly evolve from to , we see from the figure that this generalized compression ratio remains very close to four even at intermediate times. We make use of this feature for the shock detection algorithm (see appendix A.2).

Figure 5: Thermal energy density profile for the same lab frame times as in fig. 3.

In figure 5 we have plotted the thermal energy density at the same times as the number density. Unlike the density at the shock front, the thermal energy density is not expected to tend to a fixed value. The ST solution instead predicts a steep decline , which is why the final shock front thermal energy density is many orders of magnitude smaller than the initial shock front thermal energy density.

Magnetic field and particle acceleration

Figure 6: The fraction of the thermal energy that resides in the magnetic field energy for the same lab frame times as in figure 3.

We now turn to those quantities calculated in amrvac solely to aid in the construction of spectra and light curves, and that have no feedback on the dynamics. In figure 6 we have plotted . Because we assumed the number of field lines through a fluid surface element to remain frozen (see eqn. 12), the magnetic field energy density declined less rapidly than the thermal energy and as a consequence the local fraction increased. A discussion on the merit of our assumption about the magnetic field behaviour is outside the scope of this work (and from particle-in-cell simulations it can certainly be argued that it is not perfect, see e.g. Chang et al. 2008). However, our plot does show that at least it does not lead to unphysical values or strong inconsistencies. The maximum value for found in fig. 6 is 0.037 (up from 0.01 at the shock front), which is not unreasonably large and, besides, occurs far downstream in a region that will contribute negligibly to the observed flux. We emphasize that is a relative measure and that both thermal and magnetic energy densities drop steeply, both with respect to earlier times and with respect to their value at the shock front at any time. The numerical method presented in this paper for parametrizing the magnetic field energy density is quite general and can be readily modified to study different parametrizations.

Figure 7: , the fraction of electrons accelerated to a power law distribution for the same lab frame times as in fig. 3.

In figure 7 we have plotted the fraction of electrons that are accelerated to a power-law distribution. This fraction was taken to smoothly decrease from unity in the relativistic regime down to 0.1 for our simulation in the nonrelativistic regime. The rightmost profile, with the shock front arriving at 10,000 days, has down to 0.16.

Figure 8: Upper cut-off Lorentz factor of the power law distribution , normalized to 1 at the shock front, for the same lab frame times as in fig. 3.

In figure 8 we have plotted the normalized values for , the upper cut-off Lorentz factor of the power-law particle distribution. Although formally should be reset to infinity at the shock front, we picked a value corresponding to a cut-off above Hz through the entire simulation, for a fluid element heading directly towards the observer (see also section 2.2 and appendix A). In our simulation settings, this results in peak values of the order , and because these values are arbitrary as long as they are sufficiently large, we have normalized the profiles. The profiles show two things. First, they show the steep decline directly after the injection of new hot electrons. This steep decline and the orders or magnitude difference between shocked and unshocked are numerically challenging, which is why we implemented the logarithm of in our code instead (again see appendix A). Second, the width of the profile is a measure for the size of the hot region discussed in section 2.2. It can be seen from the figure that the width of the profile increases over time. The width will nevertheless remain smaller than the width of the density profile by far. In our simulations, we resolve the profile and use it to determine the local refinement level.

4 Spectra and light curves

Using the simulation data described in the previous section we have calculated spectra and light curves at various observation times and frequencies. We have saved a total number of 10,000 snapshots of the fluid profile, with seconds between consecutive snapshots, corresponding to a resolution cm. Although this resolution is of the same order as the initial shock width, it is still sufficient at the early stage because the shock initially nearly keeps up with its own radiation. The effective resolution is given by cm, which is only a factor 10 larger than the spatial resolution of cm and corresponds to a temporal resolution of seconds. It is therefore ensured that the blast wave in the initial stage is covered by over a hundred snapshots.

4.1 Expected spectral and temporal behaviour

The scaling behaviour for the critical frequencies and the flux is well known from analytical estimations assuming a homogeneous radiating slab directly behind the shock front with fluid properties determined via either the BM or ST solution (see e.g. Wijers & Galama 1999; Frail et al. 2000; Granot & Sari 2002). We summarize the interstellar medium (ISM) scalings below, with denoting the peak frequency, the synchrotron self-absorption critical frequency and the cooling break frequency. At the observer times and frequencies in this paper we find either or .

In the relativistic limit, the corresponding scalings are

(28)
(29)
(30)

for the critical frequencies. Note that now refers to observer times. The flux above both peak and self-absorption break scales as

(31)

If we get the following flux scaling below the peak break:

(32)

If we have below the self-absorption break

(33)

In the nonrelativistic limit the scalings are

(34)
(35)
(36)

for the critical frequencies. The flux above both peak and self-absorption break scales as

(37)

If we get the following flux scaling below the peak break:

(38)

If we have below the self-absorption break

(39)

The summary above shows that only the temporal behaviour of the break frequencies and fluxes is altered by the transition to the nonrelativistic regime. We therefore do not expect spectra calculated from our simulations covering the transition to differ in slope from the slopes calculated above. The light curve slope, however, may differ.

4.2 Spectra

Figure 9: Spectra at different observer times. The smooth curves show simulated spectra at different observer times: 1, 10, 100, 1,000 and 10,000 days, with later observed spectra having lower flux in the high frequency range. For comparison we have included predicted slopes at the different power law regimes after 1 day, for both (solid line) and (dashed line).

In figure 9 we have plotted spectra for a number of different observation times, ranging from 1 day to 10,000 days. For comparison we have also plotted the different power law slopes for 1 day as predicted by Granot & Sari (2002), where we have added a dependency on . We plot predictions for both and . It can be seen that the simulated spectrum still lies closer to the prediction, just as we would expect for an early time spectrum. Because of shifts in both flux level and position of the spectral breaks for different values of , the flux does not always lie in between the analytical predictions. For example, because the peak frequency for the simulation lies close to that of the prediction and flux lies also closer to but below, the resulting flux at higher frequencies ends up below both predictions.

Figure 9 proves that our method works and that the asymptotic behaviour for the spectral slopes matches the predicted slopes. For frequencies above the self-absorption break and below the cooling break this merely confirms that the synchrotron spectral function has been implemented correctly. The flux at frequencies above the cooling break however, shows the consequence of a finite and evolving upper cut-off . A slope is reproduced that matches the prediction. It has been explained above and in EW09 how this slope now arises as a product of the interplay between the hot region and the blast wave width.

At the low frequencies, where synchrotron self-absorption plays a role, the simulations also reproduce a spectral slope that corresponds to what was expected from analytical calculations. The flux level is now dictated by the radiative transfer equation through a medium that is no longer completely transparent at these frequencies. As discussed in section 2.5, the resulting flux will differ by a factor of a few from Granot & Sari (2002), due to a difference in approach when calculating the absorption coefficient from the particle distribution.

We emphasize that fig 9 covers 10 orders of magnitude in frequency, 8 orders of magnitude in flux and four orders of magnitude in observer time. As we expected from analytical calculations, the spectral slopes in the different power law regimes do not change over time. The transitions between the different regimes are smooth. An explicit calculation of the sharpnesses of the transitions will be presented in a follow-up paper.

4.3 Light curves

Figure 10: Comparison of optical (at light curves for different equations of state. The top curve has , the center curve the advanced EOS and the bottom curve . For clarity as few complications as possible are included: cooling and self-absorption are switched off, is fixed at and everywhere. Also plotted are the expected relativistic slope and nonrelativistic slope .

We will use optical light curves to illustrate the consequences of the different assumptions and model parameters. In fig. 10 we present simulated light curves for simulations that differ only in the EOS used. Electron cooling and self-absorption have been disabled, is fixed at and everywhere throughout the simulation. This allows for a clear view on both the effect of the EOS and of the transrelativistic break. The latter can be found at days for all three simulations. This is somewhat earlier than the transition time determined from the fluid flow in section 3.3.1, which we determined to be around 1290 days (the difference is due to relativistic beaming). The transition time is also later than what is usually assumed for the nonrelativistic transition by nearly a factor three.

Figure 11: Direct comparison between thermal energy density profiles for the different equations of state. The top profile has , the center curve the advanced EOS and the bottom profile has . All snapshots are taken at 515 days simulation time. The difference in radius between the blast waves is 2 percent.

The difference in flux from the different EOS assumptions can be traced to the different thermal energy profiles (and hence, for fixed , to magnetic field energies that differ with the same ratios), with having the highest . This is illustrated in fig. 11. The difference in peak thermal energy densities between the fixed adiabatic index simulations is a factor of two, as expected from the ratio . Because the flux depends on the thermal energy via the magnetic field strength and (see equations 6 and 7), the flux for is higher than that for . The light curve for the advanced equation of state lies between the two limiting cases, starting close to the curve but moving to the as the flow becomes nonrelativistic. This additional decrease in flux has the consequence that the advanced equation of state light curve will be slightly steeper in the transrelativistic phase than the fixed adiabatic index light curves.

Figure 12: Left: Comparison between complete simulation light curve (solid line) at radio frequency Hz and simulation curves where is kept fixed at 1.0 (dashed line) and 0.1 (dotted line) throughout. The complete curves start out close to but slowly evolve towards . Right: same as left, only now for optical frequency Hz. As with the radio light curves, the full curve starts near but turns to .

In figure 12 we show the effects of the detailed evolution calculation of . Aside from the full simulations, we also perform two simulations that keep fixed throughout at either 1 or 0.1 but are otherwise identical to the full simulation. At early times in the radio, before the peak frequency has passed, the curve lies above the full simulation curve, where at early times in the optical it lies below.

Figure 13: Fractional difference between complete and fixed simulation light curves. Solid line for radio, dashed line for optical.

Figure 13 shows the fractional difference between complete and fixed simulation light curves (the light curves themselves lie very close to each other on a plot using a logarithmic scale), calculated via , where is the flux. The figure shows that the late time light curves for the fixed end up below the light curves that trace the evolution of the magnetic field. This can be understood from the fact that evolving according to implies a relative rise of the magnetic field energy density relative to (but still far below) the thermal energy when the flow becomes nonrelativistic. Fixing forces the magnetic field energy density to follow . The flux spans many orders of magnitude over time.

Figure 14: Power law behaviour of the optical and radio light curves. The lines plot , assuming for every two consecutive data points the relationship . The solid line refers to the radio light curve and the dashed line to the optical light curve. The horizontal lines denote 0.5, , , , from top to bottom.

A quantitative comparison between the slopes from the radio and optical light curves for the full simulations is shown in figure 14. The horizontal lines in the plot indicate expected asymptotic values for the power law scalings. In the relativistic limit, the expected slope is before passage of and after (using ). After the cooling break passes, a further steepening to is expected. In the nonrelativistic regime the expected slopes before and after passage of the cooling break are and respectively. The plot shows that the relativistic slopes are matched very well. The radio light curve quickly tends to and after passage of the peak frequency it moves in days to , where it remains until the onset of the nonrelativistic break time. The optical light curve starts out in the intermediate regime from the passage of , with the passage of the cooling break coming too early for the light curve to settle into the pre-cooling break slope of . The post-cooling break slope is obtained instead and is again maintained until the nonrelativistic break. The light curve slopes in the nonrelativistic regime are less steep than expected. A number of factors play a role here, as discussed above. The advanced EOS leads to a steepening of the decay during the transition phase (which lasts well over 10,000 days), whereas the increase in relative to and the decrease in (leading to an increase in energy per particle) lead to less steep decay. The change in slope in the nonrelativistic regime is the result of the interplay between these different factors, with the end result being a slope less steep than expected. The final nonrelativistic slopes differ significantly from those expected from analytical models, and this has a large impact on fitting models to observational data.

5 Grb030329

In the preceding section we have systematically explored the different aspects of transrelativistic blast wave afterglows with respect to dynamics and radiation for standard values of the input parameters. We now qualitatively compare radio data for GRB030329 to simulation results using physical parameters for this GRB established by earlier authors as input. GRB030329 is one of the closest and brightest GRBs for which an afterglow was found. Because of this brightness, the afterglow could be monitored for an extended period of time at various wavelengths and after six years its radio signal is still being observed (Kamble et al., 2009). GRB030329 is a good example to use to illustrate the various aspects of the radiation code.

The redshift of GRB030329 has been determined to be (Greiner et al., 2003), which leads to a luminosity distance of cm (for a flat universe with and km s Mpc). Various authors have determined the physical properties of the GRB from analytical model fits to the data (e.g. Willingale et al. 2004; Berger et al. 2003; Sheth et al. 2003; Van der Horst et al. 2005; Huang et al. 2006; Van der Horst et al. 2008) with various assumptions for the jet structure. Here we take the physical parameters established by Van der Horst et al. (2008). From their conclusion for the jet break time and cooling frequency at this time, and assuming equipartition between accelerated particle energy and magnetic field energy (i.e. a fixed in their model), we arrive at erg (for a spherical explosion), cm, , . We assume a homogeneous medium and we set the hydrogen mass fraction in this medium to unity. Van der Horst et al. (2008) fix at unity, but we use a nonrelativistic limit . Because GRB030329 shows clear evidence of a collimated outflow, it is no longer sufficient to assume a spherical explosion. When calculating emission from a jet, we assume a hard-edged jet with opening angle 22 degrees and no lateral spreading.

Figure 15: Left: Light curves at 15 Ghz, plus data. The upper curve is calculated from a spherical explosion, the bottom from a hard-edged jet with opening angle 22 degrees, while the crosses are the data. The receding jet is clearly visible. The bottom right curve is the radiation from the counterjet alone. The transition to the nonrelativistic regime can be seen from the spherical explosion simulation at around the same time when the counterjet becomes visible. Center: Light curves at 4.8 Ghz, plus data. Right: Light curves at 1.4 Ghz, plus data.

We have plotted light curves at 15 GHz, 4.8 GHz and 1.4 GHz in figure 15, including data points from the Westerbork Synthesis Radio Telescope (WSRT) (4.8 GHz and 1.4 GHz, Van der Horst et al. 2005) for comparison and the Very Large Array (VLA) (15 GHz, Berger et al. 2003). Two things are clearly visible. First, our simulated light curves still differ strongly from the data, although largely the same input parameters have been used for the blast wave simulations as those that were derived from fitting to the dataset using an analytical model for the blast wave. The different assumptions in Van der Horst et al. (2005) account for this in part, but nevertheless this demonstrates once more the need for detailed fit prescriptions from simulations (a similar conclusion was drawn in EW09 for the ultra-relativistic case). Second, the counterjet contribution will stand out clearly for a hard-edged jet model. For now, the comparison between simulation and data is still qualitative. Newer data are available and once the simulation input parameters are fine-tuned with respect to the data as well (as opposed to estimated using an analytical fit to the data), it should be possible to address the rise of the counterjet in a more quantitative fashion.

Because of the equipartion constraint on and , both were given a relatively high value of 0.27 at the shock front. In the nonrelativistic regime, the magnetic energy density will grow relative to the thermal energy density further downstream (although both will decrease strongly in absolute value). At yrs, the last time covered by our simulation (set up to cover 10,000 days in observer time) we find that has risen to approximately at the back of the blast wave ( where the blast wave density again equals the upstream density). Even further downstream, when the density has fallen three orders of magnitude below the upstream density, peaks at . This is not unphysical, but merely an indication that magnetic fields have become dynamically important in a region of the fluid which has no consequence for the light curve.

5.1 Low frequency array light curves and resolved images

Figure 16: Left: Simulated light curve at 200 MHz for GRB030329, top curve for spherical explosion and bottom curve for hard-edged jet with opening angle 22 degrees. We have drawn the following slopes from left to right: , and . LOFAR sensitivity for 25 core and 25 remote stations after four hours of integration time is 0.273 mJy and indicated by the horizontal line. Center: Simulated light curve at 120 MHz for GRB030329. LOFAR sensitivity is 0.145 mJy. Right: Simulated light curve at 75 MHz for GRB030329. LOFAR sensitivity is 4.2 mJy, too high to be shown in the plot.

Figure 16 shows predicted light curves at the very low frequencies that can be explored in the near future by radio telescopes such as the Low Frequency Array (LOFAR), assuming four hours of integration time, 25 core stations and 25 remote stations (Nijboer & Pandey-Pommiers, 2009). GRBs are among the prime targets for LOFAR’s Transient Key Project (Fender, 2006). Most of the time all light curves lie below the self-absorption break. This, in combination with the break, a hard-edged jet model and the turnover to the nonrelativistic regime, leads to an interesting double peak structure of the light curve. First the signal rises, according to the relativistic rise in the self-absorption regime that predicts a slope of . After days a clear jet break is seen and the resulting drop in slope leads to a decreasing signal again. Around circa 150 days the critical frequency passes through the observed frequency band. The slope of the spherical explosion changes accordingly towards the predicted relativistic . Around approximately 600 days the blast wave has become nonrelativistic and the counterjet starts to contribute (but is still overwhelmed by the forward jet). The predicted nonrelativistic slope for the spherical explosion is now .

We have included LOFAR detection thresholds for four hours of integration time. These sensitivity limits are higher than those presented in Van der Horst et al. (2008), because LOFAR has been scaled down in the meantime. The spherical explosion energy is an overestimation of the actual explosion energy and the flux levels corresponding to the jet simulations lie closer to what will actually be received. However, from fig. 15 it is clear that our qualitative comparison systematically underestimates the actual flux levels. Also, the integration time used in LOFAR can easily be increased, even up to days. Fig. 16 therefore does not mean that GRB030329 will not be observable by LOFAR, but only that a larger integration time than four hours is likely required.

Figure 17: Radio images at 200 MHz. Left: after days. The intensity increases monotonically outward. The outer radius is cm. Center: after days. The intensity decreases monotonically outward. The outer radius is cm. Right: after days. A central bright ring with radius cm appears. At larger radii the intensity decreases monotonically. The outer radius is cm.

For 200 MHz we have calculated spatially resolved images as well, for spherical explosions. Three images are presented in figure 17, for three different observer times. They show three qualitatively different types of behaviour. At 15 days a limb-brightened image is observed, whereas at 240 days the image on the sky becomes limb-darkened. At 3900 days another structure is visible and a brighter ring exists within the image, at a radius of cm. This is a result of the self absorption break being different for different emitting regions of the blast wave. These images are fully consistent with predictions from Granot (2007) for the ultra-relativistic case.

6 Summary and conclusions

In this paper we present the results of detailed dynamical simulations of GRB afterglow blast waves decelerating from relativistic to nonrelativistic speeds, as well as spectra and light curves calculated from these simulations using a method first described in Van Eerten & Wijers (2009) (EW09) that we have extended to include more details of synchrotron radiation. We summarize our results and conclusions below.

We have performed, for the first time, hydrodynamical simulations of decelerating relativistic blast waves using adaptive mesh refinement techniques including a parametrisation for a shock accelerated electron distribution radiating via synchrotron radiation. From these simulated blast waves we have calculated light curves and spectra at various observer times and frequencies. An advanced equation of state was used for the dynamical simulations, with an effective adiabatic index smoothly varying between the relativistic and nonrelativistic limit. Three additional parameters were traced during hydrodynamical evolution: maximum accelerated particle Lorentz factor, magnetic field energy density and accelerated particles number density. We assumed that fewer particles were accelerated by shocks that are less relativistic. To obtain the observed flux including synchrotron self-absorption, a set of linear transfer equations were solved for beams traversing through the blast wave. This method expands upon EW09 by including self-absorption and dynamically calculated electron cooling.

We have used standard assumptions for the GRB explosion energy ( erg) and circumburst particle number density ( cm) for a homogeneous medium and particle acceleration and magnetic field parameters. By directly comparing against various analytical models and expected limiting behaviour, we draw a number of conclusions about the dynamics of our simulations:

  • We find that the transition of directly behind the shock front from the relativistic to the nonrelativistic regime occurs later than expected, around days rather than days, for the standard model parameters.

  • An analytical calculation of according to Huang et al. (1999) is found to overestimate the late time values by a factor .

  • Directly applying the Sedov-Taylor solution to late time afterglow evolution is found to overestimate the radius by a few percent and keeping the adiabatic index fixed throughout the evolution of the blast wave will lead to systematic differences of as much as ten percent.

  • The density jump across the shock may be arbitrarily high for relativistic shocks, but will be a factor of four in the nonrelativistic regime. This is known from the shock jump conditions. Our simulations show that the quantity , a combination of lab frame density and Lorentz factor directly behind the shock, will remain close to four times the unshocked density throughout the entire simulation, even though the effective adiabatic index evolves from relativistic to nonrelativistic.

  • If we assume the number of magnetic field lines through the surface of a fluid element a constant, the magnetic field energy will become relatively larger compared to the thermal energy. It will remain a small fraction however (assuming only a small amount of energy is used for magnetic field creation accross the shock). Our approach allows for different assumptions on the magnetic field energy evolution.

  • The upper cut-off Lorentz factor for the shock-accelerated relativistic power law electron distribution decreases on a distance scale much smaller than the width of the blast wave due to synchrotron losses and determines the shape of the spectrum near and above the cooling break.

Using the output from the dynamical simulations, we calculate the flux. The following general conclusions are drawn for the radiation:

  • Calculated light curves show a transition between the relativistic and nonrelativistic regime at around 1000 days in observer time, again later than expected.

  • The observed fluxes for different assumptions on the equation of state may differ by a factor of a few. This is a direct consequence of the amount of thermal energy (and therefore magnetic field energy) directly behind the shock front.

  • Implementing a changing effective adiabatic index has the consequence that the resulting light curve will slowly evolve from the relativistic limiting value to the nonrelativistic value. This transition takes tens of years in observer time and will lead to a steeper decay in the afterglow light curve than predicted by analytical models assuming a fixed index.

  • This steepening is a smaller effect than the combined effect of evolving the magnetic energy density and the accelerated particle number density. When all effects are included, the final light curve slopes differ markedly from the analytically expected values. This implies a significant complication for late time afterglow modeling.

We have applied our approach to GRB030329 as well, using physics parameters derived by Van der Horst et al. (2008) using an analytical model. It is shown that the resulting radio light curves differ up to an order of magnitude between simulation and analytical model, although this can be partly attributed to some different assumptions. Assuming a hard edged jet with an opening angle of 22 degrees, our simulated light curves show a rebrightening due to the counterjet around 1000 days. Simulated curves at radio frequencies that will be observable using LOFAR show that four hours of integration time is likely not sufficient to distinguish the signal from the noise and a larger integration time is required. Finally, spatially resolved images show a bright ring that, depending on the precise power law regime that is observed, may be located not only in the center or on the edge but also at intermediate radii within the afterglow image. This is consistent with earlier work by Granot (2007) on afterglow images in the relativistic phase.

A recent paper (Zhang & MacFadyen, 2009) has appeared discussing afterglow blast waves decelerating to nonrelativistic velocities using twodimensional simulations. The authors find that lateral expansion of a relativistic GRB jet is a very slow process and that the jet break is mostly due to the edges of the jet becoming visible. This implies the hard edged jet model that we have applied to GRB030329 is sufficient to model the jet break at days. Zhang & MacFadyen (2009) do not include synchrotron self-absorption and calculate the cooling break by assuming the cooling time throughout the entire blast wave equal to the grid time.

The approach to calculating light curves and spectra from generic fluid simulations that we present in this paper assumes that synchrotron radiation is the dominant radiative process, that particle acceleration takes place in a region far smaller than the blast wave width and that the feedback on the dynamics from the radiation is negligible. We briefly address these issues in appendix B.

7 Acknowledgements

This research was supported by NWO Vici grant 639.043.302 (RAMJW, KL), NOVA project 10.3.2.02 (HJvE). ZM performed computations on the K.U.Leuven High Performance computing cluster VIC, and acknowledges financial support from the FWO, grant G.0277.08, from the GOA/2009/009 and HPC Europa (project number: 228398). We wish to thank Atish Kamble for providing the radio data to GRB030329 for the qualitative comparison and Jonathan Granot for discussion on synchrotron self-absorption.

Appendix A Numerical implementation

a.1 Partial differential equations

amrvac was written to solve a system of coupled partial differential equations. When adding additional equations to the solver, it is therefore best to use partial differential equations. In the case of the magnetic field energy we start by rewriting equation 13 as

(40)

If we multiply this equation by and add to this the continuity equation

(41)

which we first multiply by , we obtain

(42)

This is the type of conservation equation that armvac is specialized in, and it is therefore the quantity that we calculate in amrvac.

For the evolution of the upper cut-off we follow a similar procedure. We start by simplifying equation 9 to

(43)

where , and refers to lab frame time (i.e. emission time). The Lorentz factor in the source term arises when we write the comoving time derivative in the lab frame. We can now follow a procedure similar to what we did for the magnetic energy density, but first we rewrite the equation above once more for numerical reasons. The quantity varies over many orders of magnitude in a very short time span and any quantity that depends on linearly is therefore difficult to deal with numerically. We solve this by rewriting equation 43 into

(44)

Although there still is a linear dependence on in the source term, in practice this equation provides a better starting point for amrvac. From combining with the continuity equation we get

(45)

with the quantity of interest. A similar approach to tracing the effect of cooling in the context of relativistic blast waves has also been taken by Downes et al. (2002). Although in our formalism at the shock front should be reset to infinity, and therefore to minus infinity, we just take a very low value for in order to minimize numerical diffusion. In our simulations, this arbitrarily low value corresponds to a hard cut-off of the spectrum above Hz, at frequencies sufficiently far above our observation range to be of no consequence. The ‘real’ catches up with the numerical almost instantaneously.

a.2 Shock detection method

Figure 18: Normalized profiles of (solid line and squares), (dashed line and circles) and (dotted line and triangles). The actual values at the shock front are , g cm and 0.31 respectively. The simulation time is 1400 days. In the lower left corner we zoom in on the shock front, showing exactly where we reset and .

In total, amrvac now calculates the evolution of three additional quantities: (using equation 15), (using equation 42) and (using equation 45). All three quantities get reset wherever a shock is detected. Both the reset values of and depend on the fluid variables directly behind the shock front and it is therefore important that we determine the position of the shock front as accurately as possible. Mathematically speaking, a shock is a discontinuity in the flow variables with a sudden increase in entropy across the discontinuity. In practice, however, finding a shock in a numerical approximation is more involved, both due to numerical shock diffusion and because, strictly speaking, there is a shock discontinuity across every grid cell boundary.

This has the consequence that if we try to find shocks by checking for discontinuities or for entropy jumps, we will find both shocks all over the numerical diffused shock region and at a random variety of positions where the numerical noise happens to rise above a predetermined shock treshold. This then implies that we keep on resetting the additional quantities over some region, something which is especially unwanted in the case of , given our approach where we take a fluid cell to contain a collection of electrons that have been shocked exactly at the same time and we critically rely on the size of the hot region (see section 2.2 and EW09, appendix D, for details).

Because the shocked particle number density and the magnetic field density depend directly on the fluid variables, using, for example, a jump in as a trigger, as has been done by Downes et al. (2002), is not an option in these cases either. Although it serves as an excellent indicator of the front of a shock, it will not point us to a location where we can find information on the strength of the shock, but to an arbitrarily defined position just in front of that.

In this paper we solve the issues of shock detection with two shock detection algorithms, both of them making use of the fact that directly behind the shock is four times the density just in front of the shock. For and we define the shock front to be at the peak of the Lorentz factor profile, in the region where . The numerical constant is arbitrary and could be taken closer to 4. With this method we ensure that the shock is detected at those positions where the fluid quantities are sufficiently close to their peak values, although multiple shock peaks may be detected in close proximity of each other due to numerical noise.

For it is essential that we only detect a single shock front. Here we care less about the precise fluid variable values. For the purpose of resetting we define the shock front to be at that position where crosses the value . For a single shock front, this only happens once. Although, in principle, could have been used instead of , the latter offers the significant advantage that it does not change in scale over the course of the simulation and always remains close to 4, whereas becomes arbitrarily small.

Figure 18 illustrates the use of the Lorentz factor profile peak as a shock detector. It shows that the numerical diffusion is really very small and that changes over a significantly smaller spatial scale than .

a.3 Synchrotron self-absorption

Equation 18 can also be expressed as

(46)

where

(47)

Here we have used the fact that (i.e. a slightly modified powerlaw distribution). The scaling factor () of is determined in terms of and from the requirement that the total number of accelerated electrons constitutes a fixed fraction of the available electrons. The symbol denotes the pitch angle (the angle between magnetic field and particle velocity) averaged version of the synchrotron function, a dimensionless function representing the shape of the synchrotron spectrum for a single electron in the same way represents the spectrum of a distribution of particles. The in the argument is connected to via equation 7 (see also EW09).

By changing variables from to we obtain:

(48)

where the quantities and are:

(49)

and