Contents

• Numerical simulations of radiation from blazar jets
• 1 Introduction
• 2 The model
• 3 Numerical implementation
• 4 Model input parameters
• 5 Application: Blazar 3C 279
• 6 Summary
• Appendix A: Electron injection function
• Appendix B: Broad emission line region
• References

A&A 406, 855-865 (2003)
DOI: 10.1051/0004-6361:20030794

Numerical simulations of radiation from blazar jets

R. Moderski - M. Sikora - M. Bazejowski

Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland

Received 17 May 2002 / Accepted 20 May 2003

Abstract
We present a description of our numerical code BLAZAR. This code calculates spectra and light curves of blazars during outbursts. The code is based on a model in which the non-thermal flares in blazars are produced in thin shells propagating down a conical jet with relativistic velocities. Such shells may represent layers of a shocked plasma, enclosed between the forward and reverse fronts of an internal shock. In the model adopted by us, the production of non-thermal radiation is assumed to be dominated by electrons and positrons which are accelerated directly, rather then injected by pair cascades. The code includes synchrotron emission and inverse-Compton process as the radiation mechanisms. Both synchrotron photons and external photons are included as the seed photons for Comptonization. At the present stage, the code is limited to treat the inverse Compton process only within the Thomson limit and is specialized to model radiation production in the flat spectrum radio quasars. As an example, we present the results of modeling an outburst in 3C 279 - the most extensively monitored $\gamma$-ray - bright quasar.

Key words: galaxies: active - galaxies: jets - galaxies: individual: 3C 279 - gamma rays: theory - radiation mechanisms: non-thermal

1 Introduction

One of the greatest achievements of the recently retired Compton Gamma-Ray Observatory (CGRO) is the discovery of gamma ray emission from a subclass of Active Galactic Nuclei (AGNs) known as blazars. More than 50 such sources were detected at energies above 100 MeV (Mukherjee et al. 1997). Energy constraints and $\gamma$-ray absorption arguments require that radiation emission in blazars must be beamed (Mattox et al. 1993). Indeed, it is now widely believed that the entire electromagnetic spectrum observed in these objects is dominated by non-thermal radiation produced in a jet pointing close to the line of sight (Dondi & Ghisellini 1995). Although roughly divided into two classes: flat spectrum radio quasars (FSRQ) and BL Lacs objects, all blazars share common characteristics: large-amplitude, rapid variability; smooth continuum emission in all observable bands; and high linear polarization.

Spectrum of a blazar consists of two broad components (von Montigny et al. 1995). The low energy component has a peak within IR-to-X-ray range and is usually attributed to Doppler-boosted synchrotron radiation. The high energy component peaks in the MeV-TeV energy range and is very likely produced by inverse-Compton process. Both components are highly variable, with time scales ranging from years to a fraction of a day. Analysis of $\gamma$-ray light curves seems to suggest that variability patterns of blazars are a superposition of short term flares (Magdziarz et al. 1997). The flares detected in different spectral bands appear to be correlated (Macomb et al. 1995; Wagner et al. 1995; Wehrle et al. 1998). This behavior suggests co-spatial production of the high and low energy components and indicates that significant fraction of the jet energy is dissipated in the localized events at sub-parsec distances from the center. Such events are likely to arise in internal shocks (Sikora et al. 1994; Spada et al. 2001) or at the sites of the magnetic field reconnection (Romanova & Lovelace 1992; Blackman 1996).

Both the shocks and the reconnection sites provide favorable conditions for efficient acceleration of particles. Unfortunately, the present theories of particle acceleration are still not sufficiently developed to provide quantitative predictions regarding such issues as what fraction of dissipated energy is used to accelerate particles, how much power is channeled into relativistic protons vs. relativistic electrons/positrons, and what maximum energy can a particle gain. In particular, it is not possible to deduce solely on theoretical grounds whether protons are accelerated to energies sufficient to support - via inelastic collisions with photons - the pair cascades, as is suggested by hadronic models (Mannhein & Biermann 1992), or if they are "radiatively'' passive returning all energy gained during the dissipative events back to the flow via adiabatic expansion. In the latter case, represented by the leptonic models, radiation production in blazars is totally dominated by electrons and positrons accelerated directly. There are several electron/positron acceleration mechanisms which are likely to operate in shocks (Levinson 1996; Hoshino et al. 1992; McClements et al. 1997; Shimada & Hoshino 2000; Hoshino & Shimada 2002) and reconnection sites (Larrabee et al. 2003; Lyutikov 2002), but such theoretical considerations cannot as yet predict the spectra of accelerated particles and their minimum and maximum energies.

Theories of radiation mechanisms are much better developed. They are quantitative and for a given energy distribution of the accelerated/injected particles, geometry and kinematics of the source, and external radiation field environment, one can make specific predictions regarding the radiation spectra. By confronting those theories with observations, it is possible to verify the model and determine its parameters. Using this approach, one can already exclude some radiation scenarios. This concerns e.g. models with pair cascades, which predict X-ray spectra that would be too soft as compared with observations of some FSRQ (Sikora & Madejski 2001). Therefore, for these objects, direct electron/positron acceleration mechanism is favored.

In leptonic models, both the high-energy and the low-energy components of the spectrum are produced by the same population of relativistic electrons. At low energies, the emission comes from the synchrotron process, while at higher energies, it is dominated by the inverse Compton scattering. There is still some debate about the source of seed photons for the inverse Compton process. The most obvious choice involves synchrotron photons from the low energy component. Models based on this assumption are called Self-Synchrotron-Compton (SSC) models. Originally proposed by Königl (1981) to explain production of X-rays and $\gamma$-rays (already detected by then from 3C 273 by COS-B, Swanenburg et al. 1978), currently the SSC models prove to be successful in explaining general features of BL Lac objects (Ghisellini & Maraschi 1989; Takahashi et al. 1996; Kirk et al. 1998; Tavecchio et al. 1998; Kino et al. 2002). However, in quasars - but also, in some radio selected BL Lac objects, the external radiation fields may be sufficiently dense to dominate the Compton cooling of electrons in a jet (Madejski et al. 1999). Models exploring this hypothesis are called External-Radiation-Compton (ERC) models. Several sources of photons for external radiation field have been considered: direct disc radiation (Dermer & Schlickeiser 1993); diffuse radiation from broad emission line (BEL) region (Sikora et al. 1994); infra-red radiation from hot dust (Bazejowski et al. 2000); or jet synchrotron radiation scattered back to the jet by the external gas (Ghisellini & Madau 1996).

In this paper we present a code developed to simulate non-thermal flares produced by thin shells propagating down a conical jet with relativistic speeds. Such shells approximate the geometry of relativistic plasma enclosed between the forward and reverse shock fronts formed by colliding inhomogeneities in a jet (Sikora et al. 2001; Sikora & Madejski 2002). The code is sufficiently general to treat radiation processes for any radial distribution of magnetic and external radiation fields, and for any electron/positron injection function. It includes both Comptonization of synchrotron radiation and Comptonization of external radiation, and takes into account adiabatic losses due to 2D conical expansion of the shocked plasma sheets. The main limitations of the present version of the code are that Comptonization is treated self-consistently only within the Thomson limits, and that $\gamma\gamma$-pair production is not included. Since there are observational indications that in FSRQ the $\gamma$-ray spectra have a high energy break in the 4-10 GeV band (Pohl et al. 1997; Sikora et al. 2002), while synchrotron spectra are usually steep in the UV band, the above limitations affect only marginally our models of quasars. However, they can be significant for BL Lac objects, particularly for those bright in the TeV band.

Assumptions and main features of our model, such as the equations describing evolution of the energy distribution of relativistic electrons, radiation processes, relativistic aberration and light travel effects are presented in Sect. 2. The numerical code is described in Sect. 3. Approximate analytical formulas, expressing the model input parameters as a function of observables, are derived in Sect. 4. In Sect. 5 we present results of application of our code to model the outburst of 3C 279, the best studied quasar in the $\gamma$-ray band. The work is summarized in Sect. 6.

2 The model

In our model a source of the non-thermal radiation is assumed to be a thin shell, propagating down the conical jet with a constant speed, and radiation is produced by relativistic electrons/positrons, injected into a shell within a given distance range. This picture approximates the internal shock scenario, where shocks are formed following collisions of inhomogeneities moving with different radial velocities (Sikora et al. 1994; Spada et al. 2001). Structure and dynamics of such internal shocks is in general very complex and depends on a number of parameters, such as the relative velocity of inhomogeneities, their densities, temperatures, geometry, and total masses. The shock structure can be double, single (forward or reverse), or can initially it can be double, followed by an evolution into a single shock (Daigne & Mochkovitch 1998; Bicknell & Wagner 2002). Regardless, in the case of inhomogeneities which are cold prior to their collision and which have comparable masses and comparable rest frame densities, the constant speed of the shocked plasma is a good approximation. Small thickness of the shocked plasma shell can be justified, since internal shocks are at most mildly relativistic (Sikora et al. 2001).

2.1 Electron evolution

Assuming that electron injection function and energy densities of magnetic field and of external radiation fields are uniform across the shell, one can follow evolution of electron energy distribution by solving the kinetic equation for the total population of relativistic electrons, despite the fact that each element of the conically diverging shell has its own rest frame which is different from other frames. That equation can be written in the form (Moderski et al. 2000):

$${\partial N_{\gamma} \over \partial t'} = - {\partial \over \... ... \left(N_{\gamma} {{\rm d}\gamma \over {\rm d}t'}\right) + Q ,$$ (1)

where $Q \equiv ({\rm d}N_{\gamma}/{\rm d}t')_{\rm inj}$ is the electron injection function, $\gamma$ is the random electron Lorentz factor, ${\rm d}\gamma/{\rm d}t'$ is the rate of the electron energy losses, ${\rm d}t' = {\rm d}t /\Gamma$ is the proper time, and $\Gamma$ is the Lorentz factor of the bulk motion of a shell. Within the thin shell approximation, the time derivatives can be replaced with derivatives over the radial distance r. Using a relation ${\rm d}r = \beta_{\Gamma} c \Gamma {\rm d}t'$, where $\beta_{\Gamma} = \sqrt{\Gamma^2-1}/\Gamma$, we obtain

$${\partial N_{\gamma} \over \partial r} = - {\partial \over \p... ...ma \over {\rm d}r}\right) + {Q \over c \beta_{\Gamma} \Gamma},$$ (2)

where

$${{\rm d}\gamma \over {\rm d}r} = {1 \over \beta_{\Gamma} c \G... ... {\rm d}t'\right)_{\rm rad}- {2 \over 3}{\gamma \over r} \cdot$$ (3)

The second term on the right-hand side of Eq. (3) represents the adiabatic losses due to two-dimensional (conical) expansion of the shell.

2.2 Radiative cooling rates

There are three processes contributing to the radiative energy losses of electrons:

synchrotron radiation,

$$\left({\rm d}\gamma \over {\rm d}t'\right)_{\rm S} = - {4\sigma_{\rm T} \over 3 m_{\rm e} c } u_{\rm B}' (\beta \gamma)^2 ;$$ (4)

Comptonization of synchrotron radiation,

$$\left({\rm d} \gamma \over {\rm d} t'\right)_{\rm SSC} = - {4... ...ma_{\rm T} \over 3 m_{\rm e} c } u_{\rm S}' (\beta \gamma)^2 ;$$ (5)

and Comptonization of external radiation,

$$\left({\rm d}\gamma \over {\rm d}t'\right)_{\rm ERC} = - {4 \... ...{\rm T} \over 3 m_{\rm e} c } u'_{\rm ext} (\beta \gamma)^2 ,$$ (6)

where $u'_{\rm B} = B'^2/8\pi$ is the magnetic energy density, $u'_{\rm S}$ is the energy density of the synchrotron radiation field, and  $u'_{\rm ext}$ is the energy density of the external radiation field, all as measured in the source comoving frame. B' is the strength of the magnetic field, which can be an arbitrary function of the radius.

Energy density of the synchrotron radiation produced by a thin shell is

$$u_{\rm S}' \simeq {(t'-t_0')\delta L_{\rm S}' \over \delta V'... ...j}} \simeq {\delta L_{\rm S}' \over 2 c r^2 \delta \Omega_{j}}$$ (7)

where t'-t₀' is the time that elapsed from the beginning of the shock formation, $\lambda' \simeq 2 c (t'-t_0')$ is the width of the synchrotron radiation layer at that time, and

$$\delta L'_{\rm S} = \int_{\nu_{\rm S,min}'}^{\nu_{\rm S,max}'} ~ \delta L_{{\rm S},\nu'}' {\rm d} \nu' ,$$ (8)

is the synchrotron luminosity produced in a cell of the shell, enclosed within a solid angle $\delta \Omega_{j} \ll \Omega_{j} \equiv \pi \psi_j^2$ (see Fig. 1). The monochromatic synchrotron luminosity, $\delta L_{{\rm S},\nu'}'$, is computed using Eq. (10) and integrated over the frequency range from $\nu_{\rm S,min}' \simeq \nu'_{\rm abs}$, at which the optical thickness due to the synchrotron self-absorption is equal to unity, to  $\nu_{\rm S,max}' = (2e/ 3 \pi m_{\rm e} c) \gamma_{\rm max}^2 B'$.

Energy density of the external radiation field is

 

$$u_{\rm ext}' {1 \over c} \int {I_{\rm ext}' {\rm d}\Omega'} = {1 \over c} \int {{I_{\rm ext} \over {\cal D}_{\rm in}^2} {\rm d}\Omega}$$ =

$$\simeq {\Gamma^2 \over c} \int {I_{\rm ext} (1-\beta_{\Gamma}\cos{\theta_{\rm in}})^2 {\rm d} \Omega} ~$$ (9)

where $\theta_{\rm in}$ is the angle between the trajectory of the incoming photon and the direction of the cell motion. For narrow jets this angle can be approximated by the angle between the trajectory of the incoming photon and the jet axis.

Figure 1: A schematic picture illustrating the integration procedure. A thin shell of electrons enclosed within the half-opening angle of the jet, $\psi _{j}$ is divided into cells, each $\delta \Omega _{j}$ thick. Each such cell contains  $\delta N_\gamma$ electrons. The observer is located at the angle  $\psi _{\rm obs}$ from the jet axis. Integration is performed by summation over $\theta$ contributions from cells located at different radii.

2.3 Intrinsic luminosities

2.3.1 Synchrotron radiation

The rate of the synchrotron radiation production per cell for a given electron distribution is (Chiaberge & Ghisellini 1999):

$$\delta L_{{\rm S},\nu'}'(r) = \int_{\gamma_{\rm min}}^{\gamma... ...\delta N_{\gamma}(r) F_{\rm S}(\nu',\gamma) {\rm d} \gamma ~ ,$$ (10)

where $\delta N_{\gamma} = N_{\gamma} (\delta \Omega_{j}/\Omega_{j})$, and

$$F_{\rm S}(\nu',\gamma) = {3 \sqrt{3} \over \pi} {\sigma_{\rm ... ...i \left [ K_{4/3}^2(\chi) - K_{1/3}^2(\chi) \right ] \right \}$$ (11)

is the single electron synchrotron emissivity averaged over an isotropic distribution of pitch angles (Crusius & Schlickeiser 1986), where $\chi=\nu'/(3 \gamma^2 \nu'_{\rm B})$and $\nu_{\rm B}' = e B'/(2 \pi m_{\rm e} c)$. The power $F_{\rm S}(\nu',\gamma)$ is taken with the low energy cutoff at $\nu_{\rm abs}'$.

Note that because the time scale of the electron/positron gyration in the local tangled magnetic fields is much shorter than the time scale of radiative energy losses (Dermer & Schlickeiser 1993), we can assume that the electron momentum distribution is isotropic in the cell co-moving frame. This does not imply that the synchrotron radiation field is isotropic. In general, one must transform the synchrotron radiation from neighbor cells into cell frame taking into account the change in frequency and direction of the incoming photons. However, for a thin shell and small opening angles of the jet the assumption about isotropy of the synchrotron radiation field is a reasonable approximation.

2.3.2 Synchrotron-self-Compton

Using $\delta$-function approximation, one can find (Chiang & Dermer 1999)

 

$$\delta L_{{\rm SSC},\nu'}' (r) \simeq {\sqrt 3 h c \sigma_{\rm T}... ...amma}(r) n_{{\rm S},\nu_{\rm S}'}' \nu_{\rm S}'^{-1/2} ~ {\rm d} \nu_{\rm S}' ,$$ (12)

where

$$\nu_1' = {\rm Max} \left [ \nu_{\rm abs}'; {3 \nu' \over 4 \gamma_{\rm max}^2} \right ] ,$$ (13)

and

$$\nu_2' = {\rm Min} \left [ \nu_{\rm S, max}'; {3 \nu' \over 4... ...rm min}^2} ; {3 (m_{\rm e} c^2)^2 \over 4 h^2 \nu'} \right ] ,$$ (14)

where $n_{{\rm S},\nu_{\rm S}'}'$ is the number density of the synchrotron photons per frequency. Noting that for isotropic synchrotron radiation field $n_{{\rm S},\nu'}' = u_{{\rm S},\nu'}' /h \nu'$ and

$$u_{{\rm S}, \nu'}' \simeq {\delta L_{{\rm S},\nu'}' \over 2 r^2 \delta \Omega_{j} c} ~$$ (15)

(see Eq. (8)), we finally obtain

 

$$\delta L_{{\rm SSC}, \nu'}' {\sqrt 3 \sigma_{\rm T} \over 8 r^2 \delta \Omega_{j}} {\nu'}^{1/... ...right ] \delta L_{{\rm S},\nu'}' \nu_{\rm S}'^{-3/2} ~ {\rm d} \nu_{\rm S}' ~ .$$ (16)

This is analogous to Eq. (27) of Chiang & Dermer (1999).

2.3.3 External-radiation-Compton

In the co-moving frame, the external diffuse radiation field is strongly anisotropic. Due to relativistic aberration, the external radiation appears to the jet as mostly incident from the forward direction. Within the jet, a photon scattered by a relativistic electron follows the direction of motion of that electron. Because of this, the observer will detect preferentially photons that were emitted by the electrons which during the scattering process had their momentum vector pointing at the observer. In such an approximation, the energy of the photon scattered into the observer direction is given by (Reynolds 1982):

$$\nu' \simeq \gamma^2 \nu_{\rm ext}' (1 + \cos \theta') ,$$ (17)

where $\nu_{\rm ext}' \simeq \Gamma \nu_{\rm ext}$, $h \nu_{\rm ext}$is the characteristic energy of photons in the external radiation field (assumed to have a narrow spectrum), and $\theta'$ is the angle between the bulk motion direction of a given shell segment and direction to the observer, as measured in the co-moving frame of the bulk motion. For  $\gamma \gg 1$ the rate of energy losses of an electron with a velocity vector instantaneously oriented into direction of the observer is:

$$\left \vert {\rm d} \gamma \over {\rm d} t' \right \vert_{\rm... ...over m_{\rm e} c} u_{\rm ext}' (1 + \cos \theta')^2 \gamma^2 .$$ (18)

Note that the angle-averaged value of $({\rm d} \gamma / {\rm d} t')_{\rm ERC} [\theta']$ gives Eq. (6) with $\beta \simeq 1$. Since the "blazar phenomenon'' implies observing the jet at very small angles, $\theta \ll \pi/2$, the co-moving observing angle $\theta'$ is much smaller than  $\pi - 1/\Gamma$ and then (Dermer 1995)

$$\Gamma (1 + \cos \theta') \simeq \Gamma (1 + \beta_{\Gamma} \cos \theta') = {\cal D}.$$ (19)

With that approximation, Eqs. (17) and (18) can be rewritten in the form

$$\nu' \simeq {\cal D} \gamma^2 \nu_{\rm ext} ~ ,$$ (20)

and

$$\left \vert{\rm d} \gamma \over {\rm d} t'\right \vert_{\rm E... ...m ext}' \gamma^2 \left({{\cal D}\over \Gamma }\right)^2 \cdot$$ (21)

The power per unit solid angle and frequency, emitted by a cell via the ERC process into the observer direction is, as measured in the cell co-moving frame,

$${\partial \delta L_{{\rm ERC},\nu'}'[\theta'] \over \partial ... ...} [\theta'] m_{\rm e} c^2 {{\rm d} \gamma \over {\rm d} \nu'},$$ (22)

and after combining with Eqs. (20), (21), and (9), we obtain

 

$${\partial \delta L_{{\rm ERC},\nu'}'[\theta'] \over \partial \Ome... ...nt{ (1-\beta_{\Gamma}\cos{\theta_{\rm in}})^2 ~ I_{\rm ext} ~ {\rm d} \Omega} .$$ (23)

2.4 Observed spectra and light curves

The monochromatic radiation flux observed at a given instant is

$$F_{\nu_{\rm obs}}(t_{\rm obs}) = {(1+z) \over d_{\rm L}^2} ~ ... ...nu} [t, \theta] \over \partial \Omega_{{\vec n}_{\rm obs}}} ~$$ (24)

where

$${\partial \delta L_{\nu} [t, \theta] \over \partial \Omega_{... ...[r, \theta'] \over \partial \Omega_{{\vec n}_{\rm obs}}'} ~ ,$$ (25)

$$\nu' = \nu /{\cal D} = \nu_{\rm obs}(1+z)/{\cal D} ~ ,$$ (26)

$$r = {c \beta_{\Gamma} \over 1 - \beta_{\Gamma} \cos \theta} (t-t_0) + r_0 ~ ,$$ (27)

$$\cos \theta = {\cos \theta' + \beta_{\Gamma} \over 1 + \beta_{\Gamma}\cos \theta'} ~ ,$$ (28)

$$t=t_{\rm obs}/(1+z) ~ ,$$ (29)

$${\cal D}= \Gamma (1 + \beta_{\Gamma} \cos{\theta'}) \equiv {1 \over \Gamma (1 - \beta_{\Gamma} \cos{\theta})} ~$$ (30)

and $d_{\rm L}$ is the luminosity distance. For the $\Omega_{\lambda} + \Omega_m = 1$ cosmology,

$$d_{\rm L} = {c \over H}(1+z) \int_0^z {{{\rm d}x \over [\Omega_{\lambda} + \Omega_{\rm m}(1+x)^3]^{1/2}}},$$ (31)

where H is the Hubble constant.

It should be noted here that because the synchrotron and SSC emission is isotropic in the cell co-moving frame, the co-moving power per solid angle produced by those processes is

$${\partial \delta L_{{\rm S,SSC},\nu'}'[\theta'] \over \partia... ...ec n}_{\rm obs}}'} = {\delta L_{\rm S,SSC,\nu'}' \over 4 \pi},$$ (32)

where $\delta L_{\rm S,\nu'}$ and $\delta L_{\rm SSC,\nu'}$ are given by Eqs. (10) and (12), respectively.

3 Numerical implementation

Our numerical method is similar to that used by Chiaberge & Ghisellini (1999). Electron evolution Eq. (2) is solved with the implicit difference scheme adopted from Chang & Cooper (1970). First, the uniformly spaced logarithmic energy grid is established:

$$\gamma_k = \gamma_0 \left ( {\gamma_N \over \gamma_0} \right )^{k - 1 \over n - 1} ~~~~~~ k = 0,...,n$$ (33)

with energy intervals $\Delta \gamma_k = \gamma_{k+1/2} - \gamma_{k-1/2}$. Using a definition $N^i_k = N_\gamma(\gamma_k,i \Delta r)$ we can convert Eq. (2) into a set of linear equations:

 

A_k Nⁱ⁺¹_k-1 + B_k Nⁱ⁺¹_k + C_k Nⁱ⁺¹_k+1 = Uⁱ_k , (34)

where the coefficients are

$$\begin{array}{l} A_k = 0 ~ ,\\ [3mm] B_k = 1 + \left ( \displ... ...splaystyle{{\Delta r \over \Delta \gamma_k}} ~ ,\\ \end{array}$$ (35)

and the source term is

$$U^i_k = N^i_k + {\Delta r \over c \beta_\Gamma \Gamma} Q^i_k ~ .$$ (36)

The system of Eqs. (34) yields a special case of tridiagonal matrix with the upper diagonal equal to zero. It is solved numerically by matrix decomposition, forward- and back-substitution. In practice we use tridag subroutine from Numerical Recipes (Press et al. 1992). This method of solving Eq. (2) is unstable if the sharp cutoffs are present in the electron energy distribution function. In our case, because Eq. (3) does not contain the acceleration terms, we avoid this instability by setting the upper boundary of the calculation grid to the maximum energy of the injected electrons. Another way to avoid sharp cutoffs is to introduce very steep power-law tails to mimic these cutoffs (Böttcher & Chiang 2002).

There exists an alternative way of solving Eqs. (34). If one assures that the energy grid is sufficiently wide to have $N(\gamma_n,i \Delta r) = 0$ during the entire evolution, then Eq. (34) can be solved recursively from the highest to the lowest energies using the relation

$$N^{i+1}_k = {1 \over B_k} \left ( U^i_k - C_k N^{i+1}_{k+1} \right ) ~~~~~ k = n-1, ... ,1.$$ (37)

For universality, however, we use the full tridag method. It allows the consideration of electron distributions with high energy tails without performing the calculations for unimportant grid points.

In order to calculate the electron distribution it is necessary to estimate the electron cooling function (3). This requires the knowledge of the energy density of the synchrotron radiation field which depends on the electron distribution itself. In our code this calculation is done iteratively. Initially $u'_{\rm S} = 0$ is assumed and the electron distribution is calculated using this value. Then this distribution is used to calculate synchrotron radiation density together with the frequency of synchrotron self absorption. The new value of  $u'_{\rm S}$ is used to recalculate the electron distribution. This process is repeated until convergence.

The electron distribution is then used to calculate the synchrotron luminosity (10), SSC luminosity (12), and ERC luminosity (22). For all integrations we use the Romberg's method (Press et al. 1992). The advantage of this method is that it adjusts itself, and thus only a minimal number of calculations is performed to achieve desired accuracy.

Special care must be taken when calculating light curves from Eq. (24). The radiation reaching the observer at a given time is a superposition of radiation emitted at different radii, and this in turn depends on $\theta$, which is the angle between the direction of motion of a given cell and the direction to the observer (see Eq. (27)).

One must also remember that the choice of the cell size puts a constraint on the minimum timescale that can be probed with the code. For a given cell size, $\Delta \theta$, the minimum timescale for which spatial homogeneity within the cell is maintained equals

$$\Delta t \simeq r \sin \theta_{j} \Delta \theta.$$ (38)

4 Model input parameters

The input parameters of the model are:

• the distance of the shock formation from the central source, r₀;

• the distance range of shock operation, $\Delta r_{\rm coll}$;

• the bulk Lorentz factor of the shock, $\Gamma$;

• the energy density of an external diffuse radiation field, $u_{\rm ext}$ and its characteristic frequency  $\nu_{\rm ext}$;

• the magnetic field in a shocked plasma, B';

• the electron injection function, Q;

• the angle between the line of sight and the jet axis, $\psi _{\rm obs}$;

• the half-opening angle of a jet, $\psi _{j}$.

Initial values of these parameters can be set up by using approximate analytical formulas which relate them to the observable quantities. Depending on the available observables, a different set of equations must be used. Most useful observables used to model the short term (1-10 days), $\gamma$-ray dominated outbursts in FSRQ are: the time scale of the outburst/flare, $t_{\rm fl}$ (more specifically, its growing part); radiation fluxes at the low-energy and high-energy spectral peaks, denoted respectively as  $F_{\rm LE} \equiv \nu_{\rm LE} F_{\nu_{\rm LE}}$ and  $F_{\rm HE} \equiv \nu_{\rm HE} F_{\nu_{\rm HE}}$; location of the spectral break between the X-ray and $\gamma$-ray spectral portions, $\nu_{\rm c}$; radiation flux, $\nu_{\star}F_{\nu_{\star}}$, at some $\nu_{\star} > \nu_{\rm c}$; the spectral index of the HE component at  $\nu > \nu_{\rm c}$ (as measured, e.g., by EGRET in its bandpass); and the total radiation flux of broad emission lines, $F_{\rm BEL}$ (Celotti et al. 1997). For the sake of simplicity of analytical formulas which will be used to derive the initial model parameters as a function of observables, both the light travel effects (due to finite size of the source) and the Doppler dispersion effects (due to the "conical'' expansion of a shell) are ignored. In addition, we assume that: during the outbursts the production of $\gamma$-rays is dominated by the ERC process; $B' \propto 1/r$; and electrons are accelerated in the shock at a constant rate and with a single power-law energy distribution,

$$Q=K \gamma^{-p} ~$$ (39)

for $\gamma_{\rm min}=1$ and

$$\gamma_{\rm max}= {1 \over {\cal D}} \sqrt{\nu_{\rm max}(1+z) \over \nu_{\rm BEL}}$$ (40)

(see Appendix A for justification of such an injection function).

Due to radiative and adiabatic energy losses of electrons, their energy distribution evolves with time. The strong energy dependence of the radiative energy losses of electrons causes steepening of the electron energy distribution. It takes place at $\gamma > \gamma_{\rm c}$, where  $\gamma_{\rm c}$ is the energy at which the time scale of electron energy losses,

$$t_{\rm cool}' \simeq {\gamma \over \vert {\rm d} \gamma/{\rm ... ...{\rm rad} + \vert {\rm d} \gamma/{\rm d}t'\vert_{\rm ad}} ~ ,$$ (41)

is equal to the lifetime of a shock

$$t_{\rm sh}' \simeq {t_{\rm fl} {\cal D} \over 1+z} \cdot$$ (42)

During $\gamma$-ray dominated outbursts the radiative energy losses are dominated by the ERC process, and therefore $t_{\rm cool}'=t_{\rm sh}'$gives

$$\gamma_{\rm c} = {m_{\rm e} c ~ (1+z) f(k) \over 2 \sigma_{\rm T} t_{\rm fl} u_{\rm ext}' {\cal D}} ~ ,$$ (43)

where f(k) =(1+3k)/ (2(1+k)) and $k\equiv r_0/\Delta r_{\rm coll}$.

At $\gamma =\gamma_{\rm c}$ the power law energy distribution of electrons, $N_{\gamma} \propto \gamma^{-s}$, changes the slope from s=p at  $\gamma < \gamma_{\rm c}$ (slow cooling regime) to s=p+1 at  $\gamma > \gamma_{\rm c}$ (fast cooling regime). Such a change of slope results from the fact that in the fast cooling regime ( $t_{\rm rad} < t_{\rm sh}$) the number of electrons is saturated by the radiative losses, i.e. $N_{\gamma} \sim \int_{\gamma} Q ~ {\rm d}\gamma / \vert {\rm d}\gamma/{\rm d}t'\vert \sim Q t_{\rm rad}'$. (Note that  $\gamma_{\rm c}$ is changing during the shock progression and that its particular value given by Eq. (43) is for the instant when the collision is complete.)

The break at $\gamma_{\rm c}$ is reflected in the ERC spectral component at

$$\nu_{\rm c,obs} = {\gamma_{\rm c}^2 {\cal D}^2 \nu_{\rm ext} \over 1+z} ,$$ (44)

around which the electromagnetic spectrum changes the slope by $\Delta \alpha \simeq 0.5$. We identify this break as the one observed more or less directly in FSRQ between 1 and 30 MeV. Such identification provides information about energy density of external radiation field as measured in the shock co-moving frame. Using Eqs. (43) and (44), we obtain

$$u_{\rm ext}' (r_{\rm f})= { m_{\rm e} c \over 2 \sigma_{\rm T... ...} f(k) \over t_{\rm fl} (\nu_{\rm c}/\nu_{\rm ext})^{1/2}} ~ ,$$ (45)

where $r_{\rm f} \simeq r_0 + \Delta r_{\rm coll}$. Since observations provide directly the ratio $F_{\rm HE}/F_{\rm LE}$ and since

$${F_{\rm HE} \over F_{\rm LE}} = {F_{\rm ERC} \over F_{\rm S}}... ...' \over u_{\rm B}'} \left ({{\cal D} \over \Gamma} \right)^2 ,$$ (46)

we can estimate the intensity of the magnetic field,

$$B'(r_{\rm f}) \simeq {{\cal D} \over \Gamma} ~ \sqrt{8\pi u_{\rm ext}' {F_{\rm LE} \over F_{\rm HE}}} ~ ,$$ (47)

where $u_{\rm ext}'$ is given by Eq. (45).

The short term flares (lasting typically 1-10 days) are most likely produced at distances 0.1-1.0 pc. At such distances, the external photon energy density  $u_{\rm ext}'$ is dominated by broad emission lines (see, e.g., Sikora et al. 2002). With this, and assuming spherical geometry of the BEL region, we have

$$u_{\rm ext}'(r) \simeq {4 \over 3} \Gamma^2 {(\partial L_{\rm BEL}/\partial \ln r)(r) \over 4\pi r^2 c }\cdot$$ (48)

Provided that $(\partial L_{\rm BEL}/\partial \ln r)(r)$ is known, this equation, combined with Eq. (45) and with the formula for the distance range of the shock operation,

$$\Delta r_{\rm coll} \simeq {c t_{\rm fl} \over 1 - \beta_{\Ga... ...quiv {c t_{\rm fl} \Gamma^2 \over 1+z}{{\cal D}\over \Gamma} ,$$ (49)

can be used to estimate $r_{\rm f}$ and $\Gamma$. For a broad emission line luminosity distribution approximated by (see Appendix B)

$${\partial L_{\rm BEL} \over \partial \ln r} (r) = {L_{\rm BEL... ...cr (r/r_{\rm BEL})^{-0.5} & for $r > r_{\rm BEL}$\space \cr} ~$$ (50)

where

$$r_{\rm BEL} \simeq 10^{18} (L_{\rm UV}/10^{46}~{\rm ergs ~ s}^{-1})^{0.7} ~ {\rm cm} ,$$ (51)

is the distance at which the radial distribution of the BEL luminosity has a maximum (Kaspi et al. 2000), we obtain

$$r_{\rm f} \simeq \cases{ c_{\rm r}^2 (\nu_{\rm c}/\nu_{\rm BE... ...{\rm BEL}^{1/3} & if~~~ $r_{\rm f} > r_{\rm BEL}$\space \cr} ,$$ (52)

where $\nu_{\rm BEL}$ is a characteristic frequency of the broad emission lines radiation, and

$$\Gamma = \sqrt{{(1+z) r_{\rm f} \over (1+k) c t_{\rm fl} ({\cal D}/\Gamma)}} ,$$ (53)

where $c_{\rm r} = \sigma_{\rm T} (1+z)^{1/2}/[6\pi m_{\rm e} c^3 (1+k)f(k)({\cal D}/\Gamma)]$, $L_{\rm BEL} = 4\pi d_{\rm L}^2 F_{\rm BEL}$, and $L_{\rm UV} = 4\pi d_{\rm L}^2 F_{\rm UV}$ is the luminosity of the accretion disc. The disc radiation is sometimes observed directly (during low states of FSRQ); otherwise it can be estimated from $F_{\rm UV} \simeq F_{\rm BEL}/\xi_{\rm BEL}$, assuming the canonical value of the BEL covering factor, $\xi_{\rm BEL} =0.1$.

In order to determine the parameters of the electron injection function, we use $\gamma$-ray data from the spectral band at  $\nu \gg \nu_{\rm c}$. There, deeply in the fast cooling regime,

$$p=2 \alpha_{\gamma} ~$$ (54)

and

$$N_{\gamma} \simeq \left \vert{1 \over \dot \gamma}\right \vert_{\rm tot} \int_{\gamma} Q ~ {\rm d}\gamma ~ .$$ (55)

Noting that for $\gamma$-ray dominated outbursts, $\vert \dot \gamma\vert_{\rm tot} \simeq \vert \dot \gamma\vert_{\rm ERC}$, by using approximate formula for ERC radiation production,

$$4 \pi d_{\rm L}^2 \nu_{\star} F_{\nu_{\star}} \simeq 0.5 [\ga... ...\vert_{\rm tot} m_{\rm e} c^2 \Gamma^4 ({\cal D}/\Gamma)^6 ~ ,$$ (56)

one can find that the normalization factor of the electron injection function is

$$K \simeq {8 \pi d_{\rm L}^2 [\nu_{\star}F_{\nu_{\star}}] \ove... ...rm e} c^2 \Gamma^4 ({\cal D}/\Gamma)^6} g(\alpha_{\gamma}) ~ ,$$ (57)

where $g(\alpha_{\gamma}) = (2\alpha_{\gamma}-1)\gamma_{\star}^{2(2\alpha_{\gamma}-2)}$and

$$\gamma_{\star} = {1 \over {\cal D}} \sqrt{ \nu_{\star}(1+z) \over \nu_{\rm ext}}\cdot$$ (58)

Note that for $\alpha_{\gamma} = 1$, $g(\alpha_{\gamma})= 1$.

To complete the set of the model input parameters we still need to specify the values of  ${\cal D}/\Gamma$, k and $\psi _{j}$. Since the ERC to synchrotron luminosity ratio is  $\propto ({\cal D}/\Gamma)^2$(see Eq. (46)), the view angles  $\psi_{\rm obs} \le 1/\Gamma$ are expected for strongly $\gamma$-ray dominated FSRQ. In such conditions ${\cal D}/\Gamma=1$ ranges from 2 for the observer located on the jet axis to 1 for the observer located at  $\psi_{\rm obs} = 1/\Gamma$.

The value of k is also expected to be of the order of unity, as suggested by roughly symmetrical profiles of flares (see, e.g., Sikora et al. 2001). Regarding the jet opening angle, we know they are very small (1-3 degrees) on kilo-parsec scales, but can be much larger on parsec/sub-parsec scales (see, e.g. Lobanov 1998). The specific value of $\psi _{j}$ can be determined from the value of flux in the soft/mid X-ray bands, provided that the latter is dominated by the SSC process. Since we do not follow the SSC process in our analytical approach, and because other model parameters are not very sensitive to $\psi _{j}$, the numerical simulations can be started with any $0 <\psi_{j} <1/\Gamma$ and then corrected in subsequent iterations.

In our method to calculate the model input parameters, we do not take advantage of the ratio of frequencies, $\nu_{\rm HE}$ to  $\nu_{\rm LE}$. This ratio often has been used to calculate the bulk Lorentz factor. If it were to be incorporated into our method to calculate $\Gamma$, this would give an equation for k. However, $\nu F_\nu$spectra around their maxima are quite flat and therefore, both  $\nu_{\rm HE}$ and  $\nu_{\rm LE}$ are subject to very large uncertainties. Furthermore, in FSRQ the synchrotron peak is located in the far IR and rarely is observed directly.

An additional difficulty is the fact that the amplitudes of the $\gamma$-ray flares are usually much larger than the amplitudes of the synchrotron flares. This suggests that synchrotron component is strongly diluted by the quasi-steady radiation produced at larger distances in a jet. We can deal with this case in our method by using in the Eq. (47) $\Delta F_{\rm syn}$ instead of  $F_{\rm syn}$, where  $\Delta F_{\rm syn}$ is the amount by which synchrotron flux increases during the high amplitude $\gamma$-ray flare. For the instances where the dilution is so large that $\Delta F_{\rm syn}$cannot be determined observationally, as an alternative approach, one can use the observed correlation of the soft/mid X-ray flux with the $\gamma$-ray flux. When the soft/mid X-ray flux is interpreted as produced by the SSC process, one can estimate the magnetic field intensity via an iterative process, which relies on the match of the model SSC radiation spectrum against that measured in the X-ray band.

5 Application: Blazar 3C 279

3C 279 is one of the most extensively observed blazars. Several times it was a target of multi-wavelength campaigns. The most fruitful campaign took place in the beginning of 1996 when the blazar underwent an enormous $\gamma$-ray flare. This event was monitored nearly simultaneously at many frequencies (Wehrle et al. 1998) and is analyzed below using our code.

5.1 Input parameters

3C 279 is located at a redshift z=0.538, which for the currently favored cosmology ( $\Omega_\lambda \simeq 0.7$, $\Omega_{\rm m} \simeq 0.3$ and $h \simeq 0.66$) gives the luminosity distance $d_{\rm L} \simeq 10^{28}$ cm. The data collected during the February 1996 flare (Wehrle et al. 1998; Hartman et al. 2001) give the $h\nu_{\rm c} \simeq 20$ MeV and $\nu_{\star}F_{\nu_{\star}} \simeq 2.75 \times 10^{-10}$ erg s^-1 cm^-2, where $\nu_{\star} \simeq 10^{23}$ Hz. The slope of the $\gamma$-ray spectrum during the flare in the EGRET band was $\alpha_{\gamma} \simeq 0.97$ which from Eq. (54) gives $p \simeq 1.9$. Although a reliable EGRET spectrum is available for this object, the high energy break is not directly observed, so we set $h\nu_{\rm max} = 20$ GeV. Estimate of  $F_{\rm BEL}$ is very uncertain. Analysis based on relative line intensities by Celotti et al. (1997) gives $L_{\rm BEL} = 3.4 \times 10^{44}$ erg s^-1. However, during outbursts this luminosity can be even larger (Koratkar et al. 1998). For our analysis we set $L_{\rm BEL} = 6.8 \times 10^{44}$ erg s^-1. The average photon energy of external radiation field is taken to be $h \nu_{\rm ext} = h \nu_{\rm BEL}=10$ eV. Luminosity of the accretion disc, obtained using the direct observations of the UV bump during the low state of 3C 279 (Pian et al. 1999), is $L_{\rm UV} \simeq 2.7 \times 10^{45}$ erg s^-1. Time scale of the flare is $t_{\rm fl} \simeq 1$ day (Lawson et al. 1999).

Analysis of the lightcurves of the February 1996 flare at different energies (Wehrle et al. 1998) shows that despite the very large rise of flux in the $\gamma$-ray band, the flare was almost undetectable at low energies, especially from radio to optical. This suggests that synchrotron component is significantly diluted, and synchrotron radiation produced during the flare is hidden by a more steady component. For this reason we initially set $F_{\rm syn}/F_{\rm ERC} = 10^{-2}$. The above set of observables yields (for the assumed k=1 and ${\cal D}/\Gamma = 2$) the following input parameters: $u_{\rm ext}'(r_{\rm f}) \simeq 0.21$ erg cm^-3 (Eq. (45)); $B'(r_{\rm f}) \simeq 0.46$ Gauss (Eq. (47)); $r_{\rm BEL} \simeq 4.0 \times 10^{17}$ cm (Eq. (51)); $r_{\rm f} \simeq 4.2 \times 10^{17}$ cm (Eq. (52)); $\Gamma \simeq 7.9$ (Eq. (53)); $K \simeq 3.4 \times 10^{48}$ s^-1(Eq. (57)); $\gamma_{\rm max} \simeq 3.5 \times 10^3$ (Eq. (40)). We have also set $\psi_{j} = 0.12 \simeq 1/\Gamma$.

5.2 Average spectrum

Figure 2: The average spectrum of the blazar 3C 279 during the February 1996 flare. Data points are from Hartman et al. (2001). Thick, solid line shows the averaged spectrum of our model (see text for parameters). Thin lines represent various components of the spectrum.

We first apply our model to the average spectrum of the blazar 3C 279 during the February 1996 flare. The result is presented in Fig. 2. The fine-tuning of the input parameters requires substantial change only of the normalization of electron injection function, K. We have achieved the best results setting $K = 9.5 \times 10^{49}$. This difference can easily be understood by examining Eq. (57). This quantity is very sensitive to ${\cal D}/\Gamma$. The calculation of the input parameters assumes that the source is point-like and thus that the whole radiation is beamed toward the observer. Due to the conical structure of the jet the Doppler factor decreases to the edge of the jet, and thus more electrons are required to produce the observed luminosity. The minor correction was applied also for the magnetic field. The best result yields $B'(r_{\rm f}) = 0.3$ Gauss. For the calculation we also set $\psi_{\rm obs} = 0.05$, but the averaged spectrum is not very sensitive to this value. The averaging was performed between r = r₀ and r = 3r₀.

Large discrepancy between the data and our model in the low energy component arises from the fact that the synchrotron radiation during the studied flare is substantially diluted, presumably by radiation produced at larger distances in a jet.

5.3 Light-curves

Figure 3: $\gamma$-ray (upper panel - $\nu F_\nu$ flux at 400 MeV) and X-ray (lower panel - integrated 2-10 keV flux) light-curves of the blazar 3C 279 during the February 1996 flare. Data points are from Wehrle et al. (1998). Thick, solid line shows the light-curves of our model.

Figure 3 presents the simulated light-curves together with observational data points. Two curves are shown: one is the $\nu F_\nu$ flux at 400 MeV (EGRET range) and the second is the integrated 2-20 keV flux (XTE range). For both lightcurves the steady components were added. These components were estimated from the pre-flare observations. In the case of X-ray light curve the quiescent emission is $F_{\rm q} = 1.23 \times 10^{-11}$ erg s^-1 cm^-2(Lawson et al. 1999) and contributes up to 25% of the peak flux. Thus for the X-ray band, where the radiation is presumably dominated by the SSC process, the contribution of the quiescent emission is significant. For gamma-rays we estimated the steady component from data points before 30th day of 1996. Its value is $\nu F_{\rm\nu,q} = 24.7 \times 10^{-11}$ erg s^-1 cm^-2and its contribution to the peak flux is only 12%.

The agreement between the data and our model is remarkable taking into account the relative simplicity of the model and complexity of the real phenomena as suggested by, e.g., the gamma-ray light-curve. The discrepancy at the rising part of the X-ray light curve suggests an existence of the soft X-ray precursors caused by Comptonization of the external radiation by cold electrons (Sikora & Madejski 2002).

5.4 Spectral evolution

Figure 4: Spectral index $\alpha (F_\nu \propto \nu ^{-\alpha })$ evolution (left panels) and hardness-intensity diagram (right panels) of the blazar 3C 279 during the February 1996 flare. Upper panels show the index at 400 MeV while lower panels show the index at 6 keV.

In order to demonstrate spectral changes during the flare, in Fig. 4 we present the evolution of the spectral index $\alpha$ (defined by $F_\nu \propto \nu^{-\alpha}$) at two energies: 400 MeV and 6 keV. Two plots are presented for each energy: one is the power law spectral index $\alpha$ as a function of time, the second is the so called hardness-intensity diagram showing the spectral index as a function of the flux.

Unfortunately, the comparison with observations in this case is difficult. There are no observations of the spectral index variation during the flare by EGRET. The reported slope of the $\gamma$-ray spectrum for the event was $\alpha_\gamma = 0.97$ which is in excellent agreement with the middle part of the curve in upper left corner of Fig. 4. For the XTE observations, $\alpha_x$ ranges from 0.68 to 1.26, but the errors are large, and there is a possible evidence of systematic error of 0.1 (Lawson et al. 1999). Even with the errors taken into account, the observations seem to disagree with results presented in lower left panel of Fig. 4. However, we note that in this range, the influence of the quiescent component is significant, and the proper comparison would require the analysis which is beyond the scope of this paper. The steady component may also influence the spectrum at $\gamma$-ray energies and cause the discrepancy for observations outside of the maximum of the flare.

5.5 Evolution of electron distribution

Figure 5: Evolution of the electron energy distribution during the February 1996 flare of the blazar 3C 279. Step between the curves is $\Delta r = 0.2~r_0$.

In Fig. 5 we present the time evolution of energy distribution of electrons during the flare. The evolution is followed from r = r₀ till the electrons reach the distance 3 r₀. At the very beginning of the flare, the number of all electrons increases, then the number of the most relativistic electrons starts to saturate being balanced by the radiative energy losses. After the injection of electrons stops, the high energy tail of the electron distribution decays very rapidly, while only small changes can be noticed at lowest energies.

The characteristic feature of this evolving electron spectrum is the spectral break, which moves from higher to lower values as evolution proceeds. At a distance r=2 r₀, where injection stops and the  $\gamma^2 N_{\gamma}$ reaches maximum, the approximate location of the break is $\gamma \simeq 65$. The value calculated from Eq. (43) is somewhat larger ( $\gamma_{\rm c} \simeq 105$) and this is because the approximate analytical formula for  $\gamma_{\rm c}$ does not include synchrotron and SSC energy losses. At energies lower than  $\gamma_{\rm c}$ (the slow cooling regime), the power-law energy distribution of electrons $N_{\gamma} \sim \gamma^{-s}$ has an index s=p, while for $\gamma \gg \gamma_{\rm c}$ (the fast cooling regime) the slope of the electron distribution (due to efficient radiative cooling) is steeper by unity, i.e., s=p+1. The latter effect is mostly obscured by the effect of the break at  $\gamma_{\rm max}$.

6 Summary

In this article, we have presented our numerical code BLAZAR which simulates light-curves and spectra of blazars during flares. In the code, the structure of the source responsible for the production of flares is approximated by thin uniform shells, propagating with constant speed down the conical jet. The shells are filled with relativistic electrons/positrons, which are injected within a given distance range at a constant rate. Evolution of the electron energy distribution is treated by the kinetic/continuity equation, given by Eq. (1). In that equation, the electron injection function Q is separated from the cooling term. This simplification can be justified provided time scale of acceleration of electrons is much shorter than time scale of their cooling and time scale of the shock operation. Since in FSRQ, which are main targets of the present version of the code, the cooling break is located at $\nu_{\rm c} \ll \nu_{\rm max}$, the above assumption is satisfied for all electrons except for those with  $\gamma \sim \gamma_{\rm max}$. However, this condition may not be satisfied for the X-ray selected BL Lac objects, including TeV blazars. In these objects time scale of acceleration can well be comparable to the lifetime of the source, and then Eq. (1) must be modified in the way presented by Kirk et al. (1994) (see also Kirk et al. 1998; Böttcher & Chiang 2002 for application to blazars).

Very short gyration time scale (as compared with the radiative energy losses) plus very likely presence of chaotic magnetic fields (because of MHD turbulence induced around the shock fronts) justify the assumption regarding the isotropic distribution of the electron momenta, and therefore regarding the isotropy of the synchrotron and SSC emission in the source co-moving frame. However, the isotropy approximation doesn't apply to Comptonization of external radiation. This is because due to relativistic aberration, the external diffuse radiation is seen by the source as strongly beamed from the front. Because of this, the scattered radiation is anisotropic - it is beamed into the source propagation direction already even in the source co-moving frame. This anisotropy, originally pointed out by Dermer (1995), is self-consistently treated in our code. Radiative effects of the ERC anisotropy, often ignored in other models, are particularly strong for $\psi_{\rm obs} > 1/\Gamma$.

Another advantage of our code is that it includes the adiabatic energy losses in the kinetic equation for electrons. In our model they are related to the "conical'' expansion of a shell. Inclusion of adiabatic losses is very important in calculating properly the location of the cooling break  $\nu_{\rm c}$. Those losses also strongly affect the light curves of flares observed at $\nu < \nu_{\rm c}$(Sikora et al. 2001).

As an example of application of our code, we presented in Sect. 5 the results of our modeling of the short term outburst observed in 3C 279 in February 1996 (Wehrle et al. 1998). Our results demonstrate very strong dilution of the synchrotron radiation component in this object. This effect, often ignored completely by other models, can be quite common, as suggested by much smaller amplitudes of optical flares than of $\gamma$-ray flares (Ulrich et al. 1997). For such objects, detailed fitting of the synchrotron spectrum in terms of the homogeneous/one-component model is inappropriate and can lead to very large errors in the model parameters. Our method, albeit approximate, takes into account the dilution effect.

It is worth noting that the code "BLAZAR'' was already used in several previous papers. It was applied to the study the dependence of amplitude of flares in blazars on the frequency, on the ratio  $r_0/\Delta r_{\rm coll}$, and on the time profile of an electron injection function Q(t') (Sikora et al. 2001). The code was also used to study the relative role of hot dust and broad emission line region as a source of the external seed photons for the inverse Compton process in a jet (Bazejowski et al. 2000) and was applied to demonstrate the possible unification between the MeV-blazars and GeV-blazars (Sikora et al. 2002). Finally, after some small modifications allowing to model the radiation from the shells propagating with a variable bulk Lorentz factor and supplementing the code by relevant dynamical equations, it was used to study the effect of lateral expansion of the collimated GRBs, possibly resulting in the observed breaks in their light-curves (Moderski et al. 2000).

In the present form, the code doesn't include the Klein-Nishina effects nor the $\gamma$-ray absorption and e⁺e^- pair creation. In FSRQ the Klein-Nishina effects become important at $\gamma > {\rm a\ few} \times 10^4/\Gamma$, while the absorption - at $h\nu \ge 20$ GeV, both due to interactions with external BEL photons. The high energy break at $h\nu \sim 4{-}10$ GeV, suggested by EGRET observations of FSRQ during the $\gamma$-ray flares (Pohl et al. 1997), and the steep UV spectra seem to indicate that energy distribution of electrons/positrons has a break or cutoff around $\gamma \sim {\rm a\ few} \times 10^4/\Gamma$ (see Appendix A). Hence, the KN and  $\gamma\gamma$ pair production effects are expected to be at most marginal. The above limitations, particularly the first one, are expected to be relevant only for BL Lac objects, in particular for those with strong TeV emission. In the case of our model of the 3C 279 flare $\Gamma \gamma_{\rm max}h\nu_{\rm BEL}/ m_{\rm e} c^2 \simeq 0.5$ and $hv_{\rm max} = 20$ GeV, so the effects are expected to be negligible.

Acknowledgements

This work was partially supported by the Polish Committee for Scientific Research grants no. 5 P03D 002 21 and PBZ 054/P03/2001. We are grateful for the hospitality of the Stanford Linear Accelerator Center, operated by Stanford Universty for the Department of Energy under contract no. DE-AC3-76SF00515, where some of the research described above was performed. We would like also to thank G. Madejski for his valuable comments which helped to improve the paper.

Appendix A: Electron injection function

In most FSRQ, the high energy spectra of flares can be recovered assuming a single power-law injection function. Due to cooling effect, the electron energy distribution evolves into the broken power-law, which is reflected in the spectrum of the ERC component, with the break located in the 1-30 MeV range (see Sect. 5). The low energy tails observed in several FSRQ do not show any change of the slope down to the lowest X-ray energies, and this seems to indicate that the power-law electron injection function extends down to $\gamma_{\rm min} \le$few (Tavecchio et al. 1998; Sikora et al. 2002). This is independently supported by the detection of circular polarization in radio cores of some quasars (Homan et al. 2001), provided such polarization is generated by the Faraday conversion process (Wardle et al. 1998; Ruszkowski & Begelman 2002).

It is less clear how electron injection function behaves at the highest energies and what is the maximum energy of injected electrons. EGRET data of FSRQ during flares suggest existence of a high-energy break at 4-10 GeV (Pohl et al. 1997). This high-energy break cannot be caused by the Klein-Nishina effect, because in the fast cooling regime, the luminosity of the Compton component is determined by the electron injection function and reduction of the Compton scattering cross-section in the Klein-Nishina regime can be substantially compensated by the respective increase of the electron densities. In particular, for electron injection function $\propto \gamma^{-2}$ the Klein-Nishina portion of the Compton spectrum has exactly the same spectral index, $\alpha=1$, as in the Thomson regime (Zdziarski & Krolik 1993). The high-energy break also cannot be produced by absorption of $\gamma$-rays by external photons. This is because the radiative environment is transparent for photons detected within the EGRET band (it becomes opaque only for $h\nu \ge 20$ GeV, due to absorption of such photons by optical/UV broad emission lines). Hence, the break in 4-10 GeV band, if real, must be related to the break in the electron injection function, at $\gamma \sim {\rm a\ few} \times 10^4/\Gamma$.

Is the injection function extending much above that break, but with a steeper slope as suggested by observations of the so called MeV-blazars (Sikora et al. 2002), or, there is a cutoff? One can try to answer that question by studying the spectra in the UV band, where the high energy tails of the synchrotron component dominate. However, the task is not easy, because synchrotron radiation in this band can be affected by several factors: by extinction in the host galaxy; by flattening of energy distribution of highest energy electrons due to Compton scatterings in the KN regime (Dermer & Atoyan 2002) and by  $\gamma\gamma$ pair production; and, as in the case of 3C 279, by dilution of the synchrotron component by radiation produced at other locations in a jet. Hence, the possible extension of the injection function up to energies $\gamma \gg {\rm a\ few} \times 10^4/\Gamma$ can be verified observationally in FSRQ only after the launch of the GLAST satellite in 2006.

Fortunately, the basic features of the current model of FSRQ as adopted here are not very sensitive to the presence and details of steep high energy tails of the injection function, and, in order to reproduce the spectra at the level of detail required by the current data, it is sufficient to use a single power-law injection function as a first approximation.

  
Appendix B: Broad emission line region

Reverberation mappings indicate that the BEL region in AGNs is geometrically thick, with $r_{\rm max}/r_{\rm min} > 30$(Peterson 1993). The 3D structure of the BEL region is still not fully understood, and in our code is assumed to be spherical. We approximate the radial distribution of BEL luminosity by two power-law functions which join at  $r_{\rm BEL}$ where the luminosity has a peak, i.e.

$${\partial L_{\rm BEL} \over \partial \ln r} (r) = {\partial... ...cr (r/r_{\rm BEL})^{-q_2} & for $r > r_{\rm BEL}$\space \cr} ~$$ (B.1)

where from $L_{\rm BEL} = \int (\partial L_{\rm BEL}/\partial r){\rm d}r$

$${\partial L_{\rm BEL} \over \partial \ln r} (r_{\rm BEL}) = {q_1 q_2 \over q_1+q_2} L_{\rm BEL} ~ ,$$ (B.2)

and q₁>0 and q₂>0. In the case of a smooth transition from the BEL region to the narrow emission lines region, one can expect $q_2 \sim 0.5$. The value of q₁ is very uncertain and, in order to minimize number of parameters, we adopt q₁=q₂ = 0.5.

References

 
• Bicknell, G., & Wagner, S. 2002, PASA, 19, 129 NASA ADS  
• Blackman, E. G. 1996, ApJ, 456, L87 NASA ADS  
• Bazejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107 NASA ADS  
• Böttcher, M., & Chiang, J. 2002, ApJ, 581, 127 NASA ADS  
• Chang, J. S., & Cooper, G. 1970, J. Chem. Phys., 6, 1  
• Celotti, A., Padovani, P., & Ghisellini, G. 1997, MNRAS, 286, 415 NASA ADS  
• Chiaberge, M., & Ghisellini, G. 1999, MNRAS, 306, 551 NASA ADS  
• Chiang, J., & Dermer, C. D. 1999, ApJ, 512, 699 NASA ADS  
• Crusius, A., & Schlickeiser, R. 1986, ApJ, 164, L16  
• Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275 NASA ADS  
• Dermer, C. D. 1995, ApJ, 446, 63  
• Dermer, C. D., & Atoyan, A. M. 2002, ApJ, 568, L81 NASA ADS  
• Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458 NASA ADS  
• Dondi, L., & Ghisellini, G. 1995, MNRAS, 273, 583 NASA ADS  
• Ghisellini, G., & Madeau, P. 1996, MNRAS, 280, 67 NASA ADS  
• Ghisellini, G., & Maraschi, L. 1989, ApJ, 340, 181 NASA ADS  
• Hartman, R. C., Böttcher, M., Aldering, G., et al. 2001, ApJ, 553, 683 NASA ADS  
• Homan, D. C., Attridge, J. M., & Wardle, J. F. C. 2001, ApJ, 556, 113 NASA ADS  
• Hoshino, A., Arons, J., Gallant, Y. A., & Langdon, A. B. 1992, ApJ, 390, 454 NASA ADS  
• Hoshino, A., & Shimada, N. 2002, ApJ, 572, 880 NASA ADS  
• Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533, 631 NASA ADS  
• Kino, M., Takahara, F., & Kusunose, M. 2002, ApJ, 564, 97 NASA ADS  
• Kirk, J. G., Melrose, D. B., & Priest, E. R. 1994, Plasma Astrophysics, ed. A. O. Benz, & T. J.-L. Couvoisier (Berlin: Springer)  
• Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452 NASA ADS  
• Koratkar, A., Pian, E., Urry, C. M., & Pesce, J. E. 1998, ApJ, 492, 173 NASA ADS  
• Königl, A. 1981, ApJ, 243, 700 NASA ADS  
• Larrabee, D. A., Lovelace, R. V. E., & Romanova, M. M. 2003, ApJ, 586, L72  
• Lawson, A. J., McHardy, I. M., & Marscher, A. P. 1999, MNRAS, 306, 247 NASA ADS  
• Levinson, A. 1996, MNRAS, 278, 1018 NASA ADS  
• Lobanov, A. P. 1998, A&A, 330, 79 NASA ADS  
• Lyutikov, M. 2002 [astro-ph/0211338]  
• Macomb, D. J., Akerlof, C. W., Aller, H. D., et al. 1995, ApJ, 449, L99 NASA ADS  
• Madejski, G. M., Sikora, M., Jaffe, T., Bazejowski, M., & Moderski, R. 1999, ApJ, 521, 145 NASA ADS  
• Magdziarz, P., Moderski, R., & Madejski, G. M. 1997, in Relativistic Jets in AGNs, ed. M. Ostrowski, M. Sikora, G. M. Madejski, & M. C. Begelman (Kraków: PITS), 238  
• Mannheim, K., & Biermann, P. L. 1992, A&A, 253, L21 NASA ADS  
• Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1993, ApJ, 410, 609 NASA ADS  
• McClements, K. G., Dendy, R. O., Bingham, R., Kirk, J. G., & Drury, L. O'C. 1997, MNRAS, 291, 241 NASA ADS  
• Moderski, R., Sikora, M., & Bulik, T. 2000, ApJ, 529, 151 NASA ADS  
• Mukherjee, R., Bertsch, D. L., Bloom, S. D., et al. 1997, ApJ, 490, 116 NASA ADS  
• Peterson, B. M., 1993, PASP, 105, 247 NASA ADS  
• Pian, E., Urry, C. M., Maraschi, L., et al. 1999, ApJ, 521, 112 NASA ADS  
• Pohl, M., Hartman, R. C., Jones, B. B., & Sreekumar, P. 1997, A&A, 326, 51 NASA ADS  
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C (Cambridge University Press)  
• Reynolds, S. P. 1982, ApJ, 256, 38 NASA ADS  
• Romanova, M. M., & Lovelace, R. V. E. 1992, A&A, 262, 26 NASA ADS  
• Ruszkowski, M., & Begelman, M. C. 2002, ApJ, 573, 485 NASA ADS  
• Shimada, N., & Hoshino, M. 2000, ApJ, 543, L67 NASA ADS  
• Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 NASA ADS  
• Sikora, M., Bazejowski, M., Begelman, M. C., & Moderski, R. 2001, ApJ, 554, 1 NASA ADS  
• Sikora, M., Bazejowski, M., Moderski, R., & Madejski, G. M. 2002, ApJ, 577, 78 NASA ADS  
• Sikora, M., & Madejski, G. M. 2001, in High Energy Gamma-Ray Astronomy, ed. F. A. Aharonian, & H. J. Völk, AIP Conf. Proc., 558, 275 [astro-ph/0101382]  
• Sikora, M., & Madejski, G. M. 2002, in Current High Energy Emission around Black Holes, ed. C.-H. Lee, & H.-Y. Chang (New Jersey, London, Singapore, Hong Kong: World Scientific), 132 [astro-ph/0112231]  
• Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559 NASA ADS  
• Swanenburg, B. N., Hermsen, W., Bennett, K., et al. 1978, Nature, 275, 298 NASA ADS  
• Takahashi, T., Tashiro, M., Madejski, G., et al. 1996, ApJ, 470, L89 NASA ADS  
• Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 NASA ADS  
• Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445 NASA ADS  
• von Montigny, C., Bertsch, D. L., Chiang, J., et al. 1995, ApJ, 440, 525 NASA ADS  
• Wagner, S. J., Mattox, J. R., Hopp, U., et al. 1995, ApJ, 454, L97 NASA ADS  
• Wardle, J. F. C., Homan, D. C., Ojha, R., & Roberts, D. H. 1998, Nature, 395, 457 NASA ADS  
• Wehrle, A. E., Pian, E., Urry, C. M., et al. 1998, ApJ, 497, 178 NASA ADS  
• Zdziarski, A. A., & Krolik, J. H. 1993, ApJ, 409, L33 NASA ADS