Issue 
A&A
Volume 616, August 2018



Article Number  A172  
Number of page(s)  29  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201730622  
Published online  03 September 2018 
Flaring of blazars from an analytical, timedependent model for combined synchrotron and synchrotron selfCompton radiative losses of multiple ultrarelativistic electron populations
^{1}
Universität Regensburg,
Fakultät für Mathematik,
93040 Regensburg,
Germany
email: christian.roeken@mathematik.uniregensburg.de
^{2}
RuhrUniversität Bochum,
Institut für Theoretische Physik,
Lehrstuhl IV: Weltraum und Astrophysik,
44780 Bochum,
Germany
^{3}
GeorgAugustUniversität Göttingen,
Institut für Mathematische Stochastik,
37077 Göttingen,
Germany
Received:
14
February
2017
Accepted:
22
September
2017
We present a fully analytical, timedependent leptonic onezone model that describes a simplified radiation process of multiple interacting ultrarelativistic electron populations and accounts for the flaring of GeV blazars. In this model, several monoenergetic, ultrarelativistic electron populations are successively and instantaneously injected into the emission region, that is, a magnetized plasmoid propagating along the blazar jet, and subjected to linear, timeindependent synchrotron radiative losses, which are caused by a constant magnetic field, and nonlinear, timedependent synchrotron selfCompton radiative losses in the Thomson limit. Considering a general (timedependent) multipleinjection scenario is, from a physical point of view, more realistic than the usual (timeindependent) singleinjection scenario invoked in common blazar models, as blazar jets may extend over tens of kiloparsecs and thus most likely pick up several particle populations from intermediate clouds. We analytically compute the electron number density by solving a kinetic equation using Laplace transformations and the method of matched asymptotic expansions. Moreover, we explicitly calculate the optically thin synchrotron intensity, the synchrotron selfCompton intensity in the Thomson limit, as well as the associated total fluences. In order to mimic injections of finite duration times and radiative transport, we model flares by sequences of these instantaneous injections, suitably distributed over the entire emission region. Finally, we present a parameter study for the total synchrotron and synchrotron selfCompton fluence spectral energy distributions for a generic threeinjection scenario, varying the magnetic field strength, the Doppler factor, and the initial electron energy of the first injection in realistic parameter domains, demonstrating that our model can reproduce the typical broadband behavior seen in observational data.
Key words: galaxies: active / galaxies: jets / BL Lacertae objects: general / radiation mechanisms: nonthermal / methods: analytical
© 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 lowenergy spectral component, ranging from radio to optical or Xray energies, is usually attributed to synchrotron radiation of relativistic electrons that are subjected to ambient magnetic fields. The origin of the highenergy spectral component, covering the Xray to γray regime, is still under debate. It can, for instance, be modeled by inverse Compton radiation coming from lowenergy photon fields that interact with the relativistic electrons (Böttcher 2007; Dermer et al. 1997). This process can be described either by a synchrotron selfCompton (SSC) model (see, for example, Schlickeiser & Röken 2008; Röken & Schlickeiser 2009b and references therein), where the electrons scatter their selfgenerated synchrotron photons, or by socalled 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 broadline 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 highenergy component can also be modeled via protonsynchrotron 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 2155304 (Aharonian et al. 2007, 2009) and Mrk 501 (Albert et al. 2007) in the TeV range, or Mrk 421 (Cui 2004) in the Xray domain. So far, only multizone 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 shorttime 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 onezone models, however, extreme shorttime 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 shorttime 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 onezone and multizone 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 pickup 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 selfconsistent description of the particle’s radiative cooling, the radiative transport, et cetera. Therefore, we propose a simple, but fully analytical, timedependent leptonic onezone 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 randomlyoriented, but constant, largescale 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 monoenergetic, spatially isotropically distributed electron populations, which are subjected to linear, timeindependent synchrotron radiative losses via interactions with the ambient magnetic fields and to nonlinear, timedependent 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 smallscale, 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 timedependent, relativistic kinetic equation of the volumeaveraged 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 volumeaveraged differential electron number density n = n(γ, t) (where t is the time, γ:= E_{e}∕(m_{e} c^{2}) the normalized electron energy, and [n] = cm^{−3}) of m ultrarelativistic, monoenergetic, instantaneously injected and spatially isotropically distributed electron populations in the rest frame of a nonthermal radiation source with dominant magnetic field selfgeneration and radiative loss rate L = L(γ, t) (with [L] =s^{−1}) reads (Kardashev 1962) (1)
where δ(⋅) is the Dirac distribution,q_{i} (for which [q_{i}] = cm^{−3}) the ith injection strength, γ_{i}:= E_{e,i}∕(m_{e} c^{2}) ≫ 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, timeindependent synchrotron cooling process (with a constant magnetic field B) and a nonlinear, timedependent SSC cooling process in the Thomson limit (2)
are considered. The respective cooling rate prefactors yield D_{0} = 1.3 × 10^{−9} b^{2} s^{−1} and A_{0} = 1.2 × 10^{−18} b^{2} cm^{3} s^{−1}, 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 timedependent 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^{4} 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_{e}c^{2}), 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):= γ^{2} n(γ, t) and the variable x:= 1∕γ, we can rewrite Eq. (1) in the form (3)
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 singleinjection 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_{0}. 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 piecewisedefined 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_{1} = 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 timedependent sets S_{1}, S_{2}, and S_{3} that contain indices i: 1 ≤ i ≤ j belonging to either injections in the nearinjection domain (NID) 0 ≤ G − G_{i} ≪ x_{i}, the intermediateinjection domain (IID) G − G_{i} ≈ x_{i}, or the farinjection 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 NIDIID and IIDFID 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 nexttoleadingorder terms. This results in
for all u ∈ S_{1}, v ∈ S_{2}, and w ∈ S_{3}. 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 nexttoleadingorder 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)
which depend on the actual numbers of elements of S_{1}, S_{2}, and S_{3}. 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_{1}, S_{2}, and S_{3} 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_{1} and S_{2} as well as between S_{2} and S_{3}) 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_{1} and S_{2} or of S_{2} and S_{3}. 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 nexttoleadingorder contributions only if S_{1}, or S_{2}, or S_{1} and S_{2}, or S_{3} has to be taken into account, whereas only the leadingorder terms are considered if S_{1} and S_{3}, or S_{2} and S_{3}, or S_{1} and S_{2} and S_{3} have to be employed.
2.2.1 Nearinjection domain solution
In the NID , at least the jth injection is an element of S_{1}. Further, one may – but does not necessarily need to – have injections in S_{2} and S_{3}. 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
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 secondorder 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 nonsolvable ODEs. Hence, we have to employ the simpler ODEs
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 thirdorder 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_{3} had to be fixed via the initial condition G(t = t_{j}) = G_{j} resulting in
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)
The integration constants c_{2} and c_{4},..., c_{7} 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_{1}, S_{2}, or S_{3}. 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 NIDIID or IIDFID 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 2jtuple, 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_{1}, S_{2}, and S_{3} have to be updated accordingly, meaning that the corresponding injection is moved from either S_{1} to S_{2} or from S_{2} to S_{3} and the constants involved change. For more details on the updating see Sect. 2.3 and Appendix D.
2.2.2 Intermediateinjection domain solution
In the IID , where at least the jth injection is in S_{2}, 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)
is the IID transition time associated with the transition value
The derivation of the integration constants d_{1},..., d_{5} can be found in Appendix D.
2.2.3 Farinjection domain solution
In the FID , the jth injection is, per definition, an element of S_{3}. 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)
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 firstorder geometric series approximation in the integrand of (49) (50)
Note that in the first step, we included the multiplicative identity D_{0} ∕D_{0} 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 secondorder 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 thirdorder 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 piecewisedefined solution (54)
In case only one asymptotic end exists, that is, for or , we get (55)
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_{1},..., e_{9} 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 singleinjection domain (SID), the solution branch G(t  0 ≤ t < t_{2}) 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
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)
The integration constants f_{1}, f_{2}, f_{4}, and f_{5} are fixed by the initial condition G(t = t_{1} = 0) = G_{1} = 0
whereas the integration constant f_{3} is determined by the transition condition
In the following, we use the more exact SID approximation (57)–(59). The doubleinjection domain (DID) solution branch G(t  t_{2} ≤ t < t_{3}) is given for by Sol. (26) (with j = 2 as for the other DID solutions below) if and by
if S_{2} = ∅, S_{3}≠∅ (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
if S_{1} = ∅, S_{3}≠∅. Finally, for , we apply
if either S_{1}≠∅, S_{2} = ∅ or S_{1} = ∅, S_{2}≠∅ and
if S_{1} = ∅ = S_{2}. The initial value G_{2} 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_{2}, we have to perform NIDIID updates at and/or if as well as IIDFID 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_{1}, S_{2}, and S_{3} thrice. More precisely, at the time of the second injection t = t_{2}, both the first and the second injection are contained in S_{1}, while S_{2} and S_{3} are empty. During the temporal progression toward the upper bound t_{3}, the first injection switches from S_{1} to S_{2} at , whereas the second injection continues to be in S_{1}. At , also the second injection switches over to S_{2}, leaving S_{1} empty. Last, at , the first injection switches from S_{2} to S_{3}. 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 selfCompton 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 highenergy SSC component.
3.1 Synchrotron intensity
The optically thin synchrotron intensity I_{syn.}(ϵ, t) (with [I_{syn.}] = eV s^{−1} cm^{−2} sr^{−1} and similar for the SSC intensity) for an isotropically distributed electron number density is given by (61)
is the pitchangleaveraged spectral synchrotron power of a single electron in a magnetic field of strength b, ϵ:=E_{syn.}∕(m_{e} c^{2}) is the normalized synchrotron photon energy, P_{0}:= 8.5 × 10^{23} eV s^{−1}, and ϵ_{0}:= 2.3 × 10^{−14} b (Blumenthal & Gould 1970). The CS function is discussed in detail in Appendix F. Here, we employ the approximation (63)
where a_{0}:= 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 energyintegrated 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 selfCompton 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^{2}) 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^{−25} cm^{2} 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 SSCscattered photons is ϵ_{s} ≈ 4 γ^{2} ϵ (Felten & Morrison 1966), corresponding to headon 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
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 energyintegrated SSC intensity, we find
where and . The integral is given by (65), however, with adapted integration limits.
4 Synchrotron and synchrotron selfCompton 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^{−2} sr^{−1}) is given by (69)
The Jacobian determinant yields (70)
where i: 2 ≤ i ≤ m. Substituting(70) into (69), we obtain the expression
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
Similar to the evaluation of the energyintegrated 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 selfCompton 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 NIDIID and IIDFID transition times, which are explained in more detail in Appendix C, as follows. Applying the meanvaluetheorem 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 trapezoidapproximation 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 energyintegrated 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 energyintegrated 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 nonnegative, 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 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 lefthand side and the righthand 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_{1}, 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^{5}, 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 threeinjection 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, , x_{1} ∈ {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 D_{0} = 1.3 × 10^{−9} b^{2} s^{−1}, A_{0} = 1.2 × 10^{−18} b^{2} cm^{3} 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 × 10^{5}. 
6 Summary and outlook
We introduced a fully analytical, timedependent leptonic onezone 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 volumeaveraged differential electron number density for several successively and instantaneously injected, monoenergetic, spatially isotropically distributed, interacting electron populations, which are subjected to linear, timeindependent synchrotron radiative losses and nonlinear, timedependent 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 threeinjection 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 powerlaw 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 lefthand side via integration by parts and employing the Laplace transform of R (A.1)
As the normalized initial electron energies are finite and bounded from above by γ_{i} < 1.9 × 10^{4} 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 lefthand 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 nonnegative. With the Laplace transform of K
and integration by parts as before, Eq. (A.4) reads
Since the Heaviside functions are not welldefined 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 righthand 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_{1} δ(x − x_{1}), where . Substituting this density into (A.1) by employing the relation R(x, G) = n(x, G)∕x^{2}, we get
Choosing , we can use the last term on the righthand 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 NIDIID transition point and the IID and FID errors at the IIDFID transition point , we obtain
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 seventhorder 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 rootfinding algorithm to compute numerical approximations.
One may obtain a more accurate overall approximation of Eq. (12) using the more general IID approximation
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_{2}. 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 NIDIID and IIDFID 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 NIDIID 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^{−1} 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^{−1}, 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^{−1}(G_{j})
and the area A_{t}(j) of the rightangled triangle with one cathetus given by the distance between the endpoints and the other one by the height
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 IIDFID transition time
can be computed similarly, resulting, on the one hand, in the expression
for the meanvaluetheorem 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_{2} and c_{4}, ..., c_{7} of the NID (see Sol. (26) and Sols. (31)–(33)), d_{1} , ..., d_{5} of the IID (see Sols. (38)–(41)), as well as e_{1}, ..., e_{4} and e_{6}, ..., e_{9} of the FID (see Sols. (46)–(48) and Sols. (54)–(56)) are determined. The NID integration constants c_{2}, c_{4} , c_{6} , and c_{7} are fixed via the initial condition G(t = t_{j}) = G_{j}, whereas c_{5} 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_{1}, d_{2} , d_{4} , and d_{5} are determined via the respective initial condition and d_{3} 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_{1}, e_{3} , e_{4} , e_{6} , e_{8} , and e_{9} are specified by the initial condition and the integration constants e_{2} and e_{7} 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_{1} switch over to S_{2} and/or elements of S_{2} switch over to S_{3}. 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_{4} for the case S_{2}≠∅≠S_{3} 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_{4} depends just on S_{3}, we have to consider only IIDFID updates. Hence, the conditions – and the updated constants – for all possible IIDFID transitions of the ith injection at are
Since S_{3} ∪{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
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
and on the other hand, we have
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
If, however, the (i − 1)th injection is in the IID, they result for in
Finally, when the (i − 1)th injection is already in the FID, we find in case S_{1}≠∅, S_{2} = ∅ or
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_{0} := 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 CS_{n} 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 Xrays, that is, with a lower energy limit on the order neV and an upper limit on the order keV, and z = 2 ϵ∕(3 ϵ_{0} γ^{2}) with ϵ_{0} ~ 10^{−14} b, initial electron energies γ_{i} ~ 10^{4} b^{−1∕3}, and nondimensional magnetic field strength b ~ 10^{−3} 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_{3}(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
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∕ε^{3} is Lorentzinvariant, 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^{1} and the specific incorporation of the analytical formulas. The algorithm makes heavy use of the decimal package^{2} , which is designed for high floatingpoint precision calculations. This is necessary because the range of possible values between the synchrotron and SSC cooling rate prefactors D_{0} and A_{0} , 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_{0} and A_{0}, 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 multipleinjection scenario or given by a sufficiently large value for a singleinjection 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_{1}, S_{2} , and S_{3} , 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_{2} 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_{2} < 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_{1}, S_{2} , and S_{3} is updated. The loop ends if t_{cur.} exceeds or equals either t_{2} 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
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (National Bureau of Standards) [Google Scholar]
 Aharonian F. A., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, ApJ, 664, L71 [Google Scholar]
 Aharonian F. A., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 502, 749 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Biteau, J., & Giebels, B. 2012, A&A, 548, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Blazejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Böttcher, M. 2007, Ap&SS, 309, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Böttcher, M. 2010, ArXiv eprints [1006.5048] [Google Scholar]
 Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2011, Proc. Ann. meet. French Soc. Astron. Astrophys., 555 [Google Scholar]
 Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, MNRAS, 448, 910 [NASA ADS] [CrossRef] [Google Scholar]
 Crusius, A., & Schlickeiser, R. 1988, A&A, 196, 327 [NASA ADS] [Google Scholar]
 Cui, W. 2004, ApJ, 605, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Dermer, C. D., Sturner, S. J., & Schlickeiser, R. 1997, ApJS, 109, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Eichmann, B., Schlickeiser, R., & Rhode, W. 2012, ApJ, 744, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Felten, J. E., & Morrison, P. 1966, ApJ, 146, 686 [NASA ADS] [CrossRef] [Google Scholar]
 Ghisellini, G., Tavecchio, F., Bodo, G., & Celotti, A. 2009, MNRAS, 393, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009, MNRAS, 395, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Graff, P. B., Georganopoulos, M., Perlman, E. S., & Kazanas, D. 2008, ApJ, 689, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, F. C. 1968, Phys. Rev., 167, 1159 [NASA ADS] [CrossRef] [Google Scholar]
 Kardashev, N. S. 1962, Sov. Astron., 6, 317 [NASA ADS] [Google Scholar]
 Li, H., & Kusunose, M. 2000, ApJ, 536, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Marscher, A. P., 2014, ApJ, 780, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Pohl, M., & Schlickeiser, R. 2000, A&A, 354, 395 [NASA ADS] [Google Scholar]
 Reynolds, S. P. 1982, ApJ, 256, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Röken, C., & Schlickeiser, R. 2009a, ApJ, 700, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Röken, C., & Schlickeiser, R. 2009b, A&A, 503, 309 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R. 2009, MNRAS, 398, 1483 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R., & Lerche, I. 2007, A&A, 476, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R., & Röken, C. 2008, A&A, 477, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlickeiser, R., Böttcher, M., & Menzler, U. 2010, A&A, 519, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sikora, M., Begerlman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Sikora, M., Blazejowski, M., Madejski, G. M., & Moderski, R. 2001, Proc. fourth INTEGRAL workshop, Exploring the gammaray Universe, 259 [Google Scholar]
 Sokolov, A., Marscher, A. P., & McHardy, I. M. 2004, ApJ, 613, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Weidinger, M., & Spanier, F. 2015, A&A, 573, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yan, D., & Zhang, L. 2015, MNRAS, 447, 2810 [NASA ADS] [CrossRef] [Google Scholar]
 Zacharias, M., & Schlickeiser, R. 2013, ApJ, 777, 109 [NASA ADS] [CrossRef] [Google Scholar]
The plotting algorithm is available on request via email to christian.roeken@mathematik.uniregensburg.de
For more detailed information on Python's decimal package, we refer to https://docs.python.org/2/library/decimal.html
All Figures
Fig. 1 Total synchrotron (left curves) and SSC (right curves) fluence SEDs for generic threeinjection 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, , x_{1} ∈ {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 D_{0} = 1.3 × 10^{−9} b^{2} s^{−1}, A_{0} = 1.2 × 10^{−18} b^{2} cm^{3} 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 × 10^{5}. 

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

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.