A&A 410, 101-115 (2003)
DOI: 10.1051/0004-6361:20031245

The multifrequency variability of Mrk 421

K. Katarzynski 1,2 - H. Sol2 - A. Kus 1


1 - Torun Centre for Astronomy, Nicolaus Copernicus University, ul. Gagarina 11, 87100 Torun, Poland
2 - LUTH, Observatoire de Paris-Meudon, 5 place J. Janssen, 92195 Meudon Cedex, France

Received 10 July 2002 / Accepted 25 July 2003

Abstract
We present a model of the multifrequency variability of the blazar Mrk 421. The model explains correlated variability observed from Very High Energy (VHE) gamma rays to radio frequencies. We assume that the dominant part of the stationary emission from the radio frequencies to the X-rays is generated by the synchrotron radiation of relativistic electrons ejected from the central engine. The particles move from the center of the source with relativistic velocities and form an inhomogeneous jet. We perform detailed calculations of the radiation transfer and calculate evolution of the electron energy spectrum along the jet. We explain the observed variability by the evolving synchrotron and Inverse-Compton (IC) radiation of a compact component (a blob) which travels along the jet. Two scenarios have been considered as mechanisms to generate VHE flares. The first scenario assumes that the high energy electrons, necessary for generation of the VHE flares, are injected into the jet, directly from the central engine or from an acceleration zone (e.g., a shock wave). The second scheme assumes that the high energy electrons are generated in situ by acceleration, for example by diffusive shock waves or a localized turbulence inside the jet. The particles evolve along the jet. They are cooled by the radiative processes and by the adiabatic expansion which compete with the acceleration process and the injection of high energy electrons. We present new observations we obtained in the radio domain for Mrk 421. The radio data gathered in February-April 2001 show a well defined radio outburst which corresponds to an X-ray outburst observed by RXTE-ASM and a gamma-ray flare detected by HEGRA in the TeV range. The best of our knowledge, this is the first direct observational evidence for a flare observed simultaneously in the radio range and at very high energies. Our scenario with acceleration of electrons in the middle part of the jet describes well the temporal evolution of such multispectral flare.

Key words: radiation mechanisms: non-thermal - galaxies: active - BL Lacertae objects: individual: Mrk 421

1 Introduction

Blazars are very well known for their strong variability, present from the radio frequencies across the optical wavelengths to the X-rays and in some particular cases up to VHE gamma rays. The observed variability time scales seem to depend mainly on the frequency domain. The fastest changes, of about hours or even several minutes, are observed in the X-ray and the gamma-ray ranges (e.g., Catanese et al. 1997; Pian at al. 1998; Djannati-Atai et al. 1999; Maraschi et al. 1999). Temporal scales from a few days to a few weeks or even years are observed from the radio (e.g., Hufnagel & Bregman 1992; Romero et al. 1997; Aller et al. 1999) to the optical frequencies (e.g., Tosti et al. 1998; Xie et al. 1999).

However, some fraction of the sources exhibit Intra Day Variability (IDV) even in the radio and optical ranges, with observed time scales of about one day or less (e.g., Rickett et al. 1995; Wagner et al. 1996). To explain such rapid variations it is necessary to assume very compact sources which travel towards the observer with ultra relativistic velocities (e.g., Blandford & Konigl 1979). Thus, in many cases astronomers try to explain IDV by extrinsic phenomena. There are two possible processes which can produce rapid variability of an originally quiet source. The first process is the microlensing proposed by Chang & Refsdal in 1979. The process is achromatic and should be able to produce coherent variability at any wavelength. However, detailed analyses found several arguments against such phenomena as possible explanation of IDV (see for example Wagner 1992). The second possibility is related to inhomogeneities in the electron density of the interstellar medium (ISM) on different spatial scales. Apparent variability (interstellar scintillation, ISS) is produced by diffractive (DISS) and refractive (RISS) scattering of electromagnetic waves (e.g., Rickett 1990). It is quite difficult to separate intrinsic variability from variability generated by ISS. Complex campaigns of observations performed at different bands are necessary, especially to identify RISS which may produce longer variability time scales (of about a few days or weeks, see e.g., Fiedler et al. 1987). However, correlation between variability observed at radio and optical wavelengths excludes RISS as possible source of variability. Furthermore, RISS cannot generate variability in the optical range. In 1997, Tosti et al. (1998) reported such a correlation in Mrk 421. Another argument against extrinsic origin of the variability may come from observations at frequencies higher than optical. Variability observed in X-rays or gamma rays strongly suggests that all observed activity is generated inside a source. Therefore, in our model we may safely assume that the observed variability is intrinsic and generated by the source components.

There are two basic types of models which can explain intrinsic variability. The first assumes that the observed variations are related to changes in the geometry of emitting sources (e.g., Camenzind & Krockenberger 1992; Gopal-Krishna & Wiita 1992). The most common scenario explains variability of the source as being due to blobs injected at the base of the jet. Such components move with finite angular momentum on a helical trajectory. Variability is generated mostly by the change of beaming conditions. The scenario was successfully applied in some cases (e.g., Schramm et al. 1993; Wagner et al. 1995). The second kind of models assume that variability is generated by change of emission condition inside the jet. Typical examples are the injection of fresh particles from the center of a source, or the acceleration of particles by a shock wave (e.g., Blandford & Konigl 1979; Marscher & Gear 1985; Celotti et al. 1991; Kirk et al. 1998). Depending on the physical conditions, the variability can originate at high energy bands (X-rays & gamma rays) and then spread to the radio wavelengths. However, alternative cases are also observed with activity only at the low energies (optical to radio wavelengths) or only at the high energies.

In order to further analyse variability processes in blazars and to provide tools necessary to interpret multifrequency monitoring we present here plausible scenarios for VHE flares, correlated with variability at lower frequencies. We base our models of temporal evolution on our previous study (Katarzynski et al. 2001) performed for stationary states, assuming pure leptonic modeling. This type of scenario is known to provide a good framework to explain the high energy activity of blazars (e.g., Ghisellini & Madau 1996; Dermer & Schlickeiser 1993; Sikora et al. 1994; Inoue & Takahara 1996).

2 Modeling the quiet state and the low energy emission

We assume that at least three components contribute to the stationary emission of Mrk 421 from the radio wavelengths to the optical range. The dominant part of the radiation is generated by the inner part of the jet which is seen on the radio maps as a unresolved core (e.g., Piner et al. 1999). The component radiates from the high radio frequencies ( $\nu \gtrsim 1$ GHz) to the optical or even ultraviolet ranges. Hereafter, for sake of simplicity, we call this component just the "jet''. An additional significant part of the low frequency ( $\nu \lesssim 1$ GHz) radio emission is generated in extended radio structures. These structures appear on the radio maps as several radio clouds. A third component, the thermal radiation of the host galaxy, dominates the spectrum in the optical range. However, the "jet'' mentioned above also radiates significantly at such wavelengths. It generates the optical nonthermal continuum of the radiation observed in the nucleus of the source (Ulrich et al. 1975).

We developed a model for the inner jet very similar in fundamental assumptions to the widely known jet models (e.g., Marscher 1980; Konigl 1981; Ghisellini et al. 1985; O'Dell 1988). The difference between our approach and the previous ones is that we improve calculations of the radiation transfer along the jet and we calculate evolution of the electron energy spectrum.

For sake of simplicity, we assume a constant value for the velocity of jet components ( $V_{\rm jet}$) and a constant angle to the line of sight ($\theta$). This also gives constant value of the Doppler factor along the jet:

\begin{displaymath}\delta_{\rm jet} = \frac{1}{\Gamma_{\rm jet} (1 - (V_{\rm jet}/c) \cos \theta)}
\end{displaymath} (1)

where $\Gamma_{\rm jet} = 1 / \sqrt{1 - (V_{\rm jet} / c)^2)}$ is the Lorentz factor of the jet and c is the speed of light. Thus it differs from reality, since the geometry and kinematics of the jet are more complex. A more realistic scenario may assume that $V_{\rm jet}$ decreases along the jet. However, we can also imagine a situation where $\theta$ decreases along the jet (e.g., an extended helical trajectory of the jet's components) what could stabilize the value of $\delta_{\rm jet}$.

We neglect the total inverse-Compton radiation generated within the whole structure of the jet. This type of radiation was calculated for Mrk 501 in our recent work (Katarzynski et al. 2001). These calculations show that this radiation should be observable at energies from MeV to GeV. Therefore, it cannot contribute to the emission observed in VHE gamma rays, which is analyzed in the case of this paper. Moreover, the variability time scales observed in VHE gamma rays indicate that the emission region must be at least a hundred times smaller than the total length of the jet. On the other hand, observations made by EGRET from 0.1 to 5 GeV show that radiation of Mrk 421 is relatively weak at this range ( $2.67 \times 10^{-37} \pm 4.08 \times 10^{-38}$ W m-2 Hz-1, Thompson et al. 1995). The level of the emission at this range is also variable (e.g., Macomb et al. 1995). The data on Mrk 421 as well as on other blazars at the MeV-GeV range are still quite poor. It is quite difficult to distinguish which part of the radiation comes from the regular structure of the jet and which part is generated by compact components inside the jet. Therefore, we decided to explain the whole gamma-ray emission of Mrk 421 by radiation of compact components or "blobs'' inside the jet.

In comparison with the other models, we assume that the particle density along the jet decreases relatively fast (see Sects. 2.1 and 2.2). Taking into account that the magnetic field intensity decreases as a power law function along the jets, we obtain jets where the inner parts (nearest to the black hole) are significantly more luminous than the outer parts. In addition, the relatively large density of electrons necessary to fit the observed level of radiation induces a fast decrease of the emission below GHz frequencies due to the electron's self-absorption. To explain the radio emission up to the GHz bands we have to introduce a term which describes the radiation from the extended radio structure. This additional component is assumed to be homogeneous and spherical.

The thermal emission of the host galaxy is simply reproduced by a black body spectrum. It provides an accuracy high enough for the present work and introduce only two free parameters (temperature and level of emission) to describe this radiation.

   
2.1 Geometry of the jet

In the first approach of our model we selected a conical or paraboloidal geometry of the jet. However, to simplify calculation of the radiation transfer along the jet, we approximate the geometry by a system of homogeneous cylindrical elements - slices. The slices travel from the central engine with constant velocity ( $V_{\rm jet}$) and expand adiabatically, following a conical or paraboloid geometry. The length ( $L_{\rm jet}$ (cm)) and the radius ( $R_{\rm jet}$ (cm)) of each element increase with time according to the following formula:

\begin{displaymath}L_{\rm jet} \left( t^i_{\rm jet} \right) = L^i_{\rm jet} = V_{\rm exp} t^i_{\rm jet} ,
\end{displaymath} (2)


\begin{displaymath}R_{\rm jet} \left( t^i_{\rm jet} \right) = R^i_{\rm jet} = R^...
...ft( \frac{t^i_{\rm jet}}{t^1_{\rm jet}} \right)^{r_{\rm jet}},
\end{displaymath} (3)

where $V_{\rm exp}$ (cm s-1) is the expansion velocity of the cylinder length, $R^1_{\rm jet}$ is the radius of the first cylinder, $r_{\rm jet}$ describes the jet geometry and $t^i_{\rm jet}$ are the time values defined by:

\begin{displaymath}t^1_{\rm jet} = L^1_{\rm jet} / V_{\rm exp},
\end{displaymath} (4)


\begin{displaymath}t^i_{\rm jet} = \frac{L^{i - 1}_{\rm jet} + 2 V_{\rm jet} t^{...
... jet} }{ 2 V_{\rm jet} - V_{\rm exp} },
~~~{\rm for}~~~i > 1.
\end{displaymath} (5)

This parameterization provides a correlation between the expansion rate of the slices and their velocity in space. It means that the shift of a slice in space equals half of the length increment. Therefore, a quasi-conical ( $r_{\rm jet}=1$) or paraboloidal ( $r_{\rm jet}<1$) geometry of the jet can be built by a system of moving and expanding cylindrical elements. The time moments ( $t^i_{\rm jet}$) are set to describe physical parameters calculated at the central point of each cylinder.

We assume that the shape of the jet is maintained by a continuous outflow of the slices from the central engine. The particles which reach the most distant part of the source are accumulated in the extended radio structures. The slices which form the shape are indicated by i=1...nj where nj is the total number of slices which form the jet shape.

For further simplification we assume that each slice is homogeneous with a uniform intensity of the magnetic field. Decrease of the magnetic field intensity [G] along the jet is defined by the power law function which we parameterize by:

\begin{displaymath}B(t^i_{\rm jet}) = B^i_{\rm jet} = B^1_{\rm jet} \left(\frac{t^1_{\rm jet}}{t^i_{\rm jet}}\right)^{m_{\rm jet}},
\end{displaymath} (6)

where $B^1_{\rm jet}$ is an initial magnetic field intensity. In our model we assume a relatively fast decrease of the magnetic field intensity along the jet. However, the decrease is slower than derived under the assumption of magnetic flux conservation ( $m_{\rm jet} = 2 r_{\rm jet}$). This may be explained by a turbulent structure of the magnetic field.

   
2.2 Evolution of the electron energy distribution along the jet

We describe the evolution of the electron energy spectrum $N^{\ast}_{\rm jet}$, which refers to the entire source, by the kinetic equation (e.g., Kardashev 1962):


 
$\displaystyle \frac{\partial N^{\ast}_{\rm jet} ( \gamma , t )} {\partial t}$ = $\displaystyle \frac{\partial}{\partial \gamma}~\Big\{
\Big[C^{\rm cool}_{\rm jet} (t) \gamma^2 - ( C^{\rm acc}_{\rm jet} (t)$ - $\displaystyle C^{\rm adia}_{\rm jet} (t) ) \gamma \Big]~N^{\ast}_{\rm jet} ( \gamma , t ) \Big\},$ (7)

where $C^{\rm cool}_{\rm jet}$ characterizes electron energy losses due to synchrotron and IC radiation, $C^{\rm adia}_{\rm jet}$ describes adiabatic losses due to expansion and $C^{\rm acc}_{\rm jet}$ characterizes acquisition of the electron energy in some acceleration processes (e.g., acceleration by a shock wave). The equation does not consider injection of the electrons into the emitting region nor the escape. However, in further calculations we simulate these two processes using our specific parameterization of the source geometry. The initial distribution of electron energy is given by the power law function:

 \begin{displaymath}N^{\ast}_{\rm jet}(\gamma, t^1_{\rm jet}) = K^1_{\rm jet} \le...
...ft(\gamma^{\rm cut}_{\rm jet}\right)^{-n_{\rm jet}-2} \right],
\end{displaymath} (8)

with a sharp cutoff at the highest energies ( $\gamma^{\rm cut}_{\rm jet}$). In the above formula, $ K^1_{\rm jet} $ is the initial number of electrons per cubic centimeter, $n_{\rm jet}$ describes the slope of the power law function and $\gamma^{\rm cut}_{\rm jet}$ indicates position of the cut-off in the energy spectrum. We succeeded to find the general solution of the Eq. (7), which is given in Appendix A. The methods to solve this type of equations can be found in mathematical handbooks, however, the methods described by Kardashev (1962) and by Kirk et al. (1998) are very helpful for such calculations.

To apply the solution for calculation of emission and absorption coefficients we have to convert the equations to unit volume. In the case where the evolution of a slice of volume ($V_{\rm s}$) is described by a power law function $( V_{\rm s}(t) \sim t^{2r_{\rm jet}+1}_{\rm jet})$ the conversion is given by:

\begin{displaymath}N_{\rm jet}(\gamma, t) = N^{\ast}_{\rm jet}(\gamma, t) \left( \frac{t^1_{\rm jet}}{t} \right)^{-3 d_{\rm jet}},
\end{displaymath} (9)

where $d_{\rm jet} = t \nabla {\vec v}(t) /3$ and ${\vec v}(t)$ is the velocity field, which determines expansion of $V_{\rm s}$ (e.g., Longair 1981). We assume that the total number of particles is constant during expansion of the jet slice, therefore $d_{\rm jet}=(2r_{\rm jet}+1)/3$. Note that the conversion can be replaced by an additional term $-N(\gamma, t) \nabla {\vec v}(t)$ on the left side in Eq. (7). The term for our specific parameterization of the geometry is given simply by formula: $-N(\gamma, t) 3 d_{\rm jet} / t$. However, this manner of calculating greatly complicates the procedure for analytically solving Eq. (7).

To calculate the electron energy evolution we need to specify cooling and acceleration processes and to describe decrease of the electron energy due to adiabatic expansion of the slices.

In our model of the jet radiation we selected a specific set of physical parameters (see Table 2) which implies a relatively weak radiation field energy density inside the jet. Therefore, we can consider cooling by the synchrotron process only, described by characteristic cooling time scale:

\begin{displaymath}t^{\rm cool}_{\rm jet}(\gamma,t) = \frac{1}{\gamma~C^{\rm coo...
...et} (t) = \frac{4 \sigma_{\rm T} U_{\rm B}(t)}{3 m_{\rm e} c},
\end{displaymath} (10)

where $U_{\rm B} = B^2_{\rm jet}(t) / (8 \pi)$ is the magnetic field energy density.

We consider acceleration processes which have an acceleration rate which is energy independent and decreasing along the jet, described by a characteristic acceleration time:

\begin{displaymath}t^{\rm acc}_{\rm jet}(t) = \frac{1}{C^{\rm acc}_{\rm jet} (t)...
..._{\rm jet} \left(\frac{t^1_{\rm jet}}{t}\right)^{a_{\rm jet}},
\end{displaymath} (11)

with $A^1_{\rm jet} = 1/t^1_{\rm acc_{\rm jet}}$ where $t^1_{\rm acc_{\rm jet}}$ is initial characteristic acceleration time. Note that we do not specify which acceleration scenario among the variety of models discussed in the literature (see for a short review Longair 1981) could provide the best physical conditions required by the model. Detailed description of the acceleration would introduce here too many additional parameters which would be difficult to constrain from observational data. Thus we have selected a very simple, although somewhat arbitrary, description of the acceleration with a minimal number of free parameters.

Taking into account our parameterization of the geometry evolution, we can write the coefficient which characterizes losses of the electron energy due to the adiabatic expansion by:

\begin{displaymath}C^{\rm adia}_{\rm jet} (t) = \frac{1}{3} \nabla {\vec v}(t) = \frac{d_{\rm jet}}{t} \cdot
\end{displaymath} (12)

The solution (7) provides a function which is continuous in time. However, according to our assumption on homogeneity, we acctually use a stepped particle energy energy spectrum, uniform within each slice. The spectrum is approximated by the distribution calculated at the central point of the slice: $N^i_{\rm jet} = N_{\rm jet}(\gamma, t^i_{\rm jet})$ (cm-3).

2.3 Radiation transfer along the jet

In order to calculate the radiation transfer along the jet, we consider a case where the angle between the jet symmetry axis and the line of sight ($\theta$) is relatively small (about a few degrees). Therefore, we assume that the solution of radiation transfer along paths parallel to the axis is valid for all cases investigated in this paper.

Assuming that each slice is homogeneous with constant emission ( $j_{\rm jet}$ (erg s-1 Hz-1 sterad-1)) and absorption ( $k_{\rm jet}$ (cm-1)) coefficients, we can write the intensity (erg s-1 cm-2 Hz-1 sterad-1) of the radiation for a given element:

 \begin{displaymath}I^i_{\rm jet}(\nu') = \frac{j^i_{\rm jet}(\nu')} {k^i_{\rm je...
...eft[ 1 - \exp \left( - \tau^i_{\rm jet}(\nu') \right) \right],
\end{displaymath} (13)

where $\tau^i_{\rm jet}(\nu') = L^i_{\rm jet} k^i_{\rm jet}(\nu')$. In this particular case, $j^i_{\rm jet}(\nu')$ is the synchrotron emission coefficient, $k^i_{\rm jet}(\nu')$ is the electron self-absorption coefficient (see for example Tavecchio et al. 1998; Katarzynski et al. 2001) and $\nu'$ is the frequency (Hz) measured in the source frame. Note that hereafter, for sake of readability, we omit to write the function's frequency dependency in all formulae for the radiation intensity and the flux density.

Using a transformation to the observer frame, we can express observed flux density (erg cm-2) by:

 \begin{displaymath}F_{\rm jet} = \frac{\pi \delta_{\rm jet}^3 (1+z)}{d_{\rm l}^2...
... - \left( R^{i-1}_{\rm jet} \right)^2 \right] I^{i}_{\rm tot},
\end{displaymath} (14)

where $I_{\rm tot}$ is the total intensity of the radiation calculated along the jet, which can also include the radiation of the compact component (see Appendix B). However, in this part, we describe only radiation of the jet and use a simple version of the formula for $I_{\rm tot}$ without any compact component. In the above formula $d_{\rm l}$ is the luminosity distance (e.g., Weinberg 1982), z is the redshift, $\delta_{\rm jet}$ is the Doppler factor of the jet and $R^0_{\rm jet} = 0$. We use a Hubble constant of 65 km s-1 Mps-1and a deceleration parameter q0 = 0.5 for the calculations of flux densities.

2.4 Extended radio and thermal emissions

One motivation of our present analysis is to describe simultaneously all parts of the blazar spectrum. We want to look for correlation at different wavelengths and try to interpret multifrequency light curves from the radio bands to the gamma rays. According to the assumptions which we selected for the jet model, radiation of the jet is dominated by the initial part of the jet. The density of electrons in the initial part is relatively high. Therefore the model cannot explain emission observed at lower radio frequencies (below GHz frequencies) because of electron self-absorption. The problem of self-absorption does not appear in jet models where the radio emission is dominated by the outer part of a jet that is not absorbed. To overcome this problem in our model we assume that the low frequency radio emission is dominated by extended radio structures.

The structures can be described as clouds or lobes at the end of the jet, fed by the jet itself. For simplicity, we approximate the extended structures by a spherical, homogeneous cloud, which generates the synchrotron radiation observable in low frequency radio waves ( $\nu \lesssim 1$ GHz). The observed flux density for such a source is calculated according to:

\begin{displaymath}F_{\rm ext} = \frac{\pi \delta_{\rm ext}^3 (1+z)}{d_{\rm l}^2...
...m ext}} {\rm e}^{-\tau_{\rm ext}} [\tau_{\rm ext} + 1] \right)
\end{displaymath} (15)

where $\tau_{\rm ext} = 2 R_{\rm ext} k_{\rm ext}$, $j_{\rm ext}$ and $k_{\rm ext}$ are the synchrotron emission and the electron self-absorption coefficients respectively, $\delta_{\rm ext}$ is the Doppler factor of the source. The electron energy distribution for such a simple cloud is described by the spectrum calculated for the last jet slice. However, the particle density inside the cloud is much higher due to the accumulation of particles during the long term evolution of the source. Therefore, we define the particle energy distribution by:

 \begin{displaymath}N_{\rm ext}(\gamma) = c_{\rm ext} N_{\rm jet}(\gamma, t^{n_j})
\end{displaymath} (16)

where $c_{\rm ext}$ is a free parameter, which describes the difference in the particles' density. We assume also that the magnetic field intensity inside the cloud is equal to the magnetic field intensity inside the last jet slice ( $B_{\rm ext} = B^{n_j}_{\rm jet}$). The real situation probably involves more complex geometry and distribution of particles. However, this simple approach explains reasonably well the emission of the extended radio structure and allows a correct level of permanent low frequency radio radiation to be attained.

The thermal radiation of the host galaxy, quite clearly visible at the infrared and optical wavelengths (see Fig. 3), is simply modeled by the black body radiation:

\begin{displaymath}F_{\rm th} = F^{\rm max}_{\rm th} P(\nu, T_{\rm th}) / P(3 k_{\rm B} T_{\rm th} / h, T_{\rm th})
\end{displaymath} (17)

where $P(\nu, T)$ is the Planck's function, $F^{\rm max}_{\rm th}$ is the observed flux density for the maximum of thermal radiation, $T_{\rm th}$ is the temperature, $k_{\rm B}$ is Boltzmann's constant and h is Planck's constant. This basic approach provides a good fit of the thermal emission with only two free parameters ( $F^{\rm max}_{\rm th}, T_{\rm th}$). We introduce it to explain the observed level of radiation at optical wavelengths and to generate the optical light curves.

3 Modeling of variability

In this part, we go one step further in complexity to explain multifrequency variability of Mrk 421. The approach extends our previous blob-in-jet scenario for blazar emission (Katarzynski et al. 2001). The main assumption is that the whole observed activity is generated by a single compact component (hereafter called a "blob'') placed inside the jet. We investigate two general cases where the blob may be created by rapid injection of high energy electrons or by systematic acceleration of the electrons. Apart from the two main mechanisms, we allow also some weak compression of the local jet medium. This results in an increase of the particle density and of the magnetic field intensity. The main difference between the blob and jet medium, which surrounds the blob, appears in the maximal energy and the density of the electrons. To generate activity it is necessary to provide a density inside the blob at least one order of magnitude higher than density inside the jet medium ( $K_{\rm blob} \gg K_{\rm jet}$). To explain correctly the spectral shapes observed in X-rays and VHE gamma rays it is necessary to provide a maximal energy of the electrons in the blob two orders of magnitude higher than inside the jet medium ( $\gamma^{\rm cut}_{\rm blob} \gg \gamma^{\rm cut}_{\rm jet}$).

   
3.1 Evolution of blob geometry and observed light curves

We describe the blob by several slices, as we did for the jet. The length of blob slices evolves exactly like the length of jet slices ( $L_{\rm blob} (t^i_{\rm jet}) = L_{\rm jet} (t^i_{\rm jet})$). The radius of the blob slices is less or equal to the radius of the jet, however it evolves proportionally to the radius of jet slices:

\begin{displaymath}R_{\rm blob} ( t^i_{\rm jet} ) = R^i_{\rm blob} = R^1_{\rm bl...
...ft( \frac{t^i_{\rm jet}}{t^1_{\rm jet}} \right)^{r_{\rm jet}},
\end{displaymath} (18)

where $R^1_{\rm blob}$ is the initial radius of the blob's slice. The blob components are placed at the center of the jet slices. The velocities of jet and blob slices are equal. Therefore there is no shift between the elements. In Fig. 1 we show the evolution of the blob geometry and the parameters which we use to describe this evolution.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f1.eps}}\end{figure} Figure 1: Geometry of the jet and evolution of the blob geometry. The current position of the first blob slice, related to the shape of the jet, is given by a = 1...nj. The position of the last blob slice is given by $b = a - n_{\rm b} + 1$, where $n_{\rm b}$ is the number of blob slices. The diagram illustrates also the scheme applied for calculation of the radiation transfer. In the framework of this scheme, the structure of the jet has been divided, parallel to its symmetry axis, into three specific regions. Radiation which comes from the first region ( $F^a_{\rm in}$) is generated by the blob and by the part of the jet smaller than $R^a_{\rm blob}$. Emission which comes from two other regions ( $F^a_{\rm mid}, F^a_{\rm out}$) is generated only by parts of the jet which are larger than $R^a_{\rm blob}$. A detailed description of this scheme is included in Sect. 3.3.
Open with DEXTER

To simplify the generation of the light curves we calculate the emission in steps. The blob slices are shifted towards the end of the jet relatively to the jet shape by one slice length for every step. The time intervals for every step are defined by $\Delta t^i_{\rm jet} = t^i_{\rm jet} - t^{1}_{\rm jet}$ in the source frame and are transformed into the observer's frame by $\Delta t_{\rm obs} = \Delta t_{\rm jet} /
\delta_{\rm jet}$. The time interval grows nonlinearly along the jet. Therefore, for some typical set of parameters ( $L^1_{\rm jet} \sim 10^{15}$ cm, $R^1_{\rm jet}
\sim 10^{16}$ cm, $\delta_{\rm jet} \sim 10$) we can obtain time intervals of the order of hours for the jet base and of the order of days for middle and outer parts of the jet. This provides a good temporal coverage for the light curves generated from the VHE range to the radio frequencies.

3.2 Electron energy distribution inside the blob

The evolution of blob electrons is calculated in the same way as the evolution of electrons in the jet slices. We use our general solution (see Appendix A) with almost the same assumptions:

\begin{displaymath}N^{\ast}_{\rm blob}(\gamma, t^1_{\rm jet}) = K^1_{\rm blob} \...
...(\gamma^{\rm cut}_{\rm blob}\right)^{-n_{\rm blob}-2} \right],
\end{displaymath} (19)


\begin{displaymath}t^{\rm acc}_{\rm blob}(t) = \frac{1}{C^{\rm acc}_{\rm blob} (...
...\rm blob} \left(\frac{t^1_{\rm jet}}{t}\right)^{a_{\rm blob}},
\end{displaymath} (20)


\begin{displaymath}C^{\rm adia}_{\rm blob} = \frac{d_{\rm blob}}{t},
\end{displaymath} (21)

like in the case of jet slices. The main difference appears in the cooling conditions, where we have to consider the IC cooling in addition. This can be done simply by adding the synchrotron field energy density ($U_{\rm r}$) into the cooling coefficient:

\begin{displaymath}t^{\rm cool}_{\rm blob}(\gamma,t) = \frac{1}{\gamma~C^{\rm co...
...{4 \sigma_{\rm T} (U_{\rm B} + U_{\rm r})}{3 m_{\rm e} c}\cdot
\end{displaymath} (22)

However, the radiation field energy density must be calculated under the assumption defined by the Klein-Nishina decline (see e.g., Tavecchio et al. 1998). This introduces an additional energy dependency for the physical parameter $U_{\rm r} = U_{\rm r}(t, \gamma$) which greatly complicates the solution of Eq. (7). To avoid this problem we parameterize $U_{\rm r}$ by $U_{\rm B}$:

 \begin{displaymath}U_{\rm r} = U_{\rm B} / \eta,
\end{displaymath} (23)

where the parameter $\eta$ is (hopefully) deduced from observations (see Sect. 4.2).

In addition to the blob synchrotron radiation field, there are two other radiation fields inside the blob, which may play a role in the cooling of blob electrons by IC scattering. However, we do not consider them, as justified below.

The first additional radiation field is the synchrotron radiation of the structure of the jet which surrounds the blob. The intensity of this radiation can be estimated according to the following formula: $I^{\rm rad}_{\rm jet} \sim j_{\rm jet} (R_{\rm jet} - R_{\rm blob})$ (e.g., Ghisellini et al. 1985). For a power law distribution of the electron energy spectrum ( $N(\gamma) \sim \gamma^{-n}$), the synchrotron emission coefficient can be approximated by $j \sim
K B^{1+(n-1)/2} (\nu')^{-(n-1)/2}$. Therefore, assuming $B_{\rm blob}
\gtrsim B_{\rm jet},~R_{\rm blob} \lesssim R_{\rm jet}$ and $K_{\rm blob} \gg K_{\rm jet}$, we can say that $I^{\rm rad}_{\rm jet}$ is negligible in comparison to the radiation field generated by the blob itself.

The second radiation field, which may be scattered by the electrons inside the blob, is the thermal radiation field of the accretion disk and of the medium which surrounds the active center. It is quite difficult to estimate the intensity of this radiation inside the blob. The intensity depends on the total luminosity of the nucleus of the source ( $L_{\rm nuc}$) and on the distance from the center of the source to the position of the blob ( $d_{\rm c \to b}$). To estimate the possible influence of this radiation field we use the method described by Inoue & Takahara (1996). In this method, the radiation field is approximated by the radiation of a blackbody at a given temperature ( $T_{\rm ext}$). For the estimation we assume $\tau L_{\rm nuc} = 10^{41}$ erg s-1 (a detailed description of this parameter is given in the work of Inoue & Takahara), and $T_{\rm ext} = 5~ \times~ 10^4$ K (to obtain a thermal radiation with maximum in the UV range). The calculation shows that the intensity of such a radiation field is negligible in the blob frame if $d_{\rm c \to b} \gtrsim 0.1$ pc. Therefore, keeping to a set of the parameters which reasonably describe the central engine of the source, it is possible to neglect this type of radiation field.

   
3.3 Radiation transfer including emission of the blob

To calculate the observed spectrum in the case where the blob is placed inside the jet, we have to consider absorption of the blob radiation by the jet itself. The synchrotron radiation of the blob is absorbed by the jet electrons. The VHE IC photons emitted by the blob may produce pairs by interaction with jet synchrotron photons ( $\epsilon_{\rm com} + \epsilon_{\rm syn} \to {\rm e}^+ + {\rm e}^-$). In addition to these two main processes, we can also expect absorption of the jet synchrotron radiation by electrons of the blob, especially in the case where $R^i_{\rm blob}$ is comparable with $R^i_{\rm jet}$. The calculation of the radiation transfer is more complex in this case and it becomes necessary to divide the jet into three different regions.

First we calculate the contribution to the observed flux density generated within a part of the source smaller or equal to  $R^a_{\rm blob}$:

\begin{displaymath}F^{a}_{\rm in} = \frac{\pi \delta_{\rm jet}^3 (1+z)}{d_{\rm l}^2} \sum^{a}_{i=1} f^i,
\end{displaymath} (24)


\begin{displaymath}f^x = \left\{\begin{array}{ll}
\left[(R^{x}_{\rm blob})^2 - (...
...}_{\rm tot}
\end{array},
& h^x = l^x + 1
\end{array}\right.,
\end{displaymath} (25)

where lx is the index number of the jet slice with maximal radius just smaller than $R^x_{\rm blob}$, and hx is the number of jet's slice with minimal radius just bigger than $R^x_{\rm blob}$. Note that for some cases $R^x_{\rm blob} < R^1_{\rm jet}$ therefore lx = hx = 1.

The second contribution to the total radiation comes from the medium part of the source larger than $R^a_{\rm blob}$ but smaller than $R^{h^x}_{\rm jet}$ given by:

\begin{displaymath}F^{a}_{\rm mid} = \frac{\pi \delta_{\rm jet}^3 (1+z)}{d_{\rm ...
...{\rm jet})^2 - (R^{a}_{\rm blob})^2 \right] I^{h^x}_{\rm tot}.
\end{displaymath} (26)

The third part of radiation ( $F^a_{\rm out}$) is generated by the outer part of the jet, which is bigger than $R^{h^x}_{\rm jet}$. This part of the radiation is calculated using a modified version of the Eq. (14), where the sum starts from hx + 1 instead of 1.

Note that for the middle and outer parts of the jet it is necessary to eliminate the radiation of the blob in calculations of $I_{\rm tot}$. These two parts generate radiation at a radius bigger than the current radius of the blob. Finally the total flux density is given by the sum of partial contributions:

\begin{displaymath}F^a_{\rm tot} = F^a_{\rm in} + F^a_{\rm mid} + F^a_{\rm out}.
\end{displaymath} (27)

All above formulae can be used directly for calculations of the total synchrotron spectrum, assuming that the synchrotron emission coefficient and the electron self-absorption coefficient for the blob are calculated in the same way as the jet coefficients. The magnetic field intensity inside the blob may be higher than inside the jet but it evolves in the same way. Therefore, we define it by:

\begin{displaymath}B_{\rm blob}(t) = c_{B} B_{\rm jet}(t).
\end{displaymath} (28)

With regard to the total IC radiation, according to our previous assumptions, we introduce some simplifications. All terms which appear in the formula for $I_{\rm tot}$ which concern emission coefficients for the jet ( $I_{\rm jet}$) are equal to zero. On the other hand, we consider absorption of the VHE gamma rays due to generation of the pairs. Thus, all terms with the absorption coefficients for the jet ( $k_{\rm jet}$) remain. The IC emission coefficient for the blob is calculated from the formulae given by Inoue & Takahara (1996) and Jones (1968). Here we calculate this coefficient under the assumption that only the local synchrotron radiation field generated within blob slices is scattered. The radiation field intensity for a given blob slice ( x = b ... a) is defined by:

\begin{displaymath}I^x_{\rm rad} = I^x_{\rm d} + I^x_{\rm blob} + I^x_{\rm u},
\end{displaymath} (29)


\begin{displaymath}I^x_{\rm d} = \sum_{i=b-1}^{x-1} I^{i-1}_{\rm d} \exp(-\tau^i_{\rm blob}) + I^i_{\rm blob},
\end{displaymath} (30)


\begin{displaymath}I^x_{\rm u} = \sum_{i=a+1}^{x+1} I^{i+1}_{\rm u} \exp(-\tau^i_{\rm blob}) + I^i_{\rm blob}
\end{displaymath} (31)

where $\tau^i_{\rm blob} = L^i_{\rm blob} k^i_{\rm blob}$ and the synchrotron intensity generated within the blob slice ( $I_{\rm blob}$) can be obtained by analogy to the jet slice (see Eq. (13)). Note that $I^x_{\rm d} = 0$ for x < b and $I^x_{\rm u} = 0$ for x > a. The optical depth for pair production is calculated according to the prescriptions given by Coppi & Blandford (1990) and Inoue & Takahara (1996). To simplify modeling, we calculate only absorption of the VHE gamma rays caused by this process. We do not analyze evolution and radiation of the pairs generated in this process. We perform these calculations to check how significantly the phenomenon can influence the radiation mechanism.

Finally we consider absorption of the IC spectrum by the Intergalactic Infrared Background (IIB) due to the pair generation process. We calculate this additional absorption using formulae provided by Malkan & Stecker (1998) selecting the weaker of the two proposed absorption conditions.

4 Application to Mrk 421

We apply our model to observations of Mrk 421. This object is well known as a variable source from the radio frequencies up to the VHE gamma rays. In this paper, we focus on the analysis and modeling of the radio observations, and compare them to high energy data in the perspective of multispectral investigation of high energy flares.

4.1 Observational data

First we apply our model to new observations performed at radio frequencies. The main set of data were obtained by the Torun Centre for Astronomy (TCfA) - Department of Radio Astronomy in Poland. The observations were carried out at the 32 m radio telescope using cooled ( $T_{\rm sys} \simeq 30$ K) beam switching receiving system at 5 GHz. For single determination of the flux density, we performed observation of Mrk 421 ( $F_{5 {\rm Ghz}} \sim 0.6$-0.7 Jy) and of the calibration source 3C 286 ( $F_{5 {\rm Ghz}} \sim 7.4$ Jy, Ott et al. 1994). We used the on-source/off-source method (e.g., Venturi et al. 2001) for observations. With this method, we performed a series of measurements for each source. In each series we measured the signal increment due to the source radiation eight times and the signal increment generated by the source of calibration noise, which appears in the receiving system, four times. The increments of the signal were calculated relative to the background sky radiation. This background was estimated at the distance of 3.5 half power beam width (HPBW) from the source. The accuracy of the measurements was determined from statistical computations which give a precision of about 20-30 mJy for observations obtained under the best weather conditions. The precision of the observations under poor weather conditions (cloudy sky and humidity of about 95%) is about 100 mJy. The degraded precision is related to the fluctuations of flux density generated during the propagation of radio waves through an unstable atmosphere.

The Torun data were gathered during two campaigns of radio observations. The first campaign was performed in January-April 2000. The second campaign of observations was carried out in January-March 2001. For the latter, in addition to observations at 5 GHz in Torun, we obtained data at 1.4 and 2.7 GHz from the Decimetric Radiotelescope of the Nancay Observatory in France.

Previous monitoring of the radio variability of Mrk 421, performed during several campaigns (e.g. Aller et al. 1985; Reich et al. 1993; Tosti et al. 1998; Aller et al. 1999; Venturi et al. 2001), showed that the level of the source radio emission varies on time scales of years or months. However, the temporal resolution of the previous campaigns was typically longer than a few weeks. Here, we present results of a radio monitoring with temporal resolution of about a few days. The data show two activity events with a duration time of about 30 and 15 days, for observations made in 2000 and 2001 respectively. Such relatively fast activities of Mrk 421 were probably observed also in 1992 (Reich et al. 1993) and in 1997 (Venturi et al. 2001). However, the temporal resolution of these observations was too low to precisely characterize the activity events.

We compared the results obtained from our radio campaigns with the X-ray light curves obtained from the RXTE-ASM experiment (2-10 keV) and VHE gamma-ray light curves from HEGRA telescope (Aharonian et al. 2002). We used also spectral shapes obtained in the high energy bands by RXTE (3-20 keV) and HEGRA (above 1 TeV) instruments (Krawczynski et al. 2001) to constrain physical parameters of the source.

In order to fix the level of the quiet emission of Mrk 421 we collected all available observations from the NED database. Note that the data are gathered by different types of optical telescopes, and thus the discrepancy in optical fluxes (see Fig. 3), is not related only to the variability, but is mainly due to the choice of different apertures for magnitude estimates.

   
4.2 Evaluation of physical parameters

To estimate the physical parameters which describe the source we use the method described by Katarzynski et al. (2001), derived from calculations originally proposed by Tavecchio et al. (1998). The method provides a relation between the Doppler factor and the magnetic field intensity based on the observed spectral features of the high energy emission of blazars. The main assumptions of this method are that the emitting region is spherical and homogeneous, with a uniform magnetic field. The relation $B-\delta$ can be obtained from the position of spectral breaks and peaks ( $\nu_{s_{\rm b}}, \nu_{s_{\rm p}}
\nu_{c_{\rm p}}$)[*] in the synchrotron and the IC spectrum, or from the total values of the synchrotron ( $F^t_{\rm s}$) and the IC ( $F^t_{\rm c}$) emission. To calculate total values of the emission, it is necessary to determine the slopes of the spectra ( $\alpha_1,
\alpha_2$) and values of the fluxes at the break and peak frequencies ( $F_{\rm s}(\nu_{s_{\rm b}}), F_{\rm s}(\nu_{s_{\rm p}}), F_{\rm c}(\nu_{c_{\rm p}}$)). In order to take into account absorption of the VHE radiation by IIB, the method considers two different positions of the peak of IC radiation ( $\nu_{c_{\rm p,min}}, \nu_{c_{\rm p,max}}$). Therefore, instead of two functional dependences we obtain a region of allowed physical parameters.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f2.eps}}\end{figure} Figure 2: Estimate of the Doppler factor value and magnetic field intensity obtained from the high energy emission of Mrk 421. Upper picture a) shows spectral shapes assumed for the estimate of the parameters B and $\delta $. Lower picture b) shows allowed ranges of the parameters values. The ranges are calculated for two different sizes of emitting region (R). The dot indicates values of B and $\delta $ selected for modeling in the frame of this work. The high energy data come from the work of Krawczynski et al. (2001).
Open with DEXTER

Estimates of the parameters are deduced from a high energy state gathered almost simultaneously by the RXTE and HEGRA experiments. We use also observations gathered a few hours later by RXTE for which there are no simultaneous observations from HEGRA. Therefore, we assume that $\nu_{c_{\rm p}}
F_{\rm c}(\nu_{c_{\rm p}})$ for the second state was only slightly higher in proportion to the second observation made by RXTE. The data and the spectral shapes are shown in Fig. 2a. Detailed values of the parameters deduced from the spectral properties are given in Table 1.

The range of allowed values of B and $\delta $ depends strongly on the size of the emitting region (R). We present two different results derived for $R = 4.5~ \times~ 10^{16}$ cm and $R = 3~ \times~ 10^{16}$ cm (Fig. 2b) to illustrate this effect. There is no possibility yet to definitely fix the size of the emitting region. The size can be estimated only from the variability time scales ( $R \leq c~t_{\rm var} \delta / (1+z)$), which vary from hours to days. Therefore we cannot precisely determine the physical parameters which describe the source. However, testing our variability scenarios (see Sect. 4.4), we found out that relatively weak value of the magnetic field intensity ( $B \sim 0.01$ G) allows the multifrequency emission of Mrk 421 to be explained better. Thus, finally we can select Doppler factors of about 20 and magnetic field intensities of about 0.02 G. According to the estimates, these values correspond to a source radius of about $4.5 \times 10^{16}$ cm which gives a reasonable variability time scale $t_{\rm var} \gtrsim 20$ h.

From the above estimate we also derive the ratio ($\eta$, Eq. (23)) between radiation field energy density ($U_{\rm r}$) and magnetic field energy density ($U_{\rm B}$) from the relation:

\begin{displaymath}\frac{F^{t}_{\rm s}}{F^{t}_{\rm c}} = \frac{U_{\rm B}}{U_{\rm r}} = \eta
\end{displaymath} (32)

(e.g., Rybicki & Lightman 1979) which gives $\eta \sim 3$-4.

The values of the main parameters estimated in such a way will be now used as initial parameters for modeling the jet emission. In our scenario, the observed high energy activity is generated at the base of the jet. Therefore, physical conditions at the base should not significantly differ from those estimated from the high energy emission.

Table 1: Quantities used to describe approximate spectral shapes shown in Fig. 2a. These spectral shapes were used to constrain the allowed parameter range.

4.3 Quiet emission

In the blob-in-jet model, most of the low energy emission of Mrk 421 is generated by the synchrotron radiation of the jet. The emission of the jet is almost constant on relatively short time scales (weeks and months) but may vary on time scales of years. However, we neglect here modeling of such long term variations.

The magnetic field intensity inside the external component has been fixed equal to the magnetic field intensity inside the last jet slice. Thus, we obtain a quite low value of this parameter ( $ B_{\rm ext} \sim 10^{-3}$ G). A relatively small value of the Doppler factor for such a component ( $\delta_{\rm ext} \leq 5$) appears realistic. Therefore according to the formula:

 \begin{displaymath}\nu_{\rm syn} = \frac{\delta_{\rm ext}}{1+z} \nu'_{\rm syn} \...
...{4}{3} \gamma^2 \left( \frac{e B_{\rm ext}}{2 \pi m c} \right)
\end{displaymath} (33)

we derive a relatively large value of the maximal electron energy ( $\gamma
\sim 10^4$) to explain the low frequency radio radiation. The electrons are systematically cooled during their evolution along the jet. Therefore, to explain the existence of such relatively high energy electrons inside the external component we allow continuous particle acceleration during their propagation along the jet. Thus, arriving at the end of the jet, particles have enough energy to generate the low frequency radio emission.

Illustration of this model of the quiet emission is shown in Fig. 3. The corresponding physical parameters are given in Table 2 (2001). This quiet state will be used in the following paragraph as a reference for non-active state while modeling the variability observed in 2001. The radio flux density observed at 5 GHz in 2001 ($\sim$600 mJy) was smaller than emission observed at the same frequency in 2000 ($\sim$725 mJy). This is probably related to a more complex evolution of the jet on time scales of years, which we neglected here. To take this fact into account, we simply change the particle density inside the jet ( $ K^1_{\rm jet} $) between the two epochs. Note that because of the specific parameterization of the electron energy distribution for the extended component (Eq. (16)) we have to change also the particle density inside the extended component ( $c_{\rm ext}$) to keep a constant level of radiation.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f3.eps}}\end{figure} Figure 3: The modeling of quiet emission of the jet (solid line) extended radio lobes (dotted-dashed line) and thermal emission (long-dashed line). The short-dashed lines show local synchrotron emission of every tenth jet slices. The observations were obtained from the NED database.
Open with DEXTER

   
4.4 Test of variability scenarios

In this section we investigate two general scenarios which can explain variability of Mrk 421 and other blazars from VHE gamma rays to the radio frequencies. The first approach assumes rapid injection of fresh particles. The second scenario explains activity by systematic acceleration of the electrons. In this particular case we simulate these processes at the base of the jet. We apply our tests to high energy observations obtained by the RXTE and HEGRA experiments (see Sect. 4.2). We use values of the Doppler factor and the magnetic field intensity estimated from these observations and we are try to explain observed spectral shapes as well as possible in framework of the proposed scenarios.

Table 2: Physical parameters used for the modeling of quiet state emission of Mrk 421. In the modeling of the radio observations gathered in 2001 we use almost the same set of parameters as for the monitoring in 2000. Only two parameters ( $K^1_{\rm jet}, c_{\rm ext}$) have slightly different values, which are given in the column labelled 2001. The parameters indicated by ($\star $) are estimated from the high energy emission of Mrk 421 shown in Fig. 2.

Equation (7), which describes evolution of the electron energy, does not describe injection of particles into any emitting region. Thus, to simulate the injection process we use our specific parameterization of the jet geometry. We introduce at the beginning of the jet new slices which are characterized by a maximal energy of the electrons significantly higher ( $\gamma^{\rm cut}_{\rm blob} \sim
10^6$-107) than average maximal energy of the electrons at the jet base ( $\gamma^{\rm cut}_{\rm jet} \sim 10^4$-105). Such high energies are necessary to generate high energy flares in leptonic models. To increase the efficiency of the IC scattering, we increase density of the particles inside the injected blob slices ( $K^1_{\rm blob}\sim 10^2$ cm-3) relatively to the density inside the first jet slice ( $K^1_{\rm jet} \sim 10$ cm-3). Such injection of a few new slices generates a compact and energetic blob at the base of the jet. Emission of such component causes growth of source brightness as long as new slices appear at the beginning of the jet. In this particular case, we generate the blob by injection of sixteen slices, therefore the injection time is equal to $t^{16}_{\rm d} = 18$ days in the blob frame (about 22 hours in the observer's frame).

The particles inside blob slices are cooled due to the adiabatic expansion and by synchrotron and the IC radiative losses. Thus, the brightness of the radiation quickly decreases when the injection stops. The light curves generated with this approach are shown in Fig. 4. The scenario can explain well flares observed in X-rays and gamma rays. However, on the other hand, the activity generated at lower frequencies is quite weak. The physical parameters used for this scenario are given in Table 3 in the column labeled inj.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f4.eps}}\end{figure} Figure 4: Modeling the multifrequency variability of a blazar. Examples of light curves obtained when the variability is due to injection of fresh high energy electrons at the beginning of the jet.
Open with DEXTER


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f5.eps}}\end{figure} Figure 5: Example of light curves generated by in situ acceleration of electrons at the base of the jet.
Open with DEXTER

In the second scenario we assume that the high energy electrons are produced by acceleration within some compact region at the base of the jet. We call this region the "blob'', but it can be a shock wave or a localized turbulence which accelerates electrons. We describe the emitting region as in the injection scenario, by a system of moving and evolving slices. The electron energy distribution for blob slices which appear at the beginning of the jet is very similar to the distribution selected for the jet ( $K^1_{\rm blob} = K^1_{\rm jet},~\gamma^{\rm cut}_{\rm blob} =
\gamma^{\rm cut}_{\rm jet}$). We assume in the modeling that the acceleration rate is initially very efficient ( $A^1_{\rm blob} / A^1_{\rm jet} \sim 25,~A^1_{\rm blob} =
6.8 \times 10^{-6}$ 1/s $ \to t^1_{\rm acc_{\rm blob}} = 1.47 \times 10^5$ s) however, it decreases rapidly during the evolution process ( $a_{\rm blob} = 4.5$). The efficiency of the acceleration process can be estimated from the maximal energy of the electrons ( $\gamma^{\rm max}_{\rm blob}$) which is necessary to explain the maximal observed frequency in the X-ray range ( $\nu^{\rm max}_{\rm X}$). Assuming that $\nu^{\rm max}_{\rm X} = 10^{19}$ Hz (see Fig. 2) and using formula (33), where we use values of the magnetic field intensity (B = 0.02 G) and the Doppler factor ( $\delta = 20$) estimated from the high energy emission, we obtain $\gamma^{\rm max}_{\rm blob} = 2.6 \times 10^6$. The characteristic cooling time due to synchrotron radiation for such electrons is of about $7.3 \times 10^5$(s). Introducing also losses due to the IC scattering ($\eta = 3$) we obtain $t^{\rm cool}_{\rm blob}(\gamma^{\rm max}_{\rm blob}, t^1_{\rm jet}) = 5.5 \times 10^5$ (s). Therefore, to accelerate electrons to required energy $\gamma^{\rm max}_{\rm blob}$we have to assume a process characterized by acceleration time at least equal to $t^{\rm cool}_{\rm blob}(\gamma^{\rm max}_{\rm blob})$. Detailed calculations shows that the best explanation of the observed spectral shapes can be obtained for an initial acceleration time more than three times shorter than the initial cooling time. This is related to the decrease of the cooling efficiency ( $C^{\rm cool}_{\rm blob} \sim t^{-1}$), which is slower than the decrease of the acceleration efficiency ( $C^{\rm acc}_{\rm blob}
\sim t^{-4.5}$). The characteristic time value which describes the initial losses of the particle energy due to the adiabatic expansion ( $t^1_{\rm adia_{\rm blob}} = 1/C^{\rm adia}_{\rm blob}(t^1_{\rm jet})$) is about one order of magnitude longer ( $t^1_{\rm adia_{\rm blob}} = 3.6 \times 10^6$ (s)) than the initial cooling time. Therefore, we do not consider this effect in the estimations above.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f6.eps}}\end{figure} Figure 6: Temporal evolution of the spectral energy distribution of Mrk 421 from the radio frequencies up to TeV. The solid lines show model of the total emission for the lowest of the two high energy states observed by RXTE in February 2000. The short-dashed lines show a fit for the highest state. The dotted lines show the synchrotron emission of the blob. Both fits are obtained under the assumption that the high energy emission is generated by acceleration of electrons, and the letter fit is deduced from the evolution of the former in the framework of this acceleration scenario. Observations gathered by the HEGRA telescope are simultaneous to the lowest of the two states of RXTE. The high energy data come from the work of Krawczynski et al. (2001). The bold-solid line shows the synchrotron emission of the jet, the long-dashed line shows the thermal emission of the host galaxy, and the dashed-dotted line shows the radiation of the extended radio structure. These three last contributions remain basically constant on the time scale of the active event.
Open with DEXTER

The brightness of the high energy radiation grows as long as the acceleration is efficient enough to overcome energy losses due to adiabatic expansion and radiative cooling. In this particular case, losses related to the radiative processes become strongest than the acceleration process after period of about 1.5 days measured in the observer's frame. Losses related to the adiabatic expansion become strongest than the acceleration process after about 2.5 days. When acceleration become too weak to compensate for the losses, the intensity of the emission starts to decrease. However, in this second scenario, the acceleration may compensate the losses of electron energy for a relatively longer time. Thus, significant activity is also generated at the radio frequencies. The time delay between the maximum of activity at the high energies and the maximum of the radio activity results from absorption of the radio photons by electrons. Initially, the emitting region is optically thick at the radio wavelengths but it expands during evolution and finally appears optically thin at the radio frequencies. Moreover, absorption of the blob radiation by jet electrons appears to be important in this case. Even if the blob is optically thin, the radio radiation of this component may be significantly absorbed by electrons inside the jet. The light curves generated in this scenario of particle acceleration are shown in Fig. 5. The physical parameters for this modeling are given in Table 3, in column labelled acc.

The Spectral Energy Distributions (SED) shown in Fig. 6 are generated by acceleration of particles. Testing this approach we find that we have to use relatively low values of the magnetic field intensity ( $B \sim 0.02$ G) to explain the evolution of observed spectral shapes. For larger values of B it is necessary to assume more efficient acceleration conditions, to overcome cooling of the electrons due to synchrotron and IC radiation. Such cooling processes affect only the high energy part of the electrons spectrum. Therefore, using more efficient acceleration processes with an energy independent acceleration ratio we obtain the maximal required energy of electrons. On the other hand, this significantly increases the energy of low and middle energy electrons and generates a strong emission from low energy (a few keV) X-rays to the radio frequencies, which is not observed.

Note that for set of physical parameters which we use in the model, absorption of VHE gamma rays due to pair production inside the blob and inside the jet appears negligible.

Table 3: Physical parameters used for modeling of variability. The upper indices 1/1/50/50 indicate the number of the jet slice where activity starts.

In order to correctly fit the SED, we have to select a slightly steeper slope for the initial electron energy distribution inside the blob than that inside the jet ( $n_{\rm jet} \ne n_{\rm blob}$). This can be due to some previous energy dependent acceleration efficiency, which may induce different slopes. However, energy dependent acceleration greatly complicates the solving of Eq. (7) and is not supported by our formula (A.1). So we set slightly different slopes in the blob and in the jet as initial conditions. This difference between the slopes is responsible for the drop in the optical emission which is seen at the very beginning of the activity (Fig. 5c). In this particular case we assume that $K^1_{\rm blob} = K^1_{\rm jet}$ and $R^1_{\rm blob} = R^1_{\rm jet}$. Therefore, at the very beginning of the activity, the emission of the first blob slices in the optical bands is slightly weaker than the emission of the first jet slices. This quickly changes when acceleration processes modify the electron spectrum inside the blob slices. Note that this effect is specific to this case, with its particular set of the parameters, and does not always appear.

The position of the blob inside the jet for the SED presented in Fig. 6 is indicated by a=19 and a=22 for the low and the high state respectively. This corresponds to a temporal shift of about 6 hours, well in agreement with the observations.

It is also possible to fit the observed levels of the emission using the injection scenario. However, in this case we cannot reproduce the precise evolution of the SED because the X-ray radiation grows slightly faster than the VHE radiation during the injection process (see Figs. 4c, d).

Although we have fitted the SED observed by RXTE and HEGRA instruments, we have also analyzed the light curves gathered at that time by the ASM-RXTE experiment. There is confirmation of the high energy activity, but it seems that the activity was possibly more complex than just one single flare, with a few shorter time scale flaring events superposed on the main one. However, the data available did not allow a detailed analysis of this complex pattern on short temporal scales. Such an effect of superposition of several short events on a main longer flare was already found in our previous paper in the case of Mrk 501 (Katarzynski et al. 2001).

4.5 Application to the observations

We apply our model to the two activity events of Mrk 421 observed in radio during our monitoring. The first event was in 2000 (around MJD 51621) and had apparently no high energy counterpart. The second one was in 2001 (around MJD 51970) and was observed also at high energies.

We analyze first the data obtained in 2001 because we have more observations gathered at that time (Fig. 7). Looking at the light curves registered at radio frequencies, we can see a significant increment of the emission at 5 GHz (of about 300 mJy) and 2.7 GHz (of about 500 mJy) but no activity at 1.4 GHz. This agrees qualitatively with the theoretical investigations performed in the previous section (Fig. 5a) in case of particle acceleration. Moreover, simultaneously to the radio flare, we noticed that the active event was observed in X-rays, as shown by the light curve gathered by RXTE-ASM (Fig. 5c). The activity was also observed at VHE gamma rays by HEGRA telescope (Fig. 5g, Aharonian et al. 2002). The absence of any time delay between the X-rays and the gamma rays and our highest radio frequencies, suggests that the emitting region was from the beginning of its evolution partially thin for the radio radiation and also absorption of the blob radiation by jet's electrons was not efficient.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f7.eps}}\end{figure} Figure 7: The multifrequency variability of Mrk 421 observed from the radio frequencies by Nancay and Torun Radio Astronomy Observatories to the X-rays by RXTE-ASM experiment in February-March 2001. Theoretical light curves are generated under the assumption that activity originates in the middle part of the jet by compression and acceleration of particles. Plot g) shows the diurnal integral fluxes observed above 1 TeV by the HEGRA experiment (Aharonian et al. 2002). The simultaneity between the VHE gamma-ray flare detected by HEGRA and our radio event appears better than a fraction of a day.
Open with DEXTER

To simulate such an emitting region we place it initially in the middle part of the jet (i=50) and apply the acceleration scenario. Apart from the acceleration we also introduce some weak compression of the local jet medium which increases the particle density ( $K^{50}_{\rm blob}/K^{50}_{\rm jet} = 2.5$) and the magnetic field intensity (cB = 1.5). According to our specific parameterization of the jet geometry, expansion of jet slices in the middle part of the source is much slower than at the beginning of the source. Therefore the acceleration cannot be efficiently compensated by adiabatic losses and decrease of the particle density if the density decreases according to the expansion of a slice volume ( $d_{\rm jet}=(2r_{\rm jet}+1)/3$). To overcome this problem we assume that the density of particles inside the emitting region decreases faster than the decrease due to the slice expansion ( $d_{\rm blob} > d_{\rm jet}$). This implies that there is some additional escape of particles from the emitting region. We describe the escape by a parameter: $e_{\rm blob} = 3(d_{\rm blob} - d_{\rm jet})$, which is equivalent to additional right hand side term: $-N(\gamma, t) e_{\rm blob}/t$, in Eq. (7). The escape process can be described also by a characteristic escape time ( $t_{\rm esc_{\rm blob}}(t) = t / e_{\rm blob} \to t^i_{\rm esc_{\rm blob}} = t^i_{\rm jet}/e_{\rm blob}$) which defines how long particles remain within the source boundaries. The value of the escape time can be estimated from the linear size of a source. In principle, estimation of the escape time for single slice should be performed according to the relation: $t^{i}_{\rm esc-est} \gtrsim L^i_{\rm blob} / c$. However, this approach is not precise in the case, where the blob is created by several slices. The particles which escape from a given slice may pass to neighbouring slices and vice versa. We do not simulate such detailed movements of the particles inside the blob in our model. However, we consider this effect in estimations of the escape time, where we use total length of the blob: $t^i_{\rm esc-est} \gtrsim (\sum_{j=i-n_{\rm b}+1}^{i} L^j_{\rm blob}) / c$. By analogy to the escape time we can define another characteristic time which describes the decrease of the particle density due to the expansion of a slice: $t^i_{\rm exp_{\rm blob}} = t^i_{\rm jet}/(2r_{\rm jet}+1)$. Now, we can define a time which describes escape and expansion at once: $t^{i}_{\rm ee_{\rm blob}} = t^i_{\rm jet} /
(2r_{\rm jet} + 1 + e_{\rm blob}) = t^i_{\rm jet} / (3 d_{\rm blob})$.

We compare the above defined parameters with assumed acceleration time. However, the blob's slices are expanding during evolution, therefore to compare characteristic times, we calculate an average value ($\bar{t}$) for each of them: $\bar{t}^{~i}_{xx_{\rm blob}} = t^{i-n_{\rm b}+1}_{xx_{\rm blob}} + (t^{i}_{xx_{\rm blob}} -
t^{i-n_{\rm b}+1}_{xx_{\rm blob}}) / 2$. In addition we compare parameters used for the model presented in the current subsection, with the parameters applied for the model shown in Sect. 4.4 (no escaping). We perform our comparisons for the beginning of the bob evolution which means $i=n_{\rm b}$ for the test shown in Sect. 4.4 and $i=50+n_{\rm b}-1$ for the model discussed here. First, we compare $\bar{t}_{\rm exp}$ with $\bar{t}_{\rm acc}$ for the test performed in Sect. 4.4 and we obtain: $\bar{t}^{n_{\rm b}}_{\rm exp_{\rm blob}} / \bar{t}^{~n_{\rm b}}_{\rm acc_{\rm blob}} = 2.89$, while for calculations presented here this ratio is 7.79. This means that the decrease of the particle density in previous calculations was significantly faster than in the current modeling. Next, we compare average acceleration time with time, which describes the decrease of the density due expansion and escape for the calculations presented here. The comparison gives the ratio: $\bar{t}^{~(50-n_{\rm b}+1)}_{\rm ee_{\rm blob}} / \bar{t}^{~(50-n_{\rm b}+1)}_{\rm acc_{\rm blob}} = 3.58$, very similar to the ratio obtained for calculations performed in Sect. 4.4. To obtain this value and to obtain the best fit of the observed activity we applied an average escape time ( $\bar{t}^{~(50-n_{\rm b}+1)}_{\rm esc_{\rm blob}}$) equal to $5 \times 10^6$ (s) while estimated time ( $t^{50-n_{\rm b}+1}_{\rm esc-est}$) is of about $3.2 \times 10^6$ (s).

The results of the modelling and estimations performed above show that we obtained the best fits for parameters where the ratio of average initial times which describes expansion or expansion and escape to the average initial acceleration time is about three. In the model of Sect. 4.4, a specific parameterization of the geometry provides relatively fast expansion of the blob. Therefore we do not need to assume additional escape. In the model presented in the present subsection we placed the blob in middle part of the jet. Thus, expansion of the blob's slices is relatively slow and it appears necessary to assume some escape of the particles from the emitting region. However, the value of the assumed escape time seems to be in good agreement with the estimated value of this parameter. Detailed values of the physical parameters of this model are given in Table 3 in the column labelled "2001''.

Such a model developed only based on radio and X-ray properties provides a globally good fit to the radio and X-ray light curves. However, it appears difficult to reproduce the exact level of the radiation observed at 2.7 GHz, which suggests that the real mechanism of radio emission is somewhat more complex than proposed here. The light curves computed in the optical bands predict a growth of the source brightness from 0.2 (V band) to 0.4 (U band) magnitude. The model provides a good explanation of the X-ray data if we assume an average level of the X-ray radiation at that time of about $8 \times 10^{-10}$ erg cm-2 s-1. Also activity in VHE gamma-ray detected by the HEGRA experiment, which appears simultaneously with the radio and the X-ray flares as illustrated in Fig. 7g is well explained. However, to explain the VHE activity we had to assume some non zero average level ($\sim$2.5 ph/(cm2 s)) of such radiation as well as for the X-ray radiation. We explain such observed non zero average levels of emission as remnants of previous activity events. However, for the sake of simplicity we do not simulate activity events which could precede the flare analyzed here.


  \begin{figure}
\par\resizebox{\hsize}{!}{\includegraphics{MS2898f8.eps}}\end{figure} Figure 8: The variability of Mrk 421 observed at 5 GHz by TCfA - Department of Radio Astronomy in March-April 2000. Full lines show the light curves predicted at three frequencies by a scenario with weak acceleration and compression, which fits the data at 5.0 GHz. Note that electron self-absorption effects can sometimes induce anti-correlation of fluxes observed at different radio frequencies.
Open with DEXTER

During our first campaign of observations, performed in January-April 2000, we obtained radio observations only at 5 GHz. There was one activity event (Fig. 8) observed in March-April (around MJD 51620). Note that this activity has been confirmed by independent observations made at IRAM at 100 GHz (Ungerechts, H., private communication). To obtain more information about the event, we analyzed data gathered by the RXTE-ASM experiment at 2-10 keV. The X-ray data do not show any activity of the source at that time. The emission is almost constant at level of $5 \times 10^{-10}$ erg cm-2 s-1. Moreover observations gathered at VHE gamma rays up to MJD 51615 by the CELESTE and CAT experiments (Holder et al. 2001), do not show significant variability. For the numerical simulations of this event, we can keep almost the same set of fundamental physical parameters (Table 3, "2000'') as for the multispectral flare observed in 2001. However, this time we have to consider an acceleration process which is not efficient enough ( $A^{50}_{\rm blob} / A^1_{\rm jet} \sim 3.6$) to generate the high energy electrons needed for the high energy activity. Therefore, instead of implementing an efficient acceleration mechanism, we set a higher compression of the local jet medium (about four times). This allows us to reproduce the observed level of the radio radiation without any high energy counterparts.

Such a model provides radio light curves very similar to the light curves generated by particle acceleration (Fig. 5a). However, this illustrates the possibility of a quite interesting effect, clearly visible in the light curves generated at 1.4 and 2.7 GHz, where the brightness of the source initially decreases. This can be explained as an absorption of the jet radiation which originates above the blob (between the central engine and the blob) by the electrons inside the blob.

As a conclusion, we can say that our observational data on Mrk 421 show for the first time in a BL Lac object the existence of a multispectral flare observed simultaneously in radio, X-ray and TeV gamma-ray bands. They emphasize the existence of at least two types of radio flares, with or without high energy counterparts. Both of them can be explained in the frame of our acceleration scenario with local compression of the jet medium. Such a scenario can generate reasonable light curves of the flares at all frequencies.

5 Summary

We have developed a model which is able to explain multispectral variability of blazars from the VHE gamma rays to the radio frequencies. We applied it to new observational data, which we gathered in the radio range. We assume that observed activity events are, in the simplest case, generated by a single compact component. The component can be generated by fast injection of energetic particles form the central engine at the base of the jet or from some acceleration zone. In alternate scenario we generate the component by systematic acceleration of the particles.

To test the above scenarios, we developed an original method for calculation of the radiation transfer along the jet. The method allows the electron self-absorption processes inside the source to be taken into account precisely. This appears mandatory if we want to generate precise light curves at the radio frequencies. With this method we are also able to calculate absorption of the VHE gamma-ray photons generated by IR synchrotron photons generated within the jet. However, for the set of physical parameters which we selected for the model, this process appears negligible here.

Within our model we calculated detailed evolution of the electron energy spectrum along the jet. The calculation shows that for the selected value of the magnetic field intensity, additional acceleration processes are necessary to explain energies of electrons in the outer parts of the source.

We investigated two alternative scenarios to account for the multifrequency emission of Mrk 421. Firstly, we assume that activity is generated by injection of particles. This approach can qualitatively explain the high energy activity but appears unable to explain the activity in radio. In the second approach, we generate activity by particle acceleration which is found to be more adapted to reproduce radio and high energy emission of Mrk 421. The results of our model are consistent with results obtained by other authors with somewhat similar approaches (Li & Kusunose 2000; Kusunose et al. 2000; Blazejowski et al. 2000; Sikora et al. 2001).

We successfully applied our model to recent observations obtained from the radio frequencies to X-rays and gamma rays. The observations gathered in February-March 2001 show a clear correlation between the high energy and the radio activity which can be well reproduced by our acceleration scenario. This is, to the best of our knowledge, the first detection of a multispectral flare in a BL Lac object with well-defined activity seen simultaneously in radio and in X-rays and TeV gamma rays. This analysis illustrates how monitoring of blazar variability at radio frequencies introduce a significant amount of informations about activity of the sources and may be very important for a better understanding of high energy blazar activity.

Acknowledgements

We thank D. Horns, H. Krawczynski and F. Aharonian for information and data obtained by the RXTE and HEGRA experiments. We are grateful to J. M. Martin for help with observations and in reduction of Nancay data and to M. Punch for improvement of the article.

This article made use of observations gathered by RXTE/ASM experiment, provided by HEASARC, a service of NASA/Goddard Space Fight Center. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

The project was partially supported by the "Jumelage'' for French-Polish collaboration, from the French Minister of Research, Education and Technology.

Appendix A: Solution of the kinetic equation

The general solution of the Eq. (7) for the initial conditions assumed in Eq. (8) is given by:

 \begin{displaymath}N^{\ast}_{\rm jet/blob}(\gamma, t) = \left[K^1 S_1(\gamma, t)...
...jet/blob}} -
K^2 S_1(\gamma, t)^{2} \right]
S_2(\gamma, t),
\end{displaymath} (A.1)

where

\begin{displaymath}S_1(\gamma, t) = [I_1(t) - \gamma I_2(t)]~\gamma,
\end{displaymath} (A.2)


                  $\displaystyle S_2(\gamma, t)$ = $\displaystyle \exp \Bigg[- \int^{t}_{t^{1}_{\rm jet}} \Big\{ \gamma C^{aa}(t_x) I_2(t_x))$  
  - $\displaystyle \gamma C^{aa}(t_x) I_2(t) + C^{aa}(t_x) I_1(t)$  
  - $\displaystyle 2 I_1 (t_x) C^{\rm cool}(t_x) \Big\}$ (A.3)
  / $\displaystyle \Big\{ \gamma I_2 (t_x) - \gamma I_2 (t) + I_1(t) \Big\} {\rm d} t_x \Bigg],$  


\begin{displaymath}I_1(x) = \exp \left[\int_{t^{1}_{\rm jet}}^{x} C^{aa}(y) {\rm d}y\right],
\end{displaymath} (A.4)


\begin{displaymath}I_2(x) = \int_{t^{1}_{\rm jet}}^{x} C^{\rm cool}(y) I_1(y) {\rm d}y,
\end{displaymath} (A.5)

and $K^1 = K^1_{\rm jet/blob}$, $K^2 = K^1_{\rm jet/blob}(\gamma^{\rm cut}_{\rm jet/blob})^{-2-n_{\rm jet/blob}}$, $C^{\rm cool}(t) = C^{\rm cool}_{\rm jet/blob}(t)$, $C^{aa}(t) = C^{\rm acc}_{\rm jet/blob}(t) - C^{\rm adia}_{\rm jet/blob}(t)$. Note that similar but less complex solutions were studied in detail by Kardashev (1962). We used the Kardashev's solutions to test our more general result.

Appendix B: Radiation transfer along the jet

Radiation transfer along the jet for the case where an emitting blob is placed inside the jet is described by:

\begin{displaymath}I^{j}_{\rm tot} = \sum^l_{i=j} I^{i-1}_{\rm tot} \exp \left(-\tau^{i}_{x} \right) + I^i_{x},
\end{displaymath} (B.1)

where

\begin{displaymath}I^i_x = \left\{\begin{array}{ll}
I^i_{\rm blob}, & \mbox{$b ...
...& \mbox{$i < b ~~ {\rm or} ~~ i > a$ }\\
\end{array}\right.,
\end{displaymath} (B.2)


\begin{displaymath}\tau^i_x = \left\{\begin{array}{ll}
L^i_{\rm blob} k^i_{\rm ...
...t}, & \mbox{$i < b~~{\rm or}~~i > a$ }\\
\end{array}\right.,
\end{displaymath} (B.3)

where $I^j_{\rm tot}$ is calculated under the assumption that $I^j_{\rm tot} = 0$ for i-1 < j. Definitions of the parameters a and b are given in Sect. 3.1.

References



Copyright ESO 2003