© ESO 2018

1 Introduction

Blazars are among the most energetic phenomena in nature, representing the most extreme type in the class of active galactic nuclei (Urry & Padovani 1995). They feature relativistic jets that extend over tens of kiloparsecs and are directed toward the general direction of Earth. Observations of their radiation emission show very high luminosities, rapid variabilities, and high polarizations. Moreover, apparent superluminal characteristics can be detected along the first few parsecs of the jets. The main components of blazar jets are magnetized plasmoids, which are assumed to arise in the Blandford−Znajek and the Blandford−Payne process (Blandford & Znajek 1977; Blandford & Payne 1982), constituting the major radiation zones. These plasmoids pick up – and interact with – particles of interstellar and intergalactic clouds along their trajectories (Pohl & Schlickeiser 2000), giving rise to the emission of a series of strong flares.

A blazar spectral energy distribution (SED) consists of two broad nonthermal radiation components in different domains. The low-energy spectral component, ranging from radio to optical or X-ray energies, is usually attributed to synchrotron radiation of relativistic electrons that are subjected to ambient magnetic fields. The origin of the high-energy spectral component, covering the X-ray to γ-ray regime, is still under debate. It can, for instance, be modeled by inverse Compton radiation coming from low-energy photon fields that interact with the relativistic electrons (Böttcher 2007; Dermer et al. 1997). This process can be described either by a synchrotron self-Compton (SSC) model (see, for example, Schlickeiser & Röken 2008; Röken & Schlickeiser 2009b and references therein), where the electrons scatter their self-generated synchrotron photons, or by so-called external Compton models like Dermer & Schlickeiser (1993), Sikora et al. (1994), Blazejowski et al. (2000), where the seed photons are generated in the accretion disk, the broad-line region, or the dust torus of the central black hole. Some models also consider ambient fields of IR radiation from diffuse, hot dust (see, for example, Sikora et al. 2001). Aside from these leptonic scenarios, the high-energy component can also be modeled via proton-synchrotron radiation or the emission of γ-rays arising from the decay of neutral pions formed in interactions of protons with ambient matter (see, for example, Böttcher et al. 2013; Weidinger & Spanier 2015; Cerruti et al. 2015 and references therein). Mixed models including both leptonic as well as hadronic processes are also considered in the literature (for example, Cerruti et al. 2011; Yan & Zhang 2015).

A major task is to properly understand and to account for the distinct variability patterns of the nonthermal blazar emission at all frequencies with different time scales ranging from years down to a few minutes, where the shortest variability time scales are usually observed for the highest energies of the spectral components, as in PKS 2155-304 (Aharonian et al. 2007, 2009) and Mrk 501 (Albert et al. 2007) in the TeV range, or Mrk 421 (Cui 2004) in the X-ray domain. So far, only multi-zone models, which feature an internal structure of the emission region with various radiation zones caused by collisions of moving and stationary shock waves, have been proposed to explain the extreme short-time variability of blazars (see, for example, Sokolov et al. 2004; Marscher 2014; Graff et al. 2008; Ghisellini et al. 2009; Giannios et al. 2009; Biteau & Giebels 2012). In the framework of one-zone models, however, extreme short-time variability can, a priori, not be realized as the duration of the injection into a plasmoid of finite size (with characteristic radius ) and the light crossing (escape) time in this region are naturally on the order , c being the speed of light in vacuum and the bulk Doppler factor (Eichmann et al. 2012). This leads to a minimum time scale for the observed flare duration, which may exceed the short-time variability scale in the minute range by several magnitudes. In particular, using the typical parameters and , we find in the observerframe. Thus, this type of model can only be used to account for variability on larger time scales, that is, from years down to hours.

Both one-zone and multi-zone models have been studied extensively over the last decades in order to explain the variability but also the SEDs of blazars, incorporating leptonic and hadronic interactions. In these studies, analytical as well as numerical approaches were used, containing, among others, various radiative loss and acceleration processes, details of radiative transport, diverse injection patterns, different cross sections for particle interactions, as well as particle decay and pair production/annihilation (see the review article Böttcher 2010, and the references therein). Thus, the literature on this matter is quite comprehensive. However, models featuring pick-up processes of any kind usually assume only a single injection of particles into the emission region as the cause of the flaring. This may be unrealistic as blazar jets, which may extend over tens of kiloparsecs, most likely intersect with several clouds, leading to multiple injections. In such an intricate situation, it is of particular interest to have a self-consistent description of the particle’s radiative cooling, the radiative transport, et cetera. Therefore, we propose a simple, but fully analytical, time-dependent leptonic one-zone model featuring multiple uniform injections ofnonlinearly interacting ultrarelativistic electron populations, which undergo combined synchrotron and SSC radiative losses (for previous works see Röken & Schlickeiser 2009a,b). More precisely, we assume that the blazar radiation emission originates in spherically shaped and fully ionized plasmoids that feature intrinsic randomly-oriented, but constant, large-scale magnetic fields and propagate ultrarelativistically along the general direction of the jet axis. These plasmoids pass through – and interact with – clouds of the interstellar and intergalactic media, successively and instantaneously picking up multiple mono-energetic, spatially isotropically distributed electron populations, which are subjected to linear, time-independent synchrotron radiative losses via interactions with the ambient magnetic fields and to nonlinear, time-dependent SSC radiative losses in the Thomson limit. This is the first time an analytical model that describes combined synchrotron and SSC radiative losses of several subsequently injected, interacting injections is presented (for a small-scale, purely numerical study on multiple injections see Li & Kusunose 2000). We point out that because the SSC cooling is a collective effect, that is, the cooling of a single electron depends on the entire ensemble within the emission region, injections of further particle populations into an already cooling system give rise to alterations of the overall cooling behavior (Röken & Schlickeiser 2009a; Zacharias & Schlickeiser 2013). Moreover, as we do employ Dirac distributions for the time profile of the source function of our kinetic equation and, further, do not consider any details of radiative transport, we mimic injections of finite duration times and radiative transport by partitioning each flare into a sequence of instantaneous injections, which are appropriately distributed over the entire emission region, using the quantities computed here.

The article is organized as follows. In Sect. 2, an approximate analytical solution of the time-dependent, relativistic kinetic equation of the volume-averaged differential electron number density is derived. Based on this solution, the optically thin synchrotron intensity, the SSC intensity in the Thomson limit, as well as the corresponding total fluences are calculated in Sects. 3 and 4. We explain how finite injection durations and radiative transport can be mimicked by multiple use of our results in Sect. 5, and show a parameter study for the total synchrotron and SSC fluence SEDs. Section 6 concludes with a summary and an outlook. Supplementary material, which is required for the computations of the electron number density and the synchrotron and SSC intensities, is given in Appendices A–G. Moreover, in Appendix H, we briefly describe the plotting algorithm employed for the creation of the fluence SEDs.

2 The relativistic kinetic equation

The kinetic equation for the volume-averaged differential electron number density n = n(γ, t) (where t is the time, γ:= E_e∕(m_e c²) the normalized electron energy, and [n] = cm⁻³) of m ultrarelativistic, mono-energetic, instantaneously injected and spatially isotropically distributed electron populations in the rest frame of a nonthermal radiation source with dominant magnetic field self-generation and radiative loss rate L = L(γ, t) (with [L] =s⁻¹) reads (Kardashev 1962) (1)

where δ(⋅) is the Dirac distribution,q_i (for which [q_i] = cm⁻³) the ith injection strength, γ_i:= E_e,i∕(m_e c²) ≫ 1 the ith normalized initial electron energy, and t_i the ith injection time for i: 1 ≤ i ≤ m. In this work, radiative losses in form of both a linear, time-independent synchrotron cooling process (with a constant magnetic field B) and a nonlinear, time-dependent SSC cooling process in the Thomson limit (2)

are considered. The respective cooling rate prefactors yield D₀ = 1.3 × 10⁻⁹ b² s⁻¹ and A₀ = 1.2 × 10⁻¹⁸ b² cm³ s⁻¹, which depend on the nondimensional magnetic field strength b:= ∥B∥∕Gauss = const. (Blumenthal & Gould 1970; Schlickeiser 2009). In the present context, the term linear refers to the fact that the synchrotron loss term does not depend on the electron number density n and therefore results in a linear contribution to the partial differential equation (PDE) (1), whereas the SSC loss term depends on an energy integral containing n, yielding a nonlinear contribution. The synchrotron loss term can, in principle, be modeled as nonlinear and time-dependent as well. This can be achieved by making an equipartition assumption between the magnetic and particle energy densities (Schlickeiser & Lerche 2007). We point out that except for TeV blazars, where Klein–Nishina effects drastically reduce the SSC cooling strength above a certain energy threshold, the dominant contribution of the SSC energy loss rate always originates in the Thomson regime (Schlickeiser 2009), which justifies our initial restriction. As a consequence, however, this limits the applicabilityof our model to at most GeV blazars and bounds the normalized initial electron energies from above by γ_i < 1.9 × 10⁴ b^−1∕3. We also note that the only accessible energy in the synchrotron and SSC cooling processes is the kinetic energy of the electrons. Thus, γ denotes the kinetic component of the normalized total energy E_tot.∕(m_ec²), that is,

Then, γ_tot. ∈ [1, ∞) implies γ ∈ [0, ∞). For this reason, the lower integration limit in (2) is zero. Furthermore, the Dirac distributions δ(γ − γ_i) and δ(t −t_i) in Eq. (1) determine a sequence of energy and time points . They are to be understood in the distributional sense, that is, each is a linear functional on the space of smooth test functions φ on with compact support

More precisely, for , one can rigorously define them as the mapping

with the integral of the Dirac distribution against a test function given by

2.1 Formal solution of the relativistic kinetic equation

In terms of the function R(γ, t):= γ² n(γ, t) and the variable x:= 1∕γ, we can rewrite Eq. (1) in the form (3)

where (4)

Defining the strictly increasing, continuous function G = G(t) via (5)

Eq. (3) becomes (6)

with G_i:= G(t_i). Although the nonlinear kinetic equation given in (1) is now transformed into a linear PDE, its solution can obviously serve as a Green’s function only for the single-injection scenario with m = 1. Consequently, in order to solve the generalized kinetic equation

for m > 1, where Q = Q(γ, t) is a more realistic source function, one cannot simply use Green’s method. Applying successive Laplace transformations with respect to x and G to Eq. (6) yields the solution (7)

where H(⋅) is the Heaviside step function

with , which has a jump discontinuity at k = k₀. A detailed derivation of this solution can be found in Appendix A. In order to determine the function G, we substitute Sol. (7) into (4) and, by using (5), obtain the ordinary differential equation (ODE) (8)

This equation can be regarded as a compact notation for the set of m piecewise-defined ODEs (9)

In Sect. 2.2, we present an approximate analytical solution for the general case of j injections, that is, for the interval G_j ≤ G < G_j+1, where j ∈{1,..., m} and G₁ = 0, G_m+1 = ∞, employing the method of matched asymptotic expansions. Having a separate analytical solution for each ODE of (9), in Sect. 2.3, these are connected successively requiring continuity at the transition points. Finally, by substituting (7), the electron number density results in (10)

2.2 Solution of the relativistic kinetic equation for { t t_j ≤t < t_j+1 with j ∈{1,..., m}}

In the interval G_j ≤ G < G_j+1, Eq. (8) gives (11)

In order to solve this equation, we introduce the time-dependent sets S₁, S₂, and S₃ that contain indices i: 1 ≤ i ≤ j belonging to either injections in the near-injection domain (NID) 0 ≤ G − G_i ≪ x_i, the intermediate-injection domain (IID) G − G_i ≈ x_i, or the far-injection domain (FID) G − G_i ≫ x_i, respectively. We may now rewrite Eq. (11) in terms of these sets by continuing their domains of validity to for the NID, for the IID, and for the FID with the NID-IID and IID-FID transition values and , where are positive constants, taking into account that the (j + 1)th injection can occur in each domain. However, for the computations below, we still use the former much less than, approximately equal, and much greater than relations when we carry out approximations. Then, Eq. (11) can be represented formally by (12)

Next, we express the NID and the FID summands by the Taylor series of the function evaluated at the point z = 0 (generalized geometric series)

and theIID summands by the Taylor series of the same function evaluated at the point z = 1

and suitably approximate these series by considering only the leading- and next-to-leading-order terms. This results in

for all u ∈ S₁, v ∈ S₂, and w ∈ S₃. To find an approximate analytical solution to Eq. (12), we require the function G to be extricable from the respective sums. This is already achieved with the approximations (13) and (14). The latter approximation (15), however, needs further modifications as G still couples to G_w in the denominator. To this end, we employ both geometric and generalized geometric series and again consider only the leading- and next-to-leading-order terms (16)

The approximation errors and optimal values for the interval parameters and are given in Appendix B. Substituting (13), (14), and (16) into the ODE (12), we obtain (17)

with the “constants” (18)

which depend on the actual numbers of elements of S₁, S₂, and S₃. Note that the explicit dependences of these constants on the total number of injections as well as on (the numbers of elements of) the sets S₁, S₂, and S₃ are suppressed in the subsequent calculations for simplicity if possible and given if necessary. We derive an approximate analytical solution of Eq. (17) by computing separate solutions for the NID, IID, and FID of the jth injection, which, in case , are glued together continuously at the transition points and (that define the transitions of the jth injection between S₁ and S₂ as well as between S₂ and S₃) corresponding to the transition times and specified in Appendix C. For , we regard only the NID and IID solutions with the proper continuous gluing at , and for , we consider solely the NID solution. Moreover, we have to update the above constants (18) each time an injection i: 1 ≤ i ≤ j crosses over from its NID to its IID at or from its IID to its FID at , as these transitions cause changes in either the numbers of elements of the sets S₁ and S₂ or of S₂ and S₃. In the following, due to the particular structure of Eq. (17) and for reasons of analytical solvability, the general strategy is to use both the leading- and next-to-leading-order contributions only if S₁, or S₂, or S₁ and S₂, or S₃ has to be taken into account, whereas only the leading-order terms are considered if S₁ and S₃, or S₂ and S₃, or S₁ and S₂ and S₃ have to be employed.

2.2.1 Near-injection domain solution

In the NID , at least the jth injection is an element of S₁. Further, one may – but does not necessarily need to – have injections in S₂ and S₃. Therefore, we distinguish the following four cases

Since Eqs. (19) and (20) as well as Eqs. (21) and (22) are identical except for the present constants, only two different kinds of ODEs have to be solved. Starting with Eqs. (19) and (20), we use separation of variables in order to obtain

which isequivalent to

As this solution converges strictly against the value , which can be lower than the NID transition value , it is not suitable for covering this domain. Including a second-order term in the ODE under consideration, that is, working with the equation (23)

where , leads to a solution in form of a tangent function, which may have several poles in the NID, rendering it useless as well. Adding further contributions to (23) yields analytically non-solvable ODEs. Hence, we have to employ the simpler ODEs

with the linear solution (26)

and . For Eqs. (21) and (22), we also apply separation of variables and extend the integrand with the multiplicative identity and the neutral element , leading to (27)

where . One way of deriving an approximate analytical solution of this transcendental equation is to determine G asymptotically for small and large arguments of the arctan function in case both asymptotic ends exist for and extrapolate these solutions up to an intermediate transition point, where they are connected requiring continuity, such that the entire domain of definition is covered. For small arguments , we use the third-order approximation , as the linear term in Eq. (27) and the linear contribution in the approximation of the arctan function cancel each other out. This leads to the asymptotic solution (28)

For large arguments , the arctan function can be well approximated by π∕2. We directly obtain an asymptotic solution of the form (29)

By means of the condition , we derive the NID transition value (30)

This value specifies the upper bound of the domain of validity of (28) and the lower bound of the domain of validity of (29). The corresponding transition time can be directly computed by substituting (30) into (27), in which the constant c₃ had to be fixed via the initial condition G(t = t_j) = G_j resulting in

This yields

We point out that in general, one does not have the ordered sequence because can in principle be larger than or equal to as well as smaller thanor equal to G_j. These cases arise when only one asymptotic end of the arctan function exists. Hence, for , the solution of Eq. (27) is approximately given by (31)

for by (32)

and for by (33)

The integration constants c₂ and c₄,..., c₇ are determined via the proper initial and transition conditions in Appendix D. In order to compute the constants , , , , and , we have to continually check duringthe evolution of whether an injection i: 1 ≤ i ≤ j belongs to S₁, S₂, or S₃. Therefore, we specify two different kinds of updates. The first kind occurs when a new injection enters the system, while the second kind is due to either NID-IID or IID-FID transitions. We start with the initial update at the time of the jth injection, verifying

Next, since all transition values and are known, they can be arranged as an ordered 2j-tuple, representing an increasing sequence. Assuming that a certain number of these values is contained in the time interval under consideration, whenever G reaches one of them, all sets S₁, S₂, and S₃ have to be updated accordingly, meaning that the corresponding injection is moved from either S₁ to S₂ or from S₂ to S₃ and the constants involved change. For more details on the updating see Sect. 2.3 and Appendix D.

2.2.2 Intermediate-injection domain solution

In the IID , where at least the jth injection is in S₂, the previousj − 1 injections can reside in their respective NIDs, IIDs, or FIDs. Thus, we have to discuss the four cases

As these coincide structurally with the ODEs (24) and (25) as well as (21) and (22) of the NID, we can directly write down their solutions. For Eqs. (34) and (35), we get (38)

whereas for Eqs. (36) and (37), we obtain in case (39)

if (40)

and if (41)

where

is the IID transition time associated with the transition value

The derivation of the integration constants d₁,..., d₅ can be found in Appendix D.

2.2.3 Far-injection domain solution

In the FID , the jth injection is, per definition, an element of S₃. As before, the previous j − 1 injections can be in their respective NIDs, IIDs, or FIDs depending on the initial injection parameters, and therefore we have to evaluate thefour ODEs

Since Eqs. (43)–(45) are also structurally identical to the ODEs (21) and (22), we can once again directly write down their solutions, yielding for (46)

for (47)

and for (48)

where

is the transition time corresponding to the transition value

An approximate analytical solution of Eq. (42) can be obtained in a similar way as for the NID ODEs (21) and (22) by first applying separation of variables, resulting in an integral equation of the form (49)

This integral equation is shown to be approximately equivalent to a specific transcendental equation for which we subsequently derive asymptotic solutions forG that are continued up to an intermediate transition point (if both asymptotic ends exist, else we only consider a continuation of the solution of the present asymptotic end over the entire domain), where they are glued together continuously. In more detail, since , we employ a first-order geometric series approximation in the integrand of (49) (50)

Note that in the first step, we included the multiplicative identity D₀ ∕D₀ and the neutral element in the numerator of the integrand, giving the splitting into two terms. Then, by means of simple algebraic manipulations, we expressed the integrand in a form suitable for the substitution of the geometric series. For reasons of computational simplicity, we dropped the small second-order contribution in the last step. The final integral in (50) is solved by an arctan function. Thus, introducing the parameter , Eq. (49) results in (51)

In order to find an approximate analytical solution of this transcendental equation, we again determine G asymptotically for both small and large arguments of the arctan function. Therefore, using the third-order approximation of the arctan function for small arguments , the associated asymptotic solution becomes

Further, because , the asymptotic solution for large arguments reads

Extending the domains of these solutions up to – and connecting them continuously at – the transition point (52)

which is derived from the transition condition and corresponds to the transition time (53)

we obtain for the piecewise-defined solution (54)

In case only one asymptotic end exists, that is, for or , we get (55)

or (56)

respectively.We remark that the transition time (53) was computed by fixing the integration constant

in Eq. (51) imposing the initial condition and substituting the transition value (52). The determination of the integration constants e₁,..., e₉ is given in Appendix D.

2.3 Solution of the relativistic kinetic equation for {t

The complete solution of Eq. (8) is derived as follows. Beginning with the single-injection domain (SID), the solution branch G(t | 0 ≤ t < t₂) is given for by Sol. (26) (in which j = 1 as for the other SID solutions below) with , for by Sol. (38) with , and for by

• Sol. (54) for ,

• Sol. (55) for ,

• Sol. (56) for

with and . We point out that in the SID, Eq. (8) can also be solved by directly applying separation of variables, resulting in (Schlickeiser et al. 2010)

where and . Using once again the method of matched asymptotic expansions with the transition time

leads for to the approximate solution (57)

for to (58)

and otherwise for to (59)

The integration constants f₁, f₂, f₄, and f₅ are fixed by the initial condition G(t = t₁ = 0) = G₁ = 0

whereas the integration constant f₃ is determined by the transition condition

In the following, we use the more exact SID approximation (57)–(59). The double-injection domain (DID) solution branch G(t | t₂ ≤ t < t₃) is given for by Sol. (26) (with j = 2 as for the other DID solutions below) if and by

• Sol. (31) for ,

• Sol. (32) for ,

• Sol. (33) for

if S₂ = ∅, S₃≠∅ (with the specific arguments of the constants and associated to the respective cases as for the other DID solutions below). For , we employ Sol. (38) if and

• Sol. (39) for ,

• Sol. (40) for ,

• Sol. (41) for

if S₁ = ∅, S₃≠∅. Finally, for , we apply

• Sol. (46) for ,

• Sol. (47) for ,

• Sol. (48) for

if either S₁≠∅, S₂ = ∅ or S₁ = ∅, S₂≠∅ and

• Sol. (54) for ,

• Sol. (55) for ,

• Sol. (56) for

if S₁ = ∅ = S₂. The initial value G₂ is fixed by requiring continuity of the SID and DID solution branches at the time of the second injection, yielding

In addition to the initial DID updating of constants at t = t₂, we have to perform NID-IID updates at and/or if as well as IID-FID updates at and/or if . Assuming, for example, that the transition times , , and are contained in this interval with the order , we have to update the initial DID sets S₁, S₂, and S₃ thrice. More precisely, at the time of the second injection t = t₂, both the first and the second injection are contained in S₁, while S₂ and S₃ are empty. During the temporal progression toward the upper bound t₃, the first injection switches from S₁ to S₂ at , whereas the second injection continues to be in S₁. At , also the second injection switches over to S₂, leaving S₁ empty. Last, at , the first injection switches from S₂ to S₃. This amounts to the following updating sequence:

Repeating this procedure for the remaining m − 2 injections results in the formal representation of G for t: 0 ≤ t < ∞ given by (60)

The SID contribution G_SID, the NID, IID, and FID contributions G_NID, G_IID and G_FID, and the initial constants G_i are stated explicitly in Appendix E. We note that here the updating of constants is suppressed for readability. It could, however, be written down explicitly similar to the updating of the integration constants presented in Appendix D.

3 Synchrotron and synchrotron self-Compton intensities

In this section, we calculate the optically thin synchrotron intensity and the SSC intensity in the Thomson limit. The optically thick component of the synchrotron intensity is not considered because it was shown in Schlickeiser & Röken (2008) that, for all frequencies and times, it provides only a small contribution to the high-energy SSC component.

3.1 Synchrotron intensity

The optically thin synchrotron intensity I_syn.(ϵ, t) (with [I_syn.] = eV s⁻¹ cm⁻² sr⁻¹ and similar for the SSC intensity) for an isotropically distributed electron number density is given by (61)

where the function (62)

is the pitch-angle-averaged spectral synchrotron power of a single electron in a magnetic field of strength b, ϵ:=E_syn.∕(m_e c²) is the normalized synchrotron photon energy, P₀:= 8.5 × 10²³ eV s⁻¹, and ϵ₀:= 2.3 × 10⁻¹⁴ b (Blumenthal & Gould 1970). The CS function is discussed in detail in Appendix F. Here, we employ the approximation (63)

where a₀:= 1.15. Substituting the electron number density (10) and the synchrotron power (62) with the CS function (63) into formula (61), we obtain for the optically thin synchrotron intensity (64)

with and the abbreviation . For comparisons with observational data or for generic case studies, that is, for fitting or plotting lightcurves, we have to compute the energy-integrated synchrotron intensity

with a lower integration limit ϵ_min. corresponding to the energy of the first data point in the fluence SED and an upper integration limit ϵ_max. defined by the last. Using (64), this quantity yields

where and . In order to solve the integral, we approximate the integrand by for and for . This is justified because the approximated CS function (63) is only adapted to the asymptotics and of the exact CS function and extrapolated in between. Moreover, since τ_i,min. and τ_i,max. depend on the time t, we have to consider the three cases where 1 ≤ τ_i,min., τ_i,min. < 1 ≤ τ_i,max., and τ_i,max. < 1. Accordingly, the integral results in (65)

3.2 Synchrotron self-Compton intensity

In the computation of the SSC intensity (66)

where P_SSC(ϵ_s, γ, t) is the SSC power of a single electron and ϵ_s:= E_s∕(m_e c²) is the normalized scattered photon energy, we have to employ the Thomson limit because the SSC radiative losses in (2) are already restricted to the Thomson regime. In this limit, the SSC power reads (Felten & Morrison 1966; Blumenthal & Gould 1970)

with the Thomson cross section σ_T = 6.65 × 10⁻²⁵ cm² and the total SSC photon energy density (having the dimension ). For ultrarelativistic electrons with γ ≫ 1 and synchrotron photon energies in the Thomson regime for which γ ϵ ≪ 1, the characteristic energy of the SSC-scattered photons is ϵ_s ≈ 4 γ² ϵ (Felten & Morrison 1966), corresponding to head-on collisions of the synchrotron photons with the electrons (Jones 1968; Blumenthal & Gould 1970; Reynolds 1982). Thus, it is justified to apply a monochromatic approximation in the total SSC photon energy density in form of a Dirac distribution that spikes at this characteristic energy

where

is the synchrotron photon number density. We remark that by assuming an isotropic, ultrarelativistic electron distribution, the synchrotron photon number density becomes inevitably isotropically distributed as well. Substituting the latter formulas into (66) and using Fubini’s theorem, we obtain (67)

in which the upper γ-integration limit arises from the restriction to the Thomson regime. With the electron number density (10) and the synchrotron intensity (64), the SSC intensity (67) yields, after performing both integrations, (68)

where . The double sum, with the index i referring to the ith synchrotron photon population and the index j to the jth electron population, accounts for all combinations of SSC scattering between the various electron and synchrotron photon populations.For the corresponding energy-integrated SSC intensity, we find

where and . The integral is given by (65), however, with adapted integration limits.

4 Synchrotron and synchrotron self-Compton fluences

We compute the total fluences associated with the synchrotron intensity (64) and the SSC intensity (68). For this purpose, we derive a general expression for the total fluence that, on the one hand, employs the function G and, on the other hand, explicitly displays the various approximate cases of the Jacobian determinant of the integration measure. For simplicity, the updating of constants is once more suppressed. With (ε, I, F) ∈{(ϵ, I_syn., F_syn.), (ϵ_s, I_SSC, F_SSC)}, the total fluence F (for which [F] = eV cm⁻² sr⁻¹) is given by (69)

The Jacobian determinant yields (70)

where i: 2 ≤ i ≤ m. Substituting(70) into (69), we obtain the expression

in which

and

with the characteristic function

In the following, for illustrative purposes, we calculate in detail both the synchrotron and the SSC NID fluence integrals for the cases . The remaining integrals can be solved in an analogous manner.

4.1 Synchrotron fluence

With I = I_syn.(ϵ, G) according to (64), the synchrotron NID fluence integral for reads

Using gives (71)

Similar to the evaluation of the energy-integrated intensities (cf. the paragraph directly above formula (65)), we approximate the integrand for by and for by . Thus, solving the integral, (71) yields

with the generalized incomplete gamma function defined for and Re (a) > 0 by (Abramowitz & Stegun 1972)

In the special case a = 1∕2, the generalized incomplete gamma function can be expressed in terms of the error function, namely .

4.2 Synchrotron self-Compton fluence

With I = I_SSC(ϵ_s, G) given in (68), the SSC NID fluence integral for becomes

An approximate analytical solution of the integral may be derived with the same method as in the case of the synchrotron fluence. However, one could also employ the methods used for the computations of the NID-IID and IID-FID transition times, which are explained in more detail in Appendix C, as follows. Applying the mean-value-theorem method, we first define the function

Next, weassert that there exists a mean value such that

Last, we choose the midpoint of the integration interval as an approximate value for ξ_il. The trapezoid-approximation method, on the other hand, results in the expression

5 Lightcurves and fluence spectral energy distributions

To obtain realistic lightcurves, we have to consider injections of finite duration time and radiative transport inside the emission region. Because, for reasons of computational feasibility, we have only employed Dirac distributions for the time profile of our source function and did not concern ourselves with any details of radiative transport, we have to mimic both by modeling each flare via a sequence of suitably distributed, instantaneous injections. In more detail, we partition the ith flare into n_i separate injections labeled by a subscript p ∈{1,..., n_i}, which are induced equidistantly over the entire emission region given by , that is, the pth injection of the ith flare occurs at the time t_i + κ_ip with

and for each of which we use the quantities derived in Sect. 3. Thus, the energy-integrated synchrotron intensity reads (72)

where for all i: 1 ≤ i ≤ m is a specific distribution satisfying the constraint (73)

and the function I_ip(ϵ, t) is, according to the original synchrotron intensity (64), defined by (74)

with . Similarly, one may construct a proper formula for the energy-integrated SSC intensity. We refrain from showing a parameter study for lightcurves because an adequate value for n_i, which is at least on the order , results in numerical computations that would exceed our available CPU time by far. For the associated total fluence SEDs, it is an entirely different matter as we can illustrate their functional shapes by the special case n_i = 1 for all i: 1 ≤ i ≤ m, in which each flare is modeled by a single injection (cf. Sect. 4), given that it yields lower and upper bounds. In the following, we first show by direct calculation that the total fluence of the synchrotron intensity defined in (72) is indeed bounded from below and above, up to specific constant factors, by the simple n_i = 1 case. To this end, we substitute into the expression for the total fluence and employ the variable G (75)

From Eq. (11), it can be deduced that the Jacobian determinant is positive and bounded from below and above by (76)

and since I_ip is non-negative, we can therefore bound the fluence (75) from below and above by

Then, applying the transformation and definingI_i:= I_i1 (cf. formula(74)), we obtain

Because the integral has the same value for fixed i and arbitrary p, we can always drop the subscript p in the argument of the integrand and the integration measure. Hence, by means of the constraint (73), we get

for the lower bound and

for the upper bound. Transforming the integration measure back to the time variable t and using the inverse of (76) leads to

Since the second sum on both the left-hand side and the right-hand side can be bounded from above by , we find (77)

where , proving the claim.

In Figs. 1a–c, we show a parameter study for the total synchrotron and SSC fluence SEDs (given in arbitrary units), always considering three injections of ultrarelativistic electron populations into the emission region, each modeling one flare according to the estimate (77). We separately vary the nondimensional magnetic field strength b, the Doppler factor of the plasmoid , and the first reciprocal normalized initial electron energy x₁, leaving the remaining free parameters fixed. Details on the values of these parameters and the specific variations can be found in the figure caption. We note that in all three subfigures, the black, solid SEDs correspond to the same set of parameters and thus serve as a reference curve in the parameter study. Varying only the magnetic field strength (Fig. 1a), toward higher values, we can see an increase in both the maximum synchrotron energy and emissivity, as expected. This increase naturally leads to a reduction in the SSC emission, as the leftover SSC energy budget of the electrons becomes lower for larger and higher energetic synchrotron emission. However, the maximum SSC energy increases inasmuch as the maximum synchrotron energy can now reach higher values that add to the SSC scattering process. Furthermore, shifting the Doppler factor toward higher values results in the typical increase in the apparent luminosity (Fig. 1b), which is well known from studies of the physical effects of relativistic beaming, that is, from aberration, the Doppler effect, time dilation, and retardation. In Fig. 1c, as we vary the reciprocal normalized initial electron energy of the first injection, toward higher values, both the maximum synchrotron energy and the maximum SSC energy increase, which hence broadens the SED curves on their right slopes according to the specific nonlinear cooling behavior and the strength of the parameter variation. We note that we could also create a plot in which we vary the strengths q_i, i ∈ {1, 2, 3}, of the differentinjections. However, we refrain from showing such a plot because increasing the injection strengths toward higher values obliges us to raise the number of grid points (that is, the number of time steps, which is currently 4 × 10⁵, amounting to a computation time of approximately four days) in direct proportion, as the injection strengths and the time variable are coupled by multiplication. This would result in computation times ranging from several months to years. We furthermore remark in passing that in order to detect patterns that are due to variations in the number of injections or due to their nonlinear coupling, a larger parameter study than is presented here is necessary. Finally, but most importantly, the resulting functional shapes of the various plots indicate that our model can reproduce the typical fluence SEDs of blazars visible in observational data.

Fig. 1 Total synchrotron (left curves) and SSC (right curves) fluence SEDs for generic three-injection scenarios with magnetic field strengths, Doppler boost factors, and reciprocal initial electron energies b ∈ {0.01, 0.1, 1}, , and for panel a, b= 1, , and for panel b, as well as b = 1, , x1 ∈ {0.6 × 10−4, 10−4, 10−3}, and for panel c. The injection times are always given by with the characteristic radius of the plasmoid , the injection strengths by , and the cooling rate prefactors D0 = 1.3 × 10−9 b2 s−1, A0 = 1.2 × 10−18 b2 cm3 s−1 depend on the specific value of the magnetic field strength. The number of grid points used in the plotting of each curve amounts to 4 × 105.

6 Summary and outlook

We introduced a fully analytical, time-dependent leptonic one-zone model for the flaring of blazars that employs combined synchrotron and SSC radiative losses of multiple interacting, ultrarelativistic electron populations. Our model assumes several injections of electrons into the emission region as the cause of the flaring, which differs from common blazar models, where only a single injection is considered. This is, from a physical point of view, more realistic since blazar jets may extend over distances on the order of tens of kiloparsecs, and it is thus most likely that there is a pick up of more than just one particle population from interstellar and intergalactic clouds. At the same time, it furthermore assumes the two radiative cooling processes to occur simultaneously, as would be the case in any physical scenario. In more detail, applying Laplace transformations and the method of matched asymptotic expansions, we derived an approximate analytical solution of the relativistic kinetic equation of the volume-averaged differential electron number density for several successively and instantaneously injected, mono-energetic, spatially isotropically distributed, interacting electron populations, which are subjected to linear, time-independent synchrotron radiative losses and nonlinear, time-dependent SSC radiative losses in the Thomson limit. Using this solution, we computed the optically thin synchrotron intensity, the SSC intensity in the Thomson limit, and the corresponding total fluences. Moreover, we mimicked finite injection durations and radiative transport by modeling flares in terms of sequences of instantaneous injections. Ultimately, we presented a parameter study for the total synchrotron and SSC fluence SEDs for a generic three-injection scenario with variations of the magnetic field strength, the Doppler factor, and the initial electron energy of the first injection, showing that our model can reproduce the characteristic broadband SED shapes seen in observational data. We point out that the SSC radiative loss term considered here is strictly valid only in the Thomson regime, limiting the applicability of the model to at most GeV blazars. Nonetheless, it can be generalized to describe TeV blazars by using the full Klein–Nishina cross section in the SSC energy loss rate. This leads to a model for which similar yet technically more involved methods apply. Furthermore, in order to make our simple analytical model more realistic, terms accounting for spatial diffusion and for electron escape could be added to the kinetic equation. Also, more elaborate source functions, for example, with a power-law energy dependence, a time dependence in form of rectangular functions for finite injection durations, and with a proper spatial dependence may be considered. However, judging from the complexity of our more elementary analysis, this is most likely only possible via direct numerical evaluation of the associated kinetic equation.

Acknowledgements

We are grateful to Horst Fichtner, Reinhard Schlickeiser, and Michael Zacharias for useful discussions and comments. Furthermore, we thank the anonymous referee for helpful and constructive suggestions.

Appendix A Laplace transformation method

We derive the solution R(x, G) of the PDE (6) using a composition of two Laplace transformations. First, applying a Laplace transformation with respect to the reciprocal electron energy x

to Eq. (6) gives

Evaluating the second integral on the left-hand side via integration by parts and employing the Laplace transform of R (A.1)

we obtain

and, hence, (A.2)

As the normalized initial electron energies are finite and bounded from above by γ_i < 1.9 × 10⁴ b^−1∕3 for all i: 1 ≤ i ≤ m (due to the restriction to the Thomson regime) and only radiation loss processes are considered, we know that the electron number density has support

Thus, we can write n(γ, t) = H(γ_max. − γ) n(γ, t), yielding for the second term on the left-hand side of Eq. (A.2)

where x_max.:= 1∕γ_max.. Accordingly, we find (A.3)

Second, applying a Laplace transformation with respect to the function G

to Eq. (A.3) results in (A.4)

This Laplace transformation implies that G is non-negative. With the Laplace transform of K

and integration by parts as before, Eq. (A.4) reads

Since the Heaviside functions are not well-defined at the jump discontinuities at G = G_i for i: 1≤ i ≤ m, this equation reduces to (A.5)

containing an unspecified contribution in the last term on the right-hand side, which is handled as follows. At G = 0, no energy losses have yet occurred. Therefore, the electron number density is of the form n(x, 0) = n₁ δ(x − x₁), where . Substituting this density into (A.1) by employing the relation R(x, G) = n(x, G)∕x², we get

Choosing , we can use the last term on the right-hand side of Eq. (A.5) in order to compensate for this function, yielding

We point out that both M and the sum are positive. Hence, s > −w, which allows us to apply the inverse Laplace transformations with respect to the variables s and w to M. This gives the solution of Eq. (6)

Appendix B Approximation errors and optimal domains

The errors of the NID, IID, and FID approximations (13), (14), and (16) yield

Evaluating the NID and IID errors at the NID-IID transition point and the IID and FID errors at the IID-FID transition point , we obtain

and

We derive optimal values for the interval parameters and by requiring that the total error at each transition point becomes minimal. To this end, we first impose the sufficient conditions (B.1)

which result in the solution and an equation for the zeros of a seventh-order polynomial in given by (B.2)

respectively. We verify that the total error has a minimum at by means of the second derivative test

From the Abel–Ruffini theorem, we know that there exists no general analytical solution to Eq. (B.2). Hence, one has to employ a root-finding algorithm to compute numerical approximations.

One may obtain a more accurate overall approximation of Eq. (12) using the more general IID approximation

with errors

where the unspecified expansion point a_v constitutes another degree of freedom. We note that our original IID approximation (14) corresponds to expansion points that were set to a_v = 1 for all v ∈ S₂. Then, one can determine the optimal values for all a_v such that the total errors

become minimal by making the first and second derivate tests

under the constraints (B.1) (and the associated second derivate tests). To guarantee that the minima do not coincide with expansion points at infinity, one has to further impose a constraint that limits the potential values for the expansion points to a suitable finite interval. We finally remark that, due to the complexity of the resulting equations, this procedure isfeasible only numerically. Moreover, by numerical trial and error with standard blazar parameter values, we found out that the special case (14) yields an adequate approximation close to the optimal values.

Appendix C NID-IID and IID-FID transition times

In the following, we present two different methods for the determination of the transition times and . To this end, it is advantageous to start from the integral equation representation of the ODE (11) evaluated at the respective transition point or . In the former NID-IID case, we obtain (C.1)

The first method applies the first mean value theorem for integration. Thus, as , there exists a point such that the transition time (C.1) can be written as

Becausethe mean value theorem is merely an existence theorem, we have to approximate the value of p_j by means of an additional input. Since J⁻¹ is strictly increasing, apossible choice for p_j is the midpoint of the integration interval. The second method employs a trapezoid approximation. Again due to the strictly increasing functional shape of J⁻¹, the integralin (C.1) can be approximated by the area A(j) of a trapezoid, which is computed as the sum of the area A_r(j) of the rectangle defined by the distance between the endpoints G_j and of the integration interval and the height J⁻¹(G_j)

and the area A_t(j) of the right-angled triangle with one cathetus given by the distance between the endpoints and the other one by the height

This leads to the formula

for the transition time. To obtain a more accurate approximation, one could use additional supporting points in the interval , giving rise to a finer trapezoid decomposition of the integral. We note that in the present study, the trapezoid method yields a better approximation of the transition time. The IID-FID transition time

can be computed similarly, resulting, on the one hand, in the expression

for the mean-value-theorem method and, on the other hand, in the expression

for the trapezoid method.

Appendix D Integration constants – initial and transition conditions and updating

The integration constants c₂ and c₄, ..., c₇ of the NID (see Sol. (26) and Sols. (31)–(33)), d₁ , ..., d₅ of the IID (see Sols. (38)–(41)), as well as e₁, ..., e₄ and e₆, ..., e₉ of the FID (see Sols. (46)–(48) and Sols. (54)–(56)) are determined. The NID integration constants c₂, c₄ , c₆ , and c₇ are fixed via the initial condition G(t = t_j) = G_j, whereas c₅ is fixed via the transition condition . The initial condition provides continuity of G at the transition from the (j − 1)th injection domain to the jth injection domain at t = t_j and the transition condition between the nonlinear and linear NID solution branches at . Similarly, the IID integration constants d₁, d₂ , d₄ , and d₅ are determined via the respective initial condition and d₃ via the transition condition , where the former condition yields continuity between the NID and IID solution branches at and the latter between the nonlinear and linear IID solution branches at . Last, the FID integration constants e₁, e₃ , e₄ , e₆ , e₈ , and e₉ are specified by the initial condition and the integration constants e₂ and e₇ by the transition conditions and , respectively. These conditions guarantee continuity between the IID and FID solution branches at and between the nonlinear and linear FID solution branches at and . Furthermore, the integration constants have to be updated if for i ∈{1, ..., j − 1}, as elements of S₁ switch over to S₂ and/or elements of S₂ switch over to S₃. Because these updates cause shifts in the NID, IID, and FID transition times , , , and toward higher values, we have to account for updating conditions that cover transitions to other solution branches. Below, the determination of the NID integration constant c₄ for the case S₂≠∅≠S₃ is shown in detail. The remaining integration constants are computed accordingly. Applying the initial condition to the first branch of Sol. (31), we obtain

As c₄ depends just on S₃, we have to consider only IID-FID updates. Hence, the conditions – and the updated constants – for all possible IID-FID transitions of the ith injection at are

• first branch Sol. (31) → first branch Sol. (31)

• first branch Sol. (31) → Sol. (32)

Since S₃ ∪{i}≠∅ and , it is obvious that transitions from the first branch of (31) to (26), from the first to the second branch of (31), and from the first branch of (31) to (33) are not possible.

Appendix E Components of G(t | 0 ≤ t < ∞) and initial values G_i

We specify the expressions for the SID, NID, IID, and FID components G_SID, G_NID, G_IID, and G_FID of Sol. (60). First, the SID component simply yields

The more elaborate NID component is given by

where

and

with k ∈{1, 2, 3}, is the characteristic function. For the IID and FID components, we obtain similar expressions. On the one hand, we have

where

and on the other hand, we have

where

as well as

Furthermore, we determine the initial values G_i for all i : 3 ≤ i ≤ m by requiring continuity between the solution branches of the (i − 1)th and ith injection domains. If the ith injection enters the system while the (i − 1)th injection is still in the NID, the initial values become in case

and in case

If, however, the (i − 1)th injection is in the IID, they result for in

whereas for , they yield

Finally, when the (i − 1)th injection is already in the FID, we find in case S₁≠∅, S₂ = ∅ or

and in case S₁ = ∅ = S₂

Appendix F CS function

For , the CS function is defined by (Crusius & Schlickeiser 1988) (F.1)

where K_a is the modified Bessel function and W_{a, b} denotes the Whittaker function (Abramowitz & Stegun 1972). On account of the degree of complexity of (F.1), an approximate function that is adapted to its asymptotics

where a₀ := 1.15, is usually employed. Standard approximations are therefore given by

Figure F.1 shows the absolute values of the relative deviations of these approximations with respect to (F.1) in percent, that is,

Fig. F.1 Absolute values of the relative deviations of the approximate functions CSn for all n∈{1, 2, 3} with respect to the exact CS function in percent.

Since the CS function is primarily used for the computation of the synchrotron intensity, where the spectrum covers the energy range from radio waves to X-rays, that is, with a lower energy limit on the order neV and an upper limit on the order keV, and z = 2 ϵ∕(3 ϵ₀ γ²) with ϵ₀ ~ 10⁻¹⁴ b, initial electron energies γ_i ~ 10⁴ b^−1∕3, and nondimensional magnetic field strength b ~ 10⁻³ up to b ~ 10, we have to consider values z ≪ 1 up to z ≫ 1. Thus, the only suitable choice is the approximate function CS₃(z), which coincides with both asymptotic ends of (F.1) and has the smallest overall deviation.

Appendix G Relativistic beaming

Accounting for relativistic beaming, the energy ε ∈{ϵ, ϵ_s} and the time t, which are defined in the plasmoid rest frame, are related to the corresponding quantities in an observer frame (denoted with an asterisk) by

where

is the boost factor, the Lorentz factor of the plasmoid, β := v∕c, and θ^⋆ the angle between the jet axis and the line of sight of the observer. For blazars, one can assume that θ^⋆ ↘ 0 and thus

Furthermore, because the ratio I∕ε³ is Lorentz-invariant, meaning that

one directly finds that the intensity and the fluence transform as

Appendix H Plotting algorithm for the fluence spectral energy distributions

The numerical implementation of G, the synchrotron and SSC intensities, and the corresponding total fluence SEDs is carried out with Python. Here, we describe the functionality of the algorithm¹ and the specific incorporation of the analytical formulas. The algorithm makes heavy use of the decimal package² , which is designed for high floating-point precision calculations. This is necessary because the range of possible values between the synchrotron and SSC cooling rate prefactors D₀ and A₀ , the injection strengths q_i , and the reciprocal initial electron energies x_i spans several orders of magnitude, which may lead to a loss of accuracy in expressions where these parameters come up. This problem occurs in its most severe form in the evaluation of the formulas for the transition times , , , and , yielding deviations from the expected values by more than 50% with standard floating point precision. Even with higher precision, however, it is preferable to avoid the evaluation of the transition times altogether. Thus, we compute G on a grid, using fixed time steps. This makes it possible to read out the values of G at the grid points and directly compare them to the values of , , , and in order to determine the actual solution branch without referring to the transition times. The algorithm begins with the definitions of the free parameters, namely the injection times t_i , the injection strengths q_i , and the reciprocal initial electron energies x_i for i : 1≤ i ≤ m, the nondimensional magnetic field strength b, the synchrotron and SSC cooling rate prefactors D₀ and A₀, the Lorentz boost of the plasmoid, and the upper time boundary of the grid t_end. The value oft_end is chosen as one and a half times the injection time of the final injection for a multiple-injection scenario or given by a sufficiently large value for a single-injection scenario. In a realistic setting, however, t_end corresponds to the end of the observation time. The time grid is set up homogeneously and linearly, and the number of grid points can be chosen arbitrarily. All computations are performed in the plasmoid rest frame. For the later evaluation of the fluence SEDs, the relevant quantities are transformed into the observer frame (see Appendix G). Ordered lists of the initial and transition values G_i, , , , , , , and , as well as a list keeping track of the elements of the sets S₁, S₂ , and S₃ , and a list of the various solution branches are implemented. These lists are constantly updated during runtime. The algorithm constructs G incrementally in two loops. The first loop covers G from the time of the first injection t = 0 to the time of the second injection t = t₂ or – in a scenario with only a single injection – to the upper time boundary t_end using the solutions (57)–(59). The second loop computes G from the time of the second injection to the time t_end in case t₂ < t_end employing the solutions (26), (31)–(33), (38)–(41), (46)–(48), and (54)–(56). During each step, the current time t_cur. is incremented by a fixed value and G_cur. := G(t_cur.) is determined according to the proper solution branch, which is automatically selected via the abovementioned lists. Moreover, at each grid point, the analytical expressions for G are glued together continuously. In more detail, after initializing the values of , , , , as well as and at t = 0, the first loop starts its iteration, if , in the first SID solution branch of (57). At each step, it evaluates G at the current time grid point and checks whether G_cur. exceeds or equals , that is, whether G needs to be expressed by the second SID solution branch of (57). If necessary, G is altered accordingly and connected to the previous solution branch continuously. (In case or , the first loop makes use of the solutions (58) or (59), respectively.) Also, if G_cur. becomes larger than or equal to or , the list for S₁, S₂ , and S₃ is updated. The loop ends if t_cur. exceeds or equals either t₂ or t_end. In the latter case, the numerical constructionof G is completed. The second loop constructs G in a similar way as the first loop, but now more cases, which arise from the more elaborate structure of the analytical multiple injection solution, have to be taken into account. Once the second loop ends, the values of G are known at every point of the grid. This allows us to evaluate the synchrotron and SSC intensities at each grid point by simply substituting these values into the corresponding analytical formulas (64) and (68). The associated total fluences are approximated by sums of the areas of rectangles, each of which is defined by the intensity at the left grid point and the size of the time step. Alternatively, one could implement the approximate analytical formulas derived in Sect. 4. The total synchrotron and SSC fluence SEDs are then computed as the product of the respective energy and fluence. Since the number of grid points can be increased arbitrarily, the precision of these computations is limited only by the machine accuracy and the available CPU time.

References

All Figures

Fig. 1 Total synchrotron (left curves) and SSC (right curves) fluence SEDs for generic three-injection scenarios with magnetic field strengths, Doppler boost factors, and reciprocal initial electron energies b ∈ {0.01, 0.1, 1}, , and for panel a, b= 1, , and for panel b, as well as b = 1, , x1 ∈ {0.6 × 10−4, 10−4, 10−3}, and for panel c. The injection times are always given by with the characteristic radius of the plasmoid , the injection strengths by , and the cooling rate prefactors D0 = 1.3 × 10−9 b2 s−1, A0 = 1.2 × 10−18 b2 cm3 s−1 depend on the specific value of the magnetic field strength. The number of grid points used in the plotting of each curve amounts to 4 × 105.

In the text

	Fig. F.1 Absolute values of the relative deviations of the approximate functions CSn for all n∈{1, 2, 3} with respect to the exact CS function in percent.
In the text