© ESO, 2010

1. Introduction

Blazars, a subclass of active galactic nuclei, show variability that is observed across the electro-magnetic spectrum and can, in some cases, be as fast as a few minutes (Aharonian et al. 2007; Albert et al. 2007). The working blazar scenario is that we observe their radiation coming from a jet that is directed, within a small angle, to our line of sight. The observed spectrum, which is clearly nonthermal in origin, shows a characteristic double peak in a νF_ν plot, and is produced from a population of relativistic particles, presumably accelerated by shock waves within the jet. While there are still open questions concerning the origin of the high-energy gamma radiation, there is a consensus that the lower frequency emission, which usually extends from radio to X-rays, comes from electron synchrotron radiation. The observed variability can then be attributed directly to conditions in the acceleration region (for a recent review of the various blazar models see Böttcher 2010).

The detailed X-ray observations of many blazars reveal that these objects show a rich and complex structure in their flaring behavior. Especially interesting are the trends in the lags between the soft and hard X-ray bands. The first observations, which were detailed enough to show time structure in various energy bands (Takahashi et al. 1996), revealed the blazar Mrk 421 as having a soft-lag flare, meaning the soft X-rays lagged behind the hard ones. This was explained as synchrotron cooling: lower energy electrons cool with a slower rate than higher energy ones and therefore soft X-rays appear after the harder ones.

Subsequent observations, however, showed cases of the opposite trend (Takahashi et al. 2000). This could not be explained by synchrotron cooling, but it was attributed to particle acceleration. As fresh particles are accelerated to high energies, they radiate first soft and then later hard photons; i.e. hard lags behind soft.

Long, uninterrupted observations (Brinkmann et al. 2005) have revealed that Mrk 421 shows both kinds of lags, with typical timescales on the order of 10³ s. The same mixed behavior has also been verified for the same source by Tramacere et al. (2009). Finally, blazar 1ES 1218+304 showed a hard-lag pattern during a giant flare (Sato et al. 2008).

While several time-dependent models have been put forward to explain time variations (Mastichiadis & Kirk 1997; Li & Kusunose 2000; Böttcher & Chiang 2002), hard lag behavior needs a particle acceleration scheme to be taken explicitly into account (Kirk et al. 1998; Kusunose et al. 2000). Kirk et al. (1998, hereafter KRM) have shown that both soft and hard lags can be obtained during a flare once a detailed scheme for particle acceleration is considered. Adopting a two-zone model in which particles are energized in an acceleration zone by the first-order Fermi process and then escape and radiate in a radiation zone, they argue that flares produced by an increase in the number of the injected particles can produce soft lags if the radiating electrons are in the energy regime where cooling is faster than acceleration and hard lags when the two timescales are comparable.

More recently, two-zone models that include SSC losses have been developed in (Graff et al. 2008) and (Weidinger & Spanier 2010). In particular, Weidinger & Spanier include both first and second-order Fermi acceleration in their acceleration zone. A detailed study of the two types of acceleration and of their effect on SSC spectra, but in a one-zone model, can be found in Katarzyński et al. (2006).

In the present paper we further investigate the expected variability within the two-zone model by letting its key parameters vary in time. While we keep the same assumptions as in KRM, we add radiation coming from the acceleration zone to the photon spectrum. In Sect. 2 we present the basic points of the model and then examine some interesting cases in Sect. 3. In Sect. 4 we apply our results to 1ES 1218+304 and discuss them in Sect. 5.

2. The model

We employ the two-zone acceleration model as developed by KRM. According to this, electrons are accelerated by a nonrelativistic shock wave that propagates along a cylindrical jet, and then they escape in a wider region behind the shock where they radiate the bulk of their energy. Following the KRM notation, we refer to the region around the advancing shock as the “acceleration zone” (AZ) and to the escape region as the “radiation zone” (RZ). The total size of the source is assumed to have a finite extent.

2.1. Acceleration zone

The evolution of the electron distribution function (EDF) in the AZ, N(γ,t), is given by the continuity equation $\frac{\mathit{\partial N}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial \gamma }}\mathrm{\left(}\mathit{N}\mathit{\gamma ̇}\mathrm{\right)}\mathrm{+}\frac{\mathit{N}}{{\mathit{t}}_{\mathrm{esc}}}\mathrm{=}\mathit{Q\delta }\mathrm{(}\mathit{\gamma }\mathrm{-}{\mathit{\gamma }}_{\mathrm{inj}}\mathrm{)}\mathit{.}$(1)The energy of each particle in this region changes with the rate $\mathit{\gamma ̇}\mathrm{=}\frac{\mathit{\gamma }}{{\mathit{t}}_{\mathrm{acc}}}\mathrm{-}\mathit{\alpha }{\mathit{\gamma }}^{\mathrm{2}}\mathit{,}$(2)where the first term describes first-order Fermi-type acceleration with the energy-independent timescale t_acc, and the second describes the energy losses due to synchrotron radiation in the magnetic field B with $\mathit{\alpha }\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{T}}{\mathit{B}}^{\mathrm{2}}}{\mathrm{6}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}\mathit{c}}\mathrm{·}$(3)Particles are injected in the AZ at some energy γ_injm_ec² with the rate Q and escape from it to the RZ with the energy-independent rate .

In our model we assume that variability can be produced by the change of one or more of the parameters t_acc, t_esc, Q and B in the AZ and thus treat them as time-dependent in Eq. (1). We define the dimensionless functions that include the time dependency through the relations t_acc(t) = t_a,0f_a(t), t_esc(t) = t_e,0f_e(t), Q(t) = Q₀f_q(t), and B(t) = B₀f_B(t), where t_a,0, t_e,0, Q₀, and B₀ are the unperturbed values of the parameters.

The solution of Eq. (1) can be obtained semi-analytically with the method of characteristics. This method requires the solution of the initial value problem $\frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{\gamma ̇}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}\mathit{,}\mathit{\gamma }\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}{\mathit{\gamma }}_{\mathrm{inj}}\mathit{,}$(4)which can be written in the form $\frac{{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathit{\gamma }}\mathrm{-}\frac{\mathit{\psi }\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{\gamma }}_{\max }}\mathrm{=}\frac{{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}}{{\mathit{\gamma }}_{\mathrm{inj}}}\mathrm{-}\frac{\mathit{\psi }\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}}{{\mathit{\gamma }}_{\max }}\mathit{,}$(5)where ${\mathit{\gamma }}_{\max }\mathrm{=}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}{\mathrm{)}}^{-1}\mathit{,}$(6)and the functions ϕ_a, ψ are given by the relations ${\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\exp \left({\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}\mathit{/}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}}\frac{\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}}{{\mathit{f}}_{\mathrm{a}}\mathrm{\left(}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}}\right)$(7)and $\mathit{\psi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}\mathit{/}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}{\mathit{f}}_{\mathrm{B}}^{\mathrm{2}}\mathrm{\left(}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}$(8)The solution of Eq. (1) can then be written in the generic closed form $\mathit{N}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{Q}\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}{\mathit{\gamma }}_{\mathrm{inj}}^{\mathrm{2}}}{{\mathit{\gamma }}^{\mathrm{2}}\mathit{\gamma ̇}\mathrm{\left(}{\mathit{\gamma }}_{\mathrm{inj}}\mathit{,}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}}\frac{{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}}{\left(\frac{{\mathit{\varphi }}_{\mathrm{e}}\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{\right)}}{{\mathit{\varphi }}_{\mathrm{e}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\right)}^{\mathit{s}\mathrm{-}\mathrm{1}}\mathit{S}\mathrm{\left(}\mathit{\gamma }\mathrm{;}{\mathit{\gamma }}_{\mathrm{inj}}\mathit{,}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(9)where s = 1 + t_a,0/t_e,0, S(x;a,b) = 1 if a ≤ x ≤ b and zero otherwise, and the function ϕ_e is given by Eq. (7) with f_a replaced by f_e. The quantity t_∗ = t_∗(γ,t) in the case where both the acceleration timescale and the magnetic field are time independent (f_a = f_B = 1) takes the analytic expression ${\mathit{t}}_{\mathrm{\ast }}\mathrm{=}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}\ln \left(\frac{{\mathit{\gamma }}_{\max }\mathit{/}{\mathit{\gamma }}_{\mathrm{inj}}\mathrm{-}\mathrm{1}}{{\mathit{\gamma }}_{\max }\mathit{/}\mathit{\gamma }\mathrm{-}\mathrm{1}}\right)\mathit{,}$(10)while in every other case it is obtained from the solution of the implicit Eq. (5). The upper limit of the EDF is obtained from Eq. (5) for t_∗ = 0 and is given by ${\mathit{\gamma }}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{inj}}{\mathit{\varphi }}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{inj}}}{{\mathit{\gamma }}_{\max }}\mathit{\psi }\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{·}$(11)In the case where all parameters are time independent, the unperturbed distribution function of the AZ is $\begin{array}{ccc}{\mathit{N}}_{\mathrm{0}}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}& \mathrm{=}& {\mathit{Q}}_{\mathrm{0}}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}{\left(\frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{inj}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\gamma }}_{\max }}\right)}^{\mathrm{1}\mathrm{-}\mathit{s}}{\mathit{\gamma }}^{\mathrm{-}\mathit{s}}{\left(\mathrm{1}\mathrm{-}\frac{\mathit{\gamma }}{{\mathit{\gamma }}_{\max }}\right)}^{\mathit{s}\mathrm{-}\mathrm{2}}\\ & & \mathrm{\times }\mathit{S}\mathrm{\left(}\mathit{\gamma }\mathrm{;}{\mathit{\gamma }}_{\mathrm{inj}}\mathit{,}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}\end{array}$(12)and its upper limit is ${\mathit{\gamma }}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\gamma }}_{\mathrm{inj}}{\mathrm{e}}^{\mathit{t}\mathit{/}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}}}{\mathrm{1}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{inj}}}{{\mathit{\gamma }}_{\max }}\mathrm{(}{\mathrm{e}}^{\mathit{t}\mathit{/}{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(13)a relation that is valid whenever f_a = f_B = 1. The time it takes for steady state to be established is a few tens t_a,0, and then the EDF is a power-law of index s up to the maximum Lorentz factor γ_max (Eq. (6)) where acceleration balances energy losses.

2.2. Radiation zone

The escaping particles enter the radiation zone where they radiate a part of their energy in the constant magnetic field B₀. Assuming that synchrotron losses dominate, we can write an equation for the evolution of their differential density n(γ,x,t) in the RZ as $\frac{\mathit{\partial n}}{\mathit{\partial t}}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial \gamma }}\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{n}}^{\mathrm{)}}\mathrm{=}\frac{\mathit{N}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}}{{\mathit{t}}_{\mathrm{esc}}}\mathit{\delta }\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{sh}}\mathrm{(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{,}$(14)where x_sh(t) = u_sht is the position of the shock and u_sh is the shock speed. The solution of Eq. (14) is easily found: $\mathit{n}\mathrm{\left(}\mathit{\gamma ,x,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{N}\mathrm{(}{\mathit{\gamma }}_{\mathrm{\ast }}\mathit{,x}\mathit{/}{\mathit{u}}_{\mathrm{sh}}\mathrm{)}}{{\mathit{u}}_{\mathrm{sh}}{\mathit{t}}_{\mathrm{esc}}\mathrm{(}\mathit{x}\mathit{/}{\mathit{u}}_{\mathrm{sh}}\mathrm{)}}{\left(\frac{{\mathit{\gamma }}_{\mathrm{\ast }}}{\mathit{\gamma }}\right)}^{\mathrm{2}}\mathit{S}\left(\mathit{x}\mathrm{;}\mathrm{0}\mathit{,}{\mathit{x}}_{\mathrm{sh}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\right)$(15)with ${\mathit{\gamma }}_{\mathrm{\ast }}\mathrm{=}{\left(\frac{\mathrm{1}}{\mathit{\gamma }}\mathrm{-}{\mathit{\alpha }}_{\mathrm{0}}\mathrm{(}\mathit{t}\mathrm{-}\mathit{x}\mathit{/}{\mathit{u}}_{\mathrm{sh}}\mathrm{)}\right)}^{-1}\mathit{.}$(16)Following KRM, we impose the additional restriction x > x_sh(t) − L on the spatial limits of Eq. (15), in order to take the finite size of the RZ, L, into account. This is better expressed through the time the shock needs to travel the RZ ${\mathit{t}}_{\mathrm{b}}\mathrm{=}\mathit{L}\mathit{/}{\mathit{u}}_{\mathrm{sh}}\mathit{.}$(17)Particles not fulfilling this restriction are assumed to escape the RZ in a region of practically zero magnetic field, so they do not contribute to the total radiation. The EDF in the RZ is then obtained by spatially integrating the function n(γ,x,t).

When all parameters are time-independent, the steady-state distribution in the RZ is a broken power law having the same energy limits as the EDF in the AZ. The power law has an index s for energies less than the breaking energy γ_br ≃ γ_maxt_a,0/t_b and (s + 1) for γ > γ_br. The break is formed because of the finite size of the source, since the particles with γ < γ_br leave the source before they cool. Radiation from these particles is ignored because of the assumption that the magnetic field declines substantially outside the source. The time it needs the EDF in the RZ to reach steady state is around t_b.

2.3. Radiation spectra

Knowing the electron distribution functions in both zones, it is easy to compute the emitted synchrotron spectrum in the source frame through the relation $\begin{array}{ccc}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& {\mathit{I}}_{\mathit{\nu }}^{\mathrm{AZ}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}{\mathit{I}}_{\mathit{\nu }}^{\mathrm{RZ}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\\ & \mathrm{=}& {\mathrm{\int }}_{\mathrm{1}}^{\mathrm{\infty }}\mathrm{d}\mathit{\gamma N}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}{\mathit{I}}_{\mathit{\gamma }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{+}{\mathrm{\int }}_{\mathrm{1}}^{\mathrm{\infty }}\mathrm{d}\mathit{\gamma }{\mathit{I}}_{\mathit{\gamma }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{xn}\mathrm{\left(}\mathit{\gamma ,x,t}\mathrm{\right)}\end{array}$(18)where I_γ(ν) is the single particle emissivity for synchrotron radiation that is given by ${\mathit{I}}_{\mathit{\gamma }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\sqrt{\mathrm{3}}{\mathit{e}}^{\mathrm{3}}\mathit{B}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathit{F}\left(\frac{\mathit{\nu }}{\mathrm{2}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{\nu }}_{\mathrm{0}}}\right)\mathit{,}$(19)with $\mathit{F}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathit{x}}^{\mathrm{2}}\left({\mathit{K}}_{\mathrm{4}\mathit{/}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}{\mathit{K}}_{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}\frac{\mathrm{3}\mathit{x}}{\mathrm{5}}\mathrm{\left(}{\mathit{K}}_{\mathrm{4}\mathit{/}\mathrm{3}}^{\mathrm{2}}\mathrm{\right(}\mathit{x}\mathrm{)}\mathrm{-}{\mathit{K}}_{\mathrm{1}\mathit{/}\mathrm{3}}^{\mathrm{2}}\mathrm{(}\mathit{x}{\mathrm{\left)}}^{\mathrm{\right)}}\right)$(20)(Crusius & Schlickeiser 1986), and the characteristic frequency is ${\mathit{\nu }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}\frac{\mathit{eB}}{{\mathit{m}}_{\mathrm{e}}\mathit{c}}\mathrm{·}$(21)In Eq. (18) one can see the two differences between our approach and the one adopted by KRM. The first is that we calculate radiation from both zones, while KRM ignore the contribution from the AZ. The second difference is that we do not treat light-travel effects inside the source. This assumption is equivalent to taking u_sh ≪ c, and in that case the shock speed is no longer a parameter of the problem. Our model thus considers the case of non-relativistic shocks. In the opposite case, light-travel effects can be important depending on the frequency (see the discussion in KRM) and then the intrinsic variations are smoothed out on the light-crossing timescale. A study of these phenomena in the internal-shock model for blazars can be found in Böttcher & Dermer (2010).

3. Flare profiles

The model that we presented in the previous section can produce a wide variety of flaring behaviors by varying one or more of its basic parameters. The temporal behavior of these is essentially another free parameter because no theory can provide this information. While a short duration impulsive change is assumed in many applications (Mastichiadis & Kirk 1997, KRM), here we adopt the Lorentz profile for the changes, which has a pulse shape and analytic expression ${\mathit{f}}_{\mathrm{L}}\mathrm{\left(}\mathit{t}\mathrm{;}{\mathit{t}}_{\mathrm{0}}\mathit{,w,n}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\frac{{\mathit{w}}^{\mathrm{2}}}{\mathrm{4}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{w}}^{\mathrm{2}}}\mathrm{·}$(22)At t = t₀ the pulse shows a maximum or a minimum, i.e. f_L(t₀;t₀,w,n) = n, depending on whether n > 1 or n < 1, while it is f_L = 1 for  | t − t₀ |  ≫ w/2. The quantity w is the full width at half maximum (or minimum, FWHM) of the pulse, since .

3.1. Change in the injection rate

A change in the injection rate can be attributed to the encounter of the shock with a higher or lower density area. The increase in the injection rate leads to an increase in the overall normalization of the EDFs, hence to the production of a flare. The case of a step-function change was examined in KRM. Here we study it for a Lorentz type change in more detail.

A few snapshots of the synchrotron spectrum when the injection rate varies according to f_q(t) = f_L(t;t₀,w,n) with w = 3t_a,0, n = 3 and t₀ ≫ t_b – so that steady state has been reached long before t = t₀ – are shown in Fig. 1. For the steady-state spectrum, the set of parameters γ_inj = 10, B = 0.5   G, t_b = 100t_a,0, Q₀ = 10⁴⁶   s^-1, and t_a,0, t_e,0 such that γ_max = 10⁶ and s = 1.8 is adopted. The maximum synchrotron frequency in the source frame, , is then a few keV. The corresponding light curves are shown in Fig. 2.

Fig. 1 Snapshots of the synchrotron spectrum in the source frame. The curve labeled “a” is the steady-state spectrum, while labels “b” to “g” are the snapshots at t = t0 + 4kta,0 for k = 2,...,7, respectively. The vertical arrows define the frequency range (νf,νl) for which the light curves of Fig. 2 are calculated.

Fig. 2 Flares that are produced by a Lorentzian change of the injection rate in the frequency range νf = 10-3νmax ≤ ν ≤ 100.6νmax = νl with step 100.4νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.

It is evident from Fig. 2 that the form of the flares is a function of frequency. To quantify this, we define the following parameters for each flare. If y(t) = I_ν(t)/I_ν(0) is the form of the light curve normalized to the steady-state luminosity, I_ν(0), then we define the time at which the flare peaks as t_pk and its peak value y_pk ≡ y(t_pk). We also define the flux doubling time during the rise of the flare, t_r, and the flux halving time during the decay, t_d, through the relations y(t_pk − t_r) = y(t_pk + t_d) = (1 + y_pk)/2. The ratio t_r/t_d is then a measure of the time symmetry of the flare, and w_fl = t_r + t_d is the FWHM of the flare. We should note here that we do not fit each flare with a specific function, but instead calculate the above quantities independently from the exact form of the flares.

In Fig. 3 we plot the quantities t_pk, t_r/t_d, w_fl, and y_pk as functions of frequency. In the first plot for ν ≲ 0.02ν_max, the time at which the flare peaks decreases with frequency, while it increases for ν ≳ 0.02ν_max. This means that for ν ≲ 0.02ν_max the low-frequency flares precede the high-frequency ones, or, in other words, the flares show a “soft lag”. The opposite happens for ν ≳ 0.02ν_max where the flares show a “hard lag”.

Fig. 3 Characteristics of the flares shown in Fig. 2.

As noted in KRM, who also find this result, soft lags can be attributed to electrons for which the acceleration timescale is much faster than the cooling timescale, therefore for all practical purposes these can be considered as pre-accelerated or injected at high energies with a ready power law (see, for example, Mastichiadis & Kirk 1997). On the other hand, hard lags are observed at higher frequencies where the two timescales are comparable, therefore the observer sees the wave of freshly accelerated particles moving to high energies.

Deeper insight can be gained by considering the exact changes that the electron populations experience. The distribution function of the AZ can be written in this case as $\mathit{N}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{0}}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}{\mathit{f}}_{\mathrm{L}}\mathrm{\left(}{\mathit{t}}_{\mathrm{\ast }}\mathrm{;}{\mathit{t}}_{\mathrm{0}}\mathit{,w,n}\mathrm{\right)}$(23)with t_∗ given by Eq. (10). This result indicates that the initial pulse propagates in the EDF from low to high energies, since its peak is located at γ = γ₁(t − t₀), as can be easily found by setting t_∗ = t₀. Moreover, this moving pulse preserves the shape of the initial pulse, and thus it is symmetric around its center and has an amplitude equal to n. If we only consider radiation from the AZ, we expect to get roughly symmetric flares of the same amplitude that exhibit hard lags.

The situation is more complicated in the radiation zone, since the EDF there results from the integration of Eq. (15) along the extent of the source. The hard-lag behavior of the AZ is also seen here, but only in high energies, where the contribution from the particles residing just before the shock dominates in the integral. At lower energies the perturbation propagates in the opposite direction, as the intrinsic variations of N(γ,t) are smoothed out by the spatial integration, and only the effect of the fading pulse in the AZ stands out.

In the second and third plots of Fig. 3, one sees that the width of the flares decreases with frequency, while the ratio t_r/t_d increases. This means that the flares become narrower and more symmetric as the frequency increases. If the rise time of the flares is connected with the acceleration timescale and the decay with the cooling timescale, then the ratio t_r/t_d will be proportional to $\frac{{\mathit{t}}_{\mathrm{acc}}}{{\mathit{t}}_{\mathrm{cool}}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{a}\mathit{,}\mathrm{0}}}{\mathrm{(}{\mathit{\alpha }}_{\mathrm{0}}\mathit{\gamma }{\mathrm{)}}^{-1}}\mathrm{=}\frac{\mathit{\gamma }}{{\mathit{\gamma }}_{\max }}\mathrm{\simeq }{\left(\frac{\mathit{\nu }}{{\mathit{\nu }}_{\max }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(24)The plot of this curve in the t_r/t_d diagram reveals that this is not the case, making the relation of t_r, t_d to the acceleration and cooling timescales more complicated.

Finally, in the last plot of Fig. 3 one sees that the amplitude of the flares increases with frequency until it saturates to a constant value that is somewhat lower than the amplitude of the electron variation. The level of saturation depends on the length of the pulse, because the longer the duration of the burst becomes, the closer the EDF gets to a steady state.

3.2. Change in the escape timescale

As a next case we examine a perturbation in the escape timescale. More precisely we assume that some physical process momentarily impedes particle escape, i.e. decreases the rate at which particles escape from the acceleration region. An inspection of the EDF (Eq. (12)) shows that an increase in the escape timescale leads to a decrease in the electron index. This will produce a flare, the form of which is shown in Fig. 4 at various frequencies. We assume that the change of the escape timescale again follows a Lorentzian variation, i.e. f_e(t) = f_L(t;t₀,w,n) with w = 3t_a,0, n = 4/3, and t₀ ≫ t_b, while the parameters of the steady-state spectrum are the same as before.

Fig. 4 Flares that are produced by a Lorentzian change in the escape timescale in the frequency range νf = 10-3νmax ≤ ν ≤ 100.6νmax = νl with step 100.4νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.

The specific features of the flares are shown in Fig. 5. As before, the flares exhibit a soft lag in low frequencies and a hard lag close to the maximum frequency. The frequency range in which hard lags are observed is much less than in the previous case, and also the flares now peak closer to the peak of the perturbation pulse, i.e. closer to t = t₀.

The shape of the flares is asymmetric in time at low frequencies, with a rise time that is shorter than the decay one (second plot of Fig. 5); however, this becomes symmetric, i.e. t_r ≃ t_d, at higher frequencies. The rising part of the curve indicates that t_r/t_d is related to t_acc/t_cool, and thus the rise of the flares is controlled by the acceleration and the decay by the cooling timescale.

The width of the low-frequency flares decreases as in the previous case, but at higher frequencies it increases. We note, however, that w_fl ≫ w, i.e. even a narrow pulse in the escape timescale produces a rather broad X-ray flare. The amplitude of the flares increases with frequency as in the previous case, but now it does not saturate to a maximum value.

The EDF of the AZ, Eq. (9), reads in this case $\begin{array}{ccc}\mathit{N}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}& \mathrm{=}& {\mathit{N}}_{\mathrm{0}}\mathrm{\left(}\mathit{\gamma ,t}\mathrm{\right)}\exp [\frac{\mathit{w}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathit{s}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{2}\sqrt{\mathit{n}}}\\ & & \mathrm{\times }\left({\tan }^{-1}\left(\mathrm{2}\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{\mathit{w}\sqrt{\mathit{n}}}\right)\mathrm{-}{\tan }^{-1}\left(\mathrm{2}\frac{{\mathit{t}}_{\mathrm{\ast }}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{\mathit{w}\sqrt{\mathit{n}}}\right)\right)]\end{array}$(25)where t_∗ is given by Eq. (10). The perturbation in the EDF propagates again from low to high energies, but is no longer pulse-shaped and is much broader. Moreover, the perturbation moves twice as fast now, since the location of its peak is given by γ = γ₁(2(t − t₀)). One notices that, qualitatively speaking, the trends in the present case are similar to the previous one. There is a simple explanation for this. The perturbation introduced causes the escape timescale to increase first and then decrease back to its original value. This in turn causes the electron spectral index to flatten and then steepen again. Therefore, the number of accelerated particles increases and then decreases to its unperturbed value. Thus this case can be considered as equivalent to the previous one which treated variations in the injection of particles at low energies. In the present case, the flattening of the EDF is partly compensated for by the particles decreasing in the RZ at the same time (see Eq. (15)).

Fig. 5 Characteristics of the flares shown in Fig. 4.

3.3. Change in the acceleration timescale

As a next example we investigate changes in the acceleration timescale and, more specifically, the variations induced when this decreases. This change has two effects on the EDFs and, therefore, on the produced spectrum. First, an inspection of the EDF (Eq. (12)) reveals that if the acceleration timescale decreases, then the electron index becomes flatter. Second, the maximum electron energy will increase – see Eq. (6). The flares produced for a Lorentzian change with w = 3t_a,0 and n = 3/4 are shown in Fig. 6.

As frequency increases, the flare shapes become more asymmetric and peak at earlier times. In this case we only get soft-lag flares. This can easily be explained since the increase of γ_max in the EDFs leads to the simultaneous increase in the maximum emitted frequency in the spectrum. This produces a large flare at high frequencies, since during pre flare this part of the spectrum was in the exponentially decaying synchrotron regime.

Fig. 6 Flares that are produced by a Lorentzian change in the acceleration timescale in the frequency range νf = 10-2.1νmax ≤ ν ≤ 100.6νmax = νl with step 100.3νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.

3.4. Change in the magnetic field

As a last case we examine variations in the magnetic field strength, which are assumed to occur in the AZ only. A decrease in the magnetic field leads to an increase in the maximum electron energy. The maximum photon frequency also increases, even if B decreases, and this produces a flare. The form of the flares produced for a Lorentzian change in the magnetic field with w = 3t_a,0 and n = 0.8 is shown in Fig. 7. A new feature that appears is a small decrease in the flux preceding the low-frequency flares. This is because the initial drop in the magnetic field leads to the decrease in the total power radiated by the electrons, and this is seen in low frequencies since the high-frequency spectrum rises due to the increase of γ_max, which compensates for this reduction.

The flares exhibit only soft lags for a similar reason to that of the previous case. The flares also become narrower and more symmetric with frequency. In total, this case has many similarities with the change in the acceleration timescale, since in both cases the maximum electron energy follows the same temporal behavior.

Fig. 7 Flares that are produced by a Lorentzian change in the magnetic field in the frequency range νf = 10-1.9νmax ≤ ν ≤ 100.2νmax = νl with step 100.3νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.

4. Application to 1ES 1218+304

1ES 1218+304 is a high-frequency peaked BL Lac object at a redshift z = 0.182. It was observed in May 2006 with Suzaku, and the analysis of the observations (Sato et al. 2008) revealed that the flares in the various X-ray energy bands exhibit a hard lag. Moreover, the flare shapes are asymmetric in time, and their amplitude increases with photon energy. These features resemble the ones found in the previous section for the change in the injection rate or in the escape timescale.

The radio to X-rays spectrum of 1ES 1218+304 is shown in Fig. 8, together with the steady-state spectrum of our model as given in Eq. (18). The transformation to the observer’s frame is made through the relations $\begin{array}{ccc}& & {\mathit{F}}_{\mathit{\nu }}^{\mathrm{obs}}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\frac{{\mathit{I}}_{\mathit{\nu }}}{\mathrm{4}\mathit{\pi }{\mathit{d}}_{\mathrm{L}}^{\mathrm{2}}}\\ & & {\mathit{\nu }}_{\mathrm{obs}}\mathrm{=}\frac{\mathit{\delta }}{\mathrm{1}\mathrm{+}\mathit{z}}\mathit{\nu }\end{array}$ where δ is the Doppler factor and d_L = 880   Mpc the luminosity distance of 1ES 1218+304, for the flat universe with cosmological parameters Ω_Λ = 0.7, Ω_m = 0.3 and H₀ = 70   kms^-1   Mpc^-1.

The maximum photon frequency in the spectrum is   ×  10¹⁸   Hz, leading to the relation $\mathit{\delta B}{\mathit{\gamma }}_{\max }^{\mathrm{2}}\mathrm{\simeq }{\mathrm{10}}^{\mathrm{12}} \mathrm{G}\mathit{.}$(28)We choose a random combination of these parameters so that Eq. (28) is satisfied and also take γ_inj = 10, t_b = 20t_a,0, and s = 1.8. The value of t_b is dictated by the frequency where the spectrum appears to break,   ×  10¹⁵   Hz, while the value of the electron index is chosen so that the photon index in radio frequencies is α_r ≃ 0.6. The parameter γ_inj has a small effect on the normalization of the spectrum, so its value is unimportant.

Fig. 8 Radio to X-rays spectrum of 1ES 1218+304 (starred symbols, data taken from Sato et al. 2008) together with the steady-state theoretical spectrum (solid line) for the parameters given in text. The contributions from the acceleration zone (dashed line) and from the radiation zone (dotted line) are also shown. The optical data can be attributed to the host galaxy of 1ES 1218+304 (Rüger et al. 2010), so are not included in the fit.

The remaining parameters can be further constrained by the variability timescale of 1ES 1218+304. From the values of t_pk(ε) given in Sato et al. for the time at which the flares peak, we deduce that the variation responsible for producing these flares is the change in the injection rate. This is because the energy range in which t_pk is rising is quite broad and rules out the case of the escape timescale change. As in Sato et al. we define the amount of hard lag simply by the difference in the peak times of the flares in the various energy bands from the peak time of the flare in the highest energy band, namely ${\mathit{\tau }}_{\mathrm{hard}}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{=}{\mathit{t}}_{\mathrm{pk}}\mathrm{\left(}{\mathit{\epsilon }}_{\max }\mathrm{\right)}\mathrm{-}{\mathit{t}}_{\mathrm{pk}}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathit{,}$(29)where is the logarithmic mean energy of the 5−10   keV energy band. The theoretical curve for τ_hard(ε) for the Lorentzian change of the injection rate with parameters w = 1.7t_a,0 and n = 1.4 is shown in the upper panel of Fig. 9 along with the observational data. From the fit, we estimate the acceleration timescale t_a,0 = 1.9  ×  10⁵   δ₁₀   s, where δ₁₀ = δ/10. When combined with Eq. (28), this relation enables us to express the magnetic field and maximum electron energy of 1ES 1218+304 as and γ_max = 1.3  ×  .

Since the Doppler factor cannot be determined, we choose δ₁₀ = 1, and the unperturbed value of the injection rate is then Q₀ = 3.2  ×  10⁴⁵   s^-1. With these parameters we can estimate the size of the source L = 1.1  ×  10¹⁶(u_sh/0.1c)   cm and the energy content of the source in electrons as ℰ_e = 3  ×  10⁴⁸   erg.

In the lower panel of Fig. 9 we plot the theoretical curve of the asymmetry of the flares, t_r/t_d, together with the observational data. Our definition of t_r and t_d differs from the one in Sato et al., but the ratio of the two quantities is the same in both cases. One sees that we reproduce the general trend of the data quite well, but we overestimate the ratio t_r/t_d in the energy interval 1−3   keV.

Fig. 9 Amount of hard lag and the ratio tr/td of 1ES 1218+304 in the various X-ray energy bands, together with the corresponding theoretical curves.

5. Conclusions

In the present paper we have examined the expected variability of blazars in the context of the two-zone acceleration model. By varying four of the key parameters of the model, we found that flaring can be induced in all cases. In all cases the simulated flares exhibit soft lags; only in two cases, i.e. when varying the injection rate of particles or their escape rate, one can find frequency regimes where hard lags appear. Thus, we find that this scheme cannot produce hard lags at much lower frequencies than the maximum synchrotron frequency. Another conclusion is that the flares tend to mimic the electron variations at high frequencies. As a rule at low frequencies, flares tend to be asymmetric, with rise times much faster than decay ones.

We were also able to make an acceptable fit to the hard-lag flare of 1ES 1218+304 (Sato et al. 2008) by varying the rate of particles injected into the acceleration process. We note that the fit could be improved if we were to vary two parameters instead of one, but this lies outside the scope of the present paper. The fitting we obtained gives parameters well within the ones accepted from blazar modeling. The acceleration timescale deduced (t_a,0 ≃ 2  ×  10⁵   δ₁₀ s) might be long for first-order Fermi acceleration (Tammi & Duffy 2009); however, we have to emphasize that from the setup of the model, both acceleration and escape timescales were assumed to be independent of the energy γ and the magnetic field strength B, therefore we did not try to connect them directly with theories of particle acceleration (Drury 1983; Blandford & Eichler 1987).

We restricted our analysis to the cases where synchrotron losses dominate and did not consider the inverse Compton scattering as a possible electron energy loss mechanism. This can be justified for sources where the energy density due to magnetic fields exceeds the one due to photons. In this case the variability patterns due to inverse Compton scattering will, in general, follow the ones due to synchrotron (Kirk & Mastichiadis 1999). On the other hand, an example of a model dealing with X- and gamma-ray variability where inverse Compton losses might be important can be found in Mastichiadis & Moraitis (2008), however, this is an one-zone acceleration model.

Despite its limitations, our model is able to produce a wide variety of flaring behaviors that could, in principle, be tested against the growing X-ray data on blazar variability. We note that recently, Garson et al. (2010) have also attributed the properties of X-ray flares to intrinsic changes in the acceleration process.

Acknowledgments

The authors would like to thank the anonymous referee whose comments helped improve the present paper. K.M. acknowledges financial support from the Greek State Scholarships Foundation (IKY).

References

All Figures

	Fig. 1 Snapshots of the synchrotron spectrum in the source frame. The curve labeled “a” is the steady-state spectrum, while labels “b” to “g” are the snapshots at t = t0 + 4kta,0 for k = 2,...,7, respectively. The vertical arrows define the frequency range (νf,νl) for which the light curves of Fig. 2 are calculated.
In the text

	Fig. 2 Flares that are produced by a Lorentzian change of the injection rate in the frequency range νf = 10-3νmax ≤ ν ≤ 100.6νmax = νl with step 100.4νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.
In the text

	Fig. 3 Characteristics of the flares shown in Fig. 2.
In the text

	Fig. 4 Flares that are produced by a Lorentzian change in the escape timescale in the frequency range νf = 10-3νmax ≤ ν ≤ 100.6νmax = νl with step 100.4νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.
In the text

	Fig. 5 Characteristics of the flares shown in Fig. 4.
In the text

	Fig. 6 Flares that are produced by a Lorentzian change in the acceleration timescale in the frequency range νf = 10-2.1νmax ≤ ν ≤ 100.6νmax = νl with step 100.3νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.
In the text

	Fig. 7 Flares that are produced by a Lorentzian change in the magnetic field in the frequency range νf = 10-1.9νmax ≤ ν ≤ 100.2νmax = νl with step 100.3νmax. Frequency increases from bottom to top, and for clarity the first and the last curves are labeled “νf” and “νl” respectively.
In the text

Fig. 8 Radio to X-rays spectrum of 1ES 1218+304 (starred symbols, data taken from Sato et al. 2008) together with the steady-state theoretical spectrum (solid line) for the parameters given in text. The contributions from the acceleration zone (dashed line) and from the radiation zone (dotted line) are also shown. The optical data can be attributed to the host galaxy of 1ES 1218+304 (Rüger et al. 2010), so are not included in the fit.

In the text

	Fig. 9 Amount of hard lag and the ratio tr/td of 1ES 1218+304 in the various X-ray energy bands, together with the corresponding theoretical curves.
In the text