{'fraction_non_alphanumeric': 0.035174448174143016, 'fraction_numerical': 0.021666158071406774, 'mean_word_length': 4.215490716180371, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Energy Injection in GRB Afterglow Models

We extend the standard fireball model, widely used to interpret gamma-ray burst (GRB) afterglow light curves, to include energy injections, and apply the model to the afterglow light curves of GRB 990510, GRB 000301C and GRB 010222. We show that discrete energy injections can cause temporal variations in the optical light curves and present fits to the light curves of GRB 000301C as an example. A continuous injection may be required to interpret other bursts such as GRB 010222. The extended model accounts reasonably well for the observations in all bands ranging from X-rays to radio wavelengths. In some cases, the radio light curves indicate that additional model ingredients may be needed.


INTRODUCTION
After a quarter of a century of modest progress, a breakthrough in gamma-ray burst (GRB) research came in 1997, when X-ray and optical afterglow emission was first detected from GRB 970228 (van Paradijs et al. 1997;Costa et al. 1997), showing a redshift of z = 0.695. This confirmed the cosmological origin of GRBs and to date, cosmological redshifts have been obtained for about 50 afterglows 2 .
Initially, afterglow light curves were often sparsely sampled and had a natural explanation in the standard fireball model (see e.g. , for a recent review). Many showed a smooth power-law decay followed by a steepening one to two days after the burst. This temporal behavior has been interpreted as being due to a relativistically expanding collimated outflow (e.g. Rhoads 1999, hereafter R99) or more recently to a structured jet viewed off axis (e.g. Zhang & Mészáros 2002;Rossi et al. 2002). The standard fireball model assumes a single energy release event with subsequent interaction with an ambient medium with constant density or a stellar wind density profile. It was also generally assumed that the electron energy distribution index, p, was independent of time and had values between 2-3, based on theoretical arguments and energy considerations (e.g. Zhang & Mészáros 2004).
It is quite remarkable that the light curves of most observed afterglows where sufficiently detailed data coverage exists, can be explained within the standard fireball model or slight improvements thereof (e.g. Panaitescu 2005;Zeh et al. 2006). This most likely reflects the underlying universality of the afterglow generation mechanism. For some bursts, however, interpretations with the unmodified standard model, have proved difficult due to bumps in the light curves or because some basic model parameters turn out to have values outside the range assumed in the model. In particular, the slope of the power law energy distribution of the radiating electrons, p, is found to be less than 2 in some cases (e.g. Panaitescu & Kumar 2001a). Dai & Cheng (2001) and Bhattacharya (2001) have explored the light curves and spectral signatures in models with p < 2. Numerical calculations of relativistic shock acceleration suggest a universal value of p ≈ 2.2 − 2.3 (e.g. Achterberg et al. 2001). Other work suggests an even higher value of p ≈ 2.6 − 2.7 (Lemoine & Revenu 2006). This has been referred to as the p-problem (Björnsson et al. 2002). Recent afterglow observations do indicate, however, that the value of p may not be universal (Sehn et al. 2005).
The p-problem arises when the standard model is being used to interpret light curves that decay much slower than suggested by the model. Such a slow decay had in fact been anticipated (Rees & Mészáros 1998;Sari & Mészáros 2000, hereafter RM98 and SM00, respectively) by the natural assumption that the mass ejected in each event (and the associated energy) is released with a range of Lorentz-factors rather than a single value. The distributed energy injection into the fireball sustains higher fluxes for longer periods of time, thus resulting in slowly decaying light curves, without requiring p < 2. The cause of the slow decay is thus shifted from the electron population to the properties of the source itself. This extension of the standard model was, however, not used to interpret actual data until later (e.g. Björnsson et al. 2002Björnsson et al. , 2004. As continuous energy injection can sustain slower light curve decay rates with p > 2, it is important to explore fully the properties and light curve signatures of such models. In addition to the above quoted references, energy injection in the context of GRBs has previously been considered by e.g. Panaitescu et al. (1998) who give an example of how the afterglow of GRB 970508 may be explained by refreshed shocks, although a spherically symmetric model may also work. Cohen & Piran (1999) explore the self-similarity of spherical blast waves with continuous energy injection, while Dai & Lu (2000 consider a model in which a pulsar is assumed to provide the additional energy, a variant of which was also considered by Zhang & Mészáros (2001). Kumar & Piran (2000) explore a toy model of the interaction of two colliding shells and estimate the flux contribution from both the forward and the reverse shocks. Then, Zhang & Mészáros (2002) consider the interaction of two colliding shells in detail and determine a criteria for such collisions to be mild or violent. Panaitescu (2005) considers additional energy injection in interpreting observations of several GRBs. In addition, Granot et al. (2003) and Huang et al. (2006) conclude that discrete energy injections are needed to properly account for the optical afterglow of GRB 030329. All adopt different degrees of approximations on various model ingredients.
In recent months, observations by the Swift satellite have shown that the early light curve behavior may not be as smooth as previously thought. The origin of these variations may be due to the transition from the prompt phase to the afterglow phase or due to a strong initial reverse shock (Panaitescu et al. 2006). The slower decaying portions of the light curves may be due to continuous energy injection (Nousek et al. 2005;Granot & Kumar 2006). In particular, GRB 050319 (Cusumano et al. 2006) and GRB 050315, have shown this kind of a behavior (see Chincarini et al. 2005, for the latter burst and further examples).
In this paper we present a detailed analysis of relativistic fireballs with either discrete or continuous energy injection. We show that this has interesting observational consequences in addition to the slower decay and provides clear signatures in the light curves. An example of these effects has already been presented by Björnsson et al. (2004), where a discussion of the polarization properties can also be found. The clearest of those signatures are bumps in the light curves due to refreshed shocks, which may in fact have been observed in e.g. GRB 021004 and GRB 030329. These signatures provide additional information on the hydrodynamic properties of the fireball. We give examples of fits to three GRB afterglows and show that even with the additional energy injection, there remains some discrepancy between the model and the data, especially at radio wavelengths. In our approach we try to keep the treatment as self-consistent as possible and improve on several of the approximations of some of the previous work.
The paper is structured as follows: In Section 2 we review the effects that energy injections have on the dynamics of the fireball and outline our approach to the problem. In section 3 we present approximations for simple modeling of the dynamical evolution. In Section 4 we discuss the ingredients in the radiative component of our model and in Section 5 we give examples of bursts that may be interpreted with this modified model. Finally, in Section 6 we discuss our results and conclude the paper.

THE DYNAMICS
We begin this section by considering energy and momentum conservation. We follow the approach of RM98 and SM00 and assume that energy is released over a period that is short compared to the duration of the afterglow. We then specialize the method to three cases: An instantaneous energy injection in Section 2.1 for easy comparison with R99 (see also Paczynski & Rhoads 1993), several discrete injection events in Section 2.2, and finally continuous injection with a power law distribution of the Lorentz factors in Section 2.3. We assume the dynamical evolution to be adiabatic in all cases.
Assume that the total ejected mass with initial Lorentz factor higher than Γ can be written as M (≥ Γ) = M 0 + ∆M , where M 0 is the mass ejected with the highest Lorentz factor Γ 0 , and ∆M = Γ Γ0 (dM/dΓ)dΓ, where dM/dΓ is a function specifying the distribution of ejected mass with Lorentz factor Γ. We prefer to use x = Γ/Γ 0 as a variable, and introduce a dimensionless function F (x), defined by F (x) = M (≥ Γ)/M 0 , in the interval 1 ≥ x ≥ Γ m /Γ 0 , where Γ m is the minimum Lorentz factor occurring in the ejected mass. We then have dF/dx = (Γ 0 /M 0 )dM/dΓ. Note that in general dF/dx ≤ 0, since in realistic scenarios we must have ∆M ≥ 0. Note also that F (1) = 1. The original fireball scenario with an instantaneous energy release is equivalent to dF/dx = 0.
The energy of a mass element dM is dE = ΓdM c 2 , where c is the velocity of light, and the total energy of the ejecta with Lorentz factor greater than Γ is therefore given by where E 0 = Γ 0 M 0 c 2 . Integrating by parts we find that The expanding shell decelerates by sweeping up the ambient medium. Assuming that the instantaneous total mass of the swept up medium is M e , the total energy of the shell can be expressed as E = Γ(γ i [M (≥ Γ) + M e ]c 2 ), where γ i is the average internal Lorentz factor of the particles in the expanding shell. 3 Using this and equation (2) above, we find that the conservation of energy may be written as with f = M e /M 0 and Note that the swept up mass ratio, f , increases as Γ decreases.
In the preceding equations we have assumed that the slower moving mass shells ejected by the central engine catch up with the decelerating shock front as soon as the Lorentz factors of the two coincide. This is a simplification, as the speed difference between the ejected mass elements will cause a time delay. The dynamical evolution is however not affected by this, but it should be kept in mind when considering reverse shocks originating from shells colliding with the shock front. A method for including this time delay will be discussed in section 2.2.
Before impact with the expanding shock front, each mass element has momentum dP = ΓdM βc, where βc is the velocity of the mass element in the burster frame. Hence the total momentum before impact is with P 0 = Γ 0 M 0 β 0 c. As the momentum of the hot shell is P = Γγ i [M (≥ Γ) + M e ]βc, momentum conservation results in where Dividing equation (6) by equation (3), we obtain an expression for the dimensionless velocity of the expanding fireball, We then find that the Lorentz factor is given by This is an implicit equation for Γ, and describes the evolution of the Lorentz factor with radius (or time), once the distribution of the injected Lorentz factors, F (x), has been specified along with the initial conditions of the fireball (e.g. Γ 0 and M 0 ). A description of the swept up mass, f , is also required (see eq. [10] below). Equation (9) is the key equation for obtaining the dynamical evolution of the fireball and has to be solved numerically in all realistic cases. Since we take the evolution to be adiabatic, the dynamics can be solved for separately and the radiative properties computed afterwards. For a complete description of the fireball evolution we must solve four differential equations (R99) Here, r is the radial coordinate in the burster frame and t b , t ′ and t are the burster frame time, the comoving frame time and the observer frame time, respectively. The solid angle into which the energy is beamed is estimated by Ω = 2π(1 − cos(θ 0 + c s t ′ /(ct b ))) where θ 0 is the initial opening angle of the collimated outflow and c s ≈ c/ √ 3 is the sound speed in the comoving frame, assumed to be relativistic. The mass density of the external medium is ρ, which in general is a function of radius.

Instantaneous Energy Injection
In the standard fireball scenario it is assumed that the energy is released instantaneously with x = 1 and F (x) = 1. We then find that xF (x) + I 1 = 1; I 2 = 0 and equation (9) reduces to which agrees with the result of Paczynski & Rhoads (1993) and is the starting point of the R99 analysis (see his eq. [4]). Some authors have based their analysis of GRB afterglows on this equation (e.g. Salmonson 2003;Bianco & Ruffini 2005). GRB 990510 is an example of a burst where the afterglow may be readily interpreted with a single instantaneous energy release event. We present our fit to the observations in Section 5.1.

Discrete Energy Injection
Assume that in addition to the initial injection M 0 , the central engine expels several shells with masses M i and Lorentz factors Γ i . Converting to dimensionless variables, x i = Γ i /Γ 0 and m i = M i /M 0 we have that where is the Heaviside step function. We then obtain, xF (x) . Note that by the assumptions of our approach discussed in the first two paragraphs of Section 2, the x i 's form a decreasing sequence. Inserting equation 12 into equation (9) we obtain Because the Heaviside function is zero until x i ≥ x, it is clear from this equation that the dynamical evolution between injections is similar to the instantaneous case as the evolution of Γ depends only on f in the same way. When a shell collides with the shock front, a reverse shock may propagate back through the shell if the collision is violent (Zhang & Mészáros 2002). Depending on the circumstances, radiation from the reverse shock may make a substantial contribution to the afterglow light (Nakar & Piran 2004). According to Zhang & Mészáros (2002), the liklihood of getting a reverse shock depends on the relative energy and speed between the shells and the width of the incoming shell. As the energy is known and we neglect the dynamics in the shell collision, only the relative speed is needed. Assuming that only the shock front slows down due to interaction with the surrounding medium and that the trailing shell travels with constant speed, shells ejected with the average Lorentz factor of the shock, Γ a , will have reached the shock front at r in time t. Replacing x with x a = Γ a /Γ 0 in equation (9) is then sufficient to account for the delay. As shown by Kumar & Piran (2000), the Lorentz factor of the incoming shell is Γ a ≈ 2Γ in a constant density environment for most of the evolution. Combined with the results from Zhang & Mészáros (2002) (see their figure 3) this indicates that for a reverse shock to result, the energy of the incoming shell must be at least twice the energy of the shock front.
We emphasize that in our numerical calculations we assume the shell collisions to be instantaneous and we neglect the effects of the interaction on the light curves. These effects are expected to be small considering that the relative Lorentz factor of the collision is Γ r ≈ 1.25 (Kumar & Piran 2000). We estimate that the shells would have to be at least several hundred light seconds thick for the effects of the interaction to be measurable with current instruments; these would also be strongly smoothed because of the equal arrival time surface (EATS).
An example of an afterglow that may be interpreted with discrete energy injections is that of GRB 021004. Björnsson et al. (2004) showed that its properties, including polarization, could be explained by using such a model. Energy injection does affect the polarization levels calculated using the model of Ghisellini & Lazzati (1999) as these are most strongly determined by the Lorentz factor. Changes in the Lorentz factor therefore manifests itself not only in the light curve but also in polarization data. An analysis of a more detailed data sets for GRB 021004 has been presented by Postigo et al. (2005). The late time afterglow and the host galaxy of this burst have been studied by Fynbo et al. (2005).
The afterglow of GRB 030329 also showed variations both in the optical and the polarization light curves. These have previously been interpreted with energy injections , although other alternatives have also been explored Huang et al. 2006). An application of our approach to that particular burst will be presented elsewhere (Guizy et al. in preparation).

Continuous Injection
Here, we follow the approach of RM98 and SM00 and assume that F (x) = x −s , with s > 1, i.e. the energy is dominated by the lowest Lorentz factors. (The case s < 1 is equivalent to instantaneous injection as the energy is then dominated by the highest Lorentz factors). Then xF (x) To evaluate I 2 , we assume that Γ 2 0 x 2 − 1 ≈ Γ 0 x, which is always a good approximation if Γ 2 m ≫ 1. We then find I 2 ≈ −Z. The dynamical equation for Γ (eq. [9]) in this case becomes Note that for Z = 0, we recover the result for the instantaneous case (eq. [11]). Also note that the second term under the square root is ≈ 0 when Γ 0 ≫ 1. When considering models with continuous injection, we numerically solve equation (16), since in all relevant cases Γ 2 m ≫ 1.

SIMPLE APPROXIMATIONS
The quantities that determine the dynamical evolution of the fireball are the dimensionless mass ratios F (x) and f , representing the injected energy and the swept up mass, respectively (see eq. [9]). F (x) is an input quantity in the model, and hence each case has to be solved separately. On the other hand, f is only dependent on the evolution of r, ρ and Ω. In general, ρ ∝ r −g , where g = 0 for constant density interstellar medium (ISM) and g = 2 for a wind like environment. The evolution of f can easily be obtained in two cases: the confined outflow case when c s t ′ ≪ θ 0 r and Ω may be considered to be a constant. In this case the dynamics is that of a spherical outflow. The second case is that of a laterally spreading outflow, when c s t ′ ≫ θ 0 r and Ω is large enough to neglect the initial opening angle. These cases were referred to as the power law regime and the exponential regime, respectively, in R99. That naming convention is however too limiting as we show below, since for a continuous energy injection there is no exponential regime.
As a prescription for f in the confined case, we evaluate it as while in the laterally spreading case we adopt (see R99) To obtain a full solution of the problem and the light curve behavior, these need to be consistently included in the calculation. Approximate solutions applicable in the instantaneous case have been discussed in detail by e.g. R99 and Salmonson (2003). We emphasize that the approximations can give useful insight into the properties of the fireball but they should, however, only be expected to be valid during limited time intervals in the afterglow evolution. This is partly due to the validity of the underlying assumptions in those limited intervals, and partly due to the influence of the EATS on the evolution (see Section 4.1). Another factor is the spectral energy distribution (SED), that is not really a set of broken power laws (see Appendix A), but has an intrinsic curvature, and so do the lightcurves. In the following we first give a brief discussion of the discrete case but then concentrate on the case of continuous energy injection.

Discrete injections
If the energy is released in several discrete events (as discussed in Section 2.2), the evolution after a shell impact can be approximated with the evolution following a single release event with the appropriate initial parameters. Of these, only the energy is relevant, being equal to the total energy that has accumulated in the shock front after the impact. Integration over the EATS then smooths out the transition.
As mentioned in Section 2.2, the relative speed between the impacting shell and the decelerating shock front is an important factor in evaluating the reverse shock emission arising from the collision. The speed of the incoming shell equals the average speed of the shock front at the time of collision, and to estimate it, we need to know how the radius of the fireball, r, evolves with time in the burster frame, t b . Using Γ ≈ Γ 0 /2f, appropriate for a single release (see R99), and equations (17) and (10), it is easy to show that in the confined case the speed of the incoming shell is Γ a ≈ (4 − g)Γ. Although this should only be valid early on in the evolution, numerical calculations indicate that it is in fact appropriate during most of the evolution.

Continuous energy injection
All results derived in this section assume that the shock speed, β ≈ 1. We first note that as f is the swept up mass in units of the initial ejected mass, it is small compared to the other terms in equation (16) for a substantial fraction of the initial fireball evolution, but will dominate at later times.
Early in the evolution of the fireball, when x is close to 1, Z < Γ 0 and we recover equation (4) in R99, and an identical initial fireball evolution. The most interesting case, however, is when Z/Γ 0 > 1. From the definition of Z we see that this occurs when x = Γ/Γ 0 < (s/(2s − 1)) 1/(s−1) , i.e. when x < 2/3 for s = 2 and for x < 3/5, when s = 3. In this case equation (9) simplifies to When both Γ 0 and f can be neglected compared to Z, we have that Γ ≈ Z/2f . Dividing by Γ 0 and solving for x = Γ/Γ 0 , we find Note that in the limit, s → 1, we recover the power law behavior of the instantaneous injection case, where x ∝ f −1/2 (R99).
Laterally spreading outflow: The treatment is similar to that of R99. The g = 2 case cannot, however, be solved analytically due to the coupling of the variables. We therefore only consider the case g = 0 in this section (constant density ISM). The general case has to be solved numerically.
We start from equation (18), divide by dt ′ /dr and integrate to obtain in agreement with equation (14) of R99 if we set s = 1. Here, c 1 is a constant of integration to be determined by the initial conditions. It becomes negligible when c s t ′ ≫ θ 0 r. The exponential behavior discussed by R99 is seen to follow from equation (24). If s = 1 (formally equivalent to the instantaneous case), then equation (24) shows that f ∝ t ′ 2 so that dt ′ /dr ∝ t ′ , and the time in the comoving frame increases exponentially with r. Therefore f also increases exponentially with r.
On the other hand, if s > 1, the behavior becomes a pure power-law, f ∝ t ′ 3(1+s)/(2+s) and t ′ ∝ r (2+s)/(s−1) which shows that the exponential regime is a special case. It follows that f ∝ r 3(1+s)/(s−1) , r ∝ t (s−1)/(5+s) and Γ ∝ t −3/(5+s) . The physical reason for this change in behavior is that the dynamical evolution is now determined by the power-law energy injection rather than the lateral spreading. We note that in general, if the incoming shell has the same opening angle as the initial θ 0 , the resulting decay will be a mixture of a power-law decay for the part of the outflow within θ 0 , and an exponential for that part of the flow that is outside of θ 0 .
By assumption, the dynamical evolution is independent of the radiative properties of the fireball. To proceed, a prescription of the radiative processes is required. We adopt the standard synchrotron emission mechanism (see e.g. , for a review), that is widely assumed to be the source of radiation in GRBs (e.g. Tavani 1996;Mészáros & Rees 1997;Sari et al. 1998). The details of the radiative component of our modeling are given in Section 4 and in Appendix A.
To explore the approximate light curve behavior, F ν ∝ t −α ν −β , we need to evaluate the characteristic frequencies ν m , ν c and the flux F νm at the former (See Appendix A for definitions of these quantities).
For the confined case, the temporal behavior of the characteristic frequencies, the F νm , and the light curve decay indices, α, agree with those previously derived for the forward shock by SM00 (see Table 1 in their paper).
For the laterally spreading case we find on the other hand F m ∝ t (3s−9)/(s+5) .
(29) The corresponding light curve exponents are given in Table 1. (Note that if we let s = 1 we recover the corresponding results for the instantaneous case). From the indices in Table 1 and the results of SM00 we see that if ν m < ν < ν c , the steepening in the light curve when the observer starts to see the entire outflow surface is ∆α = α 2 − α 1 = 24(2 + β)/((s + 5)(s + 7)), with β = (p − 1)/2. For parameters typically inferred in afterglows (p = 2.5 and β ∼ 0.75) and taking s = 2 as an example, we see that ∆α ∼ 1. When ν > ν c the result is ∆α = 24(1 + β)/((s + 5)(s + 7)), with β = p/2 in this case. For typical parameters we now find ∆α ∼ 0.85. As the steepening of the light curve is similar in the two cases, the magnitude of the break does not indicate in which spectral region the observing frequency is. Clearly, a more detailed fit to the observations is needed. Inferring p < 2 from a standard model fit may be taken as a strong indication that a continuous injection may be involved.

LIGHT CURVES AND SPECTRA
As is customary in GRB fireball models and discussed in the previous section, we assume the radiation to be of synchrotron origin. We adopt the approach of Wu et al. (2004), although their method of calculating the radiation flux by integrating the synchrotron emissivity over the electron energy distribution is however too computer intensive for our purposes. Instead, we use a broken power law approximation as in Wijers & Galama (1999), smoothly joined together similar to Granot & Sari (2002). The detailed expressions are summarized in Appendix A.
Similarly to Wu et al. (2004), we calculate the synchrotron self absorption coefficient directly and therefore the self absorption frequency is not included in our smoothly joined power law approximation. We also include inverse Compton scattering, but neglect higher orders as these are suppressed by the Klein-Nishina effect (Panaitescu & Kumar 2000). The Compton Y parameter is calculated from where N e is the normalized electron distribution and τ e is the optical depth to electron scattering. The modified cooling Lorentz factor is found from the implicit equation where, γ c , is the usual cooling Lorentz factor (see eq. A4). The Comptonized emission component is then found from the synchrotron spectrum using the scalings, ν IC m = 2γ 2 m ν m , ν IC c = 2γ 2 c,IC ν c,IC and P IC max = τ e P max (Sari & Esin 2001).

The emitting region
In the standard afterglow model the fireball is described by the self-similar formulation of Blandford & McKee (1976) (hereafter BM). Analytical solutions are however only available if the ambient medium has a density profile of the form n ∝ r −g . Some authors, Lazzati et al. 2002), have adopted these solutions in models with modest density variations. Using their approach, we are able to reproduce their results only if we assume that the shock thickness follows the BM solution, ∆ ∼ r/Γ 2 . Care must be exercised here, since in this approximation the number of radiating particles may not be estimated correctly. To clarify, let us consider a fireball in a constant density ISM. Using the Blandford & McKee (1976) solution for the density behind the shock front and integrating over the swept up volume we find that the number of particles contained behind the shock is where r is the the radius of the swept up volume, and n is the (constant) ISM particle density. This result is appropriate for a homogeneous medium but leads to over-or underestimates, depending on the sign of the density gradient, in the number of radiating particles when used for variable ISM densities. In such cases it leads to incorrect flux estimates. We also remind the reader that density variations do affect the evolution of Γ which has to be calculated self-consistently to obtain correct fluxes. Instead of using the BM solution directly, we make the simplifying assumption that the shock front is homogeneous in the comoving frame. Using particle conservation and the jump conditions we find that the shock thickness in the burster frame is ∆ = M e 2πm p r 2 (1 − cos θ)Γn ′ , appropriate for all density profiles. Here, m p , is the proton mass and, n ′ = 4Γn, is the density in the comoving frame. Note that when (1 − cos θ) = 2, and M e = 4πr 3 n/3, we recover the BM shock thickness. Although apparently a crude estimate, Granot et al. (1999) showed that the difference between the calculated flux levels in such a thin shell approximation and the full solution are small. In all our flux calculations, we follow the evolution of the EATS. The total flux arriving at the observer at a given time, t, is obtained by integrating over the EATS for a given observer frequency, ν, (Granot et al. 1999): where, P ′ , is the radiation power density in the comoving frame (see also appendix A), ν ′ is the frequency in the comoving frame and is the time of emission in the burster frame. Also, d l , is the luminosity distance which depends on the cosmology. In our numerical calculations we assume a standard cosmological model with Ω Λ = 0.7, Ω m = 0.3 and H 0 = 70 km s −1 Mpc −1 .

Energy injection effects on radiation properties
To understand the effects energy injections have on the radiation properties of the afterglow it is best to start by considering the effects of a single injection. Since we neglect the dynamics of shell collisions, the Lorentz factor of the leading shell increases instantaneously as soon as another shell catches up with it. The subsequent evolution will then quickly relax to a solution which is similar to the one obtained if the total energy of both shells was injected instantaneously. As the energy of the shock front increases instantaneously, so does its luminosity. However, the increase in flux received by the observer will not be instantaneous but gradual because the integration over the EATS smooths out the transition. The timescale of this gradual brightening depends highly on both the Lorentz factor and the jet opening angle . Energy injection in wider and slower moving jets will cause the flux to change on longer timescales and thus creates smoother bumps in the light curves. Shell interactions considered in more detail than in this work, would also contribute to the smoothening.

EXAMPLES OF DATA MODELING
In this Section we apply our model to observations and give three examples of how well our numerical results fit afterglow data. We use reduced χ 2 deviation to assess the goodness of the fit. As our model extends to the entire electromagnetic spectrum, an un-weighted χ 2 is not adequate for our purposes. The number of points in the optical light curves often exceeds that in other wave bands by a factor of ten or more which tends to decrease the quality of the fit at the radio and X-ray wavelengths. To correct for this, we weight the contribution from each band with the square root of the number of data points in that band. We then normalize the result so that if the contribution from each point to the un-weighted χ 2 is unity, the weighted χ 2 equals the total number of data points. The function we use to check the goodness of the fit is then where N dp is the total number of data points, N lc is the number of light curves, N i is the number of data points for light curve i, N par is the number of parameters, F i,j is the flux in light curve i, measured at time t j with the error ∆F i,j , and F i (t j ) is the flux from our model light curve i at time t j . Dividing by the number of degrees of freedom (d.o.f.) in the first factor of the equation provides the reduced χ 2 which we denote by χ 2 d.o.f. . The second factor is the normalization, and 1/ √ N i in the sum over light-curves (lc) is the weight. Note that the fit is done for log(F ) rather than F , as this is equivalent to a liner fit in magnitudes.
Fitting the model calculations to the measurements is nontrivial since the number of model parameters can be large and the goodness of fit function, χ 2 d.o.f , has multiple local minima. This makes local gradient based search methods unsuitable for our purposes. Instead we employ a global search method known as the Evolution Strategy (Runarsson & Yao 2000). The essence of the algorithm is as follows: Initially, a population of randomly generated λ search points (parameter sets) are generated within the specified parameter bounds and the initial variance of the normally distributed search distribution is set. These λ points are then evaluated using the χ 2 d.o.f function and sorted. The best µ of those are then used as parent search points and each one of them generates λ/µ new points within their Gaussian search distribution. However, before doing so the search distribution itself, or rather its standard deviation, is varied using the log-normal distribution. In this manner the most favorable search distribution is adapted to the function landscape resulting in faster search. Initially the search distribution reaches the entire search space. The evolutionary process is repeated for a number of iterations until the search distribution's standard deviation is below a pre-specified limit.
In our fits, we did not allow the g-parameter (the ambient density profile) to vary, but assumed it to be fixed at either g = 0 (homogeneous ISM) or at g = 2 (wind). In the following sections, we set g = 0, as this provides a better fit in all three cases considered here. In addition, the angle that the line of sight makes with the jet axis (the viewing angle) is assumed to be zero. This is because it mainly affects the polarization properties of the afterglow that are not considered here.
The complex landscape of the goodness function, caused by the non linearity and complexity of the underlying model, makes it difficult to find a unique parameter set that is a well defined global minimum of the goodness function. Degeneracy of the parameters can cause several minima to be considered as candidates for the global minimum. The 1 − σ errors we quote in our parameter determinations reflect the statistical uncertainty of the model parameters as determined from the measurements and do not include estimates accounting for uncertainties in the model (see also Panaitescu 2001, for a similar discussion). The results of our fits for three GRB afterglows are given in Table 2. Several comments on the results are in order. As the afterglow light curves are weakly dependent on Γ 0 , it is not constrained by the data. We only quote lower limits on it. The electron energy slope, p, is constrained both by the spectral slopes and the afterglow slopes and is the most constrained model parameter. In the case of GRB 000301C we only quote an upper limit on ǫ B , as it is poorly constrained, mostly due to the lack of X-ray data. Finally, we note that the 1 − σ errors quoted in the table are generally not symmetric around the best fit values.
The afterglow light curves of GRB 990510 fitted with a model with an instantaneous energy release. The multi-band fit is reasonably good, resulting in χ 2 d.o.f. = 2.6 with the model parameters given in table 2. The X-ray data is from Kuulkers et al. (2000), the optical data from Harrison et al. (1999); Stanek et al. (1999) and the radio data from Harrison et al. (1999). The optical data is corrected for Galactic extinction of E(B − V ) = 0.2 (Schlegel et al. 1998). 2000). This afterglow is a good example of smooth lightcurves and can be interpreted within the standard fireball model (Panaitescu & Kumar 2001b). It is therefore an ideal burst to test the accuracy of our model and to compare its results to those obtained by others. Figure 1 shows the afterglow light curves in radio, optical and X-rays. The fit resulted in χ 2 d.o.f. = 2.6, quite an acceptable fit given the simplifying assumptions for the underlying dynamics and radiation and the fact we are fitting the entire electromagnetic spectrum with a single model. The best fit parameters are given in table 2. Most of them are similar to the parameters found by Panaitescu & Kumar (2001b), within a factor of two or three. The order of magnitude difference in ǫ e , the fixed fraction of the total internal energy assumed to be in the electron population, stems from our different definition of γ m . Although not shown here, the slope of the light curves vary continuously from α = 0.6 at 0.1 day to about α = 2.2 at 30 days. A rapid steepening sets in at 1 day after the burst but there is no sharp transition from a pre-to a post-break slope. We note that for this burst, θ ≈ 1/Γ at 0.7 days, just before the steepening of the light curve.
GRB 990510 was the first burst where a polarized afterglow was detected (Covino et al. 1999). However, as there was only one positive detection it does not constrain the model very much, although it further strengthens the case for a synchrotron origin of the afterglow emission. This polarization measurement is easily accounted for by an off-axis viewing angle. It was not included in the fit presented here.

GRB 000301C
The afterglow of GRB 000301C is well time resolved at optical and radio wavelengths and was observed in the near infrared and at mm wavelengths although no X-ray observations have been published. A short lived brightening was d.o.f. = 3.8 using all available data and is a reasonably good fit. The energy injection allows for both the sharp break in the optical/nIR bands and the smooth decline of the radio light-curve. The bump around 4 days seen in the R band causes some of the contribution to χ 2 d.o.f. along with the radio data. It rises and decays too sharply to be fitted with an energy injection. The model parameters are given in Table 2. The dotted curves show the best fit model without additional energy injection, resulting in χ 2 d.o.f. = 11. See text for further discussion. Optical and near infrared data are from Jensen et al. (2001); Rhoads & Fruchter (2001); Bhargavi & Cowsik (2000); Masetti et al. (2000) and Sagar et al. (2000). Radio and mm data are from Berger et al. (2000); Smith et al. (2001) and Frail et al. (2003). All optical data is corrected for Galactic extinction of E(B − V ) = 0.053 (Schlegel et al. 1998) and intrinsic SMC like extinction of A B = 0.1 using the profile of Pei (1992). A host contribution of 40 ± 20 µJy has been subtracted from the 8.46 GHz data. observed in the optical light curves at about 3.5 days followed by a steep achromatic break around 7 days after the burst (Sagar et al. 2000;Berger et al. 2000). It has been suggested that this bump may be caused by micro-lensing (Garnavich et al. 2000). Their χ 2 value of the fit is reduced by half by introducing the lens effect.
Here, we fit the afterglow light curves with our model containing one additional energy injection, E 1 at t 1 = 1.4 days post-burst. This is shown in figure 2, where a host contribution of 40 ± 20 µJy at 8.5 GHz has been subtracted from the data. The data in the figure is corrected for Galactic extinction of E(B − V ) = 0.053 (Schlegel et al. 1998), and an intrinsic SMC like extinction of A B = 0.1 using the profile of Pei (1992). This value for the intrinsic extinction was found by Jensen et al. (2001) which fitted a power law spectrum with various extinction curves to the measured SEDs.
Fitting the data with a model with just an initial instantaneous energy release results in χ 2 d.o.f. = 11 (dotted curves in fig. 2). Adding the energy injection E 1 , reduces this by a a factor of three to χ 2 d.o.f. = 3.8, and gives the parameters presented in table 2. Most of the reduction in the χ 2 d.o.f is due to the optical light curves, as the energy injection allows for a fairly constant emission at 8.5 GHz and the sharp steepening of the optical light curves at approximately 6-7 days (solid curves in fig. 2). The remaining contribution to the χ 2 d.o.f is due to the bump at 3.5 days in the optical light curves, that rises too sharply to be fitted with an injection event. This is because of the EATS and the homogeneity of the emitting region in our model, making the flux increase too smooth to fit the sharp rise and decline observed. This sharp bump may be due to an injection that is local in the sense that it only affects a small patch of the radiating shock front. We also note that the best fit is obtained by having the energy injection at t 1 = 1.4, i.e. just before the first optical data points, rather than at times just before the light curve bump. Such an early injection does fit the flux levels in the various bands better than an injection that is constrained to lie within the time window of the observations.
Our model parameter values agree with those of Panaitescu (2001) within a factor of 2 or 3, except for quantities that are defined differently (such as ǫ e ). Given the differences in approach and the different level of details in his work and ours, these differences are rather minor.
Finally, although energy injection that is uniform across the forward shock front is unable to account for the bump in the lightcurves at 3.5 days, other explanations may apply besides micro-lensing. These could be density inhomogeneities in the ambient medium or an inhomogeneous injection that energizes the emitting surface only partly (patchy shell).

GRB 010222
GRB 010222 had well time resolved afterglow light curves at most wavelengths and can be interpreted with the standard broken power law in the optical bands. The shallow slopes of the light curves require a hard electron energy distribution (p < 2), if interpreted within the standard model. An unknown host extinction may influence the SED slope that in turn may lead to an inconsistent determination of p, between the spectra and the light-curves.  Sagar et al. (2001) and Stanek et al. (2001). X-ray date from in't Zand et al. (2001) and radio and millimeter data from Frail et al. (2003) and references therein. The dotted curves show for comparison a fit with only instantaneous injection, resulting in χ 2 d.o.f. = 6. The increase in χ 2 d.o.f. is mostly due to contributions from the radio light curve. Björnsson et al. (2002) found that a continuous energy injection in an otherwise standard fireball model could explain both the light curves and SEDs without a hard electron energy distribution. We follow their suggestion and use continuous energy injection model to fit the afterglow data. The result is shown in figure 3 where the data points have been corrected for Galactic extinction of E(B −V ) = 0.018 (Schlegel et al. 1998) and intrinsic SMC like extinction of A B = 0.1 using the profile of Pei (1992), that was also adopted by Galama et al. (2003). The host of GRB 010222 is a bright starburst galaxy (Frail et al. 2002) and where available, the host contribution has been subtracted from the light curves.
The fit resulted in χ 2 d.o.f. = 4.6 with the model parameters presented in table 2. The fit is not particularly good, mainly because the model gives too high a flux at the radio wavelengths. Both the optical and the X-ray data are well accounted for. We note that the value of Γ m is rather high (table 2) and the continuous energy injection ceases at approximately 0.1 day, which along with the jet geometry (Γ ≈ 1/θ at about 3 days), accounts for the steepening in the light curves. If we fit the data without the continuous energy injection, we obtain χ 2 d.o.f = 5.8 instead of 4.6, with the additional contribution mainly due to a worse fit at radio wavelengths. The model parameters, however, are similar in both cases; a narrow and energetic jet with a low magnetic field strength in a low density environment. The continuous energy injection thus better accounts for the radio light curves but otherwise the evolution in both cases is similar. We are therefore able to fit the optical data reasonably well within the standard model and still have p > 2.

DISCUSSION
We have presented a detailed description of a relativistic fireball afterglow model with discrete and continuous energy injection. We have given explicit examples of fits to observations for both types of injections and shown how they can provide a better fit to afterglow measurements than the standard model with single initial energy release.
Although not a direct subject matter of the present paper, we note that even polarization measurements can in some cases be accounted for with the extended model (Björnsson et al. 2004). However, interpreting afterglow polarimetry is not straightforward. First, polarization measurements are difficult as the polarization level in general seems to be small in afterglows. The data may therefore be sensitive to polarized contributions from other sources such as the host galaxy, which however is most likely to be constant. Second, our model is sensitive to the geometry of the outflow and the assumed homogeneity of the emission region. The expanding shock front is unlikely to be homogeneous during its entire evolution, especially with the added energy injections. This is because an injected shell will not expand laterally in exactly the same way as the original shock front. The transverse size of the incoming shell will then be smaller and its energy deposition in the shock front will be concentrated in the contact zone between the two. Accounting for the effects of the EATS, could then lead to a rise in the level of polarization shortly after the injection and even abrupt changes in the polarization angle if the center of the colliding shell does not coincide with the geometrical center of the front.
The origin of the additional shells is most likely the same as the origin of the shells producing the gamma ray burst in internal shocks. It is highly unlikely that the matter ejected from the central engine always collides and collects into one shell in the gamma ray burst event. The remaining shells will then be the ones to refresh the external shock later on in the evolution of the fireball. This might even be the case in most bursts, but goes unnoticed in many cases since the energy of the colliding shells must be of the same order or higher as the total accumulated energy of the front if their interaction is to produce a measurable brightening.
When fitting GRB afterglow light curves with the standard fireball model, it is quite common to use power law approximations with a light curve break due to jet structure. The slope of the light curve, α, before and after the break is then used to determine the electron energy distribution index, p, which is then compared to the value of p inferred from the SED of the burst (see e.g. Sari et al. 1998). These results must, however, be carefully considered, as the value of α depends on s as well (see table 1), and may not be smoothly varying between two fixed asymptotes, although that seems to be the general case. This is especially true when the characteristic frequencies of the synchrotron spectrum, ν m and ν c , are near the observing frequency, where the SED does not follow a power law.
The light curve break time has traditionally been used to infer the jet opening angle from the assumption that it occurs when θ ∼ 1/Γ. This can lead to erroneous angle estimates, because energy injections can delay the appearance of the light curve break. An example of this is GRB 021004 (Björnsson et al. 2004). As a result, the opening angle will be overestimated and the energy corrected for beaming may thus not be correct (Johannesson et al. 2006). Studies based on such energy estimates (e.g. Frail et al. 2001;Bloom et al. 2003;Ghirlanda et al. 2004), should include error analysis that accounts for this 'bias'.
We thank P. Jakobsson for critical comments on the manuscript and T. P. Runarsson for helpful discussion on the fitting procedures. We thank the anonymous referee for detailed and constructive criticism. This work was supported in part by a Special Grant from the Icelandic Research Council, by the University of Iceland Research Fund and by the Graduate Research Fund of the Icelandic Centre for Research. results from these approximations with the more accurate calculations using the approach from Wu et al. (2004), we find the following p-dependence of the exponents and coefficients in the fast cooling regime to be κ 1 = 2.37 − 0.3p (A8) κ 2 = 14.7 − 8.68p + 1.4p 2 (A9) φ p = 1.89 − 0.935p + 0.17p 2 (A10) χ p = 0.06 + 0.28p.
The difference between these approximations and the more accurate calculations is generally less than 10% unless ν c is close to ν m . Then the difference can be up to 40%. However, for most parameter sets the characteristic frequencies are close, only very early in the afterglow evolution and this is therefore not a concern in the later evolution. Even then, the difference is mostly smoothed out because of the integration over the equal arrival time surface. We find the difference between the two solutions to be less than 10% for most realistic cases. We also remind the reader that the electron energy distribution is an approximation and little is known about the detailed microphysics of the shocked material.