Abstract

We examine the applicability of the stochastic electron acceleration to two high synchrotron peaked blazars, Mrk 421 and Mrk 501, assuming synchrotron self-Compton emission of gamma-rays. Our model considers an emitting region moving at relativistic speed, where non-thermal electrons are accelerated and attain a steady-state energy spectrum together with the photons they emit. The kinetic equations of the electrons and photons are solved numerically, given a stationary wavenumber spectrum of the magnetohydrodynamic (MHD) disturbances, which are responsible for the electron acceleration and escape. Our simple formulation appears to reproduce the two well-sampled, long-term averaged photon spectra. In order to fit the model to the emission component from the radio to the X-ray bands, we need both a steeper wave spectral index than the Kolmogorov spectrum and efficient particle escape. Although the model provides a natural explanation for the high-energy cutoff of the electron energy distribution, the derived physical parameters raise a problem with an energy budget if the MHD waves with the Alfvén velocity are assumed to be the acceleration agent.

1 INTRODUCTION

Blazars, one of the classes of active galactic nuclei (AGNs), have highly collimated jets with relativistic speed pointing close to our line of sight (Urry & Padovani 1995; Marscher 2010). Characteristics of X-ray variability of TeV emitting blazars (Kataoka 2008) with time-scale of ≲ 1 d measured in the observer frame suggest that the emitting region of the AGN jets that we see in blazars (so-called blazar region) has the length-scale of its structure extending up to ∼1016 cm in its comoving frame when it is estimated in this energy band. Since the electromagnetic radiation from that region is observed dominantly by virtue of the relativistic beaming effect, blazars are good candidates to investigate the energy conversion process of the AGN jets. Photon spectral energy distributions (SEDs) of blazars extend from the radio to the gamma-ray bands (Fossati et al. 1998) and are dominated by non-thermal emission. This demonstrates that non-thermal particles are accelerated in the blazar region.

Blazars are categorized based on their SED peak frequency of the component that extends from the radio up to the X-ray bands in the νFν representation (Abdo et al. 2010). For high synchrotron peaked blazars (HSPs), whose peak frequency is in the UV or X-ray band, the emission spectrum has been primarily interpreted as the synchrotron self-Compton (SSC) process: the synchrotron emission by the non-thermal relativistic electrons (and positrons) below the X-ray band and the inverse Compton (IC) scattering of these electrons off the synchrotron photons for the gamma-ray band (Maraschi, Ghisellini & Celotti 1994; Sikora & Madejski 2001; Dermer & Lott 2012).

In modelling observed SEDs, most of the studies assume either a broken power-law form with arbitrary indices for the energy distribution of the radiating non-thermal electrons (e.g. Tavecchio et al. 2001) or the energy distribution calculated by solving the kinetic equation under injection of a power-law spectrum, energy losses, and escape (or adiabatic cooling) near and downstream of a shock (e.g. Mastichiadis & Kirk 1997; Kusunose, Takahara & Li 2000). Here, we mention three implications on the shape of the electron energy distribution (EED) from such studies. (i) When a broken power-law EED is applied in the one-zone approximation, the change of the spectral index below and above the break energy tends to be larger than 1 (Tagliaferri et al. 2008; Tavecchio et al. 2010; Zhang et al. 2012). This differs from a prediction of the simple one-zone model in which the break is caused by the competition between radiation cooling and escape. (ii) It is known that the broken power-law form of an EED is not always adequate to reproduce observed well-sampled SEDs. Purely phenomenologically, adding an additional break to the EED (i.e. two breaks in total) can improve the fit (Kataoka et al. 1999; Aharonian et al. 2009; Abdo et al. 2011a,b). The physical origin of these breaks (or bending of the EED) is not clear. (iii) If a population of non-thermal electrons form an SED from the radio band (∼1010–1012 Hz) to a higher frequency (e.g. the synchrotron peak), a hard EED is required. Kataoka et al. (2000) reproduced the SED of this frequency range with the index of the EED of −1.35 for PKS 2155−304, and Kino, Takahara & Kusunose (2002) −1.4, −1.6, and −1.8 for PKS 2155−304, Mrk 421, and Mrk 501, respectively. It has been shown that in principle such hard electron indices can be attained by the first-order Fermi acceleration in the test particle approximation, although the results are sensitive to conditions near the shock (e.g. inclination angle of the background magnetic field, the shock velocity, and the nature of particle scattering upstream and downstream of the shock; e.g. Bednarz & Ostrowski 1996; Niemiec & Ostrowski 2004; Lemoine & Revenu 2006). Note that since it is uncertain whether an EED of a single population is reflected in an SED from the radio band, it is usual not to intend to fit this range (e.g. Krawczynski et al. 2004; Kusunose & Takahara 2008; Rani et al. 2013).

There is also an implication on the acceleration time-scale of electrons. (iv) If the balance between synchrotron cooling and acceleration time determines the high-energy cutoff energy of an EED, the acceleration time of electrons with Lorentz factor γ, normalized by the gyroperiod of the same electrons, is given as |$\tau _{\rm acc}(\gamma _{\rm cut}) = 3e/(\sigma _{\rm T}B \gamma _{\rm cut}^2)$|⁠, where e is the charge of an electron, σT the Thomson cross-section, B the magnetic field, and γcut the electron Lorentz factor at the cutoff. Adopting B = 0.1 G and γcut = 105, which are derived for emission regions of some TeV blazars with the conventional one-zone SSC model (Kino et al. 2002), one finds τacc ∼ 106 (Inoue & Takahara 1996). If it is assumed that only the first-order Fermi mechanism at a single shock works in the blazar region, both such a long time to accelerate the cutoff electrons and the hard spectra (iii) imply that conditions around the shock would presumably be very restricted, when extrapolated from the test particle simulations (e.g. Bednarz & Ostrowski 1996; Niemiec & Ostrowski 2004; Lemoine & Revenu 2006). This issue may give an important implication on electron acceleration processes, though the spatial dependence of the magnetic field and the IC dominance on electron cooling can change the value of τacccut) to some extent and that one could not evaluate τacccut) well when only a weak constraint on the cutoff energy is obtained (e.g. Abdo et al. 2011a,b).

Instead of interpreting the implications for the non-thermal EED mentioned above in terms of the first-order Fermi mechanism, stochastic acceleration (SA), which is described as diffusion in momentum space, can be invoked to be qualitatively consistent with them. Because the models including the SA, energy loss processes, and escape from a finite acceleration region can form various shapes of an EED (Becker, Le & Dermer 2006; Stawarz & Petrosian 2008), the SA has the possibility to account for the above features (i)–(iv) simultaneously. In addition, since the first-order Fermi mechanism necessitates scattering of particles, in the situation where it really acts slowly as mentioned in (iv), it may be accompanied by the SA effectively.

The emission models including the SA have been applied to recent observations of blazars in many papers (e.g. Katarzyński et al. 2006; Massaro et al. 2006; Ushio et al. 2010; Weidinger, Rüger & Spanier 2010; Lefa, Aharonian & Rieger 2011; Tramacere, Massaro & Taylor 2011; Yan, Zeng & Zhang 2012; Cao & Wang 2013). They often focus on either limited portions of observed spectra (e.g. near the cutoffs of the synchrotron and the IC components) or relatively hard flaring spectra, and in some cases add the SA to power-law acceleration models or employ multizone models. We also made a time-dependent simulation with a radially structured jet model including the SA in Asano et al. (2014). In general, hard flaring spectra are compatible with a feature of the SA. A possibility of the SA for a slowly variable component is implied with the observation of Mrk 421 in the X-ray band by Ushio et al. (2009).

In this paper, we focus on the issues on the shape of EED, which have been raised mainly from one-zone steady-state models which assume a (broken) power-law EED. In order to investigate the validity of the SA to these issues, we apply the one-zone steady-state model employing a simple form of the SA and SSC emission to well-sampled long-term averaged emission of blazars. The observations compared with the model are the multiwavelength campaigns on Mrk 501 and 421, reported in Abdo et al. (2011a,b).

The paper is organized as follows. We describe our model in Section 2 and the procedure of application to the observations in Section 3. The results are shown in Section 4, followed by the discussion in Section 5. Section 6 is a summary.

2 MODEL

We assume that a finite region where electrons are continuously injected into the SA (hereafter referred to as blob) is formed intermittently, propagates on a jet, and then reaches a steady state, forming the blazar region. We suppose that one of the blobs dominates an observed overall SED at a given time on average (Sikora et al. 1997), and then calculate the synchrotron and the SSC emission from the blob by solving the kinetic equations of non-thermal electrons and photons inside the blob.

The blob moves relativistically with Lorentz factor Γ in the laboratory frame at an angle Γ−1 to our line of sight. The Doppler factor of the blob is approximated as Γ. The shape of the blob is taken as cylindrical, characterized by the size longitudinal to the direction of its motion and that normal to it. Denoting them as R and R, respectively, in the blob comoving frame, one can write the volume of the blob as
(1)
The two length-scales relate to the formation of the blob in the jet. Neither acceleration/deceleration nor adiabatic change of the blob is considered.
For simplicity, we assume that the accelerated electrons are always isotropically and homogeneously distributed in the blob frame and that the time evolution of their phase-space distribution is described as the pure momentum diffusion equation (Hall & Sturrock 1967; Blandford & Eichler 1987; Ostrowski & Siemieniec-Oziębło 1997). Taking some other effects into consideration, we calculate the steady-state EED by solving the kinetic equation in the blob frame given by (Stawarz & Petrosian 2008)
(2)
where the accelerated electrons are relativistic enough that their momentum is approximated as γmec, me is the electron mass, c the speed of light, ne(γ) the number density of the accelerated electrons per unit γ, De the diffusion coefficient of the electrons in γ space due to the SA, |$\dot{\gamma }_{\rm rad}$| the rate of the energy change due to the synchrotron and SSC processes, tesc, e the characteristic time for the electrons to escape from the blob, and qinj(γ) the electron injection rate into the acceleration process per unit volume per unit γ. All the quantities in equations (2) and (3) below are defined in the blob comoving frame. The energy distribution of the photons emitted by the non-thermal electrons is isotropic and homogeneous in the blob frame. The photon kinetic equation is given by (Li & Kusunose 2000)
(3)
where ϵ is the photon energy in unit of mec2, nph(ϵ) the number density of photons per unit ϵ, the first three terms on the right-hand side are the photon production and absorption rate by the electrons per unit volume per unit ϵ, and tesc, ph the photon escape time from the blob. The particular form of each term in equations (2) and (3) is described in the following.
The diffusion of the electrons in γ space can be caused by the scattering off magnetohydrodynamic (MHD) disturbances. Assuming that the disturbances are excited uniformly throughout the blob and are a superposition of small-amplitude waves with the phase speed (in unit of c) of βw = vw/c, we can estimate the diffusion coefficient De for fast electrons as (Longair 1992; Ostrowski & Siemieniec-Oziębło 1997; Schlickeiser 2002)
(4)
where v( ≫ vw) is the electron speed (≈c at present), rg = rg(γ) = γmec2/eB the gyroradius of the electrons, B the root mean square of the magnetic field randomly oriented at scales larger than kinetic ones, δBeff the disturbed magnetic field effectively contributing to the scattering of the electrons with Lorentz factor γ with δBeff/B being sufficiently less than unity, and we consider that the gyroresonance interaction between the disturbances and the electrons at the wavenumber |$k \sim r_{\rm g}^{-1}$|⁠. The factor (rg/v)−1Beff/B)2 represents the pitch-angle diffusion coefficient under such situation (Wentzel 1974). In order to accelerate electrons to higher energies, a wave spectrum distributed in a wide k range is needed. Assuming that the wave spectrum is stationary with a power-law spectrum |$\delta B_k^2 \propto k^{-q}$|⁠, containing larger amount of energy at smaller wavenumber, that is, q > 1, we write |$\delta B_{\rm eff}^2$| as
(5)
where λmax is the longest wavelength of the wave spectrum and ζ the ratio of the energy density of the disturbed magnetic field to the mean magnetic field at the scale of λmax. By definition, ζ is less than unity. From equations (4) and (5), one sees De ∝ γq.

The escape term describes spatial transport of the particles in the blob. For simplicity, we take the average length-scale the particles propagate until escape, denoted by Resc in the blob frame, as |${\rm min}\lbrace R_\perp , R_\Vert \rbrace$|⁠. We adopt the spatial diffusion time |$t_{\rm esc,e}= t_{\rm esc,e}(\gamma ) = R_{\rm esc}^2/\kappa \propto \gamma ^{q-2}$| as the electron escape time, where κ = lc/3 = rgc(BBeff)2/9 is the spatial diffusion coefficient along the magnetic field, and l ∝ γ2 − q the mean-free path of the electrons (Blandford & Eichler 1987; Longair 1992). Since treating the spatial transport as diffusion is improper for electrons with l ≳ Resc (in particular for higher energy electrons if q < 2), we set tesc, e = Resc/c for electrons with l > Resc in our calculation. As the photon escape time, we adopt the light-crossing time Resc/c.

The rate of the energy change of the accelerated electrons due to radiation can be written as |$\dot{\gamma }_{\rm rad} = \dot{\gamma }_{\rm syn} + \dot{\gamma }_{\rm IC}$|⁠, where |$\dot{\gamma }_{\rm syn}$| and |$\dot{\gamma }_{\rm IC}$| are the rate due to the synchrotron and the IC emission, respectively. The corresponding photon source terms in equation (3) are |$\dot{n}_{\rm syn}$| and |$\dot{n}_{\rm IC}$|⁠, respectively. Target photons of the IC scattering are only the photons originally emitted through the synchrotron process by the non-thermal electrons. Synchrotron self-absorption (SSA), which is important for the radio band spectra in the observer frame, is taken into consideration only as a sink term of photons, denoted by |$\dot{n}_{\rm SSA}$|⁠. Electron–positron pair production/annihilation is not important for the parameter values we will choose and is also neglected.

The synchrotron loss rate |$\dot{\gamma }_{\rm syn}$|⁠, the corresponding photon production rate |$\dot{n}_{\rm syn}(\epsilon )$|⁠, and the rate of SSA |$\dot{n}_{\rm SSA}(\epsilon ) = - c \, \alpha _{\rm SSA} \, n_{\rm ph}(\epsilon )$| are calculated following Rybicki & Lightman (1979) and utilizing a numerical approximation given by Finke, Dermer & Böttcher (2008), where αSSA is the SSA coefficient. IC photon production rate |$\dot{n}_{\rm IC}(\epsilon )$| (Li & Kusunose 2000) and the corresponding energy change rate of the electrons |$\dot{\gamma }_{\rm IC}$| are calculated with an approximate form of the IC scattering rate, following Jones (1968) and Blumenthal & Gould (1970).

Without treating thermal population of electrons, we just inject electrons to the acceleration process monoenergetically at a constant rate, qinj(γ) ∝ δ(γ − γinj), where δ is the delta function, and γinj the injection Lorentz factor fixed at 10 in this paper. Denoting the total electron injection luminosity as Linj, we relate qinj with Linjinjmec2V.

An observed SED of photons radiated isotropically in the blob rest frame is calculated with the average photon escape rate ϵnph(ϵ)V/tesc, ph, in the νFν representation, by
(6)
where dL is the luminosity distance, computed assuming the cosmological parameters of H0 = 71 km s−1Mpc−1, Ωm0 = 0.27, and ΩΛ0 = 0.73. Photon energy in the observer frame relates to ϵ as hν = Γϵmec2/(1 + z), where h is the Planck constant, and z the redshift. Steady-state energy spectra of the electrons and photons are computed numerically by solving equations (2) and (3) (Press et al. 1992; Park & Petrosian 1996). We tested our numerical code with some analytical solutions (Stawarz & Petrosian 2008; Mertsch 2011) and found that they are in good agreement.

3 FITTING PROCEDURE

To make comparison between observations and the model, we give eight physical parameters: Γ, B, Resc, V, Linj, βw, q, and |$\zeta \lambda _{\rm max}^{1-q}$|⁠, where the last one is included in equation (5). Noting that λmax cannot be constrained from photon observations, we do not consider the value of ζ or λmax separately.

The steady-state EED, except for its normalization, is determined by the balance among acceleration, cooling, and escape; when considering only a synchrotron component, one needs to pay attention to only three quantities which consist of some of the above parameters: q, synchrotron cooling efficiency ηsyn = ηsyn(γ) = tacc/tsyn, and electron escape efficiency ηesc = ηesc(γ) = tacc/tesc, e, where tacc = tacc(γ) = γ2/De ∝ γ2 − q, tesc, e = tesc, e(γ) ∝ γq − 2, and |$t_{\rm syn}= \gamma /\dot{\gamma }_{\rm syn} \propto \gamma ^{-1}$| (Stawarz & Petrosian 2008).

From equation (2), in the energy range where the radiation cooling can be neglected (ηsyn ≪ 1), the acceleration forms a power-law EED of ne(γ) ∝ γ1 − q at the steady state in the energy range of ηesc ≪ 1, while in the range of ηesc ≳ 1, the escape gradually softens the EED and forms a non-power-law one if q < 2. Such effect of particle escape has the potential to be a natural explanation of the shape of EEDs implied by observations mentioned in Section 1 (i–iii). Since the synchrotron component extends from the radio to the X-ray bands for HSPs, in the one-zone model one needs to choose the parameter region making electrons with γ sufficiently lower than γsyn satisfy ηesc > 1, where γsyn is the Lorentz factor defined by ηsyn = 1. Namely, we consider the case that the electron escape effect yields a curved EED. Note that the case of q ≥ 2 is unsuitable for reproducing the well-sampled, gradually softening synchrotron spectra because from equation (2) a power-law EED with the index of d ln ne(γ)/d ln γ = 1/2 − (9/4 + ηesc)1/2 is formed in the energy range of ηsyn ≪ 1 for q = 2 (Stawarz & Petrosian 2008), where ηesc is constant, and because for q > 2 higher-energy electrons have shorter mean-free path, in other words, shorter tacc and longer tesc, e. Hence, we consider the range of 1 < q < 2 coupled with one of the assumptions accompanying equation (5).

We show in Fig. 1 with black curves the variation of the synchrotron spectrum with the degree of the electron escape efficiency. The black thick solid curve fits the observed spectrum of Mrk 421 with q = 1.85, ηesc(γ), and ηsyn(γ) adopted in the left-hand panel of Fig. 2 and displayed in the upper panel of Fig. 3. The plotted observational data are the synchrotron component of Mrk 421 (see Abdo et al. 2011b for details). The black thin dashed (dotted) curve represents the case of higher (lower) ηesc(γ) by a factor of 2, compared to the best-fitting curve. The red thin solid curve in the same figure represents the case of q = 5/3 with the value of ηesc(γ) for the electrons responsible for the synchrotron flux at ≈1014 Hz fixed. This curve (q = 5/3) has stronger γ dependence of the escape efficiency than the best-fitting curve (q = 1.85) because ηesc ∝ γ2(2 − q), so that the EED for the former makes a narrower spectral peak. As can be recognized from the black thin dashed and dotted curves, if one makes the escape effect less efficient for the red thin solid curve, the peak frequency becomes higher (up to which the synchrotron cooling becomes effective), and simultaneously the EED and so the SED become harder, and vice versa. This means that the value of q needs to be closer to 2, that is, weaker γ dependence of ηesc(γ) than at least the case of q = 5/3 in order to obtain SEDs broad enough to fit with the observed ones. The adjustment of both the degree and the γ dependence of the escape efficiency is needed to explain observed SEDs within this simple one-zone model.

Synchrotron spectra calculated for the different value of q, the different degree of escape efficiency ηesc(γ), and that of synchrotron cooling efficiency ηsyn(γ). The black thick solid curve is for the same value of q, ηesc, and ηsyn as adopted in the left-hand panel of Fig. 2. The black thin dashed (dotted) curve represents the case of higher (lower) ηesc by a factor of 2, compared to the best-fitting one. The blue thick dashed (dotted) curve represents the case of higher (lower) ηsyn by a factor of 2, compared to the best-fitting one. The red thin solid curve represents the case of q = 5/3, which is varied from q = 1.85 adopted for the best-fitting one. The plotted observational data are the synchrotron component of Mrk 421 from Abdo et al. (2011b) and NED.
Figure 1.

Synchrotron spectra calculated for the different value of q, the different degree of escape efficiency ηesc(γ), and that of synchrotron cooling efficiency ηsyn(γ). The black thick solid curve is for the same value of q, ηesc, and ηsyn as adopted in the left-hand panel of Fig. 2. The black thin dashed (dotted) curve represents the case of higher (lower) ηesc by a factor of 2, compared to the best-fitting one. The blue thick dashed (dotted) curve represents the case of higher (lower) ηsyn by a factor of 2, compared to the best-fitting one. The red thin solid curve represents the case of q = 5/3, which is varied from q = 1.85 adopted for the best-fitting one. The plotted observational data are the synchrotron component of Mrk 421 from Abdo et al. (2011b) and NED.

Comparison between the model SEDs (solid curve) and the observational ones (blue circle). The left-hand panel is for Mrk 421, and the right-hand panel is for Mrk 501, where the observational data are taken from Abdo et al. (2011b) and Abdo et al. (2011a), respectively. Non-simultaneous archival data taken from NED are plotted with the grey squares below 1015 Hz, for reference.
Figure 2.

Comparison between the model SEDs (solid curve) and the observational ones (blue circle). The left-hand panel is for Mrk 421, and the right-hand panel is for Mrk 501, where the observational data are taken from Abdo et al. (2011b) and Abdo et al. (2011a), respectively. Non-simultaneous archival data taken from NED are plotted with the grey squares below 1015 Hz, for reference.

Upper panel: the characteristic times of acceleration, escape, and cooling for non-thermal electrons, corresponding to the left-hand panel of Fig. 2. The red thick solid line represents the acceleration time tacc, the green dashed line the electron escape time tesc, e, the blue thin solid curve the radiation cooling time at the steady state, $\gamma /\left| \dot{\gamma }_{\rm rad} \right|$, and the cyan dotted line the synchrotron cooling time tsyn. Lower panel: Electron energy distribution of the non-thermal electrons, multiplied by γ2. The electrons are injected at γ = 10 continuously at a constant rate.
Figure 3.

Upper panel: the characteristic times of acceleration, escape, and cooling for non-thermal electrons, corresponding to the left-hand panel of Fig. 2. The red thick solid line represents the acceleration time tacc, the green dashed line the electron escape time tesc, e, the blue thin solid curve the radiation cooling time at the steady state, |$\gamma /\left| \dot{\gamma }_{\rm rad} \right|$|⁠, and the cyan dotted line the synchrotron cooling time tsyn. Lower panel: Electron energy distribution of the non-thermal electrons, multiplied by γ2. The electrons are injected at γ = 10 continuously at a constant rate.

The blue thick dashed curve and dotted one in Fig. 1, respectively, represent the case of higher and lower ηsyn(γ) by a factor of 2, compared to the black thick solid one. These curves show that although the softening of the SEDs is caused by escape effect, it is insufficient to reproduce the high-energy cutoff, and radiative cooling is therefore required in this case.

We adjust the model parameters to fit the model to the full observed SED including an SSC component, using a parameter set obtained from the fitting of the synchrotron component as the initial one and considering only the parameter region satisfying the following two conditions. (I) First, the characteristic synchrotron frequency of the electrons satisfying l = Resc is sufficiently high not to affect the observed SED. We consider only the case that the electrons with Lorentz factor 10 times larger than that corresponding to the maximum frequency of the observed synchrotron component satisfy l < Resc. This choice means that we have the restriction |$\beta _{\rm w}< \eta _{\rm esc}^{-1/2} < 1$| at this electron Lorentz factor because ηesc = (lwResc)2 by definition and at least ηesc > 1 for electrons contributing to the flux at the frequencies around and above the synchrotron peak (as mentioned above in this section). (II) Secondly, since at least the value of (δBeff/B)2 for the electrons treated in equation (2) which are responsible for observed emissions needs to be sufficiently low, we take (δBeff/B)2 < 0.2 for the electrons with the Lorentz factor 102γsyn, which leads to |$\zeta \lambda _{\rm max}^{1-q}< 0.2 (r_{\rm g}(10^2 \gamma _{\rm syn}))^{1-q}$| from equation (5).

4 RESULTS

We apply the model described in Section 2 to the well-sampled long-term averaged SEDs of two HSPs, Mrk 421, and Mrk 501. The results of the fits to the overall observed SEDs are shown in Fig. 2. The adopted parameters are tabulated in Table 1. The observational data for Mrk 421 was taken from 2009 January 19 (MJD 54850) to June 1 (MJD 54983) (fig. 8 in Abdo et al. 2011b) and for Mrk 501 from 2009 March 15 (MJD 54905) to August 1 (MJD 55044) (figs 8 and 9 in Abdo et al. 2011a). For the gamma-ray band of Mrk 501, we adopt the data plotted in fig. 9 of Abdo et al. (2011a) with filled circles, which do not include the time interval MJD 54952–54982 for the time averaging. The excluded term has significant flux variations in the gamma-ray band (e.g. a flux up to about five times higher than the average one at TeV energy) enough to cause a discrepancy between averaged Fermi-LAT data and non-flaring VERITAS data in overlapping energy range if included for the time averaging (see Abdo et al. 2011a for more details). During the campaign Mrk 421 was in a low-activity state at all wavebands (e.g. flux variations are typically less than a factor of 2 relative to the average flux level in gamma-ray band; Abdo et al. 2011b). The data we adopted are suitable ones currently available for modelling the typical low state of these objects. The archival data plotted in Fig. 2 for reference was taken from NED.1

Table 1.

Physical parameters used in the calculation shown in Fig. 2.

ΓBRescVLinjβwq|$\zeta \lambda _{\rm max}^{1-q}$|γinj
Unit(G)(1014 cm)(1047 cm3)(1039 erg s−1)(10−14 cm1 − q)
Mrk 421350.0813.51.91.60.161.854.710
Mrk 501310.0142895150.231.920.1210
ΓBRescVLinjβwq|$\zeta \lambda _{\rm max}^{1-q}$|γinj
Unit(G)(1014 cm)(1047 cm3)(1039 erg s−1)(10−14 cm1 − q)
Mrk 421350.0813.51.91.60.161.854.710
Mrk 501310.0142895150.231.920.1210
Table 1.

Physical parameters used in the calculation shown in Fig. 2.

ΓBRescVLinjβwq|$\zeta \lambda _{\rm max}^{1-q}$|γinj
Unit(G)(1014 cm)(1047 cm3)(1039 erg s−1)(10−14 cm1 − q)
Mrk 421350.0813.51.91.60.161.854.710
Mrk 501310.0142895150.231.920.1210
ΓBRescVLinjβwq|$\zeta \lambda _{\rm max}^{1-q}$|γinj
Unit(G)(1014 cm)(1047 cm3)(1039 erg s−1)(10−14 cm1 − q)
Mrk 421350.0813.51.91.60.161.854.710
Mrk 501310.0142895150.231.920.1210

As shown in Fig. 2, the two samples appear to be fitted well. Here, we should note that the reduced chi-square values of the fits shown in Fig. 2 are 3.9 and 1.3, respectively, for Mrk 421 and Mrk 501 for the gamma-ray band (3.8 × 1022–9.6 × 1026 Hz for Mrk 421 and 3.9 × 1022–1.1 × 1027 Hz for Mrk 501),2 though they are not minimum ones. The current work is not able to specify the causes of the statistical poorness of the fits. This may possibly indicate the need for non-trivial spatial structure of the emission region and/or additional emission components such as an external Compton component in order to reproduce the actual data more accurately. As an example for the latter, the Comptonization of an external radio photon field has been proposed so as to explain somewhat but systematically underestimated flux near 1023 Hz (see section 4.2 in Asano et al. 2014). Nevertheless, we consider that the fits are sufficient for investigating the validity of SA to the issues described in Section 1.

The characteristic synchrotron frequency of the electrons with the Lorentz factor γinj is about 109 Hz for Mrk 421, and about 108 Hz for Mrk 501, measured in the observer frame. The break seen near 1010 Hz in the both panels in this figure is caused by SSA. In the upper panel of Fig. 3, the characteristic times of acceleration, escape, and cooling are displayed as a function of γ for the parameters adopted for Mrk 421. The lower panel shows the calculated EED of Mrk 421. It can be recognized from this figure that gradual softening from the radio band to higher frequency is formed by an increasing dominance of escape over acceleration with an increasing electron energy and that the cutoff at the synchrotron peak is caused by the balance between acceleration and cooling time under rapid escape. Note that not γ2ne(γ) but γ3ne(γ) mimics the synchrotron spectrum (Tramacere et al. 2011). Not only the softening by the efficient escape, but also cooling is required to well reproduce the observed shape of the cutoff in this case.

The synchrotron component of Mrk 501 has a softer and broader portion, compared to Mrk 421, as seen in Fig. 2. If one interprets such an SED from the radio band in the one-zone model, q closer to 2 (adopted value is 1.92) is needed as mentioned in Section 3. Treating the slope of the wave spectrum as a free parameter is essential to apply our simple model to various objects. At least the phenomenologies such as Kolmogorov-type q = 5/3 and Kraichnan-type q = 3/2 (Goldstein, Roberts & Matthaeus 1995) are not fit for the application here; larger q but smaller than 2 is needed. We omit the figures for the time-scales and EED for Mrk 501 because it is basically similar to that for Mrk 421.

The adopted parameter values are listed in Table 1. We call them and the ones obtained with them the canonical values. By examining the relation between a model SED and physical quantities, here we show that some of them are constrained to the canonical values with the SED fitting while others are not. First, we note that the synchrotron component is not modified except by the influence of the IC cooling when we change the parameter values with the flux normalization of the synchrotron component and the values of ηesc and ηsyn fixed as a function of synchrotron frequency measured in the observer frame, νsyn, with the corresponding electron Lorentz factor γ, i.e. we obtain
(7)
(8)
(9)
(10)
where the ratio of a quantity to the canonical value is denoted by the superscript , |$(\nu F_\nu )_{\rm syn}^*$| is the ratio of the flux normalization of the synchrotron component, Lsyn the synchrotron luminosity measured in the blob frame, uB = B2/8π the energy density of the mean magnetic field, ue that of the non-thermal electrons, and |$L_{\rm syn}^* = (u_{\rm B}\gamma u_{\rm e}V)^*$|⁠. Next, for the fixed synchrotron component, the IC component is not modified when
(11)
(12)
(13)
where usyn is the energy density of the synchrotron photons, measured in the blob frame, and |$L_{\rm syn}^* = (u_{\rm syn} V/R_{\rm esc})^*$|⁠. Note that the frequency range of the Thomson scattering regime of the IC component relative to the synchrotron component and the IC cutoff frequency roughly scale as γ2 and (γΓ), respectively. From equations (7)–(13), we obtain
(14)
(15)
(16)
(17)
(18)
Thus, B, Γ, and the total electron energy ueV are well constrained to the canonical values. Moreover, the SED is unchanged if the physical parameters of the blob are changed according to equations (15)–(18) with V. That is, a parameter set adopted for an observed SED as shown in Fig. 2 includes only one degree of freedom and can vary the parameter values within the restrictions (I) and (II) in Section 3. Although Linj, one of the fitting parameters, does not appear in the above analysis explicitly, one can rewrite equation (15) in terms of Linj. Note that all the presented figures are independent of the uncertainty and our general conclusions are unaffected.

For concreteness, we have shown in Table 1 the parameter values of Resc, V, Linj, and |$\zeta \lambda _{\rm max}^{1-q}$| by adopting the maximum allowed value of βw from the restriction (I) described in Section 3. The shown parameter sets can be scaled by V in the range of 0.04 ≤ V ≤ 1 for Mrk 421 and 0.09 ≤ V ≤ 1 for Mrk 501 with the SEDs fixed, where the upper and lower bound come from (I) and (II) in Section 3, respectively.

Energy densities of respective components, measured in the blob frame, are listed in Table 2. Although the magnitude relation between ue and uB is consistent with the conventional SSC model (e.g. Kino et al. 2002), the non-thermal electrons possess much greater energy than the magnetic field, ue/uB ≫ 1, compared to the previous results. This difference originates from relatively small Resc, that is, the importance of the escape effect as seen from equations (11) and (12). Within the uncertainty of the parameters, we can only make ue (and ue/uB) larger from the canonical values since (ueV) = 1 is required to fix the SEDs from equation (15), where V ≤ 1.

Table 2.

Energy contents in the blob for the parameters shown in Table 1: the energy density of magnetic field, uB; that of photons, uph; and that of non-thermal electrons, ue.

uBuphue
Unit(10−2 erg cm−3)(10−2 erg cm−3)(10−2 erg cm−3)
Mrk 4210.030.029
Mrk 5010.00070.0022
uBuphue
Unit(10−2 erg cm−3)(10−2 erg cm−3)(10−2 erg cm−3)
Mrk 4210.030.029
Mrk 5010.00070.0022
Table 2.

Energy contents in the blob for the parameters shown in Table 1: the energy density of magnetic field, uB; that of photons, uph; and that of non-thermal electrons, ue.

uBuphue
Unit(10−2 erg cm−3)(10−2 erg cm−3)(10−2 erg cm−3)
Mrk 4210.030.029
Mrk 5010.00070.0022
uBuphue
Unit(10−2 erg cm−3)(10−2 erg cm−3)(10−2 erg cm−3)
Mrk 4210.030.029
Mrk 5010.00070.0022

5 DISCUSSION

5.1 Shape and the number of the blobs

As seen in Table 1, a relatively small Resc( ≪ (V/π)1/3) is required in order to make the escape effect work to form the gradually softening SEDs. This means that the emission region is either disc-like or spindle like.

However, the extreme shape of the blob may be alleviated if one assumes that a number of the blobs with the same properties simultaneously contribute to an observed spectrum. The number of the blobs Nb can be easily included into the analysis in Section 4 by multiplying Lsyn in equation (8) by Nb. Then V in equations (15)–(18) is replaced with (VNb). It can be recognized from equation (1) and the equation modified from equation (16) with this replacement, |$R_{\rm esc}^* = (V N_{\rm b})^*$|⁠, that for example, when increasing Nb and decreasing V with (VNb) = 1, one can scale the ratio |$R_{\rm esc}/{\rm max}\lbrace R_\perp , R_\Vert \rbrace = {\rm min}\lbrace R_\perp , R_\Vert \rbrace /{\rm max}\lbrace R_\perp , R_\Vert \rbrace$| up to unity as |$(V N_{\rm b}^{3/2})^* = N_{\rm b}^{1/2}$| if |$R_\perp > R_\Vert = R_{\rm esc}$| or as |$(V^2 N_{\rm b}^3)^* = N_{\rm b}$| if |$R_\Vert > R_\perp = R_{\rm esc}$| without changing both the SED and the parameter values (except for Linj). This is the case of multiple blobs contributing to an SED. In the above example, we have |$R_{\rm esc}/{\rm max}\lbrace R_\perp , R_\Vert \rbrace < 1$| with Nb ≲ 103 and ≲ 102 for Mrk 421 and Mrk 501, respectively. However, a large number which leads to |$R_\perp \sim R_\Vert$| is rather unlikely in the spirit of the simple one-zone model approach.

5.2 Energy budget

Let us consider the energy balance between the waves and the accelerated particles. Since we have supposed the net energy transfer between the waves and the non-thermal electrons is from the former to the latter, the energy supply from large-scale disturbances to small-scale ones has to be rapid enough to meet the assumption of the stationary wave spectrum. We estimate the time-scale of damping of the wave energy caused by the electron acceleration, as a function of k, as
(19)
measured in the blob frame, where uw(k) is the total energy density of wave modes around k, |$\dot{\mathcal {U}}_{\rm acc} \equiv \gamma ^2 m_{\rm e} c^2n_{\rm e}(\gamma )/t_{\rm acc}(\gamma )$|⁠, and |$k \sim r_{\rm g}^{-1}$|⁠. If this is too short, waves would damp through particle acceleration and then uw(k) would be modified. For illustrative purpose, we display tdam(k) in Fig. 4 for the steady state we have obtained, supposing |$u_{\rm w}(k) \sim \delta B_{\rm eff}^2/8 \pi$|⁠. By comparing tdam(k) to 1/kvw, which represents a lower limit of the time-scale of the wave cascading, we find that tdam is smaller than 1/kvw in an intermediate range of k. This means that before the blob reaching the steady state, the assumed wave spectrum would be modified and evolve towards a different equilibrium state. The model is hence inconsistent in this case.
The time-scale of damping of MHD disturbances through acceleration of non-thermal electrons in unit of 1/kvw for the case of $u_{\rm w}(k) \sim \delta B_{\rm eff}^2/8 \pi$. The horizontal axis represents wavenumber k normalized by the inverse of the gyroradius of the electrons with the injection Lorentz factor γinj. The dotted line indicates a lower limit of the wave cascading time.
Figure 4.

The time-scale of damping of MHD disturbances through acceleration of non-thermal electrons in unit of 1/kvw for the case of |$u_{\rm w}(k) \sim \delta B_{\rm eff}^2/8 \pi$|⁠. The horizontal axis represents wavenumber k normalized by the inverse of the gyroradius of the electrons with the injection Lorentz factor γinj. The dotted line indicates a lower limit of the wave cascading time.

Under our treatment of the electron transport, this encourages us to interpret vw as the sound speed (not as the Alfvén velocity) since the wave energy can be much larger than the magnetic part in this case if the plasma beta of the emitting region is much larger than 1. It seems qualitatively reasonable because the energy of the magnetic field is much less than that of the particles in the blob. This implies that the fast mode waves may be the acceleration agent. (In this case, uw(k) and so tdam(k) become larger in proportion to the plasma beta although the cascading time is also affected.) The parameter sets shown in Table 1 have the value of vw close to the upper limit of the sound speed (i.e. |$c/\sqrt{3}$|⁠), which seems likely in the blazar region. Anyhow, more sophisticated turbulence modelling accounting for particle acceleration is needed for self-consistency, which is beyond the scope of this paper.

6 SUMMARY

Motivated by the previous implications about the energy distribution of the non-thermal electrons in the blazar region (i–iv described in Section 1), we have examined the stochastic acceleration by a weak random electromagnetic field as the acceleration mechanism of the electrons by modelling the SEDs of two HSPs, Mrk 421 and Mrk 501, from the radio to the gamma-ray bands. All the model parameters are time independent (Section 2). The long-term time averaged SEDs appear to be well fitted by the emitting blob reaching a steady state (Fig. 2). The blob can be considered to represent the jet emission region(s) mainly contributing to the observed SED on average in the long term. The goodness of the fits in gamma-ray band is, however, statistically not sufficient as mentioned in Section 4 and more refined modelling are needed to reproduce the actual data more accurately. With the current simple model, we have obtained the following results.

The model can naturally produce the broad and curved electron energy distribution, which has been suggested based on the conventional models so far (Section 1). To realize such electron energy distribution, two points are crucial in the model: a stationary, power-law spectrum of MHD waves with a steeper spectral index than the Kolmogorov spectrum at the relevant wavelengths and efficient particle escape from the acceleration region (Figs 1 and 3). The former is associated with the latter through weaker energy dependence of the escape efficiency. As for the formation of the high-energy cutoff of the synchrotron component, radiative cooling becomes important.

These results lead to a few issues with the obtained parameter values. First, the shape of the blob is very asymmetric having a relatively small side length compared with the volume of the blob to realize an efficient escape. Secondly, since emitted photons also efficiently escape, in order to have a reasonable value for the ratio of the energy density of the synchrotron photons to that of the magnetic field, which is strictly constrained by an observed spectrum, our model requires a large ratio of the energy density of the non-thermal electrons to that of the magnetic field (equation 11), which ratio is actually orders of magnitude larger than the ones hitherto adopted for the two HSPs (Section 4).

A trivial but important issue of the model is whether and how such properties phenomenologically required for the emission region can be formed in the blazar jets. If we estimate the energy density of the MHD waves responsible for the electron acceleration on the assumption that they propagate at the Alfvén velocity, the derived physical parameters raise a problem with an energy budget of MHD turbulence, i.e. the turbulence would be damped effectively by the particle acceleration (Section 5.2). The model is inconsistent in this case. The difficulty may be reduced to some extent if we assume that the waves propagate at the speed of sound in a high-beta plasma (i.e. the fast mode waves).

Recently, Asano et al. (2014) considered the emission from the electrons stochastically accelerated in the continuous and radially structured jet, where the escape effect can be negligible. In Asano et al. (2014) assuming the Kolmogorov spectrum (q = 5/3) for the turbulent magnetic field, which is different from our treatment of q, they investigated the effect of the radial evolution of the physical parameters on the electron energy distribution and observed SED. They have shown that electron spectra become soft as required even when q = 5/3 by assuming the evolution of the electron injection rate, that of the momentum diffusion coefficient, and adiabatic cooling. Although their jet model is quite a contrast to the blob model in this paper, it is also able to form the curved electron energy distribution that can explain the same SED of Mrk 421 by the stochastic acceleration. Thus, although the stochastic acceleration is a viable model, the way to realize a broad but a non-power-law electron spectrum is not unique and further work is needed to solve the problem.

The authors thank D. Paneque for giving us the νFν data shown in Abdo et al. (2011a,b). This research has made use of the NASA/IPAC Extragalactic Database (NED) operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

2

For the frequency range below the X-ray band, since the observed data show different flux levels at very close frequencies with small statistical errors, we do not evaluate chi-square values.

REFERENCES

Abdo
A. A.
et al.
ApJ
,
2010
, vol.
716
pg.
30
Abdo
A. A.
et al.
ApJ
,
2011a
, vol.
727
pg.
129
Abdo
A. A.
et al.
ApJ
,
2011b
, vol.
736
pg.
131
Aharonian
F.
et al.
ApJ
,
2009
, vol.
696
pg.
L150
Asano
K.
Takahara
F.
Kusunose
M.
Toma
K.
Kakuwa
J.
ApJ
,
2014
, vol.
780
pg.
64
Becker
P. A.
Le
T.
Dermer
C. D.
ApJ
,
2006
, vol.
647
pg.
539
Bednarz
J.
Ostrowski
M.
MNRAS
,
1996
, vol.
283
pg.
447
Blandford
R.
Eichler
D.
Phys. Rep.
,
1987
, vol.
154
pg.
1
Blumenthal
G. R.
Gould
R. J.
Rev. Mod. Phys.
,
1970
, vol.
42
pg.
237
Cao
G.
Wang
J.
PASJ
,
2013
, vol.
65
pg.
109
Dermer
C.
Lott
B.
J. Phys. Conf. Ser.
,
2012
, vol.
355
pg.
012010
Finke
J. D.
Dermer
C. D.
Böttcher
M.
ApJ
,
2008
, vol.
686
pg.
181
Fossati
G.
Maraschi
L.
Celotti
A.
Comastri
A.
Ghisellini
G.
MNRAS
,
1998
, vol.
299
pg.
433
Goldstein
M. L.
Roberts
D. A.
Matthaeus
W. H.
ARA&A
,
1995
, vol.
33
pg.
283
Hall
D. E.
Sturrock
P. A.
Phys. Fluids
,
1967
, vol.
10
pg.
2620
Inoue
S.
Takahara
F.
ApJ
,
1996
, vol.
463
pg.
555
Jones
F. C.
Phys. Rev.
,
1968
, vol.
167
pg.
1159
Kataoka
J.
Proc. Sci.
,
2008
 
Characterizing X-ray variability of blazars. SISSA, Trieste, PoS(BLAZARS2008)015
Kataoka
J.
et al.
ApJ
,
1999
, vol.
514
pg.
138
Kataoka
J.
Takahashi
T.
Makino
F.
Inoue
S.
Madejski
G. M.
Tashiro
M.
Urry
C. M.
Kubo
H.
ApJ
,
2000
, vol.
528
pg.
243
Katarzyński
K.
Ghisellini
G.
Mastichiadis
A.
Tavecchio
F.
Maraschi
L.
A&A
,
2006
, vol.
453
pg.
47
Kino
M.
Takahara
F.
Kusunose
M.
ApJ
,
2002
, vol.
564
pg.
97
Krawczynski
H.
et al.
ApJ
,
2004
, vol.
601
pg.
151
Kusunose
M.
Takahara
F.
ApJ
,
2008
, vol.
682
pg.
784
Kusunose
M.
Takahara
F.
Li
H.
ApJ
,
2000
, vol.
536
pg.
299
Lefa
E.
Aharonian
F. A.
Rieger
F. M.
ApJ
,
2011
, vol.
743
pg.
L19
Lemoine
M.
Revenu
B.
MNRAS
,
2006
, vol.
366
pg.
635
Li
H.
Kusunose
M.
ApJ
,
2000
, vol.
536
pg.
729
Longair
M. S.
High Energy Astrophysics
,
1992
Cambridge
Cambridge Univ. Press
Maraschi
L.
Ghisellini
G.
Celotti
A.
Courvoisier
T.
Blecha
A.
Proc. IAU Symp. 159, Multi-Wavelength Continuum Emission of AGN
,
1994
Dordrecht
Kluwer
pg.
221
Marscher
A. P.
Belloni
T.
Lecture Notes in Phys.
,
2010
 
Vol. 794, The Jet Paradigm, Springer, Berlin, p. 173
Massaro
E.
Tramacere
A.
Perri
M.
Giommi
P.
Tosti
G.
A&A
,
2006
, vol.
448
pg.
861
Mastichiadis
A.
Kirk
J. G.
A&A
,
1997
, vol.
320
pg.
19
Mertsch
P.
J. Cosmol. Astropart. Phys.
,
2011
, vol.
12
pg.
10
Niemiec
J.
Ostrowski
M.
ApJ
,
2004
, vol.
610
pg.
851
Ostrowski
M.
Siemieniec-Oziębło
G.
Astropart. Phys.
,
1997
, vol.
6
pg.
271
Park
B. T.
Petrosian
V.
ApJS
,
1996
, vol.
103
pg.
255
Press
W. H.
Teukolsky
S. A.
Vetterling
W. T.
Flannery
B. P.
Numerical Recipes in C. The Art of Scientific Computing
,
1992
Cambridge
Cambridge Univ. Press
Rani
B.
et al.
A&A
,
2013
, vol.
552
pg.
A11
Rybicki
G. B.
Lightman
A. P.
Radiative Processes in Astrophysics
,
1979
New York
Wiley
Schlickeiser
R.
Cosmic Ray Astrophysics
,
2002
Berlin
Springer-Verlag
Sikora
M.
Madejski
G.
Aharonian
F. A.
Völk
H. J.
AIP Conf. Proc. Vol. 558, High Energy Gamma-Ray Astronomy
,
2001
New York
Am. Inst. Phys.
pg.
275
Sikora
M.
Madejski
G.
Moderski
R.
Poutanen
J.
ApJ
,
1997
, vol.
484
pg.
108
Stawarz
Ł.
Petrosian
V.
ApJ
,
2008
, vol.
681
pg.
1725
Tagliaferri
G.
et al.
ApJ
,
2008
, vol.
679
pg.
1029
Tavecchio
F.
et al.
ApJ
,
2001
, vol.
554
pg.
725
Tavecchio
F.
Ghisellini
G.
Ghirlanda
G.
Foschini
L.
Maraschi
L.
MNRAS
,
2010
, vol.
401
pg.
1570
Tramacere
A.
Massaro
E.
Taylor
A. M.
ApJ
,
2011
, vol.
739
pg.
66
Urry
C. M.
Padovani
P.
PASP
,
1995
, vol.
107
pg.
803
Ushio
M.
et al.
ApJ
,
2009
, vol.
699
pg.
1964
Ushio
M.
et al.
ApJ
,
2010
, vol.
724
pg.
1509
Weidinger
M.
Rüger
M.
Spanier
F.
Astrophys. Space Sci. Trans.
,
2010
, vol.
6
pg.
1
Wentzel
D. G.
ARA&A
,
1974
, vol.
12
pg.
71
Yan
D.
Zeng
H.
Zhang
L.
MNRAS
,
2012
, vol.
424
pg.
2173
Zhang
J.
Liang
E.-W.
Zhang
S.-N.
Bai
J. M.
ApJ
,
2012
, vol.
752
pg.
157