STOCHASTIC ACCELERATION AND THE EVOLUTION OF SPECTRAL DISTRIBUTIONS IN SYNCHRO-SELF-COMPTON SOURCES: A SELF-CONSISTENT MODELING OF BLAZARS' FLARES

In the statistical picture, the change in energy of the particles at each acceleration step n_s is expressed as

$$\begin{equation} \gamma _{n_s} =\varepsilon _{n_s} \gamma _{n_s-1} = \gamma _{n_s-1}(1+\Delta \gamma _{n_s-1}/\gamma _{n_s-1}), \end{equation} \tag{ 1 }$$

where γ is the Lorentz factor of the particle and ε is the fractional energy gain. Here we investigate the role of fluctuations of ε on the spectral shape of the accelerated particles. With this aim in mind, we express the energy gain fluctuations as

$$\begin{equation} \varepsilon = \bar{\varepsilon } +\chi, \end{equation} \tag{ 2 }$$

where the random variable χ has a probability density function with zero mean value (〈χ〉 = 0) and variance σ²_χ, and $\bar{\varepsilon }$ represents the systematic energy gain, which we treat as a non-random variable and the probability density function of ε is defined on the range ε ⩾ 0. The particle energy at step n_s can be expressed as

$$\begin{equation} \gamma _{n_s} =\gamma _{0} \Pi _{i=1}^{n_s }\varepsilon _{i}, \end{equation} \tag{ 3 }$$

where γ₀ is the initial energy of the particle. This equation clearly shows that the final energy distribution (n(γ)  =  dN(γ)/dγ) will result from the product of the random variables ε_i. The determination of an analytic expression for the distribution resulting from the multiplication of generic random variable is not an easy task (Glen et al. 2004). Using the simplifying assumption that the particles are always accelerated, namely, the acceleration probability, P_a, is set to unity and applying the multiplicative case of the central limit theorem (e.g., Cowan 1998), it is possible to show that the particle energies will be distributed as a log-normal law

$$\begin{equation} n(\gamma)=\frac{N_0}{\gamma \sigma _{\gamma } \sqrt{(}2\pi)} \exp \big [-(\ln \gamma - \mu)^2/2\sigma _{\gamma }^2\big ], \end{equation} \tag{ 4 }$$

where N₀ is the total number of particles, μ = 〈lnγ〉, σ²_γ = σ²(lnγ). We can determine these two quantities by taking the logarithm of Equation (3),

$$\begin{eqnarray} \ln \gamma _{n_s} &=&\ln \gamma _{0} + \Sigma _{i=1}^{n_s} \ln (\bar{\varepsilon }+ \chi _i) \nonumber \\ &=& \ln (\gamma _{0}\bar{\varepsilon }^{n_s})+ \Sigma _{i=1}^{n_s} \ln \Big (1 +\frac{\chi _i}{\bar{\varepsilon }} \Big) \nonumber \\ &\approx &\ln (\gamma _{0}\bar{\varepsilon }^{n_s})+\Sigma _{i=1}^{n_s} \left (\frac{\chi _i}{\bar{\varepsilon }}- \frac{\chi _i^2}{2\bar{\varepsilon }^2 } \right), \end{eqnarray} \tag{ 5 }$$

assuming that $\chi _i/\bar{\varepsilon }$ is not large. We obtain for the two parameters μ and σ_γ

$$\begin{eqnarray} \mu &=&{{\rm ln}} (\gamma _{0}) +n_s \ln \bar{\varepsilon } + n_s\left [ \left\langle \frac{\chi }{\bar{\varepsilon }}\right\rangle - \frac{1}{2}\Big (\frac{\sigma _{\chi }}{\bar{\varepsilon }}\Big)^2 -\left\langle \frac{\chi }{2\bar{\varepsilon }}\right\rangle ^2 \right ]\nonumber \\ \sigma _{\gamma }^2 &=& {n_{\rm s}} \Big [ \Big (\frac{\sigma _{\chi }}{\bar{\varepsilon }}\Big)^2+ \Big (\frac{\sigma _{\chi }}{2\bar{\varepsilon }}\Big)^4+ 2\Big (\frac{\sigma _{\chi }}{2\bar{\varepsilon }}\left\langle \frac{\chi }{2\bar{\varepsilon }}\right\rangle \Big)^2 \Big ], \end{eqnarray} \tag{ 6 }$$

where we have ignored the covariance terms since we are assuming the energy gain at each acceleration step is independent of the one at the previous step. Remembering that 〈χ〉 = 0, σ_χ = σ_ε, and ignoring the fourth-order term, we can write

$$\begin{eqnarray} \mu &=&\ln (\gamma _{0}) + n_s\Big [ \ln \bar{\varepsilon }- \frac{1}{2}\Big (\frac{\sigma _{\varepsilon }}{\bar{\varepsilon }}\Big)^2 \Big ].\nonumber\\ \sigma _{\gamma }^2 &\approx & n_s\Big (\frac{\sigma _{\varepsilon }}{\bar{\varepsilon }}\Big)^2. \end{eqnarray} \tag{ 7 }$$

This equation shows that the variance increases linearly with the number of acceleration steps and is proportional to σ_ε². Substituting μ and σ_γ into Equation (4)

$$\begin{eqnarray} n(\gamma)&=&\frac{N_0}{\gamma \sigma _{\gamma } \sqrt{(}2\pi)}\nonumber\\ &&\times \exp \left [ \frac{{-}\big (\ln ({\gamma }/{\gamma _0}) - n_s\big [\ln \bar{\varepsilon }- \frac{1}{2} ({\sigma _{\varepsilon }}/{\bar{\varepsilon }})^2 \big ]\big)^2}{2n_s ({\sigma _{\varepsilon }}/{\bar{\varepsilon }})^2} \right ].\quad \end{eqnarray} \tag{ 8 }$$

Hereafter we will consider decimal logarithms (log ≡ log₁₀, c_e = 1/log₁₀e ≈ 2.3) to make easier a comparison of the curvature results from this paper with those presented in observational papers. Taking the logarithm of Equation (8) and substituting the parameters from Equation (8) we obtain

$$\begin{eqnarray} \log n(\gamma) &=& K - \log \gamma\nonumber\\ && - \frac{\big (c_e\log {\gamma }/{\gamma _0} - n_s\big [c_e\log \bar{\varepsilon }- \frac{1}{2} ({\sigma _{\varepsilon }}/{\bar{\varepsilon }})^2 \big ]\big)^2}{c_e2n_s ({\sigma _{\varepsilon }}/{\bar{\varepsilon }})^2},\quad \end{eqnarray} \tag{ 9 }$$

where K includes all the constant factors. This is a log-parabolic law with the curvature (second degree in log γ) coefficient given by

$$\begin{eqnarray} && r=\frac{c_e}{2n_s ({\sigma _{\varepsilon }}/{\bar{\varepsilon }})^2}. \end{eqnarray} \tag{ 10 }$$

The interesting physical insight of this equation is that the curvature of the particle energy distribution is inversely proportional to the acceleration steps (n_s) and to the variance of the energy gain (σ²_ε). In the case of P_a < 1, the distribution at step n_s will be given by the convolution of different log-normal distributions for each acceleration step, with the distribution at n_s broader than that at n_s − 1 and containing fewer particles, as already noted in Peacock (1981).

Similar results are obtained considering a constant energy gain but a fluctuating number of acceleration steps. Assuming that after a time t the probability distribution for the number of steps undergone by a particle is given by a Poisson law, it is possible to show that the energy distribution follows a log-parabola whose curvature term depends on the inverse of the mean number of steps multiplied by the duration of the acceleration process.

The above statistical description provides an intuitive link between the curvature in the energy distribution of accelerated particles and the presence of a randomization process, such as the dispersion in the energy gain or in the number of acceleration steps. However, this approach does not give a complete physical description of the processes responsible for the systematic and stochastic energy gain, ignoring other physical processes, such as the radiative cooling and injection rates, or the acceleration energy dependence, necessary to give a complete description of the particles energy distribution evolution. A physical self-consistent description of stochastic acceleration in a time-dependent fashion can be achieved through a kinetic equation approach. Employing the quasi-linear approximation with the inclusion of momentum-diffusion term (Ramaty 1979; Becker et al. 2006), the equation governing the temporal evolution of n(γ) is

$$\begin{eqnarray} \frac{\partial n(\gamma,t)}{\partial t}&=& \frac{\partial }{\partial \gamma } \Big \lbrace -[S(\gamma,t) + D_A(\gamma,t) ]n(\gamma,t)\nonumber\\ && +D_{p}(\gamma,t) \frac{\partial n(\gamma,t)}{\partial \gamma } \Big \rbrace -\frac{n(\gamma,t)}{T_{{\rm esc}}(\gamma)}+Q(\gamma,t),\qquad \end{eqnarray} \tag{ 11 }$$

where D_p(γ, t) is the momentum-diffusion coefficient, D_A(γ, t) = (2/γ)D_p(γ, t) is the average energy change term resulting from the momentum-diffusion process, and S(γ, t) = −C(γ, t) + A(γ, t) is an extra term describing systematic energy loss (C) and/or gain (A), and Q(γ, t) is the injection term. In the standard diffusive shock acceleration scenario, there are several possibilities for which one can expect that energy gain fluctuations will occur, due to the momentum-diffusion term. In particular, for the case of a turbulent magnetized medium, the advection of particles toward the shock due to pitch angle scattering may be accompanied by stochastic momentum-diffusion mechanism. In this scenario, particles embedded in a magnetic field with both an ordered (B₀) and turbulent (δB) component, exchange energy with resonant plasma waves, and the related diffusion coefficient is determined by the spectrum of the plasma waves. Following the approach of Becker et al. (2006), we describe the energy distribution W(k) in terms of the wave number k = 2π/λ with a PL

$$\begin{eqnarray} &&W(k)=\frac{\delta B(k)^{2}}{8\pi }=\frac{\delta B(k_{0})^{2}}{8\pi }\left(\frac{k}{k_{0}}\right)^{-q}, \end{eqnarray} \tag{ 12 }$$

with q = 2 for the "hard-sphere" spectrum, q = 5/3 for the Kolmogorov spectrum, and q = 3/2 for the Kraichnan spectrum, the total energy density in the fluctuations being

$$\begin{eqnarray} &&U_{\delta B}=\int _{k_{0}}^{k_{\rm max}}W(k)dk. \end{eqnarray} \tag{ 13 }$$

Under these assumptions, the momentum-diffusion coefficient reads (O'Sullivan et al. 2009)

$$\begin{eqnarray} &&D_{p}\approx \beta _{A}^2\left (\frac{\delta B}{B_0}\right)^2 \Big (\frac{\rho _g}{\lambda _{{\rm max}}}\Big)^{q-1}\frac{p^2c^2}{\rho _g c}, \end{eqnarray} \tag{ 14 }$$

where β_A = V_A/c and V_A is the Alfvén waves velocity, ρ_g = pc/qB is the Larmor radius, and λ_max is the maximum wavelength of the Alfvén waves spectrum. The acceleration time for particles with Lorentz factor γ, whose Larmor radii resonate with one particular magnetic field turbulence length scale, is dictated by the momentum-diffusion coefficient (D_p) as

$$\begin{eqnarray} &&t_{\rm acc}\approx \frac{p^{2}}{D_{p}}= \frac{\rho _{g}(\gamma _{0})}{c \beta _{A}^{2}} \left.\left(\frac{B_{0}^{2}}{\delta B^{2}}\right)\right|_{\gamma _{0}}\left(\frac{\gamma }{\gamma _{0}}\right)^{2-q}. \end{eqnarray} \tag{ 15 }$$

The spatial diffusion coefficient relates to the momentum-diffusion coefficient through the relation, D_xD_p ≈ p²β²_A (Skilling 1975), hence the escape time of the particles from the acceleration region of size R depends on the spatial diffusion coefficient through the relation

$$\begin{eqnarray} &&t_{\rm esc}\approx \frac{R^{2}}{D_{x}}\approx \frac{R^{2}}{\left(c\beta _{\rm A}\right)^{2}t_{\rm acc}}. \end{eqnarray} \tag{ 16 }$$

The coefficients in Equation (11), and their related timescales, can be expressed as a PL in terms of the Lorentz factor (γ)

$$\begin{equation} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}D_{p}(\gamma)&\!\!{=}\;\;D_{p0}\left(\frac{\gamma }{\gamma _{0}}\right)^{q},\quad t_D =\frac{1}{D_{p0}}\left(\frac{\gamma }{\gamma _{0}}\right)^{2-q} \\ D_A(\gamma)&\!\!{=}\;2D_{p0}\left(\frac{\gamma }{\gamma _{0}}\right)^{q-1},\quad t_{DA} = \frac{1}{2D_{p0}}\left(\frac{\gamma }{\gamma _{0}}\right)^{2-q}\\ A(\gamma)&\!\!{=}\;A_{p0}\gamma, t_{A} =\frac{1}{A_{0}}\\ \end{array}\right., \end{equation} \tag{ 17 }$$

where D_p0 and A₀ have the dimension of the inverse of a time. Analytical solutions of the diffusion equation for relativistic electrons have frequently been discussed in the literature since the early work by Kardashev (1962), in particular for the case of the "hard-sphere" approximation. Neglecting the S and T_esc terms in Equation (11), and using a mono-energetic and instantaneous injection (n(γ, 0) = N₀δ(γ − γ₀)), the solution of the diffusion equation is (Melrose 1969; Kardashev 1962)

$$\begin{equation} n(\gamma,t)= \frac{N_0}{\gamma \sqrt{4\pi D_{p0} t}} \exp { \left \lbrace -\frac{ [\ln (\gamma /\gamma _0)-(A_{p0} - D_{p0}) t ]^2 }{4D_{p0}t} \right \rbrace, } \end{equation} \tag{ 18 }$$

i.e., a log-parabolic distribution, whose curvature term is

$$\begin{equation} r=\frac{c_e}{4D_{p0} t } \propto \frac{1}{D_{p0}t}. \end{equation} \tag{ 19 }$$

This result is fully consistent with that found in the statistical description; indeed, Equations (18) and (8) have the same functional form in both the statistical and in the diffusion equation scenario, with t playing the role of n_s, D_p0 the role of the variance of the energy gain (σ²_ε), and A_p0 the role of $\log \bar{\varepsilon }$. Hence we can write

$$\begin{equation} D_{p0}\propto \Big (\frac{\sigma _\varepsilon }{\bar{\varepsilon }}\Big)^2. \end{equation} \tag{ 20 }$$

It is interesting to note that in the case of the "hard-sphere" approximation, the curvature term is simply dictated by the ratio of the diffusive acceleration time (t_D) to the evolution time (t).

Under the "hard-sphere" approximation (q = 2), the spatial diffusion coefficient does not depend on the particle energy, since the exponent of Equation (15) is q − 2 = 0. Only three independent parameters exist in this description: the scattering time, the escape time, and the velocity of the scatterers. The spectra are purely determined by how many scatterings have been able to occur, the velocity of the scatterer, and what fraction of the injected particles has escaped out of the acceleration region. The scattering time relates to the spatial diffusion coefficient by t_scat ≈ D_x/c. Similarly, the resulting acceleration time relates to the spatial diffusion coefficient by t_acc ≈ D_x/β²_Ac ≈ t_A/β²_A. Thus, for "hard-sphere" turbulence, the scattering and acceleration timescales are independent of the particle energy (since there is equal energy density of scatterers which particles of all energies may resonantly scatter with).

The left panel in Figure 1 shows the resulting instantaneous evolution of spectra for the "hard-sphere" turbulence. The log-parabolic shape is maintained along the entire acceleration process, as shown by the solid lines representing the fit of the MC distributions by means of the law in Equation (9) (Cowsik & Sarkar 1984). The evolution of the curvature parameter, obtained from the r in the log-parabolic fit and plotted in Figure 2 with the red dashed line, clearly shows the trend due to the momentum diffusion, in agreement with the prediction from Equation (19) (blue line in the plot) demonstrating the connection between D_p0, ${\sigma _\varepsilon }/{\bar{\varepsilon }}$, and β_A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To account for the effects of turbulent magnetic field spectra softer than the "hard-sphere" case, we also consider acceleration in Kolmogorov- and Kraichnan-type turbulence. We have to therefore include a fourth parameter in the MC simulation, in addition to the three considered above: the turbulent field spectral slope q (see Equation (12)).

The right plot in Figure 1 shows the evolution of spectra for the "Kolmogorov" turbulence case. Similar spectra were obtained by Lemoine & Pelletier (2003) and O'Sullivan et al. (2009) who integrated the trajectories of charged particles in a turbulent magnetic field embedded in a fluid. The results are compared to the quasi-linear theory results (Becker et al. 2006; solid lines in Figure 1). We can identify two phases in the temporal evolution. In the first phase, the SEDs are more symmetric and the curvature evolves as in the q = 2 case, while in the second phase they develop a low-energy PL tail. Figure 2 shows that, for the Kolmogorov (green line) and the Kraichnan (black line) turbulence, r is systematically larger than the "hard-sphere" case (red line), and that for t ≳ 2 × t_acc, r approaches an asymptotic value (r ≈ 1.2 and r ≈ 1.5 for q = 5/3 and q = 3/2, respectively) ruled by the exponential cutoff in the equilibrium distribution.

We study the evolution of n(γ) and of the curvature term in a homogeneous spherical geometry, with radius R and an entangled coherent magnetic field B and a turbulent component δB, in the two cases of impulsive and continuous injection with a quasi mono-energetic source function Q(γ_inj, t) normalized to have a fixed energy input rate:

$$\begin{equation} L_{{\rm inj}}=\frac{4}{3}\pi R^3\int \gamma _{{\rm inj}} m_ec^2 Q(\gamma _{{\rm inj}},t)d\gamma _{{\rm inj}}\; ({\rm erg \,s^{-1}}). \end{equation} \tag{ 25 }$$

In our approach, we do not distinguish the acceleration region from the radiative one and during the acceleration process we take into account both synchrotron and IC cooling. According to Equation (14), to determine the order of magnitude of D_p we assume 1 ≫ δB/B ≃ 0.1–0.01 and require Alfvén waves to be at least mildly relativistic, with β_A ≃ 0.1–0.5, and their maximum wavelength to be much smaller than the accelerator size (λ_max < R). To study the effect of IC cooling on the evolution of n(γ), we consider two different sizes of the acceleration region, a compact one (R = 5 × 10¹³ cm) and a larger one (R = 1 × 10¹⁵ cm). With this choice of accelerator size, we set λ_max ≈ 10¹² cm. We stress that the choice of λ_max constrains the accelerative upper limit through ρ_g < λ_max leading to γ_max < (λ_maxqB)/m_ec², since particles with larger ρ_g (hence larger γ) cannot resonate with shorter wavelengths. Taking into account a coherent magnetic field of the order of 0.1 G and λ_max ≈ 10¹² cm we found that the purely accelerative efficiency limits the particle energy to γ_max ≲ 10^7.5. In the left panel of Figure 3, we plot t_D, given by Equation (17), as a function of λ_max, for the case of q = 2, δB/B = 0.1, and β_A = 0.5. In this case, the acceleration time is energy independent and for λ_max ≈ 10¹² cm it will be of the order of t_D = 1/D_p0 ≈ 10⁴ s. In the case of q ≠ 2, the acceleration will have an energy dependence given by Equation (17), as shown in the right panel of Figure 3 for the case of q = 3/2. In this section, we focus on the evolution of the curvature as a function of the momentum-diffusion term, and therefore use only the accelerative contributions coming from the diffusion terms (D_p(γ), D_A(γ)), neglecting the systematic extra term A(γ). All the parameters and their numerical values are given in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Parameters' Values Adopted in the Numerical Solutions of the Diffusion Equation for the Cases Studied in Section 4

Parameter		Impulsive Inj.		Cont. Inj.
R	(cm)	5 × 1013, 1 × 1015	...	...	...
B	(G)	0.1, 1.0	...	...	...
Linj	(erg s−1)	1039	...	1037	...
q		2	3/2	2	3/2
$t_{D_0}=1/D_{P0}$	(s)	1 × 104	1 × 103	1 × 104	1 × 103
Tinj	(s)	100	...	1 × 104	...
Tesc	(R/c)	∞	...	2	...
Duration	(s)	1 × 105	...	...	...
γinj		10.0	...	10.0	...

Download table as:  ASCIITypeset image

In the left panels of Figure 4 and Figure 5, we plot the evolution of energy distribution n(γ, t) (upper panels) and of γ³n(γ, t) (lower panels) in the case of the impulsive injection without escape, for q = 2, and for two values of R: 1 × 10¹⁵ cm (Figure 4) and 5 × 10¹³ cm (Figure 5). We inject a quasi-monoenergetic electron distribution with γ_inj ≈ 10. The γ³n(γ, t) representation is useful to compare the results concerning n(γ) presented in this section, with those regarding the synchrotron emission presented in Section 5. We denote by γ_p the peak energy of n(γ) and by r the curvature evaluated by means of a log-parabolic best fit over a one decade-broad interval centered at γ_p. γ_3p and r_3p represent the peak of γ³n(γ) and its curvature, respectively. In the right panels of Figures 4 and 5, we report on the corresponding temporal evolutions of the curvatures under the effect of both momentum diffusion and cooling terms. The solid black line corresponds to t = 0.2 × t_acc, where $t_{{\rm acc}}=t_{D_0}$ is the acceleration time due to momentum diffusion. As the time increases, the diffusion term acts on the distribution by means of both D_A and D_p. The effect of the latter is to make the distribution broader.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One can distinguish three phases: in the first one the energy of particles increases and the curvature parameter decreases following a law r∝t⁻¹ in agreement with the statistical scenario of Section 2 and with the Equation (19), independent of the magnetic field strength (B = 1.0 G and B = 0.1 G) and of the source size, because the accelerative contribution dominates over the radiative losses; in the second phase, the radiation losses become relevant and the distribution approaches the equilibrium with an increase of the curvature; and in the third phase, the balance between acceleration and radiation losses is established and the curvature reaches a stable value.

The equilibrium distribution reached through stochastic acceleration, is described by a relativistic Maxwellian (Schlickeiser 1985; Stawarz & Petrosian 2008),

$$\begin{equation} n(\gamma)\propto \gamma ^2 \exp {\left [\frac{-1}{f(q,\dot{\gamma })} \Big (\frac{\gamma }{\gamma _{{\rm eq}}}\Big)^{f(q,\dot{\gamma })}\right ]}, \end{equation} \tag{ 26 }$$

where $f(q,\dot{\gamma })$ is a function depending on the exponent of the diffusion coefficient and on the cooling process and γ_eq is the Lorentz factor that satisfies the condition t_cool(γ) = t_acc(γ) and is given by

$$\begin{equation} \gamma _{{\rm eq}}=\frac{1}{t_{{\rm acc}} C_0(U_B+F_{{\rm KN}}(\gamma))}, \end{equation} \tag{ 27 }$$

with t_acc equal to the fastest acceleration timescale among t_A, t_D, and t_DA. In the case of Compton-dominated cooling we have γ_eq∝(R²/t_accB²f_KN), while in the case of strong KN regime, or in general for synchrotron-dominated cooling, we have γ_eq∝(1/t_accB²). Using a PL form for the acceleration terms, and in the case of only synchrotron losses (or any cooling process that can be expressed as a PL function of γ), it is possible to give an analytic expression of $f(q,\dot{\gamma })$ (Katarzyński et al. 2006; Stawarz & Petrosian 2008). The expectation for synchrotron and IC/TH cooling process and for q = 2 is $f(q,\dot{\gamma })=3-q=1$. The curvature resulting from a log-parabolic fit over a decade centered on γ_p is r ≈ 2.5 and r_3p ≈ 6.0 in the case of γ_3p.

We first discuss the case of R = 10¹⁵ cm (Figure 4) with only synchrotron cooling (dashed lines, left panels). In terms of behavior, we note that for the larger value of B (1.0 G; red lines, right panels), the r–t trend departs from the purely accelerative one (r∝t⁻¹; green lines, right panels) early (relative to the B = 0.1 G case; blue lines in the right panels). This happens because the synchrotron equilibrium energy (vertical dot-dashed lines, left panels) is lower in the case of B = 1.0 G. For both values of B, the final values of r are close to the synchrotron equilibrium value of ≈2.5. When IC cooling is also taken into account, the final values of the curvature in n(γ) are r ≈ 2.5 and r ≈ 0.6 for B = 0.1 G and B = 1.0 G, respectively. This difference is due to the different IC cooling regimes for the two cases. To show clearly the complexity of the transition from the TH to the KN regime, and its dependence on R and B, in Figure 6 we plot the ratio $\dot{\gamma }_{{\rm IC}}/\dot{\gamma }_{{\rm Synch.}}$ (solid lines), and U_ph/U_B (dashed lines), as a function of γ and normalized to unity, for the case of q = 2, for the final step of the temporal evolution. As long as the ratio U_ph/U_B is close to $\dot{\gamma }_{{\rm IC}}/\dot{\gamma }_{{\rm Synch.}}$, electrons cool in the full TH regime, and C(γ) = C₀γ²(U_B + U_ph). On the contrary, when the electrons radiate in the full KN regime $\dot{\gamma }_{{\rm IC}}/\dot{\gamma }_{{\rm Synch.}}\ll U_{{\rm ph}}/U_B$. In this case, due to the inefficient KN cooling regime we have $\dot{\gamma }_{{\rm Synch.}}\gg \dot{\gamma }_{{\rm IC}}$, and the cooling term is dominated by the synchrotron component: C(γ) ≈ C₀γ²U_B. In the intermediate cases, it is difficult to estimate analytically the ratio $\dot{\gamma }_{{\rm IC}}/\dot{\gamma }_{{\rm Synch.}}$.

Zoom In Zoom Out Reset image size

For B = 1.0 G, the SSC equilibrium is reached at γ ≈ 3 × 10⁴ and the SSC cooling occurs between the KN and TH regimes (see the top right panel in Figure 6), hence the value of f is different from unity, as predicted for the case of full IC/TH or synchrotron cooling. When B = 0.1 G, the equilibrium energy is γ ≈ 10⁷ and electrons are in extreme KN cooling (see the top left panel in Figure 6), synchrotron losses are much higher than those due to IC scattering, and again r reaches the previous value of ≈2.5. It is also interesting to note the difference in the trends of r–t and r_3p–t. In the latter case, the trend departs from the purely accelerative regime earlier (see Figure 4, lower right panel) since the electrons with energies close to γ_3p are more energetic than those close to γ_p, and thus have much shorter cooling times.

The results for the compact region (R = 5 × 10¹³ cm) are plotted in Figure 5. Considering that the injected electron luminosity is the same (see Table 1), we expect a different response from the IC cooling, due to the higher photon density n_ph(₀). The r evolution for the synchrotron cooling case is similar to the previous case, while for the SSC emission, both for the case of B = 1.0 G and B = 0.1 G, the final value of r is about 2.5. This is due to the larger photon density which moves the IC scattering into the TH regime also for the case of B = 0.1 G (compare bottom left to top left panels in Figure 6), hence n(γ) approaches the solution of Equation (26) in the case of f = 1.

In Figure 7, we show the temporal evolution for the case of q = 3/2 (R = 1.0 × 10¹⁵ cm, B = 0.1 G). In this case, contrary to the q = 2 case, the acceleration time t_D is energy dependent, hence we study the evolution of r as a function of t/t_D(γ_inj), where t_D(γ_inj) is the diffusive acceleration time evaluated at the injection energy γ_inj. The energy dependence of t_D affects the evolution of r, and the shape of the equilibrium distribution, indeed, the r–t and r_3p–t trends show different behavior compared to the case of q = 2. The equilibrium curvature is reached for t ≳ 1 × t_D(γ_inj), and the two equilibrium curvature values are r ≈ 1.2 and r_3p ≈ 3.0, roughly half of those found for the case of q = 2, and in agreement with the result from the MC. We note that the curvature obtained by means of a log-parabolic fit of Equation (26), for the case q = 3/2 (namely, f = 1.5), is r ≈ 3.7. Hence, both the MC and the numerical solution of the diffusion equation give a result different from that predicted by the analytical solution in Equation (26).

Zoom In Zoom Out Reset image size

The case of continuous injection (see Figure 8) is more complex. The distribution developes a low-energy PL tail, but a log-parabolic bending, driven by the diffusion, is still present at high energies, hence we evaluate the curvature only at γ_3p (i.e., the representation useful to compare it to the synchrotron SED emission). Spectral curvatures are generally milder than the impulsive injection. In the left panel of Figure 8, we plot the r–t trend both for the case of impulsive (red lines) and continuous (blue lines) injections, the curvature in the continuous injection case are systematically lower in the pre-equilibrium phases and in the acceleration-dominated stage the trend is again consistent with the "hard-sphere" approximation and statistical approaches. The slope of the electron distribution in the PL tail is ≈1.06, in good agreement with the predicted one ≈1 + t_min-acc/(2t_esc) = 1.075, consistent with the results of Katarzyński et al. (2006).

Zoom In Zoom Out Reset image size

We compute the evolution of b_s and b_c, as a function of the time, for the case of $t_{D_0}=1.5\times 10^4$ s, B = 0.1 G, and q = 2, using a temporal mesh of 2 s. We plot in the top left panel of Figure 9 the instantaneous SEDs at steps of 200 s: the solid lines represent the synchrotron and IC SEDs averaged over the full duration of the acceleration process (10⁴ s). As the time is increased, the peak energy of both the synchrotron and IC SEDs moves toward higher energies with a broadening of the spectral distribution. The corresponding evolution of curvature parameters is reported in the top right panel: b_s has a trend similar to that of the electron distribution, with b_s∝(t/t_D0)^−α and α ≃ 0.6 (for comparison the cyan solid line represents the r_3p/5 trend, as predicted by the S δ-approximation). The trend of b_c, as expected, is more complex because of the transition from TH to KN regime. For t/t_acc ≲ 0.4, it follows the same trend of b_s but with systematically lower values. For t/t_acc ≳ 0.4, when the TH–KN transition occurs, b_c increases with time, approaching toward the electron curvature r value. This transition starts for values of E_s ≈ 5 × 10⁻³ keV (ν_s ≈ 10¹⁴ Hz) and E_c ≈ 0.05 GeV (ν_c ≈ 10²² Hz); and the corresponding SEDs are plotted by blue thick-dashed lines in the left panel of Figure 9.

Zoom In Zoom Out Reset image size

In the bottom panels of Figure 9, we show the case of q = 3/2. The synchrotron curvature quickly approaches the equilibrium value of b_s ≈ 0.6, consistent with the equilibrium value r_3p ≈ 3.0 discussed in Section 4.2. In this case we do not observe the TH/KN transition in the IC curvature, since the lower values of E_s and E_c keep the IC scattering mainly in the TH regime.

The other parameter affecting the evolution of the spectral distributions is the diffusion coefficient D_p0 (see Equation (15)), which we assume to vary in the range [1.5 × 10⁴, 2.4 × 10⁵] s⁻¹, studying how the main spectral parameters change. In the top left panel of Figure 10, we plot averaged SEDs for each different value of D_p0. The top right panel shows the trend of b_c versus D_p0. As expected, for larger values of D_p0, the curvature measured at the peak energy is smaller. The trend is described by a PL with an exponent of about −0.6 for D_p0 ≲ 2 × 10⁻⁵ s⁻¹ and with an exponent of about −0.25 for D_p0 ≳ 2 × 10⁻⁵ s⁻¹. This break clearly shows the transition between the TH and KN regimes (marked by a vertical dashed line); indeed it happens for the same values of D_p0 corresponding to the TH/KN transition in both the D_p0–b_c trend and the E_c–b_c plot (occurring at E_c ≈ 1 GeV; see the bottom right panel in Figure 10). The break in the D_p0–b_s trend happens when electrons radiating at E_s enter the KN cooling region, hence, due to the lower cooling level (compared to the TH cooling regime, on the left side of the vertical dashed line), the curvature decreases.

Zoom In Zoom Out Reset image size

Blue filled circles in the top left panel represent the peak positions for both SED components. For the synchrotron component, according to Equation (29), the exponent α in the case of n(γ₃) = const, should be 1.5, while the results of the computations give α = 0.6. This difference is due to the fact that we inject in the mono-energetic initial distribution always the same total power that corresponds to the same number of particles. When the peak energy increases the distribution becomes broader, implying that the same total number of particles is spread over a larger energy interval and the number of particles contributing to the synchrotron peak emission decreases. Consequently, the S_s–E_s trend gets softer compared to the predicted value of 1.5.

We verified quantitatively this effect by computing the trend n(γ_3p) versus γ²_3p, and found a PL relation with an exponent of about 0.98, in nice agreement with the difference between the exponent of 1.5 and that resulting in our simulations. In the bottom panels of Figure 10, we plot b_s versus E_s (left) and b_c versus E_c (right). The S_c–E_c relation can be fitted by a PL (orange line, top left panel in Figure 10) with the same exponent of the E_s–S_s relation, as long as the IC scattering, at E_c and above, happens in TH regime. When the KN suppression becomes relevant (green line, top left panel in Figure 10), the exponent is larger and is close to unity.

The synchrotron trend (the bottom left panel in Figure 10) clearly shows the expected anti-correlation between the peak energy and the spectral curvature, which is well fit by the function given in Equation (32), with a = 0.68, not very different from 0.6, obtained for the δ-function approximation of the synchrotron emission, and assuming that n(γ) has a purely log-parabolic shape. A simple PL fit of the same points returns an exponent −0.14.

We also investigate the role of q on the spectral evolution, setting its variation range to [3/2, 2], i.e., from the Kraichnan to the "hard-sphere" case. The relations between the spectral parameters are very similar to those found in the previous case with S_s∝E^0.6_s and S_c∝E^0.9_c. Also in this case, the synchrotron component follows the expectation with a lower curvature for harder turbulence spectra, and the IC trend shows the transition from TH to KN regime. The PL best fit of S(E_s) versus E_s gives a = 0.88, larger than that obtained for the case of D_p0. In fact, for values of q lower than 2, corresponding to less turbulence and hence diffusion, the curvature gets higher values and the peak energy lower values, compared to the "hard-sphere" case. The PL fit for b_s versus E_s returns an exponent of −0.16, practically coincident with the previous one, indicating that the average properties of these parameters are the same in both the q and D_p0 cases.

The magnetic field B drives the radiative losses which affect the evolution of the spectral parameters. In Section 4, we showed that different cooling conditions, and the transition from TH to KN, can determine very different values of γ_eq for the same acceleration conditions. Assuming that the acceleration timescale is independent of the magnetic field, Equation (27) shows that γ_eq∝1/B², implying that, as long as B is small enough to result in γ_3p ≪ γ_eq, the evolution of n(γ) around the peak value is dominated by the acceleration terms, while for values of B resulting in γ_3p ≳ γ_eq the evolution obtains a notable contribution due to cooling. In the top left panel of Figure 11, we plot the averaged SEDs. According to Equations (28) and (29), the synchrotron peak value should scale as S_s∝(E_s)². Indeed, for values of B ≲ 0.2 G we obtain an exponent equal to 2.04, very close to the value found with the δ-approximation. For higher values of the magnetic field, E_s is anti-correlated with S_s. This behavior represents a cooling signature due to the decreasing value γ_eq for increasing B values, with γ_eq getting closer to γ_3p. This is confirmed both by the shape of the synchrotron SEDs and by the b_s–B plot (top right panel in Figure 11). Indeed, S SEDs for B ≳ 0.2 G exhibit an exponential decay, meaning that the distributions have reached, or are close to reaching, the equilibrium energy. Consistently with the S shape evolution, the b_s–B relation shows an almost stable value of b_s for B ≲ 0.2 G and an increasing trend for B ≳ 0.2 G. This change, in both the S_s–E_s and b_s–B trends, is interesting and can provide a useful phenomenological tool for understanding the evolution of non-thermal sources. Another interesting feature is shown in the S_c–E_c plot: for B ≲ 0.2 G the IC peak energy is practically constant, as expected in the KN limit from the kinematical limit relating the scattered photons energy to that of the electrons: hν_IC ≈ γm_ec². In fact, photons at energies ≈E_c are produced in the KN regime and for B ≲ 0.2 G the electron peak energy γ_3p is constant, so E_c must also be constant. For B ≳ 0.2 G, γ_3p decreases because of cooling, and, accordingly, E_c also decreases. This is another interesting test that can provide a probe for B-driven flares evolving to the KN regime. The E_s–b_s plot in the bottom left panel of Figure 11 confirms the cooling signature discussed above, showing b_s uncorrelated with E_s as long as γ_3p ≪ γ_eq, and an increasing value of b_s with E_s almost stable, when γ_3p ≳ γ_eq.

Zoom In Zoom Out Reset image size

The E_s–b_s trend, and in particular the anticorrelation between these two observables parameters, is the strongest signature of a stochastic component in the acceleration.

In Figure 12, we report the scatter plot in the E_s–b_s plane for the six considered sources. The left panel reports the results obtained by changing the value of D_p0: the green dashed lines describe the trend resulting from a log-parabolic fit of the synchrotron SED over a decade in energy centered on E_s; the purple lines represent the same trend obtained by fitting a log-parabola in the fixed spectral window [0.5, 100.0] keV. Both these trends are compatible with the data and track the predicted anticorrelation between E_s and b_s. Purple data, however, give a better description, hinting that the "window" effect could be a real bias. Each of the three lines was computed for a different value of the magnetic field. It is remarkable that the variation of a single parameter, D_p0, can describe the observed behavior. The dispersion in the data is relevant and can be related to the variation of B (as partially recovered by numerical computation), or by different values of the beaming factor, R, and L_inj, during different flares, and for different objects.

Zoom In Zoom Out Reset image size

The dot-dashed thick line represents the best fit of the observed data by means of Equation (32), and returns a value of a ≈ 0.6, as expected from theoretical predictions for the case of the δ-approximation, and pure log-parabolic electron distribution. This fitted line is also compatible with the numerical trend shown by the purple lines. Note that the observed curvature values are in the range [0.1, 0.5], corresponding to r_3p ∼ [0.5, 3.0]. According to the results presented in Section 4.2, the expected equilibrium curvature in the synchrotron emission, in the full KN or TH regime, and for q = 2, should be of r_3p ≈ 6.0 and of r_3p ≈ 5.0 in the intermediate regime. In the case of q = 3/2, the equilibrium curvature should be r_3p ∼ 3.0. This is perhaps an interesting hint that, both in the flaring and the quiescent states, for q = 2, the distribution is always far from equilibrium. In the case of q = 3/2, only for E_s ≲ 1.5 keV is the curvature compatible with the equilibrium (r_3p ≃ 3.0, corresponding to b_s ∼ 0.6). For larger values of E_s, we find again curvature well below the equilibrium value. These results provide a good constraint on the values of the magnetic field B ≲ 0.1 G.

The q-driven trend (right panel) is also compatible with the data, but for values of E_s ≲ 1 keV, the D_p0-driven case seems to describe better the observed behavior, but any firm conclusion is not possible because of the dispersion of the data.

As a last benchmark for the stochastic acceleration model, we reproduce the observed correlation between E_s and S_s, which follows naturally from the variations of D_p0 and q. Considering that the redshifts of the six considered HBL objects are different, we prefer to use their peak luminosity L_s = S_s4πD²_L, where D_L is the luminosity distance.⁴ To account for the different jet power of sources, we considered two data subsets, and we assumed L_inj = 5 × 10³⁹ erg s⁻¹ for the first subset (top panels of Figure 13) and L_inj = 5 × 10³⁸ for the second (bottom panels of Figure 13). In the left panels of Figure 13, we report the D_p0-driven trend and in the right panels we show the q-driven trend. Solid lines represent the trend obtained by deriving L_s from the log-parabolic best fit of the numerically computed SEDs, centered on E_s; dashed lines are the trends obtained by fitting the numerical results in the fixed energy window [0.5, 100] keV.

Zoom In Zoom Out Reset image size

Both results give a good description of the observed data and their shapes are similar. Solid lines follow well a PL with an exponent of about 0.6, while the windowed trends (dashed lines) show a break around 1 keV and the exponent below this energy turns to about 1.5. A similar break at the same energy can be noticed in the points of Mrk 421 in the E_s–S_s plot presented by Tramacere et al. (2009), who found an exponent of ∼1.1 and of ∼0.4 below and above 1 keV, respectively. This could again be an indication that the observed values are actually affected by the bias.