© ESO, 2012

1. Introduction

Active galactic nuclei (AGN) are astrophysical sources emitting broad-band electromagnetic radiation that can be observed from radio wavelengths to very high energies (VHE, E ≳ 100 GeV). The current framework for understanding the emission from these objects depicts them as composed of a disk, feeding a super massive black hole (SMBH), possibly with a jet on each side of the accretion plane. The orientation of the system with regard to the observer’s line of sight is one of the key parameters for unifying the various subclasses of AGN (Urry & Padovani 1995). The class of blazars for which the jet is probably directed within a few degrees toward the observer is of particular interest for the study of non-thermal processes occurring in astrophysical relativistic flow.

The VHE emission of blazars, such as the BL Lac object Mrk 421 or the flat spectrum radio quasar PKS 1222+21, can sometimes exhibit doubling times of the order of ten minutes (Gaidos et al. 1996; Aleksić et al. 2011) while variability on the minute scale is observed in light curves of BL Lac objects, such as Mrk 501 (Albert et al. 2007) or PKS 2155–304 (Aharonian et al. 2007). The high quality of the VHE data on the latter object has enabled the characterization of various statistical properties of the emission, such as a skewed flux distribution, interpreted as a log-normal distribution, and a linear relation between the sample rms of the flux and the sample flux (H.E.S.S. Collaboration, Abramowski et al. 2010).

These properties were initially studied in the X-ray emission of the galactic black hole binary Cygnus X-1 (Uttley & McHardy 2001; Uttley et al. 2005) and are now observed in other accreting objects, such as the blazar BL Lac (Giebels & Degrange 2009), or non-aligned AGNs, such as the Seyfert galaxies NGC 4051 (McHardy et al. 2004) and IRAS 13224-3809 (Gaskell 2004).

They are sometimes interpreted as arising from multiplicative processes, originating e.g. in the accretion disk. The product of random variables whose logarithms have finite moments is indeed asymptotically log-normal. Furthermore, flux log-normality implies a linear relation between the flux and its rms (see Appendix D of Uttley et al. 2005).

The statistical properties of non-thermal light curves are key ingredients for answering the crucial question of the origin of variability in AGNs: do the variations of the flux come from the jet itself, i.e. the medium that conveys the accelerated particles, or from the disk that could modulate the jet emission. In accretion models such as developed by Lyubarskii (1997), inward propagating fluctuations in the disk from different radii cumulate in a multiplicative way (see e.g. Arévalo & Uttley 2006), which could explain the log-normality and the linear rms-flux relation.

However, variability on time scales shorter than the light-crossing time T_G = 2GM/c³ are difficult to achieve (Narayan & Piran 2012) in disks. Transient emission on time scales t_var ~ T_G/50 are easier to locate in a jet with a bulk Lorentz factor Γ (unless generated by a companion system, as e.g. in Rieger & Volpe 2010). Phenomena occurring on time scales T_G then appear at t_var ~ T_G/δ to the observer, where δ is the Doppler factor of the jet¹.

Values of the Doppler factor δ ≥ 50 are not uncommon for VHE blazars and permit both minute-time-scale variations (Aharonian et al. 2007) and an optical thinness to the emitted VHE γ-rays (Begelman et al. 2008). Ghisellini & Tavecchio (2008) propose a needle-in-a-jet model, where fast minijets closely aligned with the line of sight could explain the variations. Giannios et al. (2009) have investigated a jets-in-a-jet scenario, involving reconnection and a Poynting-dominated flow. Narayan & Piran (2012) considered a turbulent process, modelled by a random motion of minijets inside the jet. Models involving minijets have also been invoked to explain the fast variations in the emission of radio galaxies (Giannios et al. 2010), though the time scales of variations observed at VHE are much longer than for blazars, approximately a day for M 87 (Abramowski et al. 2012).

These models are additive (the flux is given by the sum of the contributions of several regions) and, provided there are finite moments for the individual components, the central limit theorem (CLT) can be applied. But the Gaussian flux distribution expected from the CLT is not skewed, as is the observed one.

We seek a way out the jet-or-disk origin dilemma by studying the implications of the observed linearity between the sample rms and the flux (Sect. 2.1). We discuss the properties of Pareto-distributed fluxes in Sect. 2.2 and discuss the relevance of the CLT for these random variables. We show in Sect. 3 that Pareto distributions are naturally generated by the minijets-in-a-jet statistical model and discuss the consequences of this model in Sect. 4.

2. Additive processes and Pareto distributions

Modelling of VHE blazar variability should reproduce the short time scale of the variations and the statistical properties of the flux. While the stochastic nature of the flux is understood if the emission is modulated by the disk of the SMBH, the observed time scales are difficult to explain. Inversely, a jet origin of the variability is often rejected because the models are additive and should not be able to reproduce the skewness of the flux nor the linear rms-flux relation. This argument is studied in the next section.

2.1. Additive and multiplicative processes

Within additive models, the observed flux is the sum of the contributions of several (and potentially many) regions. If the components are modelled as independent, identically distributed random variables (hereafter iid) with finite moments, their sum should follow a normal distribution according to the CLT, assuming that the number of regions is large enough – typically more than a few tens. A natural outcome of such a process would then be symmetric flux distributions, in disagreement with the high-flux tails frequently observed in light curves. Moreover, the addition of incoherent variations in various components should result in a flux being independent of its variance. A translation into statistical terms would be that the sample mean² and the sample rms of a distribution are independent if and only if the underlying distribution is a Gaussian (see the discussion in the Appendix D of Uttley et al. 2005).

A contrario, multiplicative processes naturally generate log-normality, with a characteristic tailed and skewed distribution. If we let the quantity φ be the product of a large number of iid variables, then the logarithm of φ is a sum of iid quantities and, according to the CLT, log φ follows a normal distribution, i.e. φ follows a log-normal distribution (see e.g. the simulations of Ioka & Nakamura 2002, with only three multiplicative variables).

The rms-flux relation can be explained as a consequence of a log-normal distribution of the flux (Uttley et al. 2005), but the reciprocal is not true. Indeed, let the observed flux φ be a function of an underlying random variable x, so that φ = f(x); for a log-normal distribution, f is the exponential function and x is normally distributed. A small fluctuation of x around x₀, δx, results in a small fluctuation of the flux φ around f(x₀), δφ, and the variance of the flux is $\mathit{\delta }{\mathit{\phi }}^{\mathrm{2}}\mathrm{=}{\left[\frac{\mathit{\partial f}}{\mathit{\partial x}}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathrm{\right)}\right]}^{\mathrm{2}}\mathit{\delta }{\mathit{x}}^{\mathrm{2}}\mathrm{·}$(1)With a sample flux proportional to f(x₀), a linear flux-rms relation is equivalent to $\mathit{f}\mathrm{(}{\mathit{x}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{\propto }{\left[\frac{\mathit{\partial f}}{\mathit{\partial x}}\mathrm{\left(}{\mathit{x}}_{\mathrm{0}}\mathrm{\right)}\right]}^{\mathrm{2}}\mathit{\delta }{\mathit{x}}^{\mathrm{2}}\mathit{.}$(2)Equation (2) is one of the definitions of the exponential function, so the sample rms is proportional to the sample flux if and only if the flux is the exponential of an underlying variable, as is the case for a log-normally distributed flux. Log-normality thus implies linearity between the rms and the flux³, as does any distribution that arises from the exponential of an underlying variable.

2.2. Pareto distributions and the central limit theorem

We now consider the class of distributions that, instead of arising from the exponential of normal distributions, arise from the exponential of exponential distributions exp(− αx). These are known as Pareto distributions and are characterized by a probability density function (PDF) that follows a power-law function of index 1 + α: ${\mathit{f}}_{\mathit{Y}}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\frac{\mathit{\alpha }}{{\mathit{y}}^{\mathrm{1}\mathrm{+}\mathit{\alpha }}} \mathrm{for}\mathit{y}\mathit{>}\mathrm{1.}$(3)Pareto distributions are used in seismology, in finance, and in biology and various examples of applications such as “volcanic eruptions, solar-flares, lightning strikes, river networks, forest fires, extinctions of biological species, war casualties, Internet traffic, stock returns, insurance pay-offs”, can be found in the literature (see, e.g., Newman 2005; Zaliapin et al. 2005, and reference therein).

The CLT cannot be applied to Pareto distributed iid, mostly because, for α ≤ 2, the variance of each random variable cannot be defined. A generalized central limit theorem has been established for such a high-tailed distribution (see, e.g., Voit 2005), and it states that the sum of Pareto-distributed iid converge to α-stable distributions, a class of probability distribution with fascinating properties.

Only a handful of α-stable distributions can be expressed in terms of elementary functions. Some special cases are symmetric, such as the Gaussian distribution, which corresponds to the limit α = 2, but a large skewness is generally seen, e.g., with the Landau distribution (α = 1), used by particle physicists (Groom & Klein 2000) to describe the distribution of the energy loss of a charged particle passing through matter. In particular, for 0 < α < 2 the PDF of an α-stable distribution asymptotically follows a power law of index 1 + α (Zaliapin et al. 2005; Nolan 2012).

Fig. 1 Schematic view of the geometry. The left side corresponds to the observer frame ℛobs, in which the jet is tilted by angle θ from the line of sight (along xobs) and is boosted by a Lorentz factor Γ (velocity Σ) along ${\mathit{x}}_{\mathrm{jet}}^{\mathrm{\prime }}$. The minijet is defined by its Lorentz factor γ (velocity β), in the jet frame ℛjet (top right). The orientation of the minijet along ${\mathit{x}}_{\mathrm{minijet}}^{\mathrm{\prime }}$ compared to the jet axis in its rest frame, xjet, is defined by the spherical angles ψ and ϕ. The emission is assumed isotropic in the minijet frame ℛminijet (bottom right).

The convergence of the CLT is achieved with a few tens of “regular”⁴ components, but sums of Pareto iid variables tend to the asymptotic α-stable distributions only for large numbers of components, typically 10⁴ (Zaliapin et al. 2005), which prevents the PDF of the sum from being analytically derived for an intermediate number of Paretian iid quantities. Even though the log-normal distribution is not an attractor for the sum of such iid variables, skewed distributions that closely resemble the log-normal can be expected, as is shown in Sect. 3.3.

Equations (1) and (2) prove that the sample rms is proportional to the sample flux for a single random variable as long as it is the exponential of an underlying variable. This is the case for Pareto distributions, and it will be shown in Sect. 3.3 that this relation holds for the sum of a large number of Paretian iid, certainly because the PDF of the sum follows the same power-law behaviour for high fluxes.

3. A minijets-in-a-jet statistical model

Power laws are generally seen as the signature of a scale invariant process and a wide variety of models can certainly produce such flux distributions. In the next sections, a simple kinematic model that naturally generates Paretian fluxes is developed.

The relativistic enhancement of the flux, known as Doppler boosting, is computed for an emitting region randomly oriented in a medium. For the sake of generality, both the emitting region and the medium, called minijet and jet respectively, have a relativistic motion. Beamed astrophysical objects thus fit in this general framework, as well as non-beamed objects, when the Lorentz factor of the medium is set to unity.

3.1. Doppler factor of a minijet randomly oriented in a jet

We consider a minijet, e.g. a blob of plasma, that moves randomly in a jet and is characterized by its Lorentz factor γ and associated velocity $\mathit{\beta }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\gamma }}^{-2}}$, within the jet frame. The orientation of the minijet in the jet frame can be defined with the conventional angles in spherical coordinates, which will be called ψ and ϕ. The jet is defined by a Lorentz factor Γ, associated velocity Σ, within the observer frame and its axis is misaligned from the line of sight by an angle θ. A schematic representation of these parameters is shown in Fig. 1.

To derive the Doppler factor δ of the emitting region, the minijet, one only needs to determine the energy measured in the observer frame as a function of the energy in the minijet frame, where the emission is assumed to be isotropic. We consider a photon emitted in the minijet frame with energy E_minijet and momentum p_minijet =  { p_x  minijet,p_y  minijet,p_z  minijet } . A series of transformations must be applied to derive the energy in the observer frame: the tilt of the jet, a boost, the tilt of the minijet in the jet, defined by two angles, and finally the boost of the minijet. Only photons propagating along the x_obs direction, i.e. along the line of sight, are considered because these are the only ones observed. We can thus set the transverse momentum to zero, p_y  obs = p_z  obs = 0, while the longitudinal momentum is p_x  obs = E, where the speed of light is set to unity and where E is the photon energy in the observer frame.

These equalities and the frame transformations are summarized in Eq. (4), where c and s denotes the cosine and sine functions: $\begin{array}{c}\\ \left[\begin{array}{}\\ \mathit{E}\\ \mathit{E}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]\mathrm{=}\left[\begin{array}{}\\ \mathrm{1}\\ {\mathrm{c}}_{\mathit{\theta }}{\mathrm{s}}_{\mathit{\theta }}\\ \mathrm{-}{\mathrm{s}}_{\mathit{\theta }}{\mathrm{c}}_{\mathit{\theta }}\\ \mathrm{1}\end{array}\right]\left[\begin{array}{}\\ \mathrm{\Gamma }\mathrm{\Gamma \Sigma }\\ \mathrm{\Gamma \Sigma }\mathrm{\Gamma }\\ \mathrm{1}\\ \mathrm{1}\end{array}\right]\left[\begin{array}{}\\ \mathrm{1}\\ \mathrm{1}\\ {\mathrm{c}}_{\mathit{\varphi }}{\mathrm{s}}_{\mathit{\varphi }}\\ \mathrm{-}{\mathrm{s}}_{\mathit{\varphi }}{\mathrm{c}}_{\mathit{\varphi }}\end{array}\right]\\ \mathrm{\times }\left[\begin{array}{}\\ \mathrm{1}\\ {\mathrm{c}}_{\mathit{\psi }}{\mathrm{s}}_{\mathit{\psi }}\\ \mathrm{-}{\mathrm{s}}_{\mathit{\psi }}{\mathrm{c}}_{\mathit{\psi }}\\ \mathrm{1}\end{array}\right]\left[\begin{array}{}\\ \mathit{\gamma }\mathit{\gamma \beta }\\ \mathit{\gamma \beta }\mathit{\gamma }\\ \mathrm{1}\\ \mathrm{1}\end{array}\right]\left[\begin{array}{}\\ {\mathit{E}}_{\mathrm{minijet}}\\ {\mathit{p}}_{\mathit{x}\mathrm{minijet}}\\ {\mathit{p}}_{\mathit{y}\mathrm{minijet}}\\ {\mathit{p}}_{\mathit{z}\mathrm{minijet}}\end{array}\right]\mathit{.}\end{array}$(4)The inverse of the Doppler factor, i.e. the ratio between the emitted and the observed energies, is derived by inverting Eq. (4). The time-like component of the vectorial equality then reads as δ^-1E = E_minijet, with δ given in Eq. (5): $\mathit{\delta }\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{\Sigma }\mathit{\beta }{\mathrm{c}}_{\mathit{\psi }}\mathrm{-}\mathrm{(}\mathrm{\Sigma }\mathrm{+}\mathit{\beta }{\mathrm{c}}_{\mathit{\psi }}\mathrm{\left)}{\mathrm{c}}_{\mathit{\theta }}\mathrm{\right)}\mathrm{+}\mathit{\gamma \beta }{\mathrm{s}}_{\mathit{\theta }}{\mathrm{s}}_{\mathit{\psi }}{\mathrm{c}}_{\mathit{\varphi }}}\mathrm{·}$(5)For a jet closely aligned with the line of sight, the extrema of the Doppler factor are $\begin{array}{c}\\ \mathit{\delta }\mathrm{\le }\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Sigma }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}}\mathrm{~}\mathrm{4}\mathrm{\Gamma }\mathit{\gamma }\\ \mathit{\delta }\mathrm{\ge }\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Sigma }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}}\mathrm{~}\mathrm{\Gamma }\mathit{/}\mathit{\gamma }\end{array}$(6)where the approximated expressions correspond to the ultra-relativistic limit.

It is worth noting that for reasonable Lorentz factors such as Γ = 5 and γ = 5, the maximum Doppler factor obtained for a minijet and a jet aligned with the line of sight is enhanced by a factor 2γ = 10 compared to isotropic emission in the jet frame (γ = 1), and Doppler factors as high as 10² are reached.

Figure 2 shows the Doppler factor as a function of the misalignment from the line of sight θ. The Lorentz factors of the jet and of the minijet are fixed to a value of five and θ is normalized to 1/Γ, the opening angle of the jet.

Fig. 2 Doppler factor of a minijet with a Lorentz factor γ = 5 moving in a jet of Lorentz factor Γ = 5 as a function of the misalignment between the jet and the line of sight θ, normalized to the opening angle of the jet 1/Γ. The maximum, average, and minimum Doppler factor are compared to the usual Doppler factor (solid line), derived for an isotropic emission in the jet frame (γ = 1).

The maximum and minimum values of the Doppler factor are represented with the long and short dashed lines. The average Doppler factor is derived with μ = cosψ uniformly distributed between [–1, 1] and ϕ uniformly distributed between [0, 2π]. This corresponds to an isotropic distribution of the minijet orientation in the jet frame. One can compute the average Doppler factor in the blazar case, i.e. for θ = 0 $\begin{array}{ccc}\mathit{\delta ̂}\mathrm{0}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Sigma }\mathrm{)}}{\mathrm{\int }}_{-1}^{\mathrm{1}}\mathrm{d}\mathit{\mu }\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\mathit{\beta \mu }}\\ \mathit{\delta ̂}\mathrm{0}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Sigma }\mathrm{)}\mathit{\beta }}\ln \left(\frac{\mathrm{1}\mathrm{+}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\right)\mathrm{~}\frac{\mathrm{4}\mathrm{\Gamma }}{\mathit{\gamma }}\ln \mathrm{2}\mathit{\gamma }\mathit{.}\end{array}$For γ ≳ 4, the reference Doppler factor (γ = 1), represented with the solid line in Fig. 2, is larger than the average Doppler factor, a comparison that can be extended only to quantities that are linearly dependent on the Doppler factor. The impact on the flux intensity, which is roughly a quartic function of the Doppler factor, is discussed in the following.

3.2. Flux of a single minijet

Let I(E) be the flux intensity at energy E, the quantity I(E)/E³ is a Lorentz invariant (e.g., Rybicki & Lightman 1979), and the Doppler boost δ of the energy results in an enhancement δ³ of the intensity or of the flux (see, e.g., Urry & Padovani 1995). One finds the usual enhancement of the bolometric luminosity by a factor δ⁴ after integration over the energy domain.

The non-thermal spectra of high-energy astrophysical sources can generally be approximated, at least locally, by power-law functions. We characterize the intensity of the source by the index s (photon index s + 1), i.e. $\mathit{I}\mathrm{\left(}{\mathit{E}}_{\mathrm{minijet}}\mathrm{\right)}\mathrm{\propto }{\mathit{E}}_{\mathrm{minijet}}^{\mathrm{-}\mathit{s}}$, where E_minijet is the photon energy in the source rest frame. The intensity measured by the observer is, in the observer’s frame, $\mathit{I}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathit{I}\mathrm{\left(}{\mathit{E}}_{\mathrm{minijet}}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathit{I}\mathrm{(}\mathit{E}\mathit{/}\mathit{\delta }\mathrm{)}\mathrm{\propto }{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}\mathit{E}\mathit{/}\mathit{\delta }{\mathrm{)}}^{\mathrm{-}\mathit{s}}\mathrm{\propto }{\mathit{\delta }}^{\mathrm{3}\mathrm{+}\mathit{s}}{\mathit{E}}^{\mathrm{-}\mathit{s}}\mathit{.}$(9)The distribution of the intensity can then be derived from Eq. (5), imposing a distribution for the underlying random variables. We assume in this model that the relevant physical parameter is the orientation of the minijet(s) in the jet frame. Imposing isotropy in this frame corresponds to having ψ and ϕ angles homogeneously populate the unit sphere, i.e. μ = cosψ uniformly distributed between [–1, 1] and ϕ uniformly distributed between [0, 2π]. For the sake of simplicity, we derive the probability density function of the intensity in the blazar case. In this case, combining the definition of the Doppler factor in Eq. (5) with Eq. (9) yields $\mathit{I}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{\propto }{\left[\mathit{\gamma }\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Sigma }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta \mu }\mathrm{)}\right]}^{\mathrm{-}\mathrm{3}\mathrm{-}\mathit{s}}{\mathit{E}}^{\mathrm{-}\mathit{s}}\mathrm{\equiv }\mathrm{\left(}\mathrm{4}\mathrm{\Gamma }\mathit{\gamma }{\mathrm{)}}^{\mathrm{3}\mathrm{+}\mathit{s}}\mathit{g}\mathrm{\right(}\mathit{\mu }\mathrm{)}{\mathit{E}}^{\mathrm{-}\mathit{s}}\mathit{.}$(10)The factor (4Γγ)^3+s scales the function g so that $\mathit{g}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{\Sigma }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\beta }\mathrm{)}}{\mathrm{4}}\mathrm{\times }\frac{\mathrm{1}\mathrm{-}\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta \mu }}\right)}^{\mathrm{3}\mathrm{+}\mathit{s}}\mathrm{\le }\mathrm{1}\mathrm{with}\mathit{\mu }\mathrm{\in }\mathrm{[}\mathrm{-}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{]}\mathit{.}$(11)We call I_N = g(μ) the intensity normalized to its maximum. The PDF of I_N, f_ℐ(I_N), is linked to the PDF of μ = cosψ, with μ ∈  [− 1;1]  via ${\mathit{f}}_{\mathrm{ℐ}}\mathrm{\left(}{\mathit{I}}_{\mathit{N}}\mathrm{\right)}\mathrm{=}\left|\frac{\mathit{\partial }{\mathit{g}}^{-1}\mathrm{\left(}{\mathit{I}}_{\mathit{N}}\mathrm{\right)}}{\mathit{\partial I}}\right|{\mathit{f}}_{\mathrm{𝒞}}\mathrm{\left(}{\mathit{g}}^{-1}\mathrm{\right(}{\mathit{I}}_{\mathit{N}}{\mathrm{\left)}}^{\mathrm{\right)}}\mathit{,}$(12)where ${\mathit{g}}^{-1}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{\times }\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}\mathrm{+}\mathrm{\Sigma }}{\mathrm{2}}\mathrm{\times }\frac{\mathrm{1}}{\mathrm{2}{\mathit{\gamma }}^{\mathrm{2}}}\mathrm{\times }{\mathit{x}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}\mathrm{+}\mathit{s}}}\right)$.

Then, the probability density function of the normalized intensity is ${\mathit{f}}_{\mathrm{ℐ}}\mathrm{\left(}{\mathit{I}}_{\mathit{N}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathrm{\Sigma }}{\mathrm{2}\mathit{\beta }}\mathrm{\times }\frac{\mathrm{1}}{\mathrm{4}{\mathit{\gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{3}\mathrm{+}\mathit{s}\mathrm{)}}\mathrm{\times }\frac{\mathrm{1}}{{\mathit{I}}_{\mathit{N}}^{\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}\mathrm{+}\mathit{s}}}}\mathrm{·}$(13)The flux emitted by a minijet randomly oriented in a jet (or in a non-relativistic medium) thus follows a Pareto distribution of index α = 1/(3 + s), where s is the spectral index of the intensity. The Pareto distribution holds for a non-beamed object (Γ = 1 and Σ = 0), as derived independently by Clausen-Brown & Lyutikov (2012) to model the flares of the Crab in the high-energy domain.

The Pareto distribution also holds for beamed objects misaligned from the line of sight, as shown in Fig. 3, where the distribution of the logarithm of the intensity  ∝ δ^3+s has been generated by drawing random minijet orientations. The conservation of the slope of the distribution can be understood by neglecting the term γβs_θs_ψc_ϕ in Eq. (5), in which case the previous proof holds with the inverse of the Doppler factor remaining a linear function of cosψ. The conclusions derived in the following can then be applied to both blazars, such as PKS 2155-304, and non-aligned objects, such as M 87.

Fig. 3 Distribution of the logarithm of the intensity of a minijet for various angles θ between the line of sight and the jet axis. The intensity is proportional to δ3 + s, where δ is the Doppler factor and s the spectral index. These distributions are obtained simulating random orientations of a minijet for s = 1, for a jet boost Γ = 5, and a minijet boost γ = 5, without generality loss. While the intensity decreases when the misalignment increases, the slope of the distribution remains the same as in the “blazar” case, θ = 0.

3.3. Flux of N minijets

As discussed in Sect. 2.2, Pareto-distributed random variables do not follow the central limit theorem, and their sum does not converge to a normal distribution. Assuming that the jet is composed of several randomly oriented minijets, the total emitted flux is proportional to the sum of the N independent minijet intensities, in an optically thin medium.

To illustrate this point, minijets of Lorentz factor γ = 5 are simulated for θ = 0 and a jet boost Γ = 5. These moderated values yield a single component flux spanning over eight decades, roughly between (Γ/γ)^3+s and (4Γγ)^3+s. Unless strongly suppressing the beaming (typically for γ ≲ 1.5), the dynamic range of the flux is large enough that the location of the cut-off does not affect the following results (up to at least 10⁴ additive components).

As shown in Sect. 3.2, the results for other orientations or other Lorentz factors remain identical within a multiplicative factor (i.e. a shift in log representation). The flux distributions simulated for several numbers N of minijets are shown in Fig. 4. To generate smooth distributions, 10⁸ iterations, called hereafter time steps, are performed with N ∈  { 1,10,30,10²,3 × 10²,10³,3 × 10³,10⁴ } .

Fig. 4 Distribution of the logarithm of the intensity of N independent and randomly oriented minijets. The number of minijets N increases from left to right with N ∈  { 1,10,30,102,3 × 102,103,3 × 103,104 } . Even for a large number of regions, asymmetrical, tailed distributions are obtained.

Fig. 5 Distribution of the logarithm of the flux of N minijets for N = 3 × 103 (left) and N = 104 (right), as in Fig. 4. The continuous black and grey dashed lines represent the best fit with a log-normal and normal flux distributions, respectively.

Fig. 6 Sample rms as a function of the sample flux of N = 1 minijet (top) and the sum of N = 104 minijets (bottom). Linear relations are found in both cases, with a zero x-intercept in the first case and a positive one in the second.

For a large number of minijets, typically N ≳ 10³, the distribution of the logarithm of the flux can be described with a peak, followed by a power-law tail. The distributions obtained for 3 × 10³ and 10⁴ minijets are shown Fig. 5, with a linear y-axis, together with the best-fit functions corresponding to a normal and log-normal flux.

Although the distribution of the flux of N minijets is neither normal nor log-normal, Fig. 5 illustrates how an experimental distribution, with limited statistics and a poor characterization of the high flux tail, could be interpreted as a log-normal distribution, even if arising from an additive process and not from a multiplicative one.

We have shown in Sect. 2.2 that the rms is strictly proportional to the flux if and only if the flux is the exponential of an underlying variable, which is the case for Pareto distributions. Interestingly, linearity holds when considering the sum of a large number of power-law components. To illustrate this statement, light curves are simulated with 10⁵ time steps for N = 1 and N = 10⁴ minijets. The sample mean and the sample rms are then computed in ten points wide windows. For the sake of clarity, the flux and the rms are averaged in 50 bins between the maximum flux and the minimum flux. The sample rms is plotted as function of the sample flux in Fig. 6, where the error bars correspond to the standard deviation in each bin.

The positive x-intercept in the rms-flux relation, as seen in the right-hand panel in Fig. 6, corresponds to the peak in the flux distribution shown in Fig. 4. While it could be interpreted as a constant component, this threshold flux corresponds to the value around which the average emissions of the minijets pile up within the additive model presented here.

4. Discussion and conclusion

The linear relation between the rms and the flux of individual blazars or the skewness of the distribution of the flux are interpreted by several authors as arising from multiplicative processes, favouring a variability stemming from the disk.

We are studying a minijets-in-a-jet statistical model where the variability stems from the jet itself. We first considered the enhancement of the flux due to the relativistic Doppler effect and show that an isotropic orientation of a single region in the jet frame result in a Pareto distribution of its contribution to the flux, characterized by a power-law PDF. As for log-normal variables, Pareto-distributed variables can be seen as the exponential of an underlying variable and, as such, satisfy the linear relation between the flux and its rms. These random variables do not fulfil the hypotheses under which the central limit theorem can be applied and the summation of the contributions of many regions does not yield a normal distribution of the flux. The total flux distribution thus remains highly skewed and tailed when increasing the number of contributing regions, up to at least 10⁴ regions, and could very well be interpreted as a log-normal. A noticeable difference between the attractor of the sum of Pareto distributed variables, called α-stable distributions, and the log-normal distribution is their high flux tails that follow a power law and that could be probed with long-term, finely-sampled light curves. This power-law tail is, moreover, responsible for the invariance of the linear relation between the rms and the flux when the number of minijets is increased. The study of the link between the slope of the rms-flux relation and the kinematic parameters of the mini-jets-in-a-jet statistical model is beyond the scope of this paper and will be the object of further studies.

The minijets-in-a-jet statistical model is largely inspired from the jets-in-jet model of Giannios et al. (2009), where magnetic dissipation triggers the emission of relativistic blobs of plasma. The authors reproduce the luminosity and the fast variations observed in the light curves of blazars at very high energies with a Lorentz factor of the minijet γ ~ 10 in the jet frame. Narayan & Piran (2012) refined this model and invoke relativistic turbulence to explain both the variability time scales and the scarcity of the flaring events. The random orientation of the minijet in the jet frame would be directly linked to the direction in which the reconnection region is created (magnetic dissipation) or to the wandering of the velocity vector of the minijet over the radiating time (relativistic turbulence). The isotropic orientation of the minijets in the jet frame is the most natural configuration in both scenarios.

A “minijet statistical model”⁵ has also recently been derived by Clausen-Brown & Lyutikov (2012) to explain the high-energy flares of Crab Nebula. The authors assume that magnetic reconnection triggers the emission of two-sided relativistic flows with a random orientation in the nebula frame, a particular geometry corresponding to a bulk Lorentz factor of the medium set to unity in the framework developed in this paper. Assuming that the individual components do not overlap, the authors restricted the study to a very small number of minijets, a problem that is overcome in this paper.

Unlike Clausen-Brown & Lyutikov (2012), we did not study observables such as the power spectral density or the variation time scales observed in light curves. These quantities depend on the detailed mechanism responsible for the generation of the

minijets, and their study is beyond the scope of this paper. Our model, based on a Doppler-boosted beamed emission, can certainly be generalized to other mechanisms, such as the “kinetic beaming” studied by Cerutti et al. (2012), which will probably yield similar conclusions on the statistical properties of the flux, while providing refined energy-dependent observables.

The generalization of such stochastic models, where many sub-regions are considered, growingly reduces the dominance of shock emission models in high-energy astrophysics. The springboard of such a change of paradigm points out the need for high-quality non-thermal light curves and widening studies of flux distributions, of the relation between the sample rms and the sample flux, and of the power spectral density.

Acknowledgments

The authors are grateful to the anonymous referee for the insights that have improved this article.

References

All Figures

Fig. 1 Schematic view of the geometry. The left side corresponds to the observer frame ℛobs, in which the jet is tilted by angle θ from the line of sight (along xobs) and is boosted by a Lorentz factor Γ (velocity Σ) along ${\mathit{x}}_{\mathrm{jet}}^{\mathrm{\prime }}$. The minijet is defined by its Lorentz factor γ (velocity β), in the jet frame ℛjet (top right). The orientation of the minijet along ${\mathit{x}}_{\mathrm{minijet}}^{\mathrm{\prime }}$ compared to the jet axis in its rest frame, xjet, is defined by the spherical angles ψ and ϕ. The emission is assumed isotropic in the minijet frame ℛminijet (bottom right).

In the text

	Fig. 2 Doppler factor of a minijet with a Lorentz factor γ = 5 moving in a jet of Lorentz factor Γ = 5 as a function of the misalignment between the jet and the line of sight θ, normalized to the opening angle of the jet 1/Γ. The maximum, average, and minimum Doppler factor are compared to the usual Doppler factor (solid line), derived for an isotropic emission in the jet frame (γ = 1).
In the text

Fig. 3 Distribution of the logarithm of the intensity of a minijet for various angles θ between the line of sight and the jet axis. The intensity is proportional to δ3 + s, where δ is the Doppler factor and s the spectral index. These distributions are obtained simulating random orientations of a minijet for s = 1, for a jet boost Γ = 5, and a minijet boost γ = 5, without generality loss. While the intensity decreases when the misalignment increases, the slope of the distribution remains the same as in the “blazar” case, θ = 0.

In the text

	Fig. 4 Distribution of the logarithm of the intensity of N independent and randomly oriented minijets. The number of minijets N increases from left to right with N ∈  { 1,10,30,102,3 × 102,103,3 × 103,104 } . Even for a large number of regions, asymmetrical, tailed distributions are obtained.
In the text

	Fig. 5 Distribution of the logarithm of the flux of N minijets for N = 3 × 103 (left) and N = 104 (right), as in Fig. 4. The continuous black and grey dashed lines represent the best fit with a log-normal and normal flux distributions, respectively.
In the text

	Fig. 6 Sample rms as a function of the sample flux of N = 1 minijet (top) and the sum of N = 104 minijets (bottom). Linear relations are found in both cases, with a zero x-intercept in the first case and a positive one in the second.
In the text