format: text

text:
# 1 Introduction

## 1 Introduction

The physics of high-density QCD has become an increasingly active subject of research, both from experimental and theoretical points of view. In deep inelastic scattering (DIS) the parton high density regime corresponds to the small $`x`$ region and represents the challenge of studying the interface between the perturbative and nonperturbative QCD, with the peculiar feature that this transition is taken in a kinematical region where the strong coupling constant $`\alpha _s`$ is small. By the domain of perturbative QCD we mean the region where the parton picture has been developed and the separation between the short and long distance contributions (the collinear factorization) is made possible by the use of the operator product expansion (OPE). The Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations  are the evolution equations in this kinematical region. These equations are valid at leading twist, i.e. at a large value of the photon virtuality $`Q^2`$, where a subclass of all possible Feynmann graphs are dominant and the hard coefficient function is connected to the proton by only two parton lines (For more details see e.g. Ref. ). For small values of $`Q^2`$, this picture has corrections predicted by the OPE that contribute at relative order $`𝒪(1/Q^2)`$ and beyond $`[𝒪(\frac{1}{Q^2})^n,n=2,3,\mathrm{}]`$. These are commonly called higher twist corrections.

In the limit of small values of $`x`$ $`(<10^2)`$, on the other hand, one expects to see new features inside the nucleon: the density of gluons and quarks become very high and an associated new dynamical effect is expected to stop the further growth of the structure functions. In particular, for a fixed hard scale $`Q^2\mathrm{\Lambda }_{QCD}^2`$, the OPE eventually breaks down at sufficiently small $`x`$ . Ultimately, the physics in the region of high parton densities will be described by nonperturbative methods, which is still waiting for a satisfactory solution in QCD. However, the transition from the moderate $`x`$ region towards the small $`x`$ limit may possibly be accessible in perturbation theory, and, hence, allows us to test the ideas about the onset of nonperturbative dynamics.

At this moment, there are several approaches in the literature that propose distinct evolution equations for the description of the gluon distribution in high density limit  . In general these evolution equations resum a class of higher twist diagrams dominant in this kinematical region that contain powers of the function $`\kappa (x,Q^2)\frac{3\pi ^2\alpha _sA}{2Q^2}\frac{xg(x,Q^2)}{\pi R_A^2}`$, which represents the probability of gluon-gluon interaction inside the parton cascade. Moreover, these equations match (a) the DLA limit of the DGLAP evolution equation in the limit of low parton densities $`(\kappa 0)`$; (b)the GLR equation and the Glauber-Mueller formula as first terms of the high density effects. The main differences between these approaches occurs in the limit of very large densities, where all powers of $`\kappa `$ should be resummed. Although the complete demonstration of the equivalence between these formulations in the region of large $`\kappa `$ is still an open question, some steps in this direction were given recently .

Our goal in this paper is not to demonstrate the equivalence between the distinct approaches, but to analyze the relation between the gluon distribution $`xG`$ and the $`F_2`$ structure function in the asymptotic regime of very high densities (energies). In general, a leading twist relation is used to relate these quantities and the solutions of the evolution equations used as input in the calculations. However, as discussed above, the high density approaches resum a class of higher twist diagrams, which implies that we should careful in our predictions for the quantities which would be measured. Here we study the behavior of $`F_2`$, predicted by a leading and a higher twist relation, at large densities, where analytical calculations are possible. In this regime we expect the blackness of the cross section, as predicted by Gribov many years ago . We verify that this behavior is recovered by the current high density approaches, if a leading twist relation between $`F_2`$ and $`xG`$ is used. However, a softer behavior is obtained when we consider that this relation is modified by the higher twist terms associated to the high density, demonstrating the importance of this contribution. In both cases, the Froissart boundary is not violated.

The paper is organized as follows: in the next section we briefly review the derivation of the leading and higher twist relations between $`F_2`$ and $`xG`$, as well as the black disc limit. After, in section 3, we derive the asymptotic solution of the AGL equation at fixed and running $`\alpha _s`$ and compare this solution with the predictions of the McLerran - Venugopalan high density approach. Finally, in section 4, we present our results and conclusions.

## 2 The $`F_2`$ Structure Function and the Black Disc Limit

We start from the space-time picture of the $`ep`$ processes . The deep inelastic scattering $`epe+X`$ is characterized by a large electron energy loss $`\nu `$ (in the target rest frame) and an invariant momentum transfer $`q^2Q^2`$ between the incoming and outgoing electron such that $`x=Q^2/2m_N\nu `$ is fixed. In terms of Fock states we then view the $`ep`$ scattering as follows: the electron emits a photon ($`|e>|e\gamma >`$) with $`E_\gamma =\nu `$ and $`p_{t\gamma }^2Q^2`$, after the photon splits into a $`q\overline{q}`$ ($`|e\gamma >|eq\overline{q}>`$) and typically travels a distance $`l_c1/m_Nx`$, referred as the coherence length, before interacting in the nucleon. For small $`x`$ (large $`s`$, where $`\sqrt{s}`$ is $`\gamma ^{}p`$ center-of-mass energy), the photon converts to a quark pair at a large distance before it interacts to the target. Consequently, the space-time picture of the DIS in the target rest frame can be viewed as the decay of the virtual photon at high energy (small $`x`$) into a quark-antiquark pair long before the interaction with the target. The $`q\overline{q}`$ pair subsequently interacts with the target. In the small $`x`$ region, where $`x\frac{1}{2mR}`$, the $`q\overline{q}`$ pair crosses the target with fixed transverse distance $`r_t`$ between the quarks. Following Gribov , we may write a double dispersion relation for the forward $`\gamma ^{}p`$ elastic amplitude $`A`$, related to the total cross section by the optical theorem ($`ImA=s\sigma (s,Q^2)`$), and obtain for fixed $`s`$

$`\sigma (s,Q^2)={\displaystyle \underset{q}{}}{\displaystyle \frac{dM^2}{M^2+Q^2}\frac{dM^2}{M^2+Q^2}\rho (s,M^2,M^2)\frac{1}{s}ImA_{q\overline{q}+p}(s,M^2,M^2)},`$ (1)
where $`M`$ and $`M^{}`$ are the invariant masses of the incoming and outgoing $`q\overline{q}`$ pair. If we assume that the forward $`q\overline{q}+p`$ scattering does not change the momentum of the quarks then $`A_{q\overline{q}+p}`$ is proportional to $`\delta (M^2M^2)`$, and (1) becomes

$`\sigma (s,Q^2)={\displaystyle \underset{q}{}}{\displaystyle \frac{dM^2}{(M^2+Q^2)^2}\rho (s,M^2)\sigma _{q\overline{q}+p}(s,M^2)},`$ (2)
where the spectral function $`\rho (s,M^2)`$ is the density of $`q\overline{q}`$ states, which may be expressed in terms of the $`\gamma ^{}q\overline{q}`$ matrix element . Using that $`M^2=(k_t^2+m_q^2)/[z(1z)]`$, where $`k_t`$ and $`z`$ are the transverse and longitudinal momentum components of the quark with mass $`m_q`$, we can express the integral over the mass $`M`$ of the $`q\overline{q}`$ in terms of a two-dimensional integral over $`z`$ and $`k_t`$. Instead of $`k_t`$, it is useful to work with the transverse coordinate $`r_t`$ (impact parameter representation), which is the variable Fourier conjugate to $`k_t`$, resulting

$`F_2(x,Q^2)={\displaystyle \frac{Q^2}{4\pi \alpha _{em}}}\sigma (s,Q^2)={\displaystyle \frac{Q^2}{4\pi \alpha _{em}}}{\displaystyle 𝑑z\frac{d^2r_t}{\pi }|\mathrm{\Psi }(z,r_t)|^2\sigma _{q\overline{q}}(z,r_t)},`$ (3)
where

$`|\mathrm{\Psi }(z,r_t)|^2={\displaystyle \frac{6\alpha _{em}}{(2\pi )^2}}{\displaystyle \underset{i}{\overset{n_f}{}}}e_i^2\{[z^2+(1z)^2]ϵ^2K_1(ϵr_t)^2+m_i^2K_0(ϵr_t)^2\}.`$ (4)
The photon wave function $`\mathrm{\Psi }(z,r_t)`$ is simply the Fourier transform of the matrix element for the transition $`\gamma ^{}q\overline{q}`$. Moreover, $`\alpha _{em}`$ is the electromagnetic coupling constant, $`ϵ^2=z(1z)Q^2+m_i^2`$, $`m_i`$ is the quark mass, $`n_f`$ is the number of active flavors, $`e_f^2`$ is the square of the parton charge (in units of $`e`$), $`K_{0,1}`$ are the modified Bessel functions and $`z`$ is the fraction of the photon’s light-cone momentum carried by one of the quarks of the pair. In the leading log$`(1/x)`$ approximation we can neglect the change of $`z`$ during the interaction and describe the cross section $`\sigma ^{q\overline{q}}(z,4/r_t^2)`$ as a function of the variable $`x`$. Considering only light quarks ($`i=u,d,s`$) $`F_2`$ can be expressed by

$`F_2(x,Q^2)={\displaystyle \frac{1}{4\pi ^3}}{\displaystyle \underset{u,d,s}{}}e_i^2{\displaystyle _{\frac{1}{Q^2}}^{\frac{1}{Q_0^2}}}{\displaystyle \frac{d^2r_t}{\pi r_t^4}}\sigma _{q\overline{q}}(x,r_t).`$ (5)
Using that $`\sigma _{q\overline{q}}=\frac{C_F}{C_A}(3\alpha _s(\frac{4}{r_t^2})/4)\pi ^2r_t^2xG(x,\frac{4}{r_t^2})`$ , where $`xG(x,\frac{4}{r_t^2})`$ is the nucleon gluon distribution, we get

$`F_2(x,Q^2)={\displaystyle \frac{2\alpha _s}{9\pi }}{\displaystyle _{Q_0^2}^{Q^2}}{\displaystyle \frac{dQ^2}{Q^2}}xG(x,Q^2).`$ (6)
The above expression explicit the direct dependence of the $`F_2`$ structure function with the behavior of the gluon distribution, which is related to the QCD dynamics at high energies. In the linear regime, where the high densities can be disregarded, the small $`x`$ behavior of the gluon distribution is given by the solutions of the DGLAP  and/or BFKL  evolution equations. The common feature of these equations is the steep increase of $`xG`$ as $`x`$ decreases. This steep increase cannot persist down to arbitrary low values of $`x`$ since it violates a fundamental principle of quantum theory, i.e. the unitarity. In the context of relativistic quantum field theory of the strong interactions, unitarity implies that the cross section cannot increase with increasing energy $`s`$ above $`log^2s`$: the Froissart’s theorem . In the next section we will consider the modifications in the behavior of the gluon distribution associated to the high density effects. The presence of a high parton density implies a slow growth of $`xG`$ at high energies, and consequently from Eq. (6), of the structure function. However, the Eq. (6) is a leading twist relation which will eventually breaks down when we consider the higher twist terms in the evolution.

We now consider that the relation between the structure function and the gluon distribution is modified by the particular type of higher twist terms associated with the high density corrections (See also  for a similar calculation in the infinite momentum frame). In general, the higher twist contributions should be significant at small $`x`$ and $`Q^2`$, implying large perturbative corrections to the conventional leading twist relations in this kinematical region. We estimated these corrections considering the $`s`$-channel unitarity constraint in the interaction cross section of the quark-antiquark pair with the target . In this case the structure function is given by

$`F_2(x,Q^2)={\displaystyle \frac{1}{2\pi ^3}}{\displaystyle \underset{f=u,d,s}{}}e_f^2{\displaystyle _{\frac{1}{Q^2}}^{\frac{1}{Q_0^2}}}{\displaystyle \frac{d^2r_t}{\pi r_t^4}}{\displaystyle d^2b_t\{1e^{\frac{1}{2}\sigma ^{q\overline{q}}(x,4/r_t^2)S(b_t)}\}}.`$ (7)
A similar expression was used in Ref.  for a phenomenological analysis of the $`ep`$ process, where the geometrical structure of the collision (the $`\stackrel{}{b}_t`$ dependence) was disregarded and $`xgx^\lambda (\lambda >0)`$ was assumed, resulting a very good description of the HERA data. The main point of the above expression is that this resums a large class of higher twist terms, as demonstrated in Ref. .

The use of a Gaussian parameterization for the nucleon profile function $`S(b_t)=\frac{1}{\pi R^2}e^{\frac{b^2}{R^2}}`$, where $`R`$ is the spatial gluon distribution inside the proton , simplifies the calculations. We obtain that the $`F_2`$ structure function can be written as

$`F_2(x,Q^2)={\displaystyle \frac{R^2}{2\pi ^2}}{\displaystyle \underset{u,d,s}{}}ϵ_i^2{\displaystyle _{\frac{1}{Q^2}}^{\frac{1}{Q_0^2}}}{\displaystyle \frac{d^2r_t}{\pi r_t^4}}\{C+ln(\kappa _q(x,r_t^2))+E_1(\kappa _q(x,r_t^2))\},`$ (8)
where $`\kappa _q=4/9\kappa _G=(2\alpha _s/3R^2)\pi r_t^2xG_N(x,\frac{1}{r_t^2})`$. This equation allows to estimate the high density corrections to the structure function in the DLA limit. Expanding the equation (8) for small $`\kappa _q`$, the first term (Born term) will correspond to the usual DGLAP equation in the small $`x`$ region \[Eq. (6)\].

Before we discuss the dynamics for the gluon distribution at high densities, it is important to consider a general property of the structure function (i.e. the cross section): the black disc limit. As shown before, when the coherence length $`l_c`$ considerably exceeds the diameter of a target which is at rest, the virtual photon transforms into hadron components well before the target, which implies that the small $`x`$ physics probes the interaction of various hadron wave packets with a target. The geometrical limit for the cross section of virtual photon scattering off a nucleon target follows from the assumption that a target is black for the dominant hadron components in the wave function of the virtual photon. An identical assumption is often used to deduct the Froissart limit of hadron-hadron interactions at high energies.

Under the assumption that the interaction is black, we can estimate the total cross section (2) using the following approximations: $`\rho (s,M^2)M^2`$ and $`\sigma _{q\overline{q}+p}(s,M^2)=\pi R^2`$. Consequently, we can get $`\sigma (s,M^2)R^2ln[(M_{max}^2+Q^2)/(M_{min}^2+Q^2)]`$. Using that $`M_{min}^24m_{pi}^2`$ whereas $`M_{max}^2s`$, results that

$`F_2(x,Q^2)Q^2R^2ln{\displaystyle \frac{1}{x}},`$ (9)
which is the black disc limit for the $`F_2`$ structure function. This bound, derived from a geometrical analysis, represents the maximum value allowed for the cross sections at high energies (high densities). In the next sections we will analyze the predictions of the high density approaches for the asymptotic behavior of $`F_2`$ and verify if this bound is saturated or a softer behavior is predicted.

## 3 The AGL High Density Approach

About seventeen years ago, Gribov, Levin, and Ryskin (GLR)  performed a detailed study of the high density limit of QCD in the double logarithmic approximation (DLA). They argued that the physical processes of interaction and recombination of partons become important in the parton cascade at a large value of the parton density, and that these high density corrections could be expressed in a new evolution equation \- the GLR equation. Some years ago an eikonal approach to the high density corrections was proposed in the literature by Ayala, Gay Ducati, and Levin (AGL)  which improves the GLR approach. This approach address the high density regime from the region where perturbative QCD is valid by summing corrections to the linear DGLAP evolution equation.

As we are only interested in the solutions of the AGL equation, we refer the original papers  for details in its derivation. Here we only present the AGL equation, which is given by

$`{\displaystyle \frac{^2xG(x,Q^2)}{yϵ}}={\displaystyle \frac{2Q^2R^2}{\pi ^2}}\{C+ln[\kappa _G(x,Q^2)]+E_1[\kappa _G(x,Q^2)]\},`$ (10)
where $`C`$ is the Euler constant, $`E_1`$ is the exponential function and the function $`\kappa _G`$ is defined by

$`\kappa _G(x,Q^2){\displaystyle \frac{\alpha _sN_c\pi }{2Q^2R^2}}xG(x,Q^2),`$ (11)
and represents the probability of gluon-gluon interaction inside the parton cascade. Using the above definition for $`\kappa _G`$ we can rewrite the expression (10) in a more convenient form (for fixed $`\alpha _s`$)

$`{\displaystyle \frac{^2k_G(y,ϵ)}{yϵ}}+{\displaystyle \frac{k_G(y,ϵ)}{y}}`$ $`=`$ $`{\displaystyle \frac{\alpha _sN_c}{\pi }}\{C+ln[\kappa _G(x,Q^2)]+E_1[\kappa _G(x,Q^2)]\}`$ (12)
$``$ $`F(\kappa _G).`$
Analyzing the structure of the Eq. (12) we see that it has a solution which depends only on $`y`$. In  it was shown that this solution is the asymptotic solution of the AGL equation. In this case we have that at large values of densities the asymptotic solution is given by

$`k_G^{asymp}(y)={\displaystyle \frac{\alpha _sN_c}{\pi }}y.`$ (13)
Substituting the definition of $`\kappa _G`$ \[Eq. (11)\] in the above solution, in the asymptotic regime the behavior of the gluon distribution is given by (at fixed $`\alpha _s`$)

$`xG(x,Q^2)={\displaystyle \frac{2N_cQ^2R^2}{3\pi ^2}}ln({\displaystyle \frac{1}{x}}).`$ (14)
Therefore, the gluon distribution does not saturate at small values of $`x`$, but is linearly proportional to $`lns`$ $`(s1/x)`$. However, this behavior is softer than predicted by the DGLAP equation ($`xGexp[\sqrt{ln1/x}]`$) and the BFKL equation ($`xGx^\lambda `$, $`\lambda >0`$). We obtain that the gluon distribution presents a partial saturation in its behavior.

The analyzes above was made at leading order, considering a fixed coupling constant $`\alpha _s`$. One question is the possible modifications in the asymptotic behavior associated to next-to-leading order corrections as, for instance, the running of $`\alpha _s`$. In order to answer this question we need to extend the Refs.  for running coupling constant, as we will do in the following.

Since the QCD coupling constant is given by $`\alpha _s(Q^2)=4\pi /(\beta _0ϵ)`$, where $`\beta _0=112/3n_f`$ ($`n_f`$ is the number of flavors) and $`ϵ=lnQ^2/\mathrm{\Lambda }_{QCD}^2`$, the relation between the gluon distribution and the function $`\kappa _G`$ can be expressed by

$`xG(x,Q^2)={\displaystyle \frac{2Q^2R^2}{N_c\pi }}{\displaystyle \frac{\beta _0}{4\pi }}ϵ\kappa _G(x,Q^2).`$ (15)
Using this result in Eq. (10) we obtain that the AGL equation with running $`\alpha _s`$ is given by

$`{\displaystyle \frac{^2k_G(y,ϵ)}{yϵ}}+\left({\displaystyle \frac{1}{ϵ}}+1\right){\displaystyle \frac{k_G(y,ϵ)}{y}}`$ $`=`$ $`{\displaystyle \frac{N_c\alpha _s(Q^2)}{\pi }}\{C+ln[\kappa _G(x,Q^2)]`$ (16)
$`+`$ $`E_1[\kappa _G(x,Q^2)]\}`$
$``$ $`H(\kappa _G).`$
The AGL equation at fixed $`\alpha _s`$ is a direct consequence of the above equation in the limit of large $`ϵ`$. Moreover, this equation also has a solution which depends only on $`y`$. In this case we have

$`{\displaystyle \frac{k_G^{asymp}(y,ϵ)}{y}}={\displaystyle \frac{ϵ}{1+ϵ}}H(\kappa _G)`$ (17)
with the solution

$`{\displaystyle _{k_G^{asymp}(y=y_0)}^{k_G^{asymp}(y)}}{\displaystyle \frac{d\kappa _G^{}}{H(\kappa _G^{})}}={\displaystyle \frac{ϵ}{1+ϵ}}(yy_0).`$ (18)
Consequently, the asymptotic behavior of the gluon distribution is given by

$`xG(x,Q^2)={\displaystyle \frac{ϵ}{1+ϵ}}{\displaystyle \frac{2N_cQ^2R^2}{3\pi ^2}}ln({\displaystyle \frac{1}{x}}).`$ (19)
As expected at large values of $`ϵ`$ $`(Q^2)`$ this solution reduces to the solution at fixed $`\alpha _s`$. The main difference is the behavior at small values of $`ϵ`$, where the prefactor in (19) is important. However, the partial saturation of the gluon distribution is not modified by the running of $`\alpha _s`$. It agrees with the result that the unitarity corrections are expected to be relevant before the next to leading order corrections .

Before we present our results and conclusions in the next section, let us discuss the generality of the above results for the gluon distribution by the comparison with the predictions of other high density approach. In contrast to the AGL formalism, the approach proposed by McLerran, Venugopalan and collaborators (MV-JKLW approach)  address the high density regime from the nonperturbative region by developing an effective Lagrangian approach to a very dense system. In Ref.  the authors have considered the DLA limit of the functional evolution equation derived from this approach, obtaining an equation similar, but not identical, to the AGL equation (we refer  for details). One of the main results of this work was the solution of the proposed equation for large values of densities (large $`\kappa `$), which is given by

$`xG(x,Q^2)={\displaystyle \frac{N_c(N_c1)\pi }{2}}R^2Q^2ln\left({\displaystyle \frac{1}{x}}\right).`$ (20)
The remarkable feature of this solution is the same $`Q^2`$ and $`x`$ dependence of the AGL equation in the asymptotic regime. The difference in the prefactors is a function of the distinct normalizations and approximations used in the two approaches.

## 4 The Asymptotic Behavior of $`F_2`$

The main conclusion of our analysis in the previous section was the universality of the $`x`$ and $`Q^2`$ dependence of the gluon distribution (the partial saturation) in the asymptotic regime of the high density QCD. In this section we will consider the consequences of the partial saturation in the behavior of the structure function, which is experimentally measured.

Using the solution of the AGL equation in the asymptotic regime \[Eq. (14)\] as input in the Eq. (6) we get

$`F_2(x,Q^2){\displaystyle \frac{\alpha _s}{\pi ^3}}R^2Q^2ln\left({\displaystyle \frac{1}{x}}\right).`$ (21)
We see that the partial saturation of the gluon distribution implies that the $`F_2`$ structure function does not saturate at small values of $`x`$, but is linearly proportional to $`lns`$. Basically, we verify that the AGL approach, as well as the MV-JKLW approach maximizes the prediction for the $`F_2`$ structure function in the asymptotic regime of very high densities, resulting in the black disc prediction, when the leading twist relation \[Eq. (6)\] is used.

We now consider that the relation between the structure function and the gluon distribution is given by the Eq. (8). In the asymptotic regime (large $`\kappa _q`$) we obtain

$`F_2(x,Q^2){\displaystyle \frac{R^2}{2\pi ^2}}{\displaystyle \underset{u,d,s}{}}ϵ_i^2{\displaystyle _{\frac{1}{Q^2}}^{\frac{1}{Q_0^2}}}{\displaystyle \frac{d^2r_t}{\pi r_t^4}}ln(\kappa _q(x,r_t^2)).`$ (22)
Using the asymptotic solution of the AGL equation we can determine $`\kappa _q=4/9\kappa _G`$ at large values of densities, and so

$`F_2(x,Q^2){\displaystyle \frac{R^2}{2\pi ^2}}{\displaystyle \underset{u,d,s}{}}ϵ_i^2{\displaystyle _{\frac{1}{Q^2}}^{\frac{1}{Q_0^2}}}{\displaystyle \frac{d^2r_t}{\pi r_t^4}}ln\left[{\displaystyle \frac{4\alpha _s}{3}}ln({\displaystyle \frac{1}{x}})\right].`$ (23)
Therefore, considering the contribution of the higher twist terms in the relation between the structure function and the gluon distribution, we predict the following asymptotic behavior for the structure function

$`F_2(x,Q^2){\displaystyle \frac{R^2Q^2}{3\pi ^2}}ln\left[{\displaystyle \frac{4\alpha _s}{3}}ln({\displaystyle \frac{1}{x}})\right].`$ (24)
We see that the inclusion of the higher twist terms implies a softer dependence of $`F_2`$ with the energy than obtained using the leading twist relation. In contrast with the leading twist prediction, in this case the black disc limit is not maximized. However, in both cases the structure function does not violate the Froissart boundary in the asymptotic regime of high density QCD. The demonstration of the behavior (24) using the other approaches for hdQCD is an important open question. In particular, the analyzes of the predictions for $`F_2`$ obtained in Ref.  in the infinite momentum frame for the very high density regime would be a important check of our results. Our study demonstrates the importance of the correct connection (leading or higher twist relation) between the gluon distribution and the observables when we are analyzing the interface between perturbative and nonperturbative QCD.

From the results for the structure function in the asymptotic regime we can see that this regime is characterized by the identity

$`{\displaystyle \frac{dF_2(x,Q^2)}{dlnQ^2}}=F_2(x,Q^2),`$ (25)
which is an important signature of the asymptotic regime of high density QCD. This regime should be reached for the case of an interaction with nuclei at smaller parton densities than in a nucleon, since $`\kappa _A=A^{\frac{1}{3}}\times \kappa _N`$, where $`\kappa _N`$ is given by the expression (11).

## 5 Conclusions

In this paper we have analyzed the description of the $`F_2`$ structure function in the regime of very high densities, assuming first that the leading twist relation between $`F_2`$ and $`xG`$ is valid and second that this relation is modified by the higher twist terms associated to the high density corrections. In the first case we have obtained that the corresponding $`F_2`$ structure function is linearly proportional to $`Q^2R^2ln\mathrm{\hspace{0.17em}1}/x`$, which agrees with the results predicted in the black disc limit. This result is not high density approach dependent. In the second case a softer behavior is obtained. In both cases, the $`F_2`$ structure function unitarizes and, as expected by construction, the black disc boundary is not violated in the asymptotic regime of high density QCD. We verify that when the parton density is such that the proton becomes black and the interaction probability is unity, the structure function becomes proportional to $`Q^2`$ and a soft energy dependence is predicted. As a by product, we have derived the AGL equation for running $`\alpha _s`$ and obtained a signature for the asymptotic regime of QCD.

Our main conclusion is that we should be careful to relate the gluon distribution with the observables (for instance, $`F_2`$, $`F_2^c`$, heavy quark production, …) in the high density limit, since besides the behavior of the gluon distribution, also the relation between these quantities must be modified from the usual leading twist expressions.

## Acknowledgments

VPBG thanks FAPERGS and CNPq for support. This work was partially financed by CNPq and by Programa de Apoio a Núcleos de Excelência (PRONEX), BRAZIL.


text:
# Relativistic blastwaves and synchrotron emission

## 1 Introduction

Relativistic shock fronts arise in astrophysics when a relativistic flow propagates into ambient material such as in the jets of active galactic nuclei (AGNs) and fireball models of gamma ray bursts (GRBs). These shock fronts are also sites of energetic particle acceleration which in turn leads to non-thermal emission processes such as synchrotron radiation and inverse Compton scattering. Particle acceleration at relativistic shock fronts can occur via the first order Fermi mechanism. In this picture, energetic particles of speed $`v`$ are scattered back and forward across the shock front by magnetic turbulence. Each crossing involves the particle meeting the flow head-on leading to an increase of energy. At nonrelativistic shocks, $`v>>U_{1,2}`$ where $`U_1`$ and $`U_2`$ are the upstream and downstream flow speeds in the shock rest frame respectively. With diffusive particle scattering the particle distribution is almost isotropic both upstream and downstream. This results in a power law index for the particle spectrum which only depends on the shock compression ratio, $`r=U_1/U_2`$. In relativistic flows the situation is more complicated. When $`v`$ is comparable to $`U_{1,2}`$ the timescale for isotropisation of the energetic particle population by the magnetic turbulence (which determines the degree of anisotropy of the particle distribution upstream and downstream) becomes important. Moreover even if the particle distribution were isotropic in one fluid frame it will be anisotropic when viewed in the shock or other fluid frame. As a result, unlike the case for nonrelativistic flows, the particle spectrum depends not only on the shock strength but also on the details of how the particles are scattered . This creates difficulties when studying shock acceleration in flows such as spherical, relativistic blastwaves , the energetic particle distribution is the solution of a kinetic equation which describes scattering in the local fluid rest frame.

However, there has been recent progress with ultra-relativistic shocks, $`\mathrm{\Gamma }>10`$, for which it has been shown that the index of the power-law distribution of energetic particles, which undergo isotropic scattering, is $`q=2.23`$ where $`N(E)E^q`$ is the differential number of particles with energy $`E`$ . This result has potentially important implications for fireball models of GRBs (Cavallo & Rees 1978; Goodman 1986; Paczyński 1986) in which a large amount of energy is suddenly released in the form of an optically thick electron-positron plasma. The presence of baryonic matter  results in the conversion of the initial energy into a relativistically moving bulk flow. The forward and reverse shocks will accelerate particles by the above mechanism in the presence of magnetic turbulence. Similarly, internal shocks can occur because of variability in the central engine of a GRB and these are also sites of relativistic shock acceleration. In this paper we present a model for the synchrotron emission of shock accelerated particles downstream of the forward and reverse shocks in a spherically symmetric relativistic blastwave. In contrast to the strict fireball model our initial conditions are for a gas at rest with a large amount of internal thermal energy confined to a small volume. As in the fireball case the gas will accelerate to the point where all of the initial energy is converted into bulk kinetic energy if the external medium is sufficiently tenuous. Ultimately the flow will become self-similar when the amount of swept up mass becomes significant and the flow will decelerate thereafter.

Our principle aim in this paper is to study the influence of an evolving, inhomogeneous flow on the energetic particles produced by the first order Fermi mechanism at relativistic shocks. This is of central importance in the theory of GRBs and many papers have recently approached related problems on a number of different levels. The adiabatic evolution of a GRB interacting with an external medium, and the associated emission, has been considered  where an analytic model is introduced for the Lorentz factor of the shock front and relativistic remnant. A similar hydrodynamical model has also been employed in Moderski, Sikora & Bulik  for the deceleration of the blast wave but for beamed ejecta. Such a geometry has also been discussed in a recent paper  where a full hydrodynamical simulation is carried out. In this paper we solve the full hydrodynamical equations for a spherically symmetric system which starts from rest but begins to expand relativistically, and ultimately becomes self-similar, as a result of the enormous amount of internal energy in the central engine. Moreover, particles are accelerated at all shock waves which arise in our simulation, which in this case are the forward and reverse shocks, with a spectrum determined by acceleration theory. This model will be generalised to jet type geometries in future work.

In order to achieve this we have developed a hydrodynamical code, described in Sect. 2 and Appendix A, which integrates the relativistic Euler equations reliably for Lorentz factors of up to several hundred. This code is applicable to a gas with a general equation of state in one dimension and incorporates adaptive hierarchical mesh refinement. In Sect. A.1 we present tests of our hydrodynamical code. Our model for the injection of accelerated particles into the flow, and their subsequent cooling, is presented in Sect. 3. The initial conditions are presented in Sect. 4, while in Sect. 5 we give a discussion of the hydrodynamics of blastwaves with non-zero initial radius. The results are presented and discussed in Sect. 6.

## 2 Hydrodynamics

We wish to simulate a spherically symmetric explosion which gives rise to relativistic velocities. We use the relativistic Euler system of equations in order to calculate the evolution of this system. The code used is a second-order Godunov-type scheme (e.g. van Leer 1977) which uses a linear Riemann solver for shocks, and a non-linear one for strong rarefactions. Appendix A gives the details of the code, along with some test calculations. In this section we briefly describe the equations used.

The conservation equations for inviscid relativistic hydrodynamics in spherical symmetry are

$`{\displaystyle \frac{}{t}}\left(\mathrm{\Gamma }\rho \right)+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left(r^2\mathrm{\Gamma }\rho \beta \right)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{}{t}}\left(w\mathrm{\Gamma }^2\beta \right)+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left[r^2\left(w\mathrm{\Gamma }^2\beta ^2+p\right)\right]`$ $`=`$ $`{\displaystyle \frac{2p}{r}}`$ (2)
$`{\displaystyle \frac{}{t}}\left(w\mathrm{\Gamma }^2p\right)+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left(r^2w\mathrm{\Gamma }^2\beta \right)`$ $`=`$ $`0`$ (3)
where $`\mathrm{\Gamma }`$ is the fluid Lorentz factor, $`\rho `$ is the proper density, $`\beta `$ is the velocity in $`c=1`$ units, $`w`$ is the enthalpy and $`p`$ is the proper pressure. Time, $`t`$, and distance, $`r`$, refer to the coordinates measured in the observer’s frame. We can relate the enthalpy, density and pressure by

$$w=\rho +\frac{\gamma _{}}{\gamma _{}1}p$$
(4)
If the gas is composed of $`N`$ species which each have mass $`m_i`$ and number density $`n_i`$ then, for a Synge gas ,

$$\frac{\gamma _{}}{\gamma _{}1}=\left[\underset{i=1}{\overset{N}{}}m_in_i\left(\frac{K_3(\xi _i)}{K_2(\xi _i)}1\right)\right]\left[\underset{i=1}{\overset{N}{}}\frac{m_in_i}{\xi _i}\right]^1$$
(5)
where $`\xi _i=\frac{m_i}{kT}`$ and $`K_\alpha `$ are the modified Bessel functions. In this work we assume two species - electrons and protons - which are in thermal equilibrium. The variable $`\gamma _{}`$ defined above is different to the ratio of specific heats, $`\gamma `$ which is used in the calculation of, for example, sound-speeds. With this relation we should be able to solve the evolution equations 1 to 3. However, it should be noted that in relativistic hydrodynamics there are no explicit relations giving the primitive variables $`\rho `$, $`\beta `$ and $`p`$ in terms of the conserved variables $`\mathrm{\Gamma }\rho `$, $`w\mathrm{\Gamma }^2\beta `$ and $`w\mathrm{\Gamma }^2p`$. Therefore an iterative algorithm must be developed to do this. The precise nature of this algorithm can affect how high a Lorentz factor can be simulated by the code. We have found that the code used here will simulate flows up to Lorentz factors of at least several hundred reliably.

## 3 Particle acceleration and synchrotron emission

### 3.1 The particle distribution

Relativistic shocks, in the presence of magnetic fluctuations which enable multiple shock crossings, can accelerate particles by the diffusive shock mechanism . As described in the introduction it has been shown that at ultra-relativistic shocks this produces a power law distribution of particles with index $`2.23`$, i.e. $`N(E)E^{2.23}`$, so that the energy density does not diverge at high particle energies. The acceleration timescale in relativistic flows will scale roughly with the particle mean free path in the turbulent magnetic field divided by the shock speed; although there is no published derivation of the precise acceleration timescale. In this paper we are interested in what happens to the power law population produced by highly relativistic shocks over the much longer timescale which is the size of the blastwave divided by $`c`$. Therefore, over hydrodynamical timescales, we can treat acceleration as an impulsive injection of energetic particles, with this universal power law up to arbitrarily high energies, into each fluid element which passes through a relativistic shock. We will address the role of an intrinsic cut-off in the shock accelerated population, either as a result of a finite acceleration timescale or loss processes, in a future paper. The particles subsequently suffer synchrotron and adiabatic losses from which we can calculate the emission. We need to introduce three free parameters which, when taken together with the results of our hydrodynamical simulations, will give a complete picture of the blastwave’s evolution and the associated emission of energetic particles. These parameters are

* the ratio, $`ϵ_b`$, between the magnetic field energy density and the thermal energy density,
* the fraction, $`ϵ_e`$, of the downstream thermal energy density which is converted into high energy electrons as a fluid cell is shocked and,
* the energy, $`E_{\mathrm{min}}`$, of the lowest energy particle which is produced at the shock. The differential energy spectrum is a power law in energy of index $`2.23`$ from $`E_{\mathrm{min}}`$.

Once an electron is injected at the shock it will be scattered in the local fluid frame and suffer both synchrotron and adiabatic losses. If the length scale over which the electrons are isotropised is much shorter than any other length scale of interest then the relativistic electrons will respond adiabatically to the expansion or contraction of the flow. Adiabatic losses, or indeed gains, are then described by the fact that $`p/\rho ^{1/3}`$ is constant where $`p`$ and $`\rho `$ are the particle momentum and fluid density in the local fluid frame. With $`E=pc`$ for ultra-relativistic particles the combined synchrotron and adiabatic losses are described in the comoving frame by

$`\dot{E}=\alpha B^2E^2+{\displaystyle \frac{1}{3}}{\displaystyle \frac{\dot{\rho }}{\rho }}E`$ (6)
where $`\alpha `$ is a constant and $`B`$ the magnetic field strength in the local fluid frame. Consider now a fluid element which is shocked at time $`\widehat{t}_0`$ when a power law distribution of energetic particles is injected and where $`\widehat{t}`$ is time measured in the comoving frame. This population then evolves according to the equation

$`{\displaystyle \frac{N}{\widehat{t}}}+{\displaystyle \frac{}{E}}\left(\dot{E}N\right)=Q(E,\widehat{t})`$ (7)
where $`N(E,\widehat{t})`$ is the differential number of particles of energy $`E`$ at time $`\widehat{t}`$. The losses, $`\dot{E}`$, are given by equation (6) while the injection term is $`Q(E,\widehat{t})=Q_0\delta (\widehat{t}\widehat{t}_0)H(EE_{\mathrm{min}})E^p`$ ($`H`$ being the Heviside function) which describes injection of a power law spectrum, starting at time $`\widehat{t}_0`$ with a minimum particle energy $`E_{\mathrm{min}}`$. We can solve equation (7) by finding the characteristic curves along which $`N(E,\widehat{t})dE=N(E_0,\widehat{t}_0)dE_0`$ where a particle with energy $`E_0`$ at $`\widehat{t}_0`$ cools to an energy $`E`$ at time $`\widehat{t}`$. From equation (6) we have

$`{\displaystyle \frac{d}{d\widehat{t}}}\left(E^1\rho ^{\frac{1}{3}}\right)=\alpha \rho ^{\frac{1}{3}}B^2`$ (8)
which can be solved to show that along a characteristic curve

$`E_0={\displaystyle \frac{\rho _0^{\frac{1}{3}}E}{\rho ^{\frac{1}{3}}\alpha E_{\widehat{t}_0}^{\widehat{t}}\rho ^{\frac{1}{3}}B^2𝑑\widehat{t}}}`$ (9)
$`N`$ is conserved so that

$`{\displaystyle \frac{dE_0}{dE}}={\displaystyle \frac{\rho _0^{\frac{1}{3}}\rho ^{\frac{1}{3}}}{\left(\rho ^{\frac{1}{3}}\alpha E_{\widehat{t}_0}^{\widehat{t}}\rho ^{\frac{1}{3}}B^2𝑑\widehat{t}\right)^2}}`$ (10)
and the solution to (7) becomes

$`N(E,\widehat{t})=Q_0\rho _0^{\frac{1p}{3}}\rho ^{\frac{1}{3}}E^p\left(\rho ^{\frac{1}{3}}\alpha E{\displaystyle _{\widehat{t}_0}^{\widehat{t}}}\rho ^{\frac{1}{3}}B^2𝑑\widehat{t}\right)^{p2}`$ (11)
In order to calculate the particle spectrum we need to evaluate the integral in equation (11) for each fluid element which has passed through a shock at some point. No fluid element passes through two or more shocks during these simulations; our method would have to be adapted in such a case (e.g. for overtaking internal shocks). We introduce a continuity equation for the integral quantity in equation (11). Define $`I_{t_0}^{t_0+t}\frac{\rho ^{\frac{1}{3}}B^2}{\mathrm{\Gamma }}𝑑t`$ where the factor of $`\mathrm{\Gamma }^1`$ boosts observer time to the comoving fluid time. Then the equation for the distribution of $`I`$ is simply

$$\frac{}{t}\left(\mathrm{\Gamma }\rho I\right)+\frac{1}{r^2}\frac{}{r}\left(r^2\mathrm{\Gamma }\rho I\beta \right)=B^2\rho ^{\frac{4}{3}}$$
(12)
Thus the quantity $`I`$ increases by $`\rho ^{1/3}B^2\mathrm{\Delta }t/\mathrm{\Gamma }`$ with each time-step $`\mathrm{\Delta }t`$, following the fluid element, as required. If the gradient of $`\mathrm{\Gamma }\beta `$ is ever less than a certain negative value, chosen to ensure that the fluid element has passed through a shock, then $`I`$ is set to 0 at that point. The above conservation equation can easily be solved in conjunction with the hydrodynamical equations using the same numerical methods. Note that $`I`$ is not kept at zero in unshocked fluid, but this does not affect the synchrotron emission since, in such fluid, there are no energetic particles.

### 3.2 Synchrotron emission

We approximate the emission of a single particle to be a delta function in frequency with a single particle emissivity in the fluid frame given by $`j_{\widehat{\nu }}(\gamma )=a_0\gamma ^2B^2\delta (\widehat{\nu }a_1\gamma ^2B)`$ where $`\gamma `$ is the Lorentz factor of the electron in the fluid frame, $`a_0=\sigma _Tc/6\pi `$ and $`a_1=e/2\pi m_ec`$ with $`\sigma _T`$ the Thomson cross section. With the local electron spectrum given by 11, the emissivity in the local fluid frame is then

$$\widehat{}_{\widehat{\nu }}=j_{\widehat{\nu }}(\gamma )N(\gamma ,\widehat{t})𝑑\gamma $$
(13)
The emissivity in the observer’s frame is related to that in the fluid frame by $`_\nu =D^2\widehat{}_{\widehat{\nu }}`$. The Doppler factor, $`D`$, for the fluid element with a three velocity $`\beta `$ making an angle $`\theta `$ to the line of sight is $`D=[\mathrm{\Gamma }(1\mu \beta )]^1`$ where $`\mu =\mathrm{cos}\theta `$. The observed frequency is related to the photon frequency in the fluid frame by $`\nu =D\widehat{\nu }`$. The observed intensity of radiation at an inclination $`\varphi `$ to the line between the observer and the centre of the explosion is

$`I_\nu (\varphi ,t_0)={\displaystyle _0^{\mathrm{}}}_\nu (s,\varphi ,t)𝑑s`$ (14)
where $`t_0`$ is the time of observation and $`t=t_0s/c`$ is the time of emission in the observer’s frame a distance $`s`$ away from the observer (figure 1). The observed flux density is obtained by integrating the intensity over the entire source

$`F_\nu (t_0)={\displaystyle _\mathrm{\Omega }}I_\nu (\varphi ,t_0)𝑑\mathrm{\Omega }`$ (15)
In cylindrical coordinates $`(z,r)`$ centred on the explosion, $`sLz`$ since $`\varphi 1`$ so that

$`F_\nu (t_0)={\displaystyle \frac{2\pi }{L^2}}{\displaystyle _{r_c}^{r_c}}𝑑z{\displaystyle _0^{r_c}}_\nu (z,r,t)r𝑑r`$ (16)
with $`L`$ the distance from the observer to the centre of the explosion and $`r_c`$ the size of the computational domain. If $`r_0`$ is the initial radius of the explosion then the interval of $`t_0`$ since the arrival of the first photon is $`\overline{t_0}=t_0(Lr_0)/c`$ so that the flux can be written as $`F_\nu (\overline{t_o})`$. By the change of variable $`z=r_0+c(t\overline{t_0})`$ the integration over $`z`$ can be converted to one over $`t`$, the time of emission in the observer’s frame,

$`F_\nu (\overline{t_0})={\displaystyle \frac{2\pi c}{L^2}}{\displaystyle _0^{t_{\mathrm{end}}}}𝑑t{\displaystyle _0^{r_c}}_\nu (r_0+c(t\overline{t_0}),r,t)r𝑑r`$ (17)
In the simple case where the fluid is assumed to move directly towards the observer with a constant velocity we have successfully compared our method with the exact results. The above analysis, which is essentially the same as in Komissarov & Falle , also corresponds for the latter simple case to the treatment used in Moderski et al. . In particular from those particles which have not radiated away a significant fraction of their energy we get the uncooled spectrum of $`F_\nu \nu ^{(p1)/2}`$ while the higher energy particles have a steeper, cooled spectrum of $`F_\nu \nu ^{p/2}`$. Below the minimum observed frequency the flux scales as $`F_\nu \nu ^{1/3}`$. At even lower frequencies the effect of synchrotron self-absorption will become important, although we have not explicitly calculated this, leading to a $`\nu ^2`$ part of the spectrum.

## 4 Initial conditions

In studies of the fireball model of gamma-ray bursts, the parameters used are the total mass, $`M`$, and the initial radius, $`R_0`$. In general, $`E`$, the total energy, is thought to lie somewhere in the region $`10^{51}`$-$`10^{54}`$ ergs. The ratio between the rest-mass energy in the blast (i.e. the mass of the baryonic component), and the total energy determines two critical radii. These are the radius at which the baryons have been accelerated up to their maximum velocity, $`R_c`$, and the radius at which the shell of ejected baryons have swept up their own mass in interstellar material, $`R_d`$. This latter radius is the radius at which the baryonic shell begins to decelerate. The former also depends on $`R_0`$, the initial radius of the blast. This initial radius is thought to be quite small, with $`R_010^{12}`$ cm.

In order to get all the initial energy in the blast converted into kinetic energy of the baryons, it is also required that the fireball be optically thick to pair creation until the radius of the blast exceeds $`R_c`$. If this were not the case, then the photons would escape into space before accelerating the baryons to their maximum velocity. In our case we put the initial energy of the blast into thermal pressure. This means that it is already in the kinetic energy of the baryons. Hence the latter consideration need not be taken into account in our initial conditions.

The results of two simulations are presented in this paper. Each of the simulations used the following properties:

* $`E=10^{51}`$ ergs
* Ratio of energy to mass: $`\eta =\frac{E}{Mc^2}=580`$. This value being chosen so that the Newtonian value of $`R_d`$ is twice $`R_c`$. In reality, $`R_d`$ is slightly less than $`R_c`$ due to relativistic effects.
* $`R_0=1.2\times 10^{14}`$ cm. This is much larger than the postulated physical value. However, reducing $`R_0`$ to $`10^{12}`$ cm would impose extremely severe computational costs on the simulation. This increase in the value of $`R_0`$ is not expected to change the properties of the synchrotron emission calculated as we are interested in the afterglow in this work, and not the emission from the acceleration phase.
* $`ϵ_e=0.01`$ \- this value is consistent with observed values found by, e.g. Wijers & Galama
* $`ϵ_b=0.01`$ in one simulation, and $`ϵ_b=0.1`$ in the other. This is the only way in which the two simulations differ.

We must also choose the injection energy $`E_{\mathrm{min}}`$ for the shock acceleration process. Physically this will be set by nonlinear plasma processes. The highest energy thermal particles downstream of a shock can leak into the upstream forming a beam of particles which can excite plasma turbulence. This turbulence then scatters particles back into the downstream fluid and multiple shock crossings become possible. The problem of determining the details of this process is an unsolved one in the theory of particle acceleration at relativistic shocks. However, it is clear that the value $`E_{\mathrm{min}}`$ must correspond to the high energy tail of the downstream thermal population. This is hardly surprising since the thickness of a shock wave will be set by the gyroradius of a thermal particle in a collisionless plasma and first order Fermi acceleration will only be applicable to particles with gyroradii much larger than this. Consequently we have fixed the value of $`E_{\mathrm{min}}=10m_ec^2`$ which is chosen so that it exceeds the thermal energy for all shocked fluid elements. Setting $`E_{\mathrm{min}}`$ as a factor of a few times the thermal energy at the instant a piece of fluid is shocked does not change our results appreciably.

Initially, then, there is a sphere of radius $`R_0`$ which is at rest at the origin. This sphere contains all the initial energy and mass of the blast. Outside this sphere is a gas with a density of 1 cm<sup>-3</sup>.

## 5 Hydrodynamic evolution of a spherical blastwave

In this section we discuss the qualitative evolution of a spherical blastwave with finite initial radius. This discussion applies to both relativistic and non-relativistic cases. The purpose is to elucidate, in a qualitative fashion, how the forward and the reverse shocks are formed, and how they behave, from a hydrodynamic perspective. Figure 2 shows the early-time evolution of the system, computed using the code presented here. The initial conditions used were as follows:

$`\rho `$ $`=`$ $`\{\begin{array}{cc}3\hfill & \text{on }r1\hfill \\ 1\hfill & \text{on }r>1\hfill \end{array}`$
$`u`$ $`=`$ $`0`$
$`p`$ $`=`$ $`\{\begin{array}{cc}10^3\hfill & \text{on }r1\hfill \\ 3\times 10^3\hfill & \text{on }r>1\hfill \end{array}`$
These initial conditions are similar to those in Sect. A.1.3, and are chosen so that the relevant features are clearly visible in the results.

Initially a shock is driven into the ambient medium (the “forward shock” in figure 2), while a rarefaction propagates from $`r=R_0`$ back to $`r=0`$ (the “ingoing rarefaction” in figure 2). Note that, at $`t=1.5`$, the ingoing rarefaction has still not reached the origin. The forward shock gets progressively stronger as blast material gets accelerated down the pressure gradient of the rarefaction. There is also an entropy wave (the “contact discontinuity”) which moves out from $`R_0`$ with a speed slightly lower than that of the blastwave.

In addition, a reverse shock is created almost immediately between the entropy wave and the rarefaction (see figure 2). If the system were in planar symmetry then one would not expect this reverse shock since there are only three characteristics in hydrodynamics (see, e.g. Landau & Lifschitz 1966). We get this fourth wave purely as a result of the non-planar nature of the system. It forms because, as the forward shock moves outwards it sweeps up an increasing amount of interstellar matter per unit distance. This means that the forward shock does not move outwards as fast as it would in the planar case, and the reverse shock is formed to decelerate material behind the forward shock to the appropriate speed.

When the rarefaction reaches $`r=0`$ it is reflected and enhanced. This “reflected rarefaction” can be seen clearly at $`t=3`$ in figure 2. This rarefaction will eventually catch up with the reverse shock/blastwave system (see figure 4 and Sect. 6.1). When the rarefaction catches up with the reverse shock, the ram-pressure of the material entering this shock is drastically reduced, because it has been rarefied and decelerated. Hence it can no longer support the thermal pressure between the forward and reverse shocks. As a result, the material between these two shocks expands. Since the ram-pressure of the material entering the forward shock has not changed (assuming constant external density), this necessitates the reverse shock slowing down with respect to the rest-frame of the initial blast. Indeed, it will actually begin propagating towards the origin. The speed at which it does this depends on the initial value of $`\eta `$, since this determines the strength of the initial (and so the reflected) rarefaction. For example, if $`\eta `$ is very small, then the force due to the pressure gradient will produce a relatively small increase in velocity of the material in the initial sphere. Hence the rarefaction will be weak. If, on the other hand, $`\eta `$ is large then the force due to the pressure gradient will accelerate the blast material to high velocities very quickly, producing a very strong rarefaction. For the conditions in Sect. 4 it becomes sufficiently strong to cause ejected material to flow back towards the origin also.

When the reverse shock reaches $`r=0`$ it rebounds and, as it propagates outwards again, it weakens. Once this shock is sufficiently weak to be ignored, we have reached the self-similar stage of evolution. Prior to this point the finite radius of the initial blast plays a part, and hence there is a significant length-scale in the system. If the blast was initially of zero radius then the whole process described above would happen infinitely quickly (or, equivalently, would not happen at all).

We can relate the above “hydrodynamical” picture with the one used by, for example, Rees & Mészáros , as follows. The coasting radius, $`R_c`$, is reached when all material has been accelerated down the pressure gradient of the ingoing rarefaction. This, of course, never occurs, but one can define $`R_c`$ as being the radius when “most” of the material has been accelerated.

The deceleration radius, $`R_d`$, is then the radius at which the reflected rarefaction catches up with the reverse shock. This is when the material between the forward and reverse shock expands back towards the origin, necessarily decelerating significantly as it does so.

## 6 Results

We discuss the results of the simulations in two sections. The first deals with the detailed hydrodynamic evolution of the blastwave. The second concerns the observed spectra and light-curves.

### 6.1 Hydrodynamic results

Here we show the hydrodynamic aspect of the results of our simulations. Figure 3 shows the distribution of proper density, velocity and pressure after $`1\times 10^5`$ seconds. We can see the two rarefactions mentioned in Sect. 5. The head of the reflected rarefaction lies at the location of the peak velocity, with its tail at $`r=0`$. The original rarefaction lies between the peak of velocity and the rise in the density and pressure plots at $`1\times 10^5`$ light-seconds.

The forward shock can be identified as the right-most rise in pressure and density, while the reverse shock is just to the left of this.

Figure 4 shows the same plots as figure 3, but for a time of $`7\times 10^6`$ seconds. Here we can see the reverse shock has begun to separate significantly from the forward shock. The reflected rarefaction is now the only rarefaction present in the system. The reverse shock is expanding back into this rarefaction due to the insufficient ram-pressure of the material in the reflected rarefaction to support the pressure of the material between the forward and reverse shocks.

It is possible to see oscillations in the density plot just behind the reverse shock. These oscillations arise out of numerical errors during the time that the reflected rarefaction comes close to, but is not in contact with, the reverse shock. The rarefaction is extremely strong here, and it is not surprising that numerical errors become significant. However, the influence of these errors on the overall solution, or on the calculated synchrotron emission, is negligible due to the low density of the region containing the errors.

Figure 5 shows the latter stages of the non-self-similar evolution of the system. At this stage the reverse shock has propagated all the way to $`r=0`$ and has rebounded. We can see that it has weakened greatly, and will continue to do so. Once this shock is sufficiently weak to be ignored the system will evolve in a self-similar way. We can see the entropy errors generated by the very strong reflected rarefaction remaining just ahead of the reverse shock now. These errors, although apparently significant here, do not affect the synchrotron results due to the low densities in this region, and the low relativistic boosting suffered by emission from this material. Both these effects make the influence of emission from this material negligible.

### 6.2 Spectra and light curves

In this section we present the results from the calculations of the synchrotron emission.

#### 6.2.1 Spectra

Figure 6 shows a plot of the spectrum observed 24 hours after the initial blast could have been observed for the case where $`ϵ_b=0.1`$ (hereafter referred to as the high $`ϵ_b`$ case) and for $`ϵ_b=0.01`$ (hereafter referred to as the low $`ϵ_b`$ case). It is clear that in both cases the spectrum is a broken power-law.

In the first section of the spectra, the flux goes as $`\nu ^{1/3}`$, as expected for synchrotron radiation below the peak frequency of the lowest energy electron. At lower frequencies, which are not plotted here, synchrotron self-absorption would give a flux of $`\nu ^{1/3}`$. The second part of the spectrum is, in both cases, the power law $`\nu ^{0.615}`$ indicating that this part of the spectrum is dominated by emission from electrons, with a distribution of $`E^{2.23}`$, which have not suffered significant adiabatic or synchrotron cooling.

The break from this part of the spectrum to the final, steeper, part occurs in different places in the spectra plotted. For the high $`ϵ_b`$ case, the break occurs at a lower frequency than for the low $`ϵ_b`$ case as would be expected since, in higher magnetic fields, the losses suffered by the energetic population are correspondingly higher. However, while in the final part of the spectrum the exponent predicted from simple theory would be $`p/2=1.115`$, we find the spectrum to be slightly harder. In the high $`ϵ_b`$ case, we have $`F_\nu \nu ^{1.019}`$, and in the low $`ϵ_b`$ case, $`F_\nu \nu ^{0.981}`$. The reason for the slightly harder spectra at high frequencies is the non-uniform velocity distribution in the “shell” of ejected material (see, e.g., figure 4). Material moving at high velocity towards the observer will be more heavily weighted in the spectra than material moving at lower velocity, due to relativistic beaming. Since such high velocity material occurs immediately behind the forward shock, this material will contribute more to the spectrum than material further back from the shock. The fluid just behind the shock contains an electron population which has only recently been accelerated and so emission from here will be dominated by uncooled electrons to a very high frequency. Material from further behind the shock will have emission dominated by cooled electrons down to a lower frequency. Hence, above a certain cut-off, we expect to get a mixture between emission from cooled and un-cooled electrons, with un-cooled electrons being preferentially weighted. This leads to an exponent for the spectrum lying between the uncooled value ($`\frac{1p}{2}=0.615`$) and the cooled value ($`\frac{p}{2}=1.115`$). This conclusion is further supported by the fact that the high $`ϵ_b`$ case gives a slightly steeper spectrum at these frequencies than the low $`ϵ_b`$ case. For high magnetic fields the contribution from cooled electrons will be stronger closer to the shock than in the low magnetic field case. Therefore, while we still expect the uncooled electrons to be preferentially boosted in the spectrum, the effect of the cooled electrons will be stronger for higher magnetic fields.

It is worth discussing the behaviour of the ‘critical’ frequencies in the spectra with time - i.e. where the breaks occur. The lower break, which arises due to the presence of electrons injected with the minimum energy $`E_{\mathrm{min}}`$ behaves roughly as $`t^{0.73}`$. This is the same for both the low and high $`ϵ_b`$ cases. This is a much less dramatic decrease with time than that predicted by Sari & Piran  and is due to the fact that, here, we define the minimum energy of injection a priori. This energy is chosen so that the gyroradius of an electron at this energy would be accelerated by the Fermi mechanism. The latter authors choose to define the number of electrons accelerated, and this, in conjunction with the energy in energetic particles, defines $`E_{\mathrm{min}}`$. This produces the different behaviours of the lower break in the spectra with time.

We find that the frequency of the second (upper) break goes like $`t^{1.2}`$ for the low $`ϵ_b`$ case and $`t^{1.75}`$ for the high $`ϵ_b`$ case. This is a much faster decrease with time than that predicted in Sari & Piran . This difference cannot be due to the different choice of definition of $`E_{\mathrm{min}}`$, but it could be hypothesised that it is due to the effects of adiabatic cooling. This cooling will cause the second critical frequency to decrease faster where the flow is divergent. However, simulations of this system without adiabatic cooling show that, while adiabatic cooling does have this effect on the observed spectrum, it is too small to explain the discrepancy. It is not clear what specifically causes this difference between our results and those of Sari & Piran , but it must be remembered that here we perform a much fuller treatment of the dynamics and so it is not surprising that we get some significant differences with semi-analytic approaches.

#### 6.2.2 Light curves

In this section we discuss the light curves resulting from the simulations described above. Figure 7 shows plots of the light curves at low, medium and high frequencies. These frequencies are chosen so that the points are always in the regime where $`F_\nu \nu ^{1/3}`$ (low frequencies), $`F_\nu \nu ^{\frac{1p}{2}}`$ (medium frequencies) and, for high frequencies, the frequency is always in the third section of the spectrum.

It should be noted in the following discussion that the range of times covered by these light curves is 1 to 24 hrs after the first signal of the blast reaches the observer. Initially we see an increase in all the light curves shown, with the exception of the high frequency, high $`ϵ_b`$ case, where the emission is roughly constant over the first couple of hours. This increase is much faster than that predicted by previous analytic work. For example, for the low $`ϵ_b`$ case at low frequencies, the flux goes as roughly $`t^{1.4}`$ at early times (cf Sari & Piran, 1997, where the maximum increase is $`t^{1/2}`$). This effect is not due to any problems with numerical resolution since the emission in the 1st hour is dominated by material between the forward and reverse shocks approximately $`6\times 10^6`$ seconds after the initial blast, when the two shocks are well-resolved and separated by the code.

The behaviour of the low frequency light curves is qualitatively different to the other two. It keeps increasing for the duration of the simulation. This is unsurprising as it will only begin to decrease when the critical frequency $`\nu _{\mathrm{min}}`$ corresponding to $`\mathrm{\Gamma }_{\mathrm{min}}`$ passes through the frequency of the light curve, and we have chosen this frequency so that this does not happen. It is very difficult to say that the light curve behaves anything like a power-law with breaks over this time-scale based on our results.

The medium frequency light curves initially increase, as already stated, and then fall off dramatically. Again, it is difficult to see any point at which the behaviour of the light curve is a true power-law. The high frequency light curves also begin to fall off very steeply, and again, it is difficult to see any sign of a power-law behaviour in the curves.

The lack of a clear power-law behaviour in the light curves may well be due to the restricted time-scale over which the curves are calculated (from 1 hour to 1 day). However, it is clear that, if there is a broken power-law behaviour then the breaks are smeared out, and, in addition, the final fall-off of the light-curve would seem to be much faster than previously predicted.

The unexpected behaviour of the low and medium frequency light curves is more likely to do with the different treatment of the minimum energy in the accelerated electron population. In our case, since we fix $`E_{\mathrm{min}}`$ a priori, we allow the available energy to define the number of energetic electrons. Previously it has been common to define the number of energetic electrons and hence the available energy defines $`E_{\mathrm{min}}`$. This will lead to quite different behaviours of the lower break frequency with time.

## 7 Conclusions

We have presented a model for the synchrotron emission of energetic particles downstream of relativistic, spherical shock waves. The hydrodynamical part of the problem has been solved numerically with the simpler simulations agreeing with results published elsewhere. On the other hand the particle acceleration aspect has been treated as an injection process with relativistic shocks leaving a population of energetic particles immediately downstream. The spectral index is known from semi-analytic work so that we need only specify the fraction of the downstream thermal energy which is converted into energetic particles and a lower cut-off energy. The particles subsequently lose energy by synchrotron cooling and adiabatic losses. The hydrodynamical results have captured the evolution of both the forward and reverse shocks, which are of principal interest for particle acceleration and radiative emission, as well as the rarefaction waves. The spectra emitted from our system largely agree with the simple predictions for the low energy, uncooled part of the spectrum. However, at the higher frequency, cooled, part of the electron population the relativistic boosting of the material coming straight towards the observer hardens the spectrum from the pure cooled value which comes from material further downstream and from material which is not moving directly towards the observer. Further, non-trivial behaviour has been found for the temporal variation of both the break frequencies and the light curves. Adiabatic losses coupled with integrating the spectrum over a spherical system tend to smear out the break frequencies which are predicted from simple scaling arguments.

It is clear that we need to include several other effects before bringing our results into contact with observations. A generalisation to a 2-D hydrodynamical code has already begun. The role of internal shocks will also be included along with their effect on acceleration. The energetic particles obey a transport equation, the solution of which determines the spectral index and we are working towards including such a solution in our model. More complicated radiative processes can also be included such as the inverse-Compton and synchrotron self-Compton processes. However, what we have done in this paper is to combine a realistic hydrodynamical model with results from particle acceleration theory in a relativistic flow. The computed hydrodynamical evolution, spectra and light curves will allow us to compare this model with observations in future work.

## Acknowledgments

We would like to thank Luke Drury for useful discussions on the incorporation of adiabatic losses into the electron distribution. We are grateful to an anonymous referee for comments which led to improvements in the paper. This work was supported by the TMR programme of the European Union under contract FMRX-CT98-0168.

## Appendix A Numerical Details

In this appendix we discuss the details of the code used for the calculations presented in this paper. We also present the results of three test calculations.

We employ a second order finite volume Godunov-type scheme to solve the equations 1 to 3 in spherical symmetry. Assuming that the grid cells are defined so that cell $`i`$ occupies the space $`[r_{i1/2},r_{i+1/2}]`$ then the scheme can be written as

$`\stackrel{}{U}_i^{n+1}=\stackrel{}{U}_i^n{\displaystyle \frac{3\mathrm{\Delta }t}{r_{i+1/2}^3r_{i1/2}^3}}[r_{i+1/2}^2\stackrel{}{F}_{i+1/2}^{n+1/2}`$
$`r_{i1/2}^2\stackrel{}{F}_{i1/2}^{n+1/2}]+\stackrel{}{S}_i^{n+1/2}\mathrm{\Delta }t`$ (20)
where superscripts refer to the time index and subscripts refer to the spatial index. Also,

$`\stackrel{}{U}_i^n`$ $`=`$ $`\left[(\mathrm{\Gamma }\rho ,\mathrm{\Gamma }\rho \beta ,w\mathrm{\Gamma }^2p)^T\right]_i^n`$ (21)
$`\stackrel{}{F}_{i\pm 1/2}^n`$ $`=`$ $`\left[(\mathrm{\Gamma }\rho \beta ,w\mathrm{\Gamma }^2\beta ^2+p,w\mathrm{\Gamma }^2\beta )^T\right]_{i\pm 1/2}^N`$ (22)
$`\stackrel{}{S}_i^n`$ $`=`$ $`\left[(0,{\displaystyle \frac{2p}{r}},0)^T\right]_i^n`$ (23)
It should be noted that the source term $`\stackrel{}{S}`$ should be volume-averaged over $`[r_{i1/2},r_{i+1/2}]`$ . So, for spherical symmetry, and using a first order scheme we get that the source term for the radial momentum equation is

$$S_i^n=3p_i^n\frac{r_{i+1/2}+r_{i1/2}}{r_{i+1/2}^2+r_{i+1/2}r_{i1/2}+r_{i1/2}^2}$$
(24)
For a second order scheme, which allows $`p`$ to vary linearly in a given cell the correct form for $`S_i^n`$ is

$`S_i^n`$ $`=`$ $`3p_i^n{\displaystyle \frac{r_{i+1/2}+r_{i1/2}}{r_{i+1/2}^2+r_{i+1/2}r_{i1/2}+r_{i1/2}^2}}`$ (25)
$`+g_{p}^{}{}_{i}{}^{n}\left[{\displaystyle \frac{9(r_{i+1/2}+r_{i1/2})^2(r_{i+1/2}^2+r_{i1/2}^2)}{4(r_{i+1/2}^2+r_{i+1/2}r_{i1/2}+r_{i1/2}^2)^2}}\right]`$
where $`g_{p}^{}{}_{i}{}^{n}`$ is the gradient of the pressure in cell $`i`$ at time-step $`n`$.

The fluxes $`\stackrel{}{F}`$ are calculated from the primitive variables extrapolated to the cell edges using non-linear averaging of the gradients (e.g. van Leer 1977; Downes & Ray 1998), making the code second order in space. A non-linear Riemann solver is used in the presence of strong rarefactions, and otherwise a linear Riemann solver is used in order to save computational time. The flux terms $`\stackrel{}{F}_{i\pm 1/2}^{n+1/2}`$ are calculated using a first order Godunov-scheme from the conditions at time $`n`$. This makes the code second order in time.

The details of the linear and nonlinear Riemann solvers follow closely the discussions of Falle  and Falle & Komissarov .

The code also contains an option for hierarchical grid refinement. The algorithm used for this is loosely based on that of Khokhlov . However, this option is not used in the work presented here, so we will not discuss this further.

### A.1 Tests of the hydrodynamical code

The code used in this work has been rigorously tested against, in particular, results presented in Falle & Komissarov (1996). The code reproduces their published results. Below we show three examples of these tests for completeness.

#### A.1.1 Test 1: Shock-tube problem with $`\gamma =\frac{4}{3}`$

In this test we use the following initial conditions:

$`\rho `$ $`=`$ $`1`$
$`u`$ $`=`$ $`0`$
$`p`$ $`=`$ $`\{\begin{array}{cc}10^3\hfill & \text{on }x0.5\hfill \\ 10^2\hfill & \text{on }x>0.5\hfill \end{array}`$
with gradient zero boundary conditions on $`x=0`$ and $`x=1`$. This calculation was done for a uniform grid with 400 cells. The results are shown in figure 8, along with the exact solution. The plots are of the proper density, and are at different times. It can be seen that, while the contact discontinuity becomes quite smeared (as expected, since it is a linear wave) the shock and rarefaction are treated relatively well.

In particular, the position of both the shock and the rarefaction are very well calculated, indicating that both the linear and non-linear Riemann solvers are working well.

#### A.1.2 Test 2: Shock-tube problem for a Synge gas

In this test we use the following initial conditions:

$`\rho `$ $`=`$ $`1`$
$`u`$ $`=`$ $`0`$
$`p`$ $`=`$ $`\{\begin{array}{cc}10^3\hfill & \text{on }x0.5\hfill \\ 10^2\hfill & \text{on }x>0.5\hfill \end{array}`$
with gradient zero boundary conditions on $`x=0`$ and $`x=1`$. This calculation was done for a uniform grid with 400 cells. The results are shown in figure 9, along with the exact solution. Again, it can be seen that the numerical results match the analytic results reasonably well, although numerical dissipation does affect the contact discontinuity quite significantly.

#### A.1.3 Test 3: Spherically symmetric expanding wave

This test is very similar to the simulations used in this work. Therefore it is an interesting test to run. The conditions are as follows:

$`\rho `$ $`=`$ $`10`$
$`u`$ $`=`$ $`0`$
$`p`$ $`=`$ $`\{\begin{array}{cc}10^3\hfill & \text{on }r0.5\hfill \\ 3\times 10^3\hfill & \text{on }x>0.5\hfill \end{array}`$
The grid spans $`r=0`$ to $`r=4`$ with 100 uniformly spaced cells. The results obtained are very similar to those given in Falle & Komissarov  and are shown in figure 10.

Overall, then, the code reproduces published results very well. One additional test was performed. This was a test of the code for the non-relativistic Sedov blastwave. The code successfully reproduced the behaviour of the radius of the blastwave (with $`rt^{\frac{2}{5}}`$).

It is noted in Sect. 6.1 that the simulations performed here give entropy errors in some parts of the fluid. These errors arise when the reflected rarefaction (see Sect. 5) has almost caught up with the reverse shock. In this case there is an extremely strong rarefaction, with differences of the Lorentz factor of the fluid of 30 – 40 between adjacent grid cells, just behind the reverse shock. These very extreme conditions result in negative pressures being produced (in the grid cell, not in the Riemann solver), subsequently resulting in entropy errors being apparent in the solution.

As stated in Sect. 6.1, however, these errors occur in a region of the flow where the density of emitting particles is relatively low, and the errors do not affect the simulated emission appreciably.