A&A 452, 743-749 (2006)
DOI: 10.1051/0004-6361:20035604

Channelled relativistic outflows in active galactic nuclei: analytic solutions for the evolution of particle spectra

C. Schuster ¹ - I. Lerche^1,2 - R. Schlickeiser¹ - M. Pohl³

1 - Institut für Theoretische Physik, Lehrstuhl IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum, Germany
2 - Institut für Geophysik und Geologie, Universität Leipzig, Germany
3 - Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011

Received 31 October 2003 / Accepted 4 November 2005

Abstract
Context. Calculations of highly energetic neutrino and TeV $\gamma$-ray emission from relativistic jets of Active Galactic Nuclei (AGN) are often based on a model that involves a collimated relativistic blast wave, in which the spectral evolution of energetic particles is determined by the interplay between the particle injection by sweep up of the interstellar medium, energy losses by radiation and diffusive escape.
Aims. To date such models only have been solved numerically because of the highly non-linear nature of the time-dependent equations describing particle spectra and bulk energy loss. However it is difficult to see how parameters and groupings of parameters influence the resulting numerical solutions, except through intensive investigations of many numerical simulations. Therefore analytic solutions are particularly helpful.
Methods. We provide exact mathematical solutions to a very broad class of AGN type models. By selecting different functional behaviors of parameter values, we cover a large compendium of possible situations. The comparison of the exact solutions with observational information can help to improve our understanding of the evolution of individual AGN. Exact solutions can also be used to provide controls on the appropropriateness and accuracy of numerical programs used to solve the equations.
Results. We provide an analytical description of the evolution of proton spectra according to the pick-up model. We analyze the behavior of the particle spectra in the plasma frame. The solutions are determined by two competing processes: the deceleration of the jet plasmoid and particle cooling via radiation.

Key words: methods: analytical - radiation mechanisms: non-thermal - galaxies: jets - galaxies: active - plasmas - gamma rays: theory

1 Introduction

Active galactic nuclei (AGN) have been established as powerful emitters of high-energy gamma-ray photons with energies from $\rm {MeV}$ to $\rm {TeV}$ $\gamma$-rays (for a recent review see Maraschi & Tavecchio 2005). They are also prime candidates for sources of high-energy neutrino radiation (Schuster et al. 2002; Waxman 2005) although no detections have been reported yet (Ackermann et al. 2004).

The reported $\gamma$-ray emissions are highly variable on all time scales to the observational limits of days at $\rm {GeV}$ energies and hours at $\rm {TeV}$ energies (von Montigny et al. 1995; Mukherjee et al. 1997; Catanese & Weekes 1999; Mattox et al. 1997; Gaidos et al. 1996). Combined with the huge luminosities these observations imply a relativistic Doppler amplification of the high-energy photon radiation. In many cases the relativistic bulk motion is directly observed as apparent super-luminal motion of individual emission regions in the so-called jets of VLBI observations of their radio emission (Pohl et al. 1995; Barthel et al. 1995; Piner & Kingham 1997a,b). It is therefore generally accepted that the high-energy photon (and neutrino) emission of AGNs originates in the relativistic jets of these objects. The Lorentz and Doppler factors derived from VLBI campaigns are of the order ten for general samples of AGNs (Vermeulen & Cohen 1994), but may be higher for AGNs showing prominent $\gamma$-ray emission (Homan et al. 2002, 2003).

The range of Lorentz factors present at the time of $\gamma$-ray emission, which presumably occurs before the emission region becomes visible at radio frequencies, is not known, but is likely considerably higher than ten, especially as the bulk kinetic energy of the jets is the most probable energy reservoir for particle acceleration. To accomodate all such possibilities, this paper allows the initial bulk flow Lorentz factor $\Gamma$ to be up to several hundred, although the analysis is valid for all $\Gamma^2 \gg 1$. In particular, when graphical results are presented, high initial bulk flow Lorentz factors are used to emphasize the relativistic aspects. Lower initial Lorentz factors can also be used graphically, of course, because the analytic solutions are valid as long as $\Gamma^2 \gg 1.$

How the kinetic energies of such relativistic outflows is transferred to energetic particles, and subsequently converted into high-energy photon (and neutrino) radiation, is one of the most important issues in current AGN research. Many existing radiation models of AGNs (see e.g. Böttcher & Dermer 1998; Dermer & Schlickeiser 2002; Ghisellini 2002; Sikora 1997) are unspecific on the microphysical details of this crucial point. Without detailed discussion these models often assume that a significant fraction of the energy in nonthermal baryons of the relativistic outflow is injected into nonthermal ultrarelativistic pairs and/or hadrons with power-law energy distribution functions in the comoving jet frame. This significant conversion is attributed to the scenario that the outflowing relativistic jet plasma has produced a relativistic shock wave with fully developed hydromagnetic turbulence in order to allow for efficient nonthermal diffusive particle acceleration at the collisionless shock fronts.

In a series of papers (Pohl & Schlickeiser 2000; Pohl et al. 2002; Schlickeiser et al. 2002; Schlickeiser 2003; Vainio 2004; Schlickeiser et al. 2003), the microphysical details of this energy conversion process have been investigated both for leptonic and hadronic relativistic outflows. As described in Pohl & Schlickeiser (2000) the emission region in the AGN jet is assumed to be a cloud of dense plasma (hereafter called plasmoid), which moves relativistically through the interstellar medium of the AGN host galaxy. The swept-up ambient matter is quickly isotropized in the jet frame by a relativistic two-stream instability, which provides relativistic particles in the jet without invoking any further acceleration process. The typical primary particles would be protons with an initial Lorentz factor of the order of one hundred, distributed isotropically in a jet also moving with an initial bulk Lorentz factor of the order of one hundred. In the case of hadronic outflows, inelastic proton-proton collisions with the outflow plasma would then lead to the production of TeV $\gamma$-rays by $\pi^0$-decay and high energy neutrinos via decay of charged pions ($\pi^\pm$), as well as synchrotron, inverse Compton and bremsstrahlung emission produced by secondary electrons at lower energies. Due to this efficient pick-up of interstellar particles, the outflowing cloud is decelerated and mass-loaded, so that the Lorentz factor drops quickly to the order of 10, as observed.

The time-dependent numerical calculations of Pohl & Schlickeiser (2000) indicate that the observed variability behavior and $\gamma$-ray spectra of blazars may be accounted for especially if the outflowing cloud propagates through a structured non-uniform ambient interstellar or intergalactic medium. The observed variable light curves then reflect the spatial inhomogeneities of the traversed medium because the varying ambient gas density regulates the injection rate of primary protons, electrons and hydrogen neutrals into the outflowing cloud.

There are two aspects to the physics of radiation modeling. There are available detailed scenarios for radiation production from two flow leptonic models (Henri & Pelletier 1991), or from hadronic models (Mannheim et al. 1991; Rachen 2000), but these models do not include detailed descriptions of the primary particle (proton) spectra evolution in the evolving jets. Instead such particle spectra are assumed as input and then the radiation production is calculated. What is required is the second aspect: a detailed, time-dependent, self-consistent model of primary particle evolution that can then be coupled to radiation production. Precisely that problem is the focus of this paper.

The evolution of relativistic electrons and their associated radiation has already been discussed by Pohl & Schlickeiser (2000). Here we provide analytic results of the relevant equations rather than the purely numerical results reported earlier.

One can distinguish four main parts of the relativistic pick-up process:

(1)

after rapid isotropisation in the outflow cloud due to the short radiative interaction time scales, the temporal evolution of energetic particle spectra in the comoving jet frame has to be calculated;

(2)

from the calculated particle spectra the yield of generated photons and neutrinos has to be determined in the comoving frame;

(3)

the influence of Lorentz transformations on the produced photon and neutrino radiation to an observer's rest frame (the so called Lorentz boost) has to be determined using Lorentz invariants (see the discussion in e.g. Dermer & Schlickeiser 2002);

(4)

to transform to the observer's frame one has to calculate the time behaviour of the bulk Lorentz factor from a hydrodynamical model with self-consistent mass loading induced deceleration of the cloud.

In the original model of Pohl & Schlickeiser (2000) these four parts were handled numerically. We provide predominantly analytical calculations of these four steps. These analytical calculations are of considerable value as they indicate more transparently the parameters and conditions that control the calculated photon and neutrino spectra for specific AGNs.

Not only are TeV-emitting AGN highly variable but there is a correlation between the different emission bands in these AGN spectra (especially the two bumps shape of the blazars) with the correlation between the different bands beeing important in constraining the physics of the emission mechanism. A flow through non-uniform media may account for such conditions. However this paper is concerned dominantly with exploring and obtaining the self-consistent flow and associated primary particle spectra equations. A detailed comparison with observations of the resulting model behaviors is reserved for a future paper. Once one has a self-consistent structure, it is a matter of detailed determination of all parameters in the model based on a comparison with radiation from each jet to elucidate the "best fit'' to such data. This has already been done for a few jets (Schuster 2005) using the self-consistent model developed here, but considerably more work must be undertaken for more objects.

The organization of the paper is as follows: after reviewing the basic physical features of the model for channelled relativistic outflows in Sect. 2, we show the analytical solution for the mass of the jet plasmoid in Sect. 3. This solution comes in useful for the analytical description of the energy distribution of protons in an AGN jet, derived in Sect. 4. In Sect. 5 we deal with an example of secondary particle production the resulting $\gamma$-rays. In Sect. 6 we summarize and discuss the results.

As already stated, VLBI results indicate Lorentz factors of at least order ten. The initial bulk flow models developed here allow for, but do not require, initial outflow Lorentz factors of several hundred. However a collimated flow (also called jet) can be produced by the accretion disk surrounding the central black hole (the Blandford-Payne scenario) or the central black hole itself (the Blandford-Znajek scenario). Neither of these mechanisms are able to produce such relativistic outflow (the extreme upper limit appears to be a Lorentz factor of the order of a few units. However a Lorentz factor of at least ten from radiation observations (and surely higher for $\gamma$-ray emission) already stretches this to their limits, if not beyond. Here we deal with what happens to a relativistically moving jet and its associated hadron (proton) particle spectra. The high initial bulk Lorentz factor (of order of a few hundred) used in the examples sharply emphasizes the consequences.

   
2 Basic dynamics of the model

We briefly review the assumptions of the model for a collimated hadronic outflow as given in Pohl & Schlickeiser (2000). An emission region in the jet outflow is assumed to be a plasma cloud consisting of electrons and protons with a cylindrical shape of thickness d, being small compared to the radius R. The cloud is moving with a bulk Lorentz factor $\Gamma$through the ambient medium. Viewed in the rest frame of the jet plasma cloud, the interstellar medium forms an electron-proton beam of number density $n_{\rm i} =\Gamma~ n_{\rm i}^*$, which is low compared with the plasma density $n_{\rm b}= 10^8~ n_{\rm b,8}~\rm {cm^{-3}}$ of the jet cloud.

Pohl & Schlickeiser (2000) showed that the incoming particles of the interstellar medium quickly become isotropized by self-excited $\rm Alfv\acute{e}nic$ turbulence, but retain their relative velocities with respect to the jet plasma. The jet plasma thus becomes enriched with relativistic particles, of which mainly the protons are of interest as secondary generators (e.g. $\gamma$-rays).

So the relativistic protons are injected at a rate

$$% \tilde N(\gamma,t)= \pi R^2 n_{\rm i}^{\ast} c \sqrt{\Gamma(t)^2-1}~ \delta(\gamma-\Gamma(t))$$ (1)

in the case of particle isotropization by pure electromagnetic turbulence. Note that $\gamma$ describes the Lorentz factor of the individual protons, whereas $\Gamma$ represents the bulk speed of the jet.

As a main consequence of the sweeping up of interstellar matter we deal with a system that is not stationary. Momentum conservation requires that the bulk Lorentz factor $\Gamma$ of the plasmoid decreases at a rate:

$$% {\partial \Gamma \over \partial t} = - {m_{\rm p}~ \pi~ R^{2}~ n_{\rm i}^{\ast}~ c~ (\Gamma(t) ^{2} - 1)^{3/2} / M(t)}.$$ (2)

Note the key fact that the relativistic mass

$$% M(t) = m_{\rm p} ~ V ~\left [ n_{\rm b}~ (1+ \epsilon)\; +... ...1^\infty \!\! \gamma \; n(\gamma,t) ~{\rm d} \gamma~ \right ]$$ (3)

is obtained from the particle spectrum $n(\gamma,t)$ in the jet plasma. Here we abbreviate the volume of the cylindrical disc with $V=\pi~ R^{2}~d.$ We also introduce the thermal energy $\epsilon= k_{\rm B}~T /(m_{\rm p}~c^2)\simeq 9.185$ $\times$ $10^{-14}[{\rm K}^{-1}]~ T[{\rm K}] \ll 1$, stating that the thermal particles in the jet outflow are nearly at rest.

Before the detailed discussion we comment on the physical parameters relevant for the emission regions of AGN jets:

1) from variability time scales $\delta t$ corrected for beaming ${\rm d}~\Gamma=c~ \delta t$ we estimate the size d to range from  $10^{12}~\rm {cm}$ (for sources variable on the scale of hours) to d=2.6 $\times$ $10^{13}~\rm {cm}$ (for sources variable on the scale of days) with $\Gamma=10^2 \Gamma_2$.

2) from Doppler-boosted $\gamma$-ray luminosity we estimate the density $n_{\rm p}$ of the non-thermal particles:

$$% L\simeq D^3 \langle E \rangle /\tau_{\pi}$$ (4)

with the Doppler factor $D=[\Gamma~ (1-\beta\cos\theta)]^{-1},$ which is approximately $\Gamma$ for small angles. We use the total average energy content of non-thermal particles:

$$% \langle E \rangle = n_{\rm p} \Gamma m_{\rm p} c^2\pi R^2d$$ (5)

and take the energy loss rate from Mannheim & Schlickeiser (1994)

$$% {\rm d}E/{\rm d}t= - 0.65~c~ n_{\rm b}~ \sigma_{\rm pp} (E-m_{\rm p}c^2)~ \theta [E-E_{\rm th}] \;~ {\rm erg/s},$$ (6)

where $\theta$ denotes the Heaviside function, $\sigma_{\rm pp}= 3$ $\times$ $10^{-26}~{\rm cm^2}$ the total pion cross-section and $E_{\rm th}= 1.9$ $\times$ $10^{-3}~{\rm erg}$ the threshold energy for pion production. The resulting loss time is

$$% \tau_{\pi} T/\vert{\rm d}T/{\rm d}t\vert= (E-m_{\rm p}c^2)/\vert{\rm d}E/{\rm d}t\vert$$ =

$$\simeq 1.7 \times 10^{7}~{\rm s}.$$ (7)

Inserting values for the parameters leads to:

$$% \langle E \rangle = 4.7 \times 10^{42} ~ n_{\rm p}~ \Gamma_{2}~ R_{15}^2 ~ d_{13} \;\;\; {\rm erg}.$$ (8)

Typical luminosities for blazars range from L = 2.7 $\times$ $10^{43}~{\rm erg/s}$ for Mkn 421 (Punch 1992) (with an average flux of 1.5 $\times$ $10^{-11}~{\rm [ph/cm^2/s]}$) up to L= 1.3 $\times$ $10^{44}~{\rm erg/s}$ for the blazar H 1426+428, which is located at z=0.129 (here we use the average photon flux above one TeV of 1.7 $\times$ $10^{-12}~{\rm [ph/cm^2/s]}$) (Aharonian et al. 2002) and 3C 273 with an average 7-week luminosity in the EGRET band (100 MeV-10 GeV) of L=1.7 $\times$ $10^{46}~{\rm erg/s}$ (Collmar et al. 2000). Therefore the non-thermal particles would have densities from $n_{\rm p} \simeq 100~{\rm cm^{-3}}$ for Mkn 421 and $n_{\rm p} \simeq 500~{\rm cm^{-3}}$ up to

$$% n_{\rm p}= 5.6 \times 10^4~{\rm cm^{-3}}$$ (9)

for a source like 3C 273. These values are uncorrected for absorption from background.

Our detailed modelling is performed for these values of parameters and thus relevant for AGN jets.

Assuming the background magnetic field to be uniform and directed along the x-axis and neglecting any expansion of the beam, the temporal evolution of proton spectra $N(\gamma,t)=V~n(\gamma,t)$ obeys an enhanced continuity equation

$$% {\partial N(\gamma,t) \over \partial t} + {\partial [\dot \... ...{\rm E}} + {N(\gamma,t) \over T_{\rm N}} = \tilde N(\gamma,t)$$ (10)

for continuous energy losses from elastic and inelastic scattering at a rate of

$$% -\dot \gamma = \tau_{\rm el}^{-1}~ {\gamma\over \sqrt{\gamma^2-1}} +\tau_{\rm inel}^{-1}~ {(\gamma-1)^2 \over \gamma+1}$$ (11)

with time scales $\tau_{\rm el} = 3.33$ $\times$ $10^{15}~n_{\rm b}^{-1}~\rm {s}$ and $\tau_{\rm inel}= 1.43$ $\times$ $10^{15}~n_{\rm b}^{-1}~\rm {s}$. Particle losses arising from diffusive escape are taken into account by the time scale  $T_{\rm E},$ as are losses via neutron production in $p \rightarrow n$ reactions by considering the time scale $T_{\rm N}$. For more details of the underlying model see Pohl & Schlickeiser (2000) and references therein.

The escape and diffusion of relativistic particles (protons, electrons and perhaps neutral hydrogen through charge exchange processes) is influenced by both the strength and orientation of any prevailing magnetic field in relation to the flow direction, as is known from many other astrophysical plasma kinetic situations. In particular, transverse and parallel diffusion are known to be different (Schlickeiser 2003). The influence of such a magnetic field orientation, particularly if non-uniform, on the self-consistent calculation reported here would be interesting to undertake, as would the influence of an energy-dependent diffusion tensor on the analytic solution. If, as suggested by Bogovalov & Tsinganos (1999), a non-relativistic MHD-jet is required for collimation, the assumption of a weak beam necessary for isotropisation would have to be abandoned. The latter assumption, however, appears to be corroborated by observation.

   
3 The time-dependent mass of the jet plasma

To achieve an analytical description of the particle spectrum  $n(\gamma,t)$ consider first the behavior of the time dependent mass M of the jet plasma. Transform Eq. (10) to:

$$% {\partial\Gamma \over \partial t} {\partial n \over \partia... ...rm N}} = { \lambda~} \sqrt{\Gamma^2-1}\; \delta(\gamma-\Gamma)$$ (12)

by using the fact that the bulk Lorentz factor $\Gamma$ is a monotonic and continuously decreasing function of increasing time t (cf. Eq. (2)). Note also that we deal here with the normalized density n=N/V and abbreviate the combination of physical parameters $c~ n_{\rm i}^\ast/ d$ by $\lambda$. The main astrophysical interest in the relativistic beam equations pertains to relativistic particles because these are the only particles that can produce the MeV to TeV photons. Thus, the main concern is for particles for which one can replace Eq. (11) by:

$$% - {}\dot \gamma \simeq {\gamma \over \tau} \qquad\mbox{for}\qquad\gamma \gg 1.$$ (13)

In the context of the underlying model we deal with values of the initial bulk Lorentz factor $\Gamma$ of up to several hundred. One can describe the time dependent mass  $M(\Gamma(t))$ as:

 

$$% M(\Gamma) = { n_{\rm b}(1+\epsilon) m_{\rm p} V\over \sqrt{\Gam... ...\Gamma)] \sqrt{2~Z}} -{K \over \exp(Z~[ x(\Gamma)-x(\Gamma_0)])} \Biggm ]^{-1}.$$ (14)

This result is derived in the appendix, where also the constant $K=K(\Gamma_0),$ depending only on the initial value for the bulk Lorentz factor $\Gamma_0$, is given. Here we use the abbreviation $x(\Gamma)= [{\Gamma/ \sqrt{\Gamma^2-1}}-1]$ and Z is given by a dimensionless combination of physical parameters:

$$% Z = {n_{\rm b}(1+\epsilon)~ d~ \over \tau_{\rm eff}~ c n_{\rm i}^{\ast}} \simeq {N_{\rm b}\over N_{\rm i,eff}}$$ (15)

where we recognize

$$% 1/\tau_{\rm eff}= (1/\tau + 1/T_{\rm E} +1/T_{\rm N})$$ (16)

as the effective time scale and define $N_{\rm b}= n_{\rm b} d$, and $N_{\rm i,eff}=n_{\rm i}^{\ast}~\tau_{\rm eff}c=n_{\rm i}^{\ast}~ l_{\rm eff}$ as the column depths of jet plasma and interstellar medium, respectively. With typical values for physical parameters of $n_{\rm b}= 10^8~ n_{\rm b,8}~\rm {cm^{-3}}$ for the inner density of the jet plasma, $n_{\rm i^\ast} = 10^{-1}~n_{\rm i^\ast,{-1}}~\rm {cm^{-3}}$, for the density of the surrounding medium and values for the thickness of the cylindrical plasma disk of $d = 10^{13}~d_{13}~\rm {cm},$ a wide range of values of Z is obtained.

Note especially that with the dependence of the parameters for the time scale of diffusive escape

$$% T_{\rm E} \sim d^2~ n_{\rm i}^{\ast}/ \sqrt{n_{\rm b}}$$ (17)

small values of Z are only obtainable for a thickness of about $d = 10^{13}~d_{13}~\rm {cm},$ or with large values of outer, but low inner densities.

As we will see in the next section, the temporal behavior of the mass function determines the particle distribution.

   
4 The particle distribution

The main technical problem is to solve the partial differential Eq. (10) analytically for the particle spectrum  $n(\gamma,t)$. After multiplication of Eq. (12) by $\gamma$, integration and use of the standard method of characteristics we deduce

$$% \gamma(\Gamma) = \gamma_0~ \exp\left[ {(m_{\rm p} V \lambda... ... (\tilde\Gamma^{2} - 1)^{3/2}}\; {\rm d}\tilde\Gamma\right],$$ (18)

the equation of characteristics, where we introduce $\gamma _0$ as the constant of integration. The constant is determined subject to the initial condition $\gamma(\Gamma_0) =\gamma(\Gamma_0,\gamma_0) =\gamma_0,$ so $\gamma _0$ is the Lorentz factor labeling a specific characteristic curve.

The characteristic curves can be referred to as the paths of the particles in a $\gamma$-$\Gamma$ phase space because Eq. (18) completely determines the cooling of the protons. For example a proton that starts at t=0 with an initial bulk Lorentz factor $\Gamma_0=300$ can be found after some time t with a reduced Lorentz factor $\gamma$ that coincides with the value $\gamma$ of the characteristic curve labelled $\gamma_0 = 300$ at the corresponding Lorentz factor $\Gamma$. Note also that, subject to the initial condition $n(\gamma,t_0)=0,$ no particles populate particle paths at t = 0. The characteristics may be therefore referred to as "virtual paths'' until they are enriched with particles subject to the sweep-up condition, which is determined by Eq. (1).

Figure 1: Example characteristic curves. Some particle paths intersect the line of sweep-up (dashed line of $\Gamma =\gamma$) and some (lower $\gamma _0$) never touch the sweep-up line. In between there has to be one curve that just touches the line of sweep-up. The value of Z determines for which curve and at which $\Gamma$ this osculating curve occurs. Here we use $\tau = 2$ $\times$ $10^{15}/n_{\rm b}~[1/{\rm s}]$ for the time scale of radiation losses. For a thickness, d = 3 $\times$ $10^{13}~\rm {cm}$ the value of Z= 4.7 $\times$ 105 can be obtained for example by $n_{\rm b} = 2$ $\times$ $10^8~\rm{cm^{-3}}$ and $n_{\rm i}^\ast = 0.4~\rm{cm^{-3}}$ for the number densities inside and outside the jet, respectively.

Open with DEXTER

Figure 1 gives typical characteristic curves. Additionally, the line of sweep-up of particles, $\gamma =\Gamma$, a diagonal line, is displayed. With the appropriate approximations for small and large parameters Z, the function  $\gamma(\Gamma)$, given by Eq. (18) can be integrated in closed form. For very large values of the parameter Z the function (18) reduces to:

$$% \gamma =\gamma_0~ \exp\left[ {n_{\rm b}(1+\epsilon)\over \l... ...amma^2_0-1}}-{\Gamma\over \sqrt{\Gamma^2-1}} \right) \right].$$ (19)

The remaining task, the derivation of the particle spectra of protons, $n(\gamma,t)$ can be achieved by solving Eq. (12) using the characteristics, as shown in the appendix. The result is:

 

$$% n(\gamma,t)= \sum_i {\left[\Theta(\Gamma_0-\nu_i)-\Theta(\Gamma... ...u_i)~ /[m_{\rm p} V \lambda\tau ~(\nu_i^{2} \!-\! 1)^{3/2}]}\!-\!1\right\vert},$$ (20)

where $\Theta$ is the Heaviside step-function and $\nu_i=\nu_i(\gamma_0)$ denotes the different zeros of ${\gamma(\Gamma,\gamma_0)-\Gamma,}$ that occur where the characteristic curves intersect the line of sweep up and with $\gamma(\Gamma,\gamma_0)$ given by Eq. (18) (and approximately by Eq. (19) for large Z). Equation (2) then provides the time t to bulk Lorentz factor $\Gamma$ connection through

$$% t(\Gamma) = (m_{\rm p} V \lambda)^{-1} \int_{\Gamma_0}^\Gam... ...\over \; (\tilde\Gamma^{2} - 1)^{3/2}}\; {\rm d}\tilde\Gamma,$$ (21)

so that $t_{\nu_i}$ denotes the time that passes in the rest frame of the jet plasma until the respective zero is reached. Apart from determination of the zeros of $\gamma(\Gamma,\gamma_0)-\Gamma$ the function needs no numerical support.

From Eq. (20) we can now see how the temporal behavior of the particle spectrum is entirely determined by the interplay between cooling and particle loss time scales. The particle injection by sweep-up enters via the mass $M(\nu_i)$ at the point when the relativistic particles have joined the jet plasma.

For a typical particle path (see Fig. 1) there are two zeros of  $\gamma(\Gamma,\gamma_0)-\Gamma=0$ but some characteristic curves fail completely to intersect the line of sweep-up. In between, there is one curve that will just touch this sweep-up line. In this case the particle path produces a singularity because the denominator of Eq. (20) vanishes

$$% \left\vert\; {\partial \gamma\over\partial\Gamma}\bigr\vert _{\nu_i} - 1 \;\right\vert \rightarrow 0,$$ (22)

but the singularity is integrable so that it does not influence the particle density except at an isolated point in $\Gamma$-space. This double zero is called the osculating point and varies with time because it depends directly on $\Gamma(t)$.

In Figs. 2 and 3 we display various particle distributions for some values of Z. The shape of the spectrum is altered in relation to the osculating point. So for calculations involving proton fluxes (as for example, generation of secondaries like GeV and TeV photons, due to decay of neutral pions), it is inappropriate to choose a simple power law for the injection spectra of protons.

At the osculation point of Eq. (22) ${\nu_i},$ we can write:

$$% {\gamma~ \over \tau} = {m_{\rm p} V \lambda~(\nu_i^{2} - 1... ...ff \; \; \dot \gamma = {\partial \Gamma \over \partial t}\cdot$$ (23)

Thus we deal with two competing processes. At the beginning we have a particle deceleration faster than the cooling of the particles by radiation. At that era the freshly incoming particles are the ones with the lowest energy (as shown for example in the right panel of Fig. 3). At the osculation point this behavior changes. That is the point when the deceleration slows so that radiative particle cooling becomes more effective. All particles are then cooler than the freshly swept-up particles. From Fig. 3, we can see that this is a common feature of all solutions. Only the position of the osculation point in Lorentz space $\Gamma$ changes with Z.

	Figure 2: A typical evolution of a particle distribution calculated according to Eq. (20). One can clearly trace the alteration of the shape of the distribution due to the occurrence of the osculating curve. One can also trace the formation of the second peak that happens due to double zeros along some characteristic paths. Note the associated times refer to the jet plasma frame.
Open with DEXTER

	Figure 3: The left panel investigates the solution corresponding to a large value of Z, the right panel shows a typical solution for small values of Z. Because the osculating curve touches the line of sweep-up at a small Lorentz factor, no alteration of shape is seen in the right panel.
Open with DEXTER

   
5 The evolution of high energy photon spectra

The energetic protons in the jet plasma undergo inelastic collisions and produce charged and neutral pions. The neutral pions decay into two $\gamma$-rays, while the charged pions produce high energy neutrinos.

Figure 4: Comparison of the analytical solution (28) with Monte Carlo results from Pohl & Schlickeiser (2000). Additionally, we compute the results of a numerical integration, using Eq. (27). The comparisons are nearly identical to the analytical solution, confirming that the proton distribution determines the behavior of the $F_{\nu }$ spectra. For the computation we chose physical parameters identical to those of Fig. 4 of Pohl & Schlickeiser (2000).

Open with DEXTER

From the the proton spectra derived in the previous section we now deal with the secondary products, obtained as a folding with specific source functions. The integrability depends on the form of the source function and is not simple because of the zeros of Eq. (20), that are dependent on $\gamma$ via $\nu_i(\gamma_0(\gamma,t))$. In this section we provide an alternative method for the folding using the techniques already derived.

The number of photons produced per cubic centimetre per GeV and per second can be calculated via the well-known source function:

$$% \Omega(E_{\rm s},t) = 2~n_{\rm b}~ \int\limits_1^\infty N... ...p},t)\; Q(E_{\rm s},\gamma_{\rm p})~ \;{\rm d}\gamma_{\rm p},$$ (24)

where the differential cross section ${{\rm d}\sigma(\gamma_{\rm p})/ {\rm d}\gamma_\pi}$ is included via the source function for pion production:

$$% Q(E_{\rm s},\gamma_{\rm p})= \int\limits_{\gamma_{\pi ~{\rm... ...(\gamma_{\rm p}) \over {\rm d} \gamma_\pi} {\rm d}\gamma_\pi.$$ (25)

Here we denote $f(\gamma_\pi)=1/\sqrt{\gamma_\pi^2 - 1}$ for the normalized spectrum of the two body decay in the CMS. The lower limit of integration is given by:

$$% \gamma_{\pi~ {\rm min}} = {E_{\rm s}\over m_\pi~c^2}+ {m_\pi~c^2 \over 4~ E_{\rm s}}\cdot$$ (26)

Various models have been discussed (Dermer 1986), for example scaling models (Stephens & Badhwar 1981), or excitation and decay of isobars (Stecker 1970). We use an analytical approximation of the source function based on numerical results using the Monte-Carlo particle interaction code DTUNUC (V2.2) (Möhring & Ranft 1991; Ranft et al. 1994; Ferrari et al. 1996; Engel et al. 1997), which is based on a dual parton model (Capella et al. 1994). One has the representation

$$% Q(E_{\rm s},\gamma_{\rm p}) \approx 5 \times 10^{-15} \exp[-2/\gamma_{\rm p} -0.9~(E_{\rm s}+1.2)^2],$$ (27)

where $E_{\rm s}$ is the logarithm of the energy of the secondary particle i.e. the photon. The time dependence of the emission spectra is completely determined by the changing spectra of the incoming protons (Schuster et al. 2002). The accuracy of the approximation can be estimated by a comparison with numerical results, as done in Fig. 4. As shown in the appendix, one has:

 

$$% {\Omega(E_{\rm s},\Gamma) \over 2~n_{\rm b}} = \exp \left [ -\k... ...r \exp[-\kappa~ t(\tilde\Gamma)]\; (\tilde\Gamma^2 - 1)}\; {\rm d}\tilde\Gamma,$$ (28)

where $\kappa = \left [ {1/ T_{\rm E}} + {1/ T_{\rm N}}\right]$. The integrability depends on the mass $M(\Gamma)$ function.

Figure 4 displays an example of the temporal evolution for the resulting photon spectra with a value of about Z=7.5 $\times$ 10⁶. This value has been calculated with physical parameters chosen identical to those of Fig. 4 of Pohl & Schlickeiser (2000). The number densities inside and outside the jet are $n_{\rm b} = 5$ $\times$ $10^8~\rm{cm^{-3}}$ and $n_{\rm i}^\ast = 0.2~\rm {cm^{-3}}$, respectively, and the thickness of the disk is d = 3 $\times$ $10^{13}~\rm {cm}$, whereas the initial Lorentz factor is $\Gamma_0=300$ and the red-shift z=0.5 enters with the luminosity distance of an Einstein-de Sitter model:

$$% D_{\rm L} = {2~c~ (1+z)\over H_0}~ \left[ 1-{1\over \sqrt{1+z}}\right],$$ (29)

with the Hubble constant of $H =75~{\rm km~s}^{-1}~{\rm Mpc}^{-1}$. The time scale $\tau$ for radiation losses is taken to be about 1.45 $\times$ $10^{15}/n_{\rm b}~[1/{\rm s}].$ For details of the time scales of particle losses we refer the reader to Pohl & Schlickeiser (2000). The Lorentz factors in the plasma frame are 283.62, 208.75 and 106, respectively, whereas the initial Lorentz factor is 300. For a viewing angle $\theta= 0.1^\circ$ these Lorentz factors correspond to 1 h, 10 h and 100 h in an observer frame.

Figure 4 also provides a comparison of the analytical solution (28) with two numerical results, the first calculated as the folding with a source function obtained directly from the Monte Carlo calculation of Pohl & Schlickeiser (2000) and the second being the result of the same numerical routine but fed with the analytic approximation to the source function given by Eq. (27). Figure 4 shows that the solution (28) holds for all times and that all the deviations (which are extremely minor, exept at the very highest energies abouve $10^3~\rm {GeV}$) result from the approximation (27). It is, therefore, more than adequate to use the analytic form of Eq. (20). The formalism can be used for a wider class of source functions. The requirement is:

$$% {\partial Q(E_{\rm s},\gamma) \over \partial \gamma} = C\; Q(E_{\rm s},\gamma)/\gamma^m$$ (30)

with $m \ge 1,$ because the end-point contributions (arising from integration by parts) can be neglected (in the case m=1 only the meaning of $\kappa$ would change). Thus, the formalism can also be applied to neutrino emission.

6 Discussion

This paper has developed an analytical description of the evolution of proton spectra according to the pick-up model, in which a two-stream instability leads to sweep up of ambient matter. To trace the fate of the relativistically incoming particles, we first investigated the time dependent mass of the jet plasmoid. The accretion of mass was shown to depend crucially on the depth dof the jet disk as well as the outer and inner number densities. For quantitative investigation these parameters were combined into a single dimensionless parameter Z. Small values of Z correspond to solutions with a mass accretion lasting until the jet plasmoid decelerates to low bulk Lorentz factors. For very large Z the mass loading can be neglected.

For an analytical description of particle spectra we transformed the variable $\gamma$ (the Lorentz factor of the particles) into a new variable $\gamma _0$, labeling the characteristic path of a particle, tracing the cooling of particles due to radiation losses.

This description further helped to analyze the behavior of the particle spectra in the plasma frame. The solutions are determined by two competing processes. When it takes a long time for the jet plasmoid to decelerate to low bulk Lorentz factors, the particle cooling via radiation becomes more effective. Due to the fact that particles in the jet plasma are always swept up at the Lorentz factor $\gamma =\Gamma$ one can distinguish between "young'' spectra (the deceleration is so fast that the swept-up particles are the ones with the lowest Lorentz factor $\gamma$) and "old'' spectra (the deceleration is now so slow that all particles are cooled to lower energies than the swept-up particles). Between these two extreme behaviors there is a point where both processes are equal, related to a integrable singularity in the solution. This description refers to the plasma frame only.

This paper has dealt with the structure of solutions of the relativistic outflow model, the comparison with data is presented in Schuster & Schlickeiser (2005).

The relativistic particles in the jet plasma give rise to high energy radiation due to decay of neutral pions. To investigate the effects of the proton distribution on the resulting photon spectra, we analytically approximated the results from Monte Carlo claculations for the pion source function. We confirmed that all deviations in the resulting spectra are due to the analytically approximated source function and not to the proton spectrum itself.

Acknowledgements

Partial support by the Bundesministerium für Bildung und Forschung through the DESY, grant 05 CH1PCA 6, is gratefully acknowledged.

References

• Ackermann, M., et al. 2005, Phys. Rev. D, 71, 077102 [NASA ADS] [CrossRef] (In the text)
• Aharonian, F., Akhperjanian, A., Barrio, J., et al. 2002, A&A, 384, L23 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Barthel, P. D., Conway, J. E., Myers, S. T., et al. 1995, ApJ, 444, L21 [NASA ADS] [CrossRef] (In the text)
• Böttcher, M., & Dermer, C. D. 1998, ApJ, 501, L51 [NASA ADS] [CrossRef], Gamma-ray Blazars: An Overview, International Symposium on High Energy Gamma-Ray Astronomy, AIP Conf. Proc. 745, Heidelberg, Germany, 26-30 July 2004 (In the text)
• Bogovalov, S., & Tsingnos, K. 1999, MNRAS, 305, 211 [NASA ADS] [CrossRef] (In the text)
• Capella, A., Sukhatme, U., Tan, C.-I., & Trân Tranh Vân, J. 1994, Phys Rep., 236, 227 [NASA ADS] (In the text)
• Catanese, M., & Weekes, T. C. 1999, PASP 111, 764, 1193 [NASA ADS] (In the text)
• Collmar, W., Reimer, O., Bennett, K., et al. 2000, A&A, 354, 513 [NASA ADS] (In the text)
• Dermer, C. D. 1986, A&A, 157, 223 [NASA ADS] (In the text)
• Dermer, C. D., & Schlickeiser, R. 2002, ApJ, 575, 667 [NASA ADS] [CrossRef] (In the text)
• Engel, R., Ranft, J., & Roesler, S. 1997, Phys. Rev. D, 55, 11, 6957 (In the text)
• Ferrari, A., Sala, P., R., Ranft, J., & Roesler, S. 1996, Z. Phys. C, 70, 413 [CrossRef] (In the text)
• Gaidos, J. A., Akerlof, C. W., Biller, S., et al. 1996, Nature, 383, 319 [NASA ADS] [CrossRef] (In the text)
• Ghisellini, G., Celotti, A., & Costamante, L. 2002, A&A, 386, 833 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Henri, G., & Pelletier, G. 1991, ApJ, 383, L7 [NASA ADS] [CrossRef] (In the text)
• Homan, D. C., Wardle, J. F. C., Cheung, C. C., et al. 2002, ApJ, 580, 742 [NASA ADS] [CrossRef] (In the text)
• Homan, D. C., Lister, M. L., Kellermann, K. I., et al. 2003, ApJ, 589, L9 [NASA ADS] [CrossRef] (In the text)
• Mannheim, K., Biermann, P. L., & Kruells, W. M. 1991, A&A, 251, 723 [NASA ADS] (In the text)
• Mannheim, K., & Schlickeiser, R. 1994, A&A, 286, 983 [NASA ADS] (In the text)
• Mattox, J. R., Wagner, S. J., Malkan, M., et al. 1997, ApJ, 472, 692 [NASA ADS] [CrossRef] (In the text)
• Maraschi, L., & Tavecchio, F. 2005, AIP Conf. Proc. 745, Heidelberg, Germany, 26-30 July 2004 (In the text)
• Mukherjee, R., Bertsch, D. L., Bloom, S. D., et al. 1997, ApJ, 490, 116 [NASA ADS] [CrossRef] (In the text)
• Möhring, H.-J., & Ranft, J. 1991, Z. Phys. C, 52, 643 [CrossRef] (In the text)
• Piner, B. G., & Kingham, K. A. 1997a, ApJ, 479, 684 [NASA ADS] [CrossRef] (In the text)
• Piner, B. G., & Kingham, K. A. 1997b, ApJ, 485, L61 [NASA ADS] [CrossRef]
• Pohl, M., & Schlickeiser, R. 2000, A&A, 354, 395 [NASA ADS] (In the text)
• Pohl, M., Reich, W., Krichbaum, T. P., et al. 1995, A&A, 303, 383 [NASA ADS] (In the text)
• Pohl, M., Lerche, I., & Schlickeiser, R. 2002, A&A, 383, 309 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Punch, M., Akerlof, C. W., Crawley, M. F., et al. 1992, Nature, 358, 477 [NASA ADS] [CrossRef] (In the text)
• Rachen, J. P. 2000, AIP Conf. Proc., 515, 41 (In the text)
• Ranft, J., & Roesler, S. 1994, Z. Phys. C, 62, 329 [CrossRef] (In the text)
• Schlickeiser, R. 2002, Cosmic ray astrophysics (Berlin: Springer Verlag)
• Schlickeiser, R. 2003, A&A, 410, 397 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Schlickeiser, R., Vainio, R., Böttcher, M., et al. 2002, A&A, 393, 69 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Schlickeiser, R., Vainio, R., & Pohl, M. 2003, ApJ, 596, 840 [NASA ADS] [CrossRef] (In the text)
• Schuster, C., Pohl, M., & Schlickeiser, R. 2002, A&A, 382, 829 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Schuster, C. 2004, Ph.D. Thesis
• Schuster, C., & Schlickeiser, R. 2005, AIP Conf. Proc., 745, 555 (In the text)
• Sikora, M. 1997, AIP Conf. Proc., 410, 494 (In the text)
• Siewert, M., Pohl, M., & Schlickeiser, R. 1997, A&A, 425, 405 [EDP Sciences] [NASA ADS], ApJ, 484, 108
• Stecker, F. W. 1970, Ap&SS, 6, 377 [NASA ADS] (In the text)
• Stephens, S. A., & Badhwar, G. D. 1981, Ap&SS, 76, 213 [NASA ADS] (In the text)
• Vainio, R., Pohl, M., & Schlickeiser, R. 2004, A&A, 414, 463 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Vermeulen, R. C., & Cohen, M. H. 1994, ApJ, 430, 467 [NASA ADS] [CrossRef] (In the text)
• von Montigny, C., Bertsch, D. L., Chiang, J., et al. 1995, ApJ, 440, 525 [NASA ADS] [CrossRef] (In the text)
• Waxman, E. 2005, Summary of talk presented at the Nobel Symposium 129: Neutrino Physics, Sweden 2004 [arXiv:astro-ph/0502159] (In the text)

 

  
Online Material

  
Appendix A: Derivation of the main formulae

A.1 The mass function

The procedure to obtain the mass function (14) is based on transforming Eq. (12) into an ordinary differential equation for the mass itself by multiplying Eq. (12) by $\gamma$ and subsequently performing an integration over $1\le\gamma\le\infty,$ where an integration by parts of the second term on the left hand side of Eq. (12) gives:

$$% \int \limits_1^\infty { \gamma~ \partial (\gamma n) \over \... ...gamma = \int \limits_1^\infty\!\!- \gamma ~ n \;{\rm d}\gamma,$$ (A.1)

with respect to the fact that the function  $n(\gamma,t)$ is zero at the upper and lower limits, because, as stated in Eq. (3), all particles are at least at $\gamma = 1+ k~T /(m_{\rm p}~c^2).$ There are no particles at $\gamma=\infty$, so $n(\infty,t)= 0$. The result is an ordinary differential equation for the function  $\int\gamma~ n(\gamma,\Gamma)~{\rm d}\gamma$. Adding the constant term  $n_{\rm b}(1+\epsilon)$ and comparing with Eq. (3), together with the abbreviation $I(\Gamma)= M_{\rm BW} (\Gamma)/(m_{\rm p} V)$ provides:

$$% {1\over I(\Gamma)}{Z + \Gamma \sqrt{\Gamma^2-1} \over (\Gam... ...={ (\Gamma^{2} - 1)^{-3/2} \over \tau_{\rm eff}\lambda }\cdot$$ (A.2)

This linear differential equation for $1/{I(\Gamma)}$ is solved by:

 

$$% {\tau_{\rm eff}\lambda\over I(\Gamma)} -\int\limits_{\Gamma_0} ^{\Gamma} {{\sqrt{\Gamma^2-1}\over(s^2-1)... ...[{{Zs\over {\sqrt{s^2-1}}}-{Z\Gamma\over {\sqrt{\Gamma^2-1}}}}\right] ~{\rm d}s$$ =

$$+{\tau_{\rm eff}\lambda \sqrt{\Gamma^2-1}\over n_{\rm b}~(1\!+\!\... ...\Gamma_0 \over\sqrt{\Gamma_0^2-1}} -{Z\Gamma \over{\sqrt{\Gamma^2 -1}}}\right].$$ (A.3)

Now substitute $u = s/\sqrt{(s^2-1)}$ with the advantage that one can safely set:

$$% \sqrt{u^2 - 1} \simeq \sqrt{(u - 1)~(1 + 1)}= \sqrt{2~u - 2}$$ (A.4)

because ${\it u}$ is always very close to one for $\Gamma \gg 1$. Note also that this approximation is independent of any parameter values and holds for all Z. Now the integral is solved by Eq. (14), with the constant

$$% K= \sqrt{2~x(\Gamma_0)} -{\sqrt{\pi}~ {\rm erfi} \bigl(\sqr... ...} \bigl[ Z~x(\Gamma_0) \bigl]} - {1\over \sqrt{\Gamma_0^2-1}}$$ (A.5)

depending on the initial bulk Lorentz factor $\Gamma_0$.

A.2 The particle spectrum

Abbreviate $\gamma~ N(\gamma,\Gamma)/V$ with $G(\gamma,\Gamma).$ After dividing Eq. (10) by the deceleration ${\partial\Gamma /\partial t}$ and incorporating Eq. (2), the equation for the particle spectrum reads

 

$$% {{\rm d} G \over {\rm d}\Gamma} = {(T -1/\lambda \tau) I(\Gamma... ...{3/2}}~ G(\Gamma ) -{\gamma~ I(\Gamma)~ \delta(\gamma-\Gamma) \over \Gamma^2-1}$$ (A.6)

derived by the method of characteristics. The complete solution to Eq. (A.6) is:

 

$$% G(\Gamma,\gamma_0)= \exp \left[ (T-1/\lambda \tau)\; A(\Gamma)\... ...- 1) \exp\left[ (T- 2/\lambda \tau ) A(\tilde\Gamma)\right]}{\rm d}\tilde\Gamma$$ (A.7)

with the definition:

$$% A(\Gamma) =\int\limits_{\Gamma_0}^\Gamma I(\tilde\Gamma)/(\tilde\Gamma^2-1)^{3/2}{\rm d}\tilde\Gamma$$ (A.8)

so that time t is given in terms of $\Gamma$ through $t(\Gamma) = -A(\Gamma)/\lambda$. The Dirac $\delta$-function in Eq. (A.1) is transformed to:

$$% \delta(\gamma(\Gamma)-\Gamma)= \sum_i {\delta(\Gamma-\nu_i... ...vert {\partial \gamma(\nu_i)/ \partial \Gamma}-1\right\vert}.$$ (A.9)

Here $\nu_i=\nu_i(\gamma_0)$ denotes the different zeros of  $\gamma(\Gamma,\gamma_0)-\Gamma$ which are the values where the characteristic curves intersect the line of sweep up. Equation (A.7) then reads:

$$G(\Gamma,\gamma_0)= - \gamma_0~\exp \left[ (T-1/\lambda \tau)\; A... ...\nu_i^2-1)~ \left\vert {\partial \gamma(\nu_i)/ \partial \Gamma}-1\right\vert},$$ (A.10)

where we denote the Heaviside step-function by $\Theta$. Using Eq. (18) we calculate $\partial \gamma/\partial\Gamma$, the slope of the characteristic curve evaluated at the zero  $\nu_i(\gamma_0)$:

$$% {\partial \gamma(\nu_i)\over\partial\Gamma} ={\partial \ga... ...mbda\tau\right)\over \lambda \tau~ (\nu_i^{2} - 1)^{3/2}}\cdot$$ (A.11)

Equation (20) is obtained by transforming back to the set of variables $\gamma$ and t using the inverse transfornation of Eq. (18) for the variable $\gamma$.

A.3 Secondary particles

Equation (28) is derived with a method similar to the derivation for the mass function. Multiply Eq. (12) by  $Q(E_{\rm s},\gamma_{\rm p})$ to obtain (after an integration over $1\le\gamma\le\infty)$:

$${\partial\Gamma \over \partial t} {\partial \over \partial \Gamma... ...{1/ T_{\rm E}} + {1/ T_{\rm N}}\right] = Q(\Gamma)~ \lambda~ \sqrt{\Gamma^2-1}.$$ (A.12)

The second term in Eq. (A.2) is suitable for an integration by parts:

$$% \int\limits_1^\infty {Q \over\tau} {\partial (\gamma n) \ov... ...r \tau} \; {\partial Q \over \partial \gamma}\; {\rm d}\gamma.$$ (A.13)

The end-point contributions are zero for the same reasons as given previously, and use of Eq. (27) allows one to write:

$$% - \int\limits_1^\infty {\partial Q \over \partial \gamma} ... ...2~ \int\limits_1^\infty {Qn\over \gamma\tau} \; {\rm d}\gamma.$$ (A.14)

One can safely neglect this term as long as $\gamma$ $\gg$ 1. The result is a differential equation for the function $\Omega (E_{\rm s},\Gamma)$(= $\int_1^\infty Q(E_{\rm s},\gamma)~ n(\gamma, \Gamma)~ {\rm d} \gamma$) that we then divide by  $Q_{\rm Es}$ (the part of Eq. (27) that is independent of the Lorentz factor $\gamma$), and use the abbreviation $\int_1^\infty \exp[-2/\gamma]\; n(\gamma,\Gamma)~{\rm d} \gamma = \Xi(\Gamma)$ to derive the ordinary differential equation:

$$% {\partial \Xi(\Gamma) \over \partial \Gamma}= { \kappa ~ \... ...)^{3/2}} - {\exp[-2/\Gamma]\; I(\Gamma) \over (\Gamma^2 - 1)}$$ (A.15)

with $\kappa = \left [ {1/ T_{\rm E}} + {1/ T_{\rm N}}\right]$. The solution to the linear differential Eq. (A.15) is Eq. (28).