Modelling the steady state spectral energy distribution of the BL-Lac Object PKS 2155-304 using a selfconsistent SSC model

In this paper we present a fully selfconsistent SSC model with particle acceleration due to shock and stochastic acceleration (Fermi-I and Fermi-II-Processes respectively) to model the quiescent spectral energy distribution (SED) observed from PKS 2155. The simultaneous August/September 2008 multiwavelength data of H.E.S.S., Fermi, RXTE, SWIFT and ATOM give new constraints to the high-energy peak in the SED concerning its curvature. We find that, in our model, a monoenergetic injection of electrons at $\gamma_0 = 910$ into the model region, which are accelerated by Fermi-I- and Fermi-II-processes while suffering synchrotron and inverse Compton losses, finally leads to the observed SED of PKS 2155-30.4 shown in H.E.S.S. and Fermi-LAT collaborations (2009). In contrast to other SSC models our parameters arise from the jet's microphysics and the spectrum is evolving selfconsistently from diffusion and acceleration. The $\gamma_0$-factor can be interpreted as two counterstreaming plasmas due to the motion of the blob at a bulk factor of $\Gamma = 58$ and opposed moving upstream electrons at moderate Lorentz factors with an average of $\gamma_u \approx 8$.


INTRODUCTION
Among the class of active galactic nuclei (AGN), blazars are showing a spectral energy distribution (SED) that is strongly dominated by nonthermal emission across a wide range of wavelengths, from radio waves to gamma rays, and rapid, large-amplitude variability. Presumably, these characteristics are due to a highly relativistic jet which covers a small angle to the line-of-sight, emitting the observable Doppler-boosted synchrotron and inverse Compton radiation. In SSC models the characteristic double humped spectra of blazars are explained by electrons in the jet emitting Correspondence to: M. Weidinger mweidinger@astro.uni-wuerzburg.de synchrotron radiation while being accelerated in a magnetic field, which gives the first peak in the SED. These high energy electrons upscatter the very same synchrotron photons to TeV energies due to the inverse Compton effect, resulting in a second peak in the SED. Another approach to the explanation of the double humped structure are proton initiated electromagnetic cascades (e.g., Mannheim, 1993) or external Compton models.
The key issue is to understand which physical mechanisms are leading to such SEDs. In particular that means to explain the ultrarelativistic electron spectra within the jet, which are believed to be responsible for the gamma radiation. The high peaked BL Lac objects (HBLs) as a subclass of blazars show a peak in their SED in the X-ray regime, suggesting that an inverse Compton peak should occur at correspondingly high gamma-ray energies. In fact, a large fraction of the known nearby HBLs have already been discovered with Cerenkov telescopes, such as H.E.S.S., MAGIC, and VERITAS. Since 2008 the Fermi satellite measures at these high gamma-ray energies. The energy range of the Fermi data is slightly different from the H.E.S.S. and VERITAS-Telescopes which gives new constraints to the SEDs. The first Fermi data published is from PKS 2155-30.4, a HBL at redshift z = 0.117 (luminosity distance: d L = 1.67 · 10 27 cm) (H. E. S. S. and Fermi-LAT collaborations, 2009).
We present a selfconsistent SSC model that is not only able to model the SED of PKS 2155-30.4 shown in H. E. S. S. and Fermi- LAT collaborations (2009) but also to partly explain the "ad-hoc" injected particle spectra of many SSC models. Therefore we introduce and solve the kinetic equation describing the synchrotron-self-Compton emission numerically in two different zones within the jet (see section 2). We use the exact Klein-Nishina cross section which is important at the relevant very high gamma-energies to describe the inverse Compton radiation and energy losses of the electrons. The emphasis lies on the accurate treatment of the two possible particle acceleration mechanisms (Fermi-Iand Fermi-II) which are able to produce high energy electrons as well as on the selfconsistent treatment of the radiation processes.

Model geometry
We extend the well-established SSC model by Fermi-I and Fermi-II acceleration mechanisms to a selfconsistent SSC model with two zones in a nested setup. Both regions (the acceleration-and the radiation zone) forming the blob are assumed to be spherical and homogeneous containing isotropically distributed non-thermal electrons and a randomly oriented magnetic field. The acceleration zone is assumed to be spatially significantly smaller than the surrounding radiation zone. Furthermore every electron leaving the accelerationzone enters the radiation zone. These assumptions are common place in SSC models (e.g., Kirk et al., 1998). To derive the kinetic equations describing the time evolution of n e (γ), N e (γ) (n e in the acceleration zone, N e in the radiation zone) as the differential electron densities we use the one dimensional diffusion approximation (eq. (1)) of the relativistic Vlasov equation (e.g. Schlickeiser, 2002), which is applicable due to the assumptions made above.
where f (p, t) is a particle distribution function, with the particle's momentum value p. F describes the contributing processes, such as synchrotron radiation or acceleration, in momentum space. Catastrophic particle gains and losses are considered via S(p, t).
Making use of the relativistic approximation E ≈ pc = γmc and the relation n(p, t) = 4πp 2 f (p, t) one can derive the kinetic equations governing the model.

Acceleration zone
While the blob propagates through the jet, electrons are continuously injected into the acceleration zone when considering the blob's rest frame, leading to an injection function which we assume monoenergetic and time independent. These low-to mid-energy electrons are accelerated systematically and stochastically due to Fermi-I and Fermi-II processes while suffering synchrotron and inverse Compton losses. Energy losses due to inverse Compton scattering are calculated using the full Klein-Nishina cross section, see eq.
(15). This leads to P IC (γ) given in eq. (8) with the corresponding radiation field n PH of the acceleration zone. Due to the non equilibrium of magnetic and radiative energy in the acceleration zone the energy losses via inverse Compton scattering can become quite significant and must not be neglected, also the Thomson limit is not appropriate here. The synchrotron losses are calculated using eq. (3) from Ginzburg and Syrovatskii (1969) for isotropic particle distributions P s (γ) = 1 6π with the Thomson cross section σ T . According to Schlickeiser (1984) particle acceleration via parallel shockfronts and stochastic acceleration caused by scattering at Alfvén waves leads to for the function F . With the parallel spatial diffusion coefficient K || , which is momentum independent for hard spheres and the characteristic speeds v A for the Alfvén mediated stochastic acceleration and v S for parallel shockfronts. Substituting p → γ in eq. (1) and eq. (4) according to the relativistic approximation mentioned above, one will finally find eq. (5); the kinetic equation of the acceleration zone.
where the characteristic acceleration timescale t acc is given by Eq. (6) for t acc is a direct consequence of the derivation of eq. (5) out of eq. (1) using eqn. (4) and (3). The expression in eq. (6) includes the analytical timescale for non-relativistic shock acceleration. According to Bednarz and Ostrowski (1996) and especially to Ellison et al. (1990) the acceleration timescale for parallel relativistic shock waves decreases approximately by a factor of 3. We did not take into account this behavior for it is unclear how the analytical expression looks like in that case. Secondly we are omitting the energy dependency of t acc using hard spheres for the plasma instabilities anyway. This issue is irrelevant for the modelling (for we are setting numerical values for t acc ) and the type of energy spectrum (a powerlaw) produced by Fermi-I acceleration is identical in the non-relativistic and relativistic case (Sokolov et al., 2004). But it has to be kept in mind for the interpretation of the results. The parameter a ≈ v 2 s /v 2 A determines the ratio of shock to stochastic acceleration. t esc = ηR acc /c is the characteristic timescale for electrons escaping from the acceleration region, where η is an empirical factor set to η = 10 and R acc the radius of the acceleration sphere. All escaping electrons enter the radiation zone downstream the jet. The seperation in two zones can firstly explain the injected electron spectra and secondly takes account of a much more confined shock region for Fermi-I acceleration will probably not occur in the whole blob region when considering physical sources. Our model can be compared with the model presented by Katarzyński et al. (2006). The kinetic equation (eq. 3 in their paper) is almost similar to the kinetic equation in the acceleration zone eq. 5. One major difference to our model is their sole use of stochastic acceleration. In fact their model is the limit of our model for a → 0. Additionally they limit themselves to radiation in the acceleration zone, which is useful when not taking into account shock acceleration. Besides that there are number of minor differences regarding the exact treatment of inverse Compton losses and the derivation of escape rates. Due to the small spatial extent the acceleration zone does not contribute to the SED directly, i.e. n ph (ν) is only calculated in order to determine the inverse Compton loss rate for the electrons in the acceleration zone.

Radiation zone
The electrons are not accelerated here. Thus the kinetic equation takes the simple form Electrons in the radiation zone suffer synchrotron (∝ β s γ 2 ) and inverse Compton losses (eq. (8)), other energy losses are irrelevant in jetsystems of such low electron density (see e.g. Böttcher and Chiang, 2002).
The integrals in eq. (8) are solved numerically using the full Klein-Nishina cross section for a single electron given in eq. (15). The photon energies are rewritten in terms of the electrons rest mass, i.e. hν = αmc 2 for the scattered photons and hν = α 1 mc 2 for the target photons. The integration bounds of the outer integral in eq. (8) are a direct consequence of the kinematics. The non-trivial dependency of P IC (γ) from N ph and thus of N e from eq. (8) makes the numerical treatment of the kinetic equations inevitable leading to a time resolved model. The loss rates for the electron distribution of PKS 2155 in the steady state are shown in Fig. 1, which indicates that in our case the inverse Compton losses would be slightly overestimated in the often used Thomson limit for all Lorentzfactors γ because of the dependency on the photon field P thom ∝ dννN ph γ 2 and its special shape due to the modified injection at γ 0 = 910. This would not be the case for low energetic injected electrons and the resulting powerlaw-like photon distribution. For high Lorentzfactors γ however the deviation of the Thompson limit for the inverse Compton scattering to the real Klein-Nishina treatment becomes more significant in each case. Again electrons escaping the blob are taken into account via t esc,rad = ηR rad /c with the empirical factor η set to η = 10. Both electronic PDEs are connected via the catastrophic particle loss/gain-term. The factor (R acc /R rad ) 3 ensures particle conservation.
To determine the time-dependent spectral energy distribution of blazars we solve the differential equation for the differential photon number density, obtained from the radiative transfer equation, including the corresponding terms with respect to the SSC model, To describe the synchrotron photon production rate R s in a convenient way we use the well known Melrose approxima-tion (see e.g., Pohl, 2002) with the characteristic frequency eq. (11) for an electron with a Lorentz factor of γ.
In optically thick regimes the emitted synchrotron radiation is absorbed by the emitting electrons itself. This is described by the synchrotron self absorption coefficient, (with γ c = f (ν c ) −1 ). Here we made use of the monochromatic approximation (e.g., Felten and Morrison, 1966) for the synchrotron power: In SSC models the second hump in the SED of a blazar is due to inverse Compton scattered photons by the synchrotron radiation emitting electrons themselves. Here the full Klein-Nishina cross section from Blumenthal and Gould (1970) is used to calculate the inverse Compton photon production rate.
The last catastrophic term in eq. (9) describes photons escaping from the emitting region, where is the approximate escape time, with R rad the radius of the emitting blob. The escape time is chosen to be the light crossing time of the photons. The photon lossrate due to the pair production of electrons and positrons is not taken into account for two reasons.
Firstly it is insignificant compared to the dominating synchrotron and inverse Compton processes. This is a consequence of the relatively low density . Secondly it would violate the selfconsistency of our model for positrons are not treated, hence violating energy conservation.
To compute the SEDs in our model we must shift the frame of reference from the blob to the observer. For a sphere of radius R the observed flux at distance r is With the Lorentz boosted intensity I obs ν = δ 3 I blob ν due to the bulk motion with a doppler factor δ of the blob and the Lorentz transformed, red shifted frequency ν obs = δ/(1 + z)ν. Where I blob ν is calculated from the photon unit density for homogenous spheres.

NUMERICS
In our model we numerically solve the kinetic equations forward in time in order to obtain a model SED. The downstream motion of the electrons induces the sequence of solving the acceleration zone's equation before the kinetic equation of the radiation zone in each time step. The simple Euler scheme was found adequate to do the time integration.
In the acceleration zone we had to combine the Crank-Nicholson scheme (Press, 2002) with Godunov's method to provide both correct treatment of the characteristics and stability for the derivation in γ. In the radiation zone the characteristic flows, due to the absence of acceleration, only in one direction making the Crank-Nicholson scheme sufficient. With our carefully tested code it is possible to calculate the dynamics of SEDs in a range of 20 orders of magnitude. The implemented code complies particle conservation in each zone alone and both together as well as the conservation of the total energy (i.e. of the electrons and the photons) over typical simulation times with a maximum error of O(5%). For negligible stochastic acceleration, i.e. a → ∞, and without a radiation field, i.e. no inverse Compton losses, the steady state solution for the kinetic equation yields n e,steady (γ) = C 1 γ 2 1 γ − β s t acc tacc−tesc tesc with γ max = (t acc β s ) −1 and the constant C determined by the injection function Q. The implemented numeric model converges against this solution for sufficient simulation time.
Additionally it was tested against the steady state analytical solution with Fermi-II processes given in Schlickeiser (1984) with no significant deviations. Setting t esc,rad → ∞ and neglecting inverse Compton scattering the spectral index of the powerlaw part of the electron distribution in the radiation zone is (analytically) reduced by one compared to the one in the acceleration zone, which also was confirmed by the implemented code. The inverse Compton scattering rate was confirmed against the approximate analytical results (low energetic Thompson regime and extreme Klein-Nishina limit) before implementing. Concerning the photon distributions we validated the expected spectral indices in the steady state solution for the different frequency regimes, which together with the energy conservation between electrons and photons approves the integrity of the model. A detailed description of the used numeric techniques as well as the implemented model also in context with the variability of the sources will be given in a paper yet to be published.

RESULTS
The recent Fermi data give new constraints on the gammaray peak of the HBL PKS 2155-30.4 concerning its curvature. This is leading to a deep dip between the optical/X-ray and the gamma-ray peak. We are able to model the SED of PKS 2155-30.4 with our model by setting for the monoenergetic injection into the acceleration zone. This is rather unusual but required to model the SED of PKS 2155. Such moderate but not small Lorentz factors can be explained e.g. by two counterstreaming plasmas. If the upstream electrons would be at rest, the bulk doppler factor of δ = 116 would automatically lead to γ 0 = Γ ≈ 58. Assuming speculatively that the upstream electrons moving in the opposite direction of the blob with a mean velocity of v u hence a upstream Lorentz-Factor γ u the γ 0 factor in the blob's rest frame must be calculated according to the relativistic superposition: Solving eq. (20) for our setup we find γ u ≈ 8 for the upsteam electrons which are streaming towards the blob. The numerically solved steady state electron density in the acceleration zone is shown in Fig. 2. We also show the time development for a "switched on" injection, i.e. n e (γ)| t<0 = 0 ∀ γ and Q = Q 0 δ(γ − γ 0 )ϑ(t), until the steady state is reached. In Fig. 2 it can clearly be seen tvhat accelerating electrons  Fig. 3 as arising from the injection function Q0. The corresponding intrinsic times are t = 1000 s (dashed black), t = 5000 s (dashed blue), t = 1 · 10 4 s (dashed red), t = 2 · 10 4 s (dashed green). The steady state with its rising and falling powerlaw and the exponential cutoff at γ ≈ 10 5 is reached at about t = 10 5 s.
using Fermi-I and Fermi-II processes leads to powerlaw electron distributions with an exponential cut-off often used as the ad-hoc injection function in onezone-SSC models (Chiang and Böttcher, 2002) (right side of γ 0 in Fig. 2), thus explaining them using the diffusion theory derived from plasma physics. By injecting electrons with eq. (19) and significant stochastic acceleration (i.e. a = O(1)) we are also able to produce rising electron spectra before decreasing in a powerlaw and an exponential cut-off, like introduced in Böttcher and Chiang (2002) (left side of γ 0 in Fig. 2). The Fermi-II processes are responsible for the rising power-law and exponential cut-off, whereas the ratio of t acc /t esc determines the spectral index of the power-law at γ > γ 0 . It can clearly be seen from Fig. 2 that the convergence against the steady state solution for the electron density begins relatively rapid while slowing down eventually. The simulation time when the steady state is reached corresponds to the escape time of the electrons in the acceleration zone. When concerning variability and time resolved lightcurves of blazars this is an advantage of the twozone model because the rising part in such lightcurves corresponds partially to the escape time of the acceleration zone (while the falling part is connected to the response time of the system, t rad,esc ). An acceleration zone electron density as shown by the solid black curve in Fig. 2, leads to the desired broken power-law electron spectrum in the radiation zone which finally is able to model the SED of PKS 2155-30.4 (see Fig. 3 and Fig. 4). We used the parameters in Table 1 for the model SED in Fig.  4 (black, solid line). The black dashed curve in Fig. 4 corre-   Fig. 4). The intrinsic times are t = 1 · 10 4 s (dashed black), t = 2 · 10 4 s (dashed blue), t = 5 · 10 4 s (dashed red), t = 1 · 10 5 s (dashed green), the complete steady state is reached at about t = 2 · 10 6 s, which correlates to the response time of the system due to tesc,rad.
sponds to a fit assuming a black body for the thermal contribution of the host galaxy thus the ATOM optical data is not to be taken into account for the SSC modelling. The curvature and deep dip in the model SED is a direct consequence of the rising part in the electron density of the acceleration zone. Thus it can be modeled by varying the ratio a of shock to stochastic acceleration. All the parameters in Table 1 Aharonian et al. (2005) data of H.E.S.S. a few years ago which show the same flux level as the recent data. We used the EBL studies described in Primack et al. (2005) to do the EBL deabsorption for the H.E.S.S. datapoints, a correction of the Fermi data is not necessary.
The time development of the SED due to a switched on in-jection of electrons into the acceleration zone at time t 0 = 0 is shown in Fig. 3. It can clearly be seen that the final state of the model SSC correlates with the response time of the radiation zone t esc,rad and that the convergence again begins fast and slows down rapidly at higher simulation times.

CONCLUSIONS
Our model is able to explain the injection function of many onezone-SSC models as shock and stochastic acceleration of electrons upstream the jet entering the blob while continuously suffering synchrotron losses. By introducing Fermi-II acceleration we get rid of the sharp cut-off introduced in Kardashev (1962) or Kirk et al. (1998) which probably does not occur in physical sources. Additionally we are able to model relatively complex electron densities with increasing and decreasing parts through the stochastic acceleration of electrons, only by varying the monoenergetic injection to higher γ 0 . In contrast to the ad-hoc injection of some onezone-SSC models such Lorentz factors have a physically reasonable, but highly speculative, explanation as upstream previously accelerated but already partially cooled electrons. These electrons are averagely moving in the opposite direction of the blob with a mean Lorentz factor of γ u ≈ 8 resulting, together with the motion of the blob, in γ 0 ≈ 900 for the monoenergetic injection function used in the acceleration zone of our model.
As recent data points out, these complex electron distributions are necessary to model the new constraints concerning the gamma-ray peak of blazar's SEDs if one does not simply shift the synchrotron-peak to achieve the inverse Compton spectrum (e.g., Kataoka et al., 2000). The curvature of the peak, and thus the deep dip between the two humps, is a direct consequence of the rising part in the responsible electron distribution within the blob. This constraint rules out many SSC models, which are not able to produce such electron spectra. With our model we are able to form the curvature of the gamma-ray peak and the dip by varying the influence of the Fermi-II processes. The shape and position of the synchrotron peak in the model SED is dominated by t acc and R acc , R rad as well as B. For the parameters concerning the acceleration arise from plasmaphysics considerations we gain insight into the jets microphysics while modelling observed SEDs. We have also shown that in such environments the Thomson approximation for the inverse Compton effect can not always be applied, especially when considering time resolution and hence non equilibria of the energy distribution in the blob.
Here we only introduced steady state solutions of our model, but due to the spatially relatively small acceleration region, which is at least an oder of magnitude smaller than the emitting region, this twozone-SSC model is able to selfconsistently model the rising part in the lightcurves of flaring blazars which are connected to the behavior in the acceleration zone, especially the energy transport from low to high energies. This, together with the consequences of the model geometry on the observable SEDs and lightcurves of blazars like in Sokolov et al. (2004), will be subject of a following paper.   Primack et al. (2005). The dashed black curve shows a thermal fit for the contribution of the host-galaxy. Our model SSC fit, arising from the steady state electron distribution in the radiation zone is shown in the solid black curve, a moderate energy injection at γ0 ≈ 910 into the acceleration zone together with stochastic and systematic acceleration is needed to meet the curvature of the VHE peak given via the Fermi data.