Flaring of Blazars from an Analytical, Time-dependent Model for Combined Synchrotron and Synchrotron Self-Compton Radiative Losses of Multiple Ultrarelativistic Electron Populations
Christian Röken, Florian Schuppan, Katharina Proksch, Sebastian Schöneberg
FFlaring of Blazars from an Analytical, Time-dependent Model forCombined Synchrotron and Synchrotron Self-Compton RadiativeLosses of Multiple Ultrarelativistic Electron Populations
Christian R¨oken ∗ Universit¨at Regensburg, Fakult¨at f¨ur Mathematik, 93040 Regensburg, Germany
Florian Schuppan
Ruhr-Universit¨at Bochum, Institut f¨ur Theoretische Physik,Lehrstuhl IV: Weltraum- und Astrophysik, 44780 Bochum, Germany
Katharina Proksch
Georg-August-Universit¨at G¨ottingen, Institut f¨ur Mathematische Stochastik, 37077 G¨ottingen, Germany
Sebastian Sch¨oneberg
Ruhr-Universit¨at Bochum, Institut f¨ur Theoretische Physik,Lehrstuhl IV: Weltraum- und Astrophysik, 44780 Bochum, Germany (Dated: September 2017)
ABSTRACT.
We present a fully analytical, time-dependent leptonic one-zone model that de-scribes a simplified radiation process of multiple interacting ultrarelativistic electron populations,accounting for the flaring of GeV blazars. In this model, several mono-energetic, ultrarelativisticelectron populations are successively and instantaneously injected into the emission region, i.e., amagnetized plasmoid propagating along the blazar jet, and subjected to linear, time-independentsynchrotron radiative losses, which are caused by a constant magnetic field, and nonlinear, time-dependent synchrotron self-Compton radiative losses in the Thomson limit. Considering a general(time-dependent) multiple-injection scenario is, from a physical point of view, more realistic thanthe usual (time-independent) single-injection scenario invoked in common blazar models, as blazarjets may extend over tens of kiloparsecs and, thus, most likely pick up several particle populationsfrom intermediate clouds. We analytically compute the electron number density by solving a kineticequation using Laplace transformations and the method of matched asymptotic expansions. More-over, we explicitly calculate the optically thin synchrotron intensity, the synchrotron self-Comptonintensity in the Thomson limit, as well as the associated total fluences. In order to mimic injectionsof finite duration times and radiative transport, we model flares by sequences of these instantaneousinjections, suitably distributed over the entire emission region. Finally, we present a parameterstudy for the total synchrotron and synchrotron self-Compton fluence spectral energy distributionsfor a generic three-injection scenario, varying the magnetic field strength, the Doppler factor, andthe initial electron energy of the first injection in realistic parameter domains, demonstrating thatour model can reproduce the typical broad-band behavior seen in observational data.
Contents
I. Introduction II. The Relativistic Kinetic Equation (cid:8) t ∈ R ≥ | t j ≤ t < t j +1 with j ∈ { , ..., m } (cid:9)
61. Near-injection Domain Solution 72. Intermediate-injection Domain Solution 93. Far-injection Domain Solution 10C. Solution of the Relativistic Kinetic Equation for { t ∈ R ≥ } ∗ e-mail: [email protected] a r X i v : . [ a s t r o - ph . H E ] S e p III. Synchrotron and Synchrotron Self-Compton Intensities
IV. Synchrotron and Synchrotron Self-Compton Fluences
V. Lightcurves and Fluence Spectral Energy Distributions VI. Summary and Outlook Acknowledgments A. Laplace Transformation Method B. Approximation Errors and Optimal Domains C. NID-IID and IID-FID Transition Times D. Integration Constants – Initial and Transition Conditions and Updating E. Components of G ( t | ≤ t < ∞ ) and Initial Values G i F. The CS Function G. Relativistic Beaming H. Plotting Algorithm for the Fluence Spectral Energy Distributions References I. INTRODUCTION
Blazars are among the most energetic phenomena in nature, representing the most extreme type in the class ofactive galactic nuclei [39]. They feature relativistic jets that extend over tens of kiloparsecs and are directedtoward 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 detectedalong the first few parsecs of the jets. The main components of blazar jets are magnetized plasmoids, whichare assumed to arise in the Blandford-Znajek and the Blandford-Payne process [6, 7], constituting the majorradiation zones. These plasmoids pick up – and interact with – particles of interstellar and intergalactic cloudsalong their trajectories [28], giving rise to the emission of a series of strong flares.A blazar spectral energy distribution (SED) consists of two broad non-thermal radiation components in dif-ferent domains. The low-energy spectral component, ranging from radio to optical or X-ray energies, is usuallyattributed to synchrotron radiation of relativistic electrons subjected to ambient magnetic fields. The originof the high-energy spectral component, covering the X-ray to γ -ray regime, is still under debate. It can, forinstance, be modeled by inverse Compton radiation coming from low-energy photon fields that interact with therelativistic electrons [10, 18]. This process can be described either by a synchrotron self-Compton (SSC) model(see, e.g., [30, 32] and references therein), where the electrons scatter their self-generated synchrotron photons, orby so-called external Compton models like [8, 17, 36], 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 IRradiation from diffuse, hot dust, see, e.g., [37]). Aside from these leptonic scenarios, the high-energy componentcan also be modeled via proton-synchrotron radiation or the emission of γ -rays arising from the decay of neutralpions formed in interactions of protons with ambient matter (see, e.g., [12, 14, 40] and references therein). Mixedmodels including both leptonic as well as hadronic processes are also considered in the literature (e.g., [13, 41]).A major task is to properly understand and to account for the distinct variability patterns of the non-thermalblazar emission at all frequencies with different time scales ranging from years down to a few minutes, where theshortest variability time scales are usually observed for the highest energies of the spectral components, as in PKS2155-304 [2, 3] and Mrk 501 [4] in the TeV range, or Mrk 421 [16] in the X-ray domain. So far, only multi-zonemodels, which feature an internal structure of the emission region with various radiation zones caused by collisionsof moving and stationary shock waves, have been proposed to explain the extreme short-time variability of blazars(see, e.g., [5, 21–23, 27, 38]). 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 R ) and the light crossing (escape) time in this region are naturally of the order O (2 R / ( c D )), c being the speedof light in vacuum and D the bulk Doppler factor [19]. This leads to a minimum time scale for the observedflare duration, which may exceed the short-time variability scale in the minute range by several magnitudes. Inparticular, using the typical parameters R = 10 cm and D = 10, we find 2 R / ( c D ) ≈ . II. THE RELATIVISTIC KINETIC EQUATION
The kinetic equation for the volume-averaged differential electron number density n = n ( γ, t ) (where t is thetime, γ := 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 non-thermal radiation source with dominant magnetic field self-generation and radiative loss rate L = L ( γ, t ) (with[ L ] = s − ) reads [25] ∂n∂t − ∂∂γ (cid:0) L n (cid:1) = m (cid:88) i =1 q i δ ( γ − γ i ) δ ( t − t i ) , (1)where δ ( · ) is the Dirac distribution, q i (for which [ q i ] = cm − ) the i th injection strength, γ i := E e ,i / ( m e c ) (cid:29) i th normalized initial electron energy, and t i the i th injection time for i : 1 ≤ i ≤ m . In this work, radiativelosses 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 L = L syn. + L SSC = γ (cid:18) D + A (cid:90) ∞ γ n d γ (cid:19) (2)are considered. The respective cooling rate prefactors yield D = 1 . × − b s − and A = 1 . × − b cm s − ,which depend on the nondimensional magnetic field strength b := (cid:107) B (cid:107) / Gauss = const. [9, 33]. In the presentcontext, the term linear refers to the fact that the synchrotron loss term does not depend on the electron numberdensity n , thus, resulting in a linear contribution to the partial differential equation (PDE) (1), whereas the SSCloss term depends on an energy integral containing n , yielding a nonlinear contribution. The synchrotron lossterm can, in principle, be modeled as nonlinear and time-dependent, too. This can be achieved by making anequipartition assumption between the magnetic and particle energy densities [34]. We point out that, exceptfor TeV blazars, where Klein-Nishina effects drastically reduce the SSC cooling strength above a certain energythreshold, the dominant contribution of the SSC energy loss rate always originates in the Thomson regime [33],justifying our initial restriction. But as a consequence, this limits the applicability of our model to at most GeVblazars and bounds the normalized initial electron energies from above by γ i < . × b − / . We also note thatthe 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 ), i.e., γ = E tot. − m e c m e c = γ tot. − . 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( γ i , t i ) i ∈{ ,...,m } . They are to be understood in the distributional sense, that is, each is a linear functional on thespace of smooth test functions ϕ on R with compact support C ∞ ( R ) := (cid:8) ϕ | ϕ ∈ C ∞ ( R ) , supp ϕ compact (cid:9) . More precisely, for k ∈ R , one can rigorously define them as the mapping δ : C ∞ ( R ) → R , ϕ (cid:55)→ ϕ (0)with the integral of the Dirac distribution against a test function given by (cid:90) ∞−∞ ϕ ( k ) δ ( k ) d k = ϕ (0) for all ϕ ∈ C ∞ ( R ) . A. 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 ∂R∂t + J ∂R∂x = m (cid:88) i =1 q i δ ( x − x i ) δ ( t − t i ) , (3)where J = J ( t ) := D + A (cid:90) ∞ R ( x, t ) x d x . (4)Defining the strictly increasing, continuous function G = G ( t ) via0 < d G d t := J , (5)Eq. (3) becomes ∂R∂G + ∂R∂x = m (cid:88) i =1 q i δ ( x − x i ) δ ( G − G i ) (6)with G i := G ( t i ). Albeit the nonlinear kinetic equation (1) is now transformed into a linear PDE, its solutioncan obviously serve as a Green’s function only for the single-injection scenario with m = 1. Consequently, inorder to solve the generalized kinetic equation ∂n∂t − ∂∂γ (cid:0) L n (cid:1) = Q 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 R ( x, G ) = m (cid:88) i =1 q i H ( G − G i ) δ ( x − x i − G + G i ) , (7)where H ( · ) is the Heaviside step function H ( k − k ) := (cid:40) k < k k > k with k, k ∈ R , which has a jump discontinuity at k = k . A detailed derivation of this solution can be found inAppendix A. In order to determine the function G , we substitute solution (7) into (4) and, by using (5), obtainthe ordinary differential equation (ODE)d G d t = D + A m (cid:88) i =1 q i H ( G − G i ) (cid:0) G − G i + x i (cid:1) . (8)This equation can be regarded as a compact notation for the set of m piecewise-defined ODEs d G d t = D + A q (cid:0) G − G + x (cid:1) for G ≤ G < G d G d t = D + A (cid:32) q (cid:0) G − G + x (cid:1) + q (cid:0) G − G + x (cid:1) (cid:33) for G ≤ G < G ... ... ... ...d G d t = D + A (cid:32) q (cid:0) G − G + x (cid:1) + q (cid:0) G − G + x (cid:1) + . . . + q m (cid:0) G − G m + x m (cid:1) (cid:33) for G m ≤ G < ∞ . (9)In Section II B, we present an approximate analytical solution for the general case of j injections, i.e., for theinterval G j ≤ G < G j +1 , where j ∈ { , ..., m } and G = 0, G m +1 = ∞ , employing the method of matchedasymptotic expansions. Having a separate analytical solution for each ODE of (9), in Section II C, these areconnected successively requiring continuity at the transition points. Finally, by substituting (7), the electronnumber density results in n ( γ, t ) = γ − R ( γ, t ) = γ − m (cid:88) i =1 q i H (cid:0) G ( t ) − G i (cid:1) δ (cid:18) γ − γ i − G ( t ) + G i (cid:19) = m (cid:88) i =1 q i H ( t − t i ) δ (cid:32) γ − γ i γ i (cid:0) G ( t ) − G i (cid:1) (cid:33) . (10) B. Solution of the Relativistic Kinetic Equation for (cid:8) t ∈ R ≥ | t j ≤ t < t j +1 with j ∈ { , ..., m } (cid:9) In the interval G j ≤ G < G j +1 , Eq. (8) givesd G d t = D + A j (cid:88) i =1 q i (cid:0) G − G i + x i (cid:1) . (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 (cid:28) x i , the intermediate-injection domain (IID) G − G i ≈ x i , or the far-injection domain (FID) G − G i (cid:29) x i , respectively. We may nowrewrite Eq. (11) in terms of these sets by continuing their domains of validity to G i ≤ G < min (cid:0) G j +1 , G (N → I)T ( i ) (cid:1) for the NID, G (N → I)T ( i ) ≤ G < min (cid:0) G j +1 , G (I → F)T ( i ) (cid:1) for the IID, and G (I → F)T ( i ) ≤ G < G j +1 for the FID withthe NID-IID and IID-FID transition values G (N → I)T ( i ) := G i + x i /ξ (N → I) i and G (I → F)T ( i ) := G i + ξ (I → F) i x i , where ξ (N → I) i , ξ (I → F) i > j + 1)th injection can occur in eachdomain. 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 byd G d t = D + A (cid:32) (cid:88) u ∈ S q u (cid:0) G − G u + x u (cid:1) + (cid:88) v ∈ S q v (cid:0) G − G v + x v (cid:1) + (cid:88) w ∈ S q w (cid:0) G − G w + x w (cid:1) (cid:33) . (12)Next, we express the NID and the FID summands by the Taylor series of the function (1 + z ) − evaluated atthe point z = 0 (generalized geometric series)1(1 + z ) = ∞ (cid:88) n =0 ( − n ( n + 1) z n for | z | < , and the IID summands by the Taylor series of the same function evaluated at the point z = 11(1 + z ) = ∞ (cid:88) n =0 ( − n n + 12 n +2 ( z − n for | z − | < , and suitably approximate these series by considering only the leading- and next-to-leading-order terms. Thisresults in 1( G − G u + x u ) = 1 x u ∞ (cid:88) n =0 ( n + 1) (cid:18) − G − G u x u (cid:19) n ≈ x u (cid:18) − G − G u ) x u (cid:19) (13)1( G − G v + x v ) = 1 x v ∞ (cid:88) n =0 n + 12 n +2 (cid:18) − G − G v x v (cid:19) n ≈ x v (cid:18) − G − G v x v (cid:19) (14)1( G − G w + x w ) = 1( G − G w ) ∞ (cid:88) n =0 ( n + 1) (cid:18) − x w G − G w (cid:19) n ≈ G − G w ) (cid:18) − x w G − G w (cid:19) (15)for all u ∈ S , v ∈ S , and w ∈ S . To find an approximate analytical solution to Eq. (12), we require thefunction 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 thedenominator. To this end, we employ both geometric and generalized geometric series and again consider onlythe leading- and next-to-leading-order terms1( G − G w + x w ) ≈ G ∞ (cid:88) k =0 ( k + 1) (cid:18) G w G (cid:19) k (cid:34) − x w G ∞ (cid:88) l =0 (cid:18) G w G (cid:19) l (cid:35) ≈ G (cid:18) G w − x w ) G (cid:19) . (16)The approximation errors and optimal values for the interval parameters ξ (N → I) i and ξ (I → F) i are given in AppendixB. Substituting (13), (14), and (16) into the ODE (12), we obtaind G d t ≈ C ( j ; S , S ) + C ( j ; S , S ) G + C ( j ; S ) G + C ( j ; S ) G (17)with the “constants” C ( j ; S , S ) := D + A (cid:32) (cid:88) u ∈ S q u x u (cid:20) G u x u (cid:21) + 12 (cid:88) v ∈ S q v x v (cid:20) G v x v (cid:21)(cid:33) , C ( j ; S , S ) := − A (cid:32) (cid:88) u ∈ S q u x u + 14 (cid:88) v ∈ S q v x v (cid:33) , C ( j ; S ) := A (cid:88) w ∈ S q w , and C ( j ; S ) := 2 A (cid:88) w ∈ S q w ( G w − x w ) , (18)which depend on the actual numbers of elements of S , S , and S . Note that the explicit dependences of theseconstants 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 anapproximate analytical solution of Eq. (17) by computing separate solutions for the NID, IID, and FID of the j th injection, which, in case G (N → I)T ( j ) < G (I → F)T ( j ) < G j +1 , are glued together continuously at the transitionpoints G (N → I)T ( j ) and G (I → F)T ( j ) (that define the transitions of the j th injection between S and S as well asbetween S and S ) corresponding to the transition times t (N → I)T ( j ) and t (I → F)T ( j ) specified in Appendix C. For G (N → I)T ( j ) < G j +1 ≤ G (I → F)T ( j ), we regard only the NID and IID solutions with the proper continuous gluingat G (N → I)T ( j ), and for G j +1 ≤ G (N → I)T ( j ) < G (I → F)T ( j ), we consider solely the NID solution. Moreover, we haveto update the above constants (18) each time an injection i : 1 ≤ i ≤ j crosses over from its NID to its IID at t = t (N → I)T ( i ) or from its IID to its FID at t = t (I → F)T ( i ), as these transitions cause changes in either the numbersof 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-ordercontributions only if S , or S , or S and S , or S has to be taken into account, whereas only the leading-orderterms are considered if S and S , or S and S , or S and S and S have to be employed.
1. Near-injection Domain Solution
In the NID G j ≤ G < min (cid:0) G j +1 , G (N → I)T ( j ) (cid:1) , at least the j th 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 fourcases d G d t = C ( j ; S ) + C ( j ; S ) G for S = ∅ = S (19)d G d t = C ( j ; S , S ) + C ( j ; S , S ) G for S (cid:54) = ∅ , S = ∅ (20)d G d t = C ( j ; S ) + C ( j ; S ) G for S = ∅ , S (cid:54) = ∅ (21)d G d t = C ( j ; S , S ) + C ( j ; S ) G for S (cid:54) = ∅ (cid:54) = S . (22)Since Eqs. (19) and (20) as well as Eqs. (21) and (22) are identical except for the present constants, only twodifferent kinds of ODEs have to be solved. Starting with Eqs. (19) and (20), we use separation of variables inorder to obtain (cid:90) d G C + C G = 1 C ln ( C + C G ) = t + c , c = const. ∈ R , which is equivalent to G ( t ) = 1 | C ( j ) | (cid:104) C ( j ) − exp (cid:0) −| C ( j ) | ( t + c ) (cid:1)(cid:105) for t j ≤ t < min (cid:0) t j +1 , t (N → I)T ( j ) (cid:1) . As this solution converges strictly against the value C ( j ) / | C ( j ) | , which can be smaller than the NID transitionvalue G (N → I)T ( j ), it is not suitable for covering this domain. Including a second-order term in the ODE underconsideration, that is, working with the equationd G d t = C + C G + C G , (23)where C = const. ∈ R , leads to a solution in form of a tangent function, which may have several poles in the NID,rendering it useless, too. Adding further contributions to (23) yields analytically non-solvable ODEs. Hence, wehave to employ the simpler ODEsd G d t = C ( j ; S ) for S = ∅ = S (24)d G d t = C ( j ; S , S ) for S (cid:54) = ∅ , S = ∅ (25)with the linear solution G ( t ) = C ( j ) t + c ( j ) for t j ≤ t < min (cid:0) t j +1 , t (N → I)T ( j ) (cid:1) (26)and c = const. ∈ R . For Eqs. (21) and (22), we also apply separation of variables and extend the integrandwith the multiplicative identity C / C and the neutral element C − C , leading to (cid:90) G d G C + C G = 1 C (cid:18) G − C (cid:90) d G C + C G (cid:19) = 1 C (cid:34) G − (cid:115) C C arctan (cid:32)(cid:115) C C G (cid:33)(cid:35) = t + c , (27)where c = const. ∈ R . One way of deriving an approximate analytical solution of this transcendental equation isto determine G asymptotically for small and large arguments of the arctan function (in case both asymptotic endsexist for G : G j ≤ G < min( G j +1 , G (N → I)T ( j ))) and extrapolate these solutions up to an intermediate transitionpoint, where they are connected requiring continuity, such that the entire domain of definition is covered. Forsmall arguments (cid:112) C / C G (cid:28)
1, we use the third-order approximation arctan ( z ) | z (cid:28) ≈ z − z /
3, as the linearterm in Eq. (27) and the linear contribution in the approximation of the arctan function cancel each other out.This leads to the asymptotic solution G ( t ) (cid:39) (3 C t + c ) / , c = const. ∈ R . (28)For large arguments (cid:112) C / C G (cid:29)
1, the arctan function can be well approximated by π/
2. We directly obtainan asymptotic solution of the form G ( t ) (cid:39) C t + c (cid:48) , c (cid:48) = const. ∈ R . (29)By means of the condition (cid:112) C / C G = 1, we derive the NID transition value G (N)T ( j ) = (cid:115) C ( j ) C ( j ) . (30)This value specifies the upper bound of the domain of validity of (28) and the lower bound of the domain ofvalidity of (29). The corresponding transition time can be directly computed by substituting (30) into (27), inwhich the constant c had to be fixed via the initial condition G ( t = t j ) = G j resulting in c = 1 C (cid:34) G j − (cid:115) C C arctan (cid:32)(cid:115) C C G j (cid:33)(cid:35) − t j . This yields t (N)T ( j ) = t j + 1 C ( j ) (cid:32)(cid:115) C ( j ) C ( j ) (cid:34) − π (cid:32)(cid:115) C ( j ) C ( j ) G j (cid:33)(cid:35) − G j (cid:33) . We point out that in general, one does not have the ordered sequence G j < G (N)T ( j ) < G (N → I)T ( j ) because G (N)T ( j )can in principle be larger than or equal to G (N → I)T ( j ) as well as smaller than or equal to G j . These cases arisewhen only one asymptotic end of the arctan function exists. Hence, for G j < G (N)T ( j ) < G (N → I)T ( j ), the solutionof Eq. (27) is approximately given by G ( t ) = (cid:40)(cid:0) C ( j ) t + c ( j ) (cid:1) / for t j ≤ t < min (cid:0) t j +1 , t (N)T ( j ) (cid:1) C ( j ) t + c ( j ) for min (cid:0) t j +1 , t (N)T ( j ) (cid:1) ≤ t < min (cid:0) t j +1 , t (N → I)T ( j ) (cid:1) , c , c = const. ∈ R , (31)for G j < G (N → I)T ( j ) ≤ G (N)T ( j ) by G ( t ) = (cid:0) C ( j ) t + c ( j ) (cid:1) / for t j ≤ t < min (cid:0) t j +1 , t (N → I)T ( j ) (cid:1) , c = const. ∈ R , (32)and for G (N)T ( j ) ≤ G j by G ( t ) = C ( j ) t + c ( j ) for t j ≤ t < min (cid:0) t j +1 , t (N → I)T ( j ) (cid:1) , c = const. ∈ R . (33)The integration constants c and c , ..., c are determined via the proper initial and transition conditions inAppendix D. In order to compute the constants C ( j ; S ), C ( j ; S , S ), C ( j ; S ), C ( j ; S , S ), and C ( j ; S ),we have to continually check during the evolution of G ∈ (cid:2) G j , min (cid:0) G j +1 , G (N → I)T ( j ) (cid:1)(cid:1) whether an injection i : 1 ≤ i ≤ j belongs to S , S , or S . Therefore, we specify two different kinds of updates. The first kind occurswhen 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 j th injection verifying G j − G i < x i /ξ (N → I) i for the i th injection being in S x i /ξ (N → I) i ≤ G j − G i < ξ (I → F) i x i for the i th injection being in S ξ (I → F) i x i ≤ G j − G i for the i th injection being in S . Next, since all transition values G (N → I)T ( i ) and G (I → F)T ( i ) are known, they can be arranged as an ordered 2 j -tuple representing an increasing sequence. Assuming that a certain number of these values is contained in thetime interval under consideration, whenever G reaches one of them, all sets S , S , and S have to be updatedaccordingly, i.e., the corresponding injection is moved from either S to S or from S to S and the constantsinvolved change. For more details on the updating see Section II C and Appendix D.
2. Intermediate-injection Domain Solution
In the IID G (N → I)T ( j ) ≤ G < min (cid:0) G j +1 , G (I → F)T ( j ) (cid:1) , where at least the j th injection is in S , the previous j − G d t = C ( j ; S ) for S = ∅ = S (34)d G d t = C ( j ; S , S ) for S (cid:54) = ∅ , S = ∅ (35)d G d t = C ( j ; S ) + C ( j ; S ) G for S = ∅ , S (cid:54) = ∅ (36)d G d t = C ( j ; S , S ) + C ( j ; S ) G for S (cid:54) = ∅ (cid:54) = S . (37)As these coincide structurally with the ODEs (24) and (25) as well as (21) and (22) of the NID, we can directlywrite down their solutions. For Eqs. (34) and (35), we get G ( t ) = C ( j ) t + d ( j ) for t (N → I)T ( j ) ≤ t < min (cid:0) t j +1 , t (I → F)T ( j ) (cid:1) , d = const. ∈ R , (38)0whereas for Eqs. (36) and (37), we obtain in case G (N → I)T ( j ) < G (I)T ( j ) < G (I → F)T ( j ) G ( t ) = (cid:40)(cid:0) C ( j ) t + d ( j ) (cid:1) / for t (N → I)T ( j ) ≤ t < min (cid:0) t j +1 , t (I)T ( j ) (cid:1) C ( j ) t + d ( j ) for min (cid:0) t j +1 , t (I)T ( j ) (cid:1) ≤ t < min (cid:0) t j +1 , t (I → F)T ( j ) (cid:1) , d , d = const. ∈ R , (39)if G (N → I)T ( j ) < G (I → F)T ( j ) ≤ G (I)T ( j ) G ( t ) = (cid:0) C ( j ) t + d ( j ) (cid:1) / for t (N → I)T ( j ) ≤ t < min (cid:0) t j +1 , t (I → F)T ( j ) (cid:1) , d = const. ∈ R , (40)and if G (I)T ( j ) ≤ G (N → I)T ( j ) G ( t ) = C ( j ) t + d ( j ) for t (N → I)T ( j ) ≤ t < min (cid:0) t j +1 , t (I → F)T ( j ) (cid:1) , d = const. ∈ R , (41)with the transition time t (I)T ( j ) = t (N → I)T ( j ) + 1 C ( j ) (cid:32)(cid:115) C ( j ) C ( j ) (cid:34) − π (cid:32)(cid:115) C ( j ) C ( j ) G (N → I)T ( j ) (cid:33)(cid:35) − G (N → I)T ( j ) (cid:33) associated to the transition value G (I)T ( j ) = (cid:115) C ( j ) C ( j ) . The derivation of the integration constants d , ..., d can be found in Appendix D.
3. Far-injection Domain Solution
In the FID G (I → F)T ( j ) ≤ G < G j +1 , the j th injection is, per definition, an element of S . As before, the previous j − G d t = D + C ( j ; S ) G + C ( j ; S ) G for S = ∅ = S (42)d G d t = C ( j ; S ) + C ( j ; S ) G for S (cid:54) = ∅ , S = ∅ (43)d G d t = C ( j ; S ) + C ( j ; S ) G for S = ∅ , S (cid:54) = ∅ (44)d G d t = C ( j ; S , S ) + C ( j ; S ) G for S (cid:54) = ∅ (cid:54) = S . (45)Since Eqs. (43), (44), and (45) are also structurally identical to the ODEs (21) and (22), we can once againdirectly write down their solutions, yielding for G (I → F)T ( j ) < G (F)T , ( j ) < G j +1 G ( t ) = (cid:40)(cid:0) C ( j ) t + e ( j ) (cid:1) / for t (I → F)T ( j ) ≤ t < t (F)T , ( j ) C ( j ) t + e ( j ) for t (F)T , ( j ) ≤ t < t j +1 , e , e = const. ∈ R , (46)for G (I → F)T ( j ) < G j +1 ≤ G (F)T , ( j ) G ( t ) = (cid:0) C ( j ) t + e ( j ) (cid:1) / for t (I → F)T ( j ) ≤ t < t j +1 , e = const. ∈ R , (47)and for G (F)T , ( j ) ≤ G (I → F)T ( j ) < G j +1 G ( t ) = C ( j ) t + e ( j ) for t (I → F)T ( j ) ≤ t < t j +1 , e = const. ∈ R , (48)1where t (F)T , ( j ) = t (I → F)T ( j ) + 1 C ( j ) (cid:32)(cid:115) C ( j ) C ( j ) (cid:34) − π (cid:32)(cid:115) C ( j ) C ( j ) G (I → F)T ( j ) (cid:33)(cid:35) − G (I → F)T ( j ) (cid:33) is the transition time corresponding to the transition value G (F)T , ( j ) = (cid:115) C ( j ) C ( j ) . 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 (cid:90) G d GD G + C G + C = t + e , e = const. ∈ R . (49)This integral equation is shown to be approximately equivalent to a specific transcendental equation for whichwe subsequently derive asymptotic solutions for G that are continued up to an intermediate transition point (ifboth asymptotic ends exist, else we only consider a continuation of the solution of the present asymptotic endover the entire domain), where they are glued together continuously. In more detail, since C / ( C G ) (cid:28)
1, weemploy a first-order geometric series approximation in the integrand of (49) (cid:90) G d GD G + C G + C = GD − D (cid:90) C G + C D G + C G + C d G = GD − D (cid:90) d G D G C G + C = GD − D (cid:90) d G D G C (cid:18) C C G (cid:19) − = GD − D (cid:90) d G D G C ∞ (cid:88) n =0 (cid:18) − C C G (cid:19) n ≈ GD − D (cid:90) d G D C (cid:18) G − C C (cid:19) − D C C ≈ GD − D (cid:90) d G D C (cid:18) G − C C (cid:19) . (50)Note that in the first step, we included the multiplicative identity D /D and the neutral element ( C G + C ) − ( C G + C ) in the numerator of the integrand, giving the splitting into two terms. Then, by means of simplealgebraic manipulations, we expressed the integrand in a form suitable for the substitution of the geometricseries. For reasons of computational simplicity, we dropped the small second-order contribution − D C / (4 C )in the last step. The final integral in (50) is solved by an arctan function. Thus, introducing the parameter β = β ( j ) := D / C ( j ) / (cid:0) C / ( j ) (cid:1) , Eq. (49) results in GD − C / D / arctan (cid:18) β (cid:20) C G C − (cid:21)(cid:19) = t + e . (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-orderapproximation of the arctan function for small arguments β (2 C G/ C − (cid:28)
1, the associated asymptoticsolution becomes G ( t ) (cid:39) (3 C t + e ) / + C C , e = const. ∈ R . Further, because arctan ( x ) | x (cid:29) ≈ π/
2, the asymptotic solution for large arguments β (2 C G/ C − (cid:29) G ( t ) (cid:39) D t + e (cid:48) , e (cid:48) = const. ∈ R . G (F)T , ( j ) = (cid:115) C ( j ) D + C ( j )2 C ( j ) , (52)which is derived from the transition condition β (2 C G/ C −
1) = 1 and corresponds to the transition time t (F)T , ( j ) = t (I → F)T ( j ) + 1 D (cid:32)(cid:115) C ( j ) D (cid:34) − π (cid:18) β ( j ) (cid:20) C ( j ) G (I → F)T ( j ) C ( j ) − (cid:21)(cid:19)(cid:35) − G (I → F)T ( j ) + C ( j )2 C ( j ) (cid:33) , (53)we obtain for G (I → F)T ( j ) < G (F)T , ( j ) < G j +1 the piecewise-defined solution G ( t ) = (cid:0) C ( j ) t + e ( j ) (cid:1) / + C ( j )2 C ( j ) for t (I → F)T ( j ) ≤ t < t (F)T , ( j ) D t + e ( j ) for t (F)T , ( j ) ≤ t < t j +1 , e , e = const. ∈ R . (54)In case only one asymptotic end exists, that is, for G (I → F)T ( j ) < G j +1 ≤ G (F)T , ( j ) or G (F)T , ( j ) ≤ G (I → F)T ( j ) < G j +1 ,we get G ( t ) = (cid:0) C ( j ) t + e ( j ) (cid:1) / + C ( j )2 C ( j ) for t (I → F)T ( j ) ≤ t < t j +1 , e = const. ∈ R , (55)or G ( t ) = D t + e ( j ) for t (I → F)T ( j ) ≤ t < t j +1 , e = const. ∈ R , (56)respectively. We remark that the transition time (53) was computed by fixing the integration constant e = 1 D (cid:32) G (I → F)T − (cid:115) C D arctan (cid:18) β (cid:20) C G (I → F)T C − (cid:21)(cid:19)(cid:33) − t (I → F)T in Eq. (51) imposing the initial condition G ( t = t (I → F)T ) = G (I → F)T and substituting the transition value (52).The determination of the integration constants e , ..., e is given in Appendix D. C. Solution of the Relativistic Kinetic Equation for { t ∈ R ≥ } The complete solution of Eq. (8) is derived as follows. Beginning with the single-injection domain (SID), thesolution branch G ( t | ≤ t < t ) is given for 0 ≤ t < min (cid:0) t , t (N → I)T (1) (cid:1) by Sol. (26) (in which j = 1 as for theother SID solutions below) with C = C (1; S = { } ), for t (N → I)T (1) ≤ t < min (cid:0) t , t (I → F)T (1) (cid:1) by Sol. (38) with C = C (1; S = { } ), and for t (I → F)T (1) ≤ t < t by • Sol. (54) for t (I → F)T (1) < t (F)T , (1) < t , • Sol. (55) for t (I → F)T (1) < t ≤ t (F)T , (1) , • Sol. (56) for t (F)T , (1) ≤ t (I → F)T (1) < t with C = D and C = C (1; S = { } ). We point out that in the SID, Eq. (8) can also be solved by directlyapplying separation of variables, resulting in [35] (cid:90) G D G + A q d G = G D − D (cid:90) d G D G A q = G D − ( A q ) / D / arctan (cid:32)(cid:115) D A q G (cid:33) = t + f , where G := G + x and f = const. ∈ R . Using once again the method of matched asymptotic expansions withthe transition time t (S)T = 1 D (cid:32)(cid:114) A q D (cid:34) − π (cid:32)(cid:115) D A q x (cid:33)(cid:35) − x (cid:33) < t (S)T < t to the approximate solution G ( t ) = (cid:40) (3 A q t + f ) / − x for 0 ≤ t < t (S)T D t + f for t (S)T ≤ t < t , f , f = const. ∈ R , (57)for 0 < t ≤ t (S)T to G ( t ) = (3 A q t + f ) / − x for 0 ≤ t < t , f = const. ∈ R , (58)and otherwise for t (S)T ≤ < t to G ( t ) = D t + f for 0 ≤ t < t , f = const. ∈ R . (59)The integration constants f , f , f , and f are fixed by the initial condition G ( t = t = 0) = G = 0 f = x D − ( A q ) / D / arctan (cid:32)(cid:115) D A q x (cid:33) f = f = x f = 0 , whereas the integration constant f is determined by the transition condition G (cid:0) t = t (S)T | ≤ t < t (S)T (cid:1) = G (cid:0) t = t (S)T | t (S)T ≤ t < t (cid:1) f = (cid:0) A q t (S)T + x (cid:1) / − x − D t (S)T . 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 t ≤ t < min (cid:0) t , t (N → I)T (2) (cid:1) by Sol. (26) (with j = 2 as for the otherDID solutions below) if S (cid:54) = ∅ , S = ∅ and by • Sol. (31) for t < t (N)T (2) < min (cid:0) t , t (N → I)T (2) (cid:1) , • Sol. (32) for t < min (cid:0) t , t (N → I)T (2) (cid:1) ≤ t (N)T (2) , • Sol. (33) for t (N)T (2) ≤ t < min (cid:0) t , t (N → I)T (2) (cid:1) if S = ∅ , S (cid:54) = ∅ (with the specific arguments of the constants C and C associated to the respective cases asfor the other DID solutions below). For t (N → I)T (2) ≤ t < min (cid:0) t , t (I → F)T (2) (cid:1) , we employ Sol. (38) if S (cid:54) = ∅ , S = ∅ and • Sol. (39) for t (N → I)T (2) < t (I)T (2) < min (cid:0) t , t (I → F)T (2) (cid:1) , • Sol. (40) for t (N → I)T (2) < min (cid:0) t , t (I → F)T (2) (cid:1) ≤ t (I)T (2) , • Sol. (41) for t (I)T (2) ≤ t (N → I)T (2)if S = ∅ , S (cid:54) = ∅ . Finally, for t (I → F)T (2) ≤ t < t , we apply • Sol.(46) for t (I → F)T (2) < t (F)T , (2) < t , • Sol.(47) for t (I → F)T (2) < t ≤ t (F)T , (2) , • Sol.(48) for t (F)T , (2) ≤ t (I → F)T (2) < t if either S (cid:54) = ∅ , S = ∅ or S = ∅ , S (cid:54) = ∅ and • Sol.(54) for t (I → F)T (2) < t (F)T , (2) < t , • Sol.(55) for t (I → F)T (2) < t ≤ t (F)T , (2) , • Sol.(56) for t (F)T , (2) ≤ t (I → F)T (2) < t S = ∅ = S . The initial value G is fixed by requiring continuity of the SID and DID solution branches atthe time of the second injection, yielding G = G (cid:0) t = t | t (S)T ≤ t < t ; Sol.(57) (cid:1) for 0 < t (S)T < t G (cid:0) t = t | ≤ t < t ; Sol.(58) (cid:1) for 0 < t ≤ t (S)T G (cid:0) t = t | ≤ t < t ; Sol.(59) (cid:1) for t (S)T ≤ < t . In addition to the initial DID updating of constants at t = t , we have to perform NID-IID updates at t = t (N → I)T (1) and/or t = t (N → I)T (2) if (cid:8) t (N → I)T (1) , t (N → I)T (2) (cid:9) ∩ [ t , t ) (cid:54) = ∅ as well as IID-FID updates at t = t (I → F)T (1)and/or t = t (I → F)T (2) if (cid:8) t (I → F)T (1) , t (I → F)T (2) (cid:9) ∩ [ t , t ) (cid:54) = ∅ . Assuming, for example, that the transition times t (N → I)T (1), t (I → F)T (1), and t (N → I)T (2) are contained in this interval with the order t < t (N → I)T (1) < t (N → I)T (2) 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 [32] that,for all frequencies and times, it provides only a small contribution to the high-energy SSC component. A. Synchrotron Intensity The optically thin synchrotron intensity I syn. ( (cid:15), t ) (with [ I syn. ] = eV s − cm − sr − and similar for the SSCintensity) for an isotropically distributed electron number density is given by I syn. ( (cid:15), t ) = R π (cid:90) ∞ n ( γ, t ) P ( (cid:15), γ ) d γ , (61)where the function P ( (cid:15), γ ) = P (cid:15)γ CS (cid:18) (cid:15) (cid:15) γ (cid:19) (62)is the pitch-angle-averaged spectral synchrotron power of a single electron in a magnetic field of strength b , (cid:15) := E syn. / ( m e c ) is the normalized synchrotron photon energy, P := 8 . × eV s − , and (cid:15) := 2 . × − b [9]. The CS function is discussed in detail in Appendix F. Here, we employ the approximation CS ( z ) ≈ a z / (cid:0) z / exp ( z ) (cid:1) , z ∈ R ≥ , (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 I syn. ( (cid:15), t ) = I , syn. (cid:15) / m (cid:88) i =1 q i H ( t − t i ) Y / i ( t )1 + (cid:18) (cid:15) (cid:15) (cid:19) / Y / i ( t ) exp (cid:18) (cid:15) (cid:15) Y i ( t ) (cid:19) (64)with I , syn. := 3 / R P a (cid:15) / / (2 / π ) and the abbreviation Y i ( t ) := G ( t ) − G i + x i . For comparisons withobservational data or for generic case studies, that is, for fitting or plotting lightcurves, we have to compute theenergy-integrated synchrotron intensity¯ I syn. ( t ; (cid:15) min. , (cid:15) max. ) = (cid:90) (cid:15) max. (cid:15) min. I syn. ( (cid:15), t ) d (cid:15) with a lower integration limit (cid:15) min. corresponding to the energy of the first data point in the fluence SED andan upper integration limit (cid:15) max. defined by the last. Using (64), this quantity yields¯ I syn. ( t ; (cid:15) min. , (cid:15) max. ) = (cid:18) (cid:15) (cid:19) / I , syn. m (cid:88) i =1 q i H ( t − t i ) Y i ( t ) (cid:90) τ i, max. τ i, min. (cid:101) τ / (cid:101) τ / exp ( (cid:101) τ ) d (cid:101) τ , where τ i, min. := 2 (cid:15) min. Y i / (3 (cid:15) ) and τ i, max. := 2 (cid:15) max. Y i / (3 (cid:15) ). In order to solve the integral, we approximatethe integrand by (cid:101) τ / for (cid:101) τ ≤ − (cid:101) τ ) for (cid:101) τ > 1. This is justified because the approximated CS function(63) is only adapted to the asymptotics (cid:101) τ (cid:28) (cid:101) τ (cid:29) CS function and extrapolated in between.6Moreover, since τ i, min. and τ i, max. depend on the time t , we have to consider the three cases where 1 ≤ τ i, min. , τ i, min. < ≤ τ i, max. , and τ i, max. < 1. Accordingly, the integral results in (cid:90) τ i, max. τ i, min. (cid:101) τ / (cid:101) τ / exp ( (cid:101) τ ) d (cid:101) τ ≈ H ( τ i, min. − (cid:2) exp ( − τ i, min. ) − exp ( − τ i, max. ) (cid:3) + H ( τ i, max. − H (1 − τ i, min. ) (cid:20) (cid:16) − τ / i, min. (cid:17) − exp ( − τ i, max. ) + exp ( − (cid:21) + 34 H (1 − τ i, max. ) (cid:104) τ / i, max. − τ / i, min. (cid:105) . (65) B. Synchrotron Self-Compton Intensity In the computation of the SSC intensity I SSC ( (cid:15) s , t ) = R π (cid:90) ∞ n ( γ, t ) P SSC ( (cid:15) s , γ, t ) d γ , (66)where P SSC ( (cid:15) s , γ, t ) is the SSC power of a single electron and (cid:15) s := E s / ( m e c ) is the normalized scattered photonenergy, we have to employ the Thomson limit because the SSC radiative losses in (2) are already restricted tothe Thomson regime. In this limit, the SSC power reads [9, 20] P SSC ( (cid:15) s , γ, t ) = 43 σ T c γ E ( (cid:15) s , t )with the Thomson cross section σ T = 6 . × − cm and the total SSC photon energy density E ( (cid:15) s , t ) (havingthe dimension [ E ] = eV cm − ). For ultrarelativistic electrons with γ (cid:29) γ (cid:15) (cid:28) 1, the characteristic energy of the SSC-scattered photons is (cid:15) s ≈ γ (cid:15) [20],corresponding to head-on collisions of the synchrotron photons with the electrons [9, 24, 29]. Thus, it is justifiedto apply a monochromatic approximation in the total SSC photon energy density in form of a Dirac distributionthat spikes at this characteristic energy E ( (cid:15) s , t ) = 14 π (cid:90) ∞ (cid:15) N ( (cid:15), t ) δ (cid:0) (cid:15) s − γ (cid:15) (cid:1) d (cid:15) , where N ( (cid:15), t ) = 4 π I syn. ( (cid:15), t ) c (cid:15) is the synchrotron photon number density. We remark that by assuming an isotropic, ultrarelativistic electrondistribution, the synchrotron photon number density becomes inevitably isotropically distributed, too. Substi-tuting the latter formulas into (66) and using Fubini’s theorem, we obtain I SSC ( (cid:15) s , t ) = R σ T (cid:15) s π (cid:90) ∞ I syn. ( (cid:15), t ) (cid:15) (cid:90) /(cid:15) n ( γ, t ) δ (cid:0) (cid:15) s − γ (cid:15) (cid:1) d γ d (cid:15) , (67)in which the upper γ -integration limit arises from the restriction to the Thomson regime. With the electronnumber density (10) and the synchrotron intensity (64), the SSC intensity (67) yields, after having performedboth integrations, I SSC ( (cid:15) s , t ) = I , SSC (cid:15) / m (cid:88) i,j =1 q i q j H ( t − t i ) H ( t − t j ) H (cid:16) − (cid:15) s Y j ( t ) (cid:17) (cid:2) Y i ( t ) Y j ( t ) (cid:3) / (cid:18) (cid:15) s (cid:15) (cid:19) / (cid:2) Y i ( t ) Y j ( t ) (cid:3) / exp (cid:18) (cid:15) s (cid:15) (cid:2) Y i ( t ) Y j ( t ) (cid:3) (cid:19) , (68)where I , SSC := R σ T I , syn. / (2 / π ). The double sum, with the index i referring to the i th synchrotronphoton population and the index j to the j th electron population, accounts for all combinations of SSC scattering7between the various electron and synchrotron photon populations. For the corresponding energy-integrated SSCintensity, we find ¯ I SSC ( t ; (cid:15) s , min. , (cid:15) s , max. ) = (6 (cid:15) ) / I , SSC m (cid:88) i,j =1 q i q j H ( t − t i ) H ( t − t j ) (cid:2) Y i ( t ) Y j ( t ) (cid:3) × (cid:90) min( τ ij, max. , Y i Y j / (3 (cid:15) )) τ ij, min. (cid:101) τ / (cid:101) τ / exp ( (cid:101) τ ) d (cid:101) τ , where τ ij, min. := (cid:15) s , min. ( Y i Y j ) / (6 (cid:15) ) and τ ij, max. := (cid:15) s , max. ( Y i Y j ) / (6 (cid:15) ). The integral is given by (65),however, with adapted integration limits. IV. SYNCHROTRON AND SYNCHROTRON SELF-COMPTON FLUENCES We compute the total fluences associated with the synchrotron intensity (64) and the SSC intensity (68). Forthis 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 ofthe integration measure. For simplicity, the updating of constants is once more suppressed. With ( ε, I, F ) ∈ (cid:8) ( (cid:15), I syn. , F syn. ) , ( (cid:15) s , I SSC , F SSC ) (cid:9) , the total fluence F (for which [ F ] = eV cm − sr − ) is given by F ( ε ) = (cid:90) ∞ I ( ε, t ) d t = m (cid:88) i =1 (cid:90) G i +1 G i I ( ε, G ) d t d G d G . (69)The Jacobian determinant yieldsd t d G = ( G + x ) D ( G + x ) + A q for 0 ≤ G < G C ( i ) for the domains of validity of Eqs. (24) , (25) , (34) , and (35) G C ( i ) G + C ( i ) for the domains of validity of Eqs. (21) , (22) , (36) , (37) , (43)-(45) G D G + C ( i ) G + C ( i ) for G (I → F)T ( i ) ≤ G < G i +1 and S = ∅ = S , (70)where i : 2 ≤ i ≤ m . Substituting (70) into (69), we obtain the expression F ( ε ) = (cid:90) G ( G + x ) I ( ε, G ) D ( G + x ) + A q d G + m (cid:88) i =2 (cid:34)(cid:90) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1) G i B NID ( G ; i ) I ( ε, G ) d G + (cid:90) min (cid:0) G i +1 ,G (I → F)T ( i ) (cid:1) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1) B IID ( G ; i ) I ( ε, G ) d G + (cid:90) G i +1 min (cid:0) G i +1 ,G (I → F)T ( i ) (cid:1) B FID ( G ; i ) I ( ε, G ) d G (cid:35) , in which B NID ( G ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) C ( i ; S ) + χ ( S ) (cid:2) − χ ( S ) (cid:3) C ( i ; S , S )+ (cid:2) − χ ( S ) (cid:3) χ ( S ) G C ( i ; S ) G + C ( i ; S ) + χ ( S ) χ ( S ) G C ( i ; S , S ) G + C ( i ; S ) , B IID ( G ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) C ( i ; S ) + χ ( S ) (cid:2) − χ ( S ) (cid:3) C ( i ; S , S )+ (cid:2) − χ ( S ) (cid:3) χ ( S ) G C ( i ; S ) G + C ( i ; S ) + χ ( S ) χ ( S ) G C ( i ; S , S ) G + C ( i ; S ) , and B FID ( G ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) G D G + C ( i ; S ) G + C ( i ; S ) + χ ( S ) (cid:2) − χ ( S ) (cid:3) G C ( i ; S ) G + C ( i ; S )+ (cid:2) − χ ( S ) (cid:3) χ ( S ) G C ( i ; S ) G + C ( i ; S ) + χ ( S ) χ ( S ) G C ( i ; S , S ) G + C ( i ; S )with the characteristic function χ ( S k ) := (cid:40) S k (cid:54) = ∅ S k = ∅ , k ∈ { , , } . In the following, for illustrative purposes, we calculate in detail both the synchrotron and the SSC NID fluenceintegrals for the cases S (cid:54) = ∅ , S = ∅ . The remaining integrals can be solved in an analogous manner. A. Synchrotron Fluence With I = I syn. ( (cid:15), G ) according to (64), the synchrotron NID fluence integral for S (cid:54) = ∅ , S = ∅ reads I NIDsyn. ( (cid:15) ; i ) := 1 C ( i ) (cid:90) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1) G i I syn. ( (cid:15), G ) d G = I , syn. (cid:15) / C ( i ) i (cid:88) l =1 q l (cid:90) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1) G i Y / l ( G )1 + (cid:18) (cid:15) (cid:15) (cid:19) / Y / l ( G ) exp (cid:18) (cid:15) (cid:15) Y l ( G ) (cid:19) d G . Using τ l := 2 (cid:15) Y l / (3 (cid:15) ) gives I NIDsyn. ( (cid:15) ; i ) = (3 (cid:15) ) / I , syn. / C ( i ) (cid:15) − / i (cid:88) l =1 q l (cid:90) τ l (cid:0) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1)(cid:1) τ l ( G i ) d (cid:101) τ (cid:101) τ / (cid:0) (cid:101) τ / exp ( (cid:101) τ ) (cid:1) . (71)Similar to the evaluation of the energy-integrated intensities (cf. the paragraph directly above formula (65)), weapproximate the integrand for (cid:101) τ ≤ (cid:101) τ − / and for (cid:101) τ > (cid:101) τ − / exp ( − (cid:101) τ ). Thus, solving the integral, (71)yields I NIDsyn. ( (cid:15) ; i ) ≈ (3 (cid:15) ) / I , syn. / C ( i ) (cid:15) − / i (cid:88) l =1 q l (cid:34) H (cid:0) τ l ( G i ) − (cid:1) Γ (cid:18) , τ l ( G i ) , τ l (cid:0) min (cid:0) G i +1 , G (N → I)T ( i ) (cid:1)(cid:1)(cid:19) + H (cid:16) τ l (cid:0) min (cid:0) G i +1 , G (N → I)T ( i ) (cid:1)(cid:1) − (cid:17) H (cid:0) − τ l ( G i ) (cid:1) (cid:20) Γ (cid:18) , , τ l (cid:0) min (cid:0) G i +1 , G (N → I)T ( i ) (cid:1)(cid:1)(cid:19) + 65 (cid:0) − τ / l ( G i ) (cid:1)(cid:21) + 65 H (cid:16) − τ l (cid:0) min (cid:0) G i +1 , G (N → I)T ( i ) (cid:1)(cid:1)(cid:17) (cid:104) τ / l (cid:0) min (cid:0) G i +1 , G (N → I)T ( i ) (cid:1)(cid:1) − τ / l ( G i ) (cid:105)(cid:35) with the generalized incomplete gamma function defined for a ∈ C and Re( a ) > a, y, z ) := (cid:90) zy t a − exp ( − t ) d t . a = 1 / 2, the generalized incomplete gamma function can be expressed in terms of the errorfunction, namely Γ(1 / , y, z ) = √ π (cid:2) erf( z ) − erf( y ) (cid:3) . B. Synchrotron Self-Compton Fluence With I = I SSC ( (cid:15) s , G ) given in (68), the SSC NID fluence integral for S (cid:54) = ∅ , S = ∅ becomes I NIDSSC ( (cid:15) s ; i ) := 1 C ( i ) (cid:90) min (cid:0) G i +1 ,G (N → I)T ( i ) (cid:1) G i I SSC ( (cid:15) s , G ) d G = I , SSC (cid:15) / C ( i ) i (cid:88) k,l =1 q k q l H (cid:18) (cid:15) s − G i + G l − x l (cid:19) × (cid:90) min (cid:0) G i +1 ,G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1) G i (cid:2) Y k ( G ) Y l ( G ) (cid:3) / (cid:18) (cid:15) s (cid:15) (cid:19) / (cid:2) Y k ( G ) Y l ( G ) (cid:3) / exp (cid:18) (cid:15) s (cid:15) (cid:2) Y k ( G ) Y l ( G ) (cid:3) (cid:19) d G . An approximate analytical solution of the integral may be derived with the same method as in the case of thesynchrotron fluence. However, one could also employ the methods used for the computations of the NID-IIDand IID-FID transition times, which are explained in more detail in Appendix C, as follows. Applying themean-value-theorem method, we first define the function C ∞ (cid:0) R ≥ (cid:1) (cid:51) H kl ( (cid:15) s , G ) := (cid:2) Y k ( G ) Y l ( G ) (cid:3) / (cid:18) (cid:15) s (cid:15) (cid:19) / (cid:2) Y k ( G ) Y l ( G ) (cid:3) / exp (cid:18) (cid:15) s (cid:15) (cid:2) Y k ( G ) Y l ( G ) (cid:3) (cid:19) . Next, we assert that there exists a mean value ξ il ∈ (cid:2) G i , min (cid:0) G i +1 , G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1)(cid:3) such that I NIDSSC ( (cid:15) s ; i ) = I , SSC (cid:15) / C ( i ) i (cid:88) k,l =1 q k q l H (cid:18) (cid:15) s − G i + G l − x l (cid:19) H kl ( (cid:15) s , ξ il ) × (cid:104) min (cid:0) G i +1 , G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1) − G i (cid:105) . Lastly, we choose the midpoint of the integration interval (cid:2) G i + min (cid:0) G i +1 , G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1)(cid:3) / ξ il . The trapezoid-approximation method on the other hand results in the expression I NIDSSC ( (cid:15) s ; i ) ≈ I , SSC (cid:15) / C ( i ) i (cid:88) k,l =1 q k q l H (cid:18) (cid:15) s − G i + G l − x l (cid:19) min (cid:0) G i +1 , G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1) − G i × (cid:104) H kl ( (cid:15) s , G i ) + H kl (cid:0) (cid:15) s , min (cid:0) G i +1 , G (N → I)T ( i ) , /(cid:15) s + G l − x l (cid:1)(cid:1)(cid:105) . V. LIGHTCURVES AND FLUENCE SPECTRAL ENERGY DISTRIBUTIONS To obtain realistic lightcurves, we have to consider injections of finite duration time and radiative transport insidethe emission region. But since, for reasons of computational feasibility, we have only employed Dirac distributionsfor 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 i th flare into n i separate injections labeled by a subscript p ∈ { , ..., n i } , which0are induced equidistantly over the entire emission region given by 2 R , i.e., the p th injection of the i th flareoccurs at the time t i + κ ip with κ ip := p − n i − R c , and for each of which we use the quantities derived in Section III. Thus, the energy-integrated synchrotronintensity reads I syn. ( t ; (cid:15) min. , (cid:15) max. ) = (cid:90) (cid:15) max. (cid:15) min. I syn. ( (cid:15), t ) d (cid:15) = (cid:90) (cid:15) max. (cid:15) min. m (cid:88) i =1 n i (cid:88) p =1 q i ( p ) H ( t − t i − κ ip ) I ip ( (cid:15), t ) d (cid:15) , (72)where (cid:0) q i ( p ) (cid:1) p ∈{ ,...,n i } for all i : 1 ≤ i ≤ m is a specific distribution satisfying the constraint q i = n i (cid:88) p =1 q i ( p ) (73)and the function I ip ( (cid:15), t ) is, according to the original synchrotron intensity (64), defined by I ip ( (cid:15), t ) := I , syn. (cid:15) / Y / ip ( t )1 + (cid:18) (cid:15) (cid:15) (cid:19) / Y / ip ( t ) exp (cid:18) (cid:15) (cid:15) Y ip ( t ) (cid:19) (74)with Y ip ( t ) := G ( t ) − G ( t i + κ ip ) + x i . Similarly, one may construct a proper formula for the energy-integratedSSC intensity. We refrain from showing a parameter study for lightcurves because an adequate value for n i ,which is at least of the order O (10 ), results in numerical computations that would exceed our available CPUtime by far. For the associated total fluence SEDs, it is an entirely different matter as we can illustrate theirfunctional shapes by the special case n i = 1 for all i : 1 ≤ i ≤ m , in which each flare is modeled by a singleinjection (cf. Section IV), given that it yields lower and upper bounds. In the following, we first show by directcalculation that the total fluence of the synchrotron intensity I syn. ( (cid:15), t ) defined in (72) is indeed bounded frombelow and above, up to specific constant factors, by the simple n i = 1 case. To this end, we substitute I syn. ( (cid:15), t )into the expression for the total fluence and employ the variable G F syn . ( (cid:15) ) = (cid:90) ∞ I syn . ( (cid:15), t ) d t = m (cid:88) i =1 n i (cid:88) p =1 q i ( p ) (cid:90) ∞ t i + κ ip I ip ( (cid:15), t ) d t = m (cid:88) i =1 n i (cid:88) p =1 q i ( p ) (cid:90) ∞ G ( t i + κ ip ) I ip ( (cid:15), G ) d t d G d G . (75)From Eq. (11), it can be deduced that the Jacobian determinant is positive and bounded from below and aboveby (cid:18) D + A j (cid:88) i =1 q i x i (cid:19) − ≤ d t d G ≤ D , (76)and since I ip is non-negative, we can, therefore, bound the fluence (75) from below and above by m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) − n i (cid:88) p =1 q i ( p ) (cid:90) ∞ G ( t i + κ ip ) I ip ( (cid:15), G ) d G ≤ F syn . ( (cid:15) ) ≤ D m (cid:88) i =1 n i (cid:88) p =1 q i ( p ) (cid:90) ∞ G ( t i + κ ip ) I ip ( (cid:15), G ) d G . Then, applying the transformation (cid:101) G ip := G − G ( t i + κ ip ) + G i and defining I i := I i (cf. formula (74)), weobtain m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) − n i (cid:88) p =1 q i ( p ) (cid:90) ∞ G i I i ( (cid:15), (cid:101) G ip ) d (cid:101) G ip ≤ F syn . ( (cid:15) ) ≤ D m (cid:88) i =1 n i (cid:88) p =1 q i ( p ) (cid:90) ∞ G i I i ( (cid:15), (cid:101) G ip ) d (cid:101) G ip . Because the integral has the same value for fixed i and arbitrary p , we can always drop the subscript p in theargument of the integrand and the integration measure. Hence, by means of the constraint (73), we get m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) − q i (cid:90) ∞ G i I i ( (cid:15), G ) d G = m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) − (cid:18) n i (cid:88) p =1 q i ( p ) (cid:19) (cid:90) ∞ G i I i ( (cid:15), (cid:101) G i ) d (cid:101) G i ≤ F syn . ( (cid:15) )1 -17 -15 -13 -11 -9 -7 -5 -3 -1 E [GeV] -7 -5 -3 -1 E · F [ a r b . u . ] b =0 . b =0 . b =1 (a) -17 -15 -13 -11 -9 -7 -5 -3 -1 E [GeV] -7 -5 -3 -1 E · F [ a r b . u . ] D =5 D =10 D =20 (b) -17 -15 -13 -11 -9 -7 -5 -3 -1 E [GeV] -7 -5 -3 -1 E · F [ a r b . u . ] x =6 · − x =10 − x =10 − (c) FIG. 1: Total synchrotron (left curves) and SSC (right curves) fluence SEDs for generic three-injection scenarios withmagnetic field strengths, Doppler boost factors, and reciprocal initial electron energies b ∈ { . , . , } , D = 10, and( x i ) i =1 , , = (10 − , − , − ) for (a), b = 1, D ∈ { , , } , and ( x i ) i =1 , , = (10 − , − , − ) for (b), as well as b = 1, D = 10, x ∈ { . × − , − , − } , and ( x i ) i =2 , = (10 − , − ) for (c). The injection times are always givenby ( t i ) i =1 , , = (0 , . , 3) 2 R /c with the characteristic radius of the plasmoid R = 10 cm, the injection strengthsby ( q i ) i =1 , , = (1 . × cm − , × cm − , × cm − ), and the cooling rate prefactors D = 1 . × − b s − , A = 1 . × − b cm s − depend on the specific value of the magnetic field strength. The number of grid points usedin the plotting of each curve amounts to 4 × . F syn . ( (cid:15) ) ≤ D m (cid:88) i =1 (cid:18) n i (cid:88) p =1 q i ( p ) (cid:19) (cid:90) ∞ G i I i ( (cid:15), (cid:101) G i ) d (cid:101) G i = 1 D m (cid:88) i =1 q i (cid:90) ∞ G i I i ( (cid:15), G ) d G for the upper bound. Transforming the integration measure back to the time variable t and using the inverse of(76) leads to D m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) − q i (cid:90) ∞ t i I i ( (cid:15), t ) d t ≤ F syn . ( (cid:15) ) ≤ D m (cid:88) i =1 (cid:18) D + A i (cid:88) k =1 q k x k (cid:19) q i (cid:90) ∞ t i I i ( (cid:15), t ) d t . Since the second sum on both the left-hand side and the right-hand side can be bounded from above by (cid:80) mk =1 q k /x k , we find (cid:18) A D m (cid:88) k =1 q k x k (cid:19) − F syn . ( (cid:15) ) ≤ F syn . ( (cid:15) ) ≤ (cid:18) A D m (cid:88) k =1 q k x k (cid:19) F syn . ( (cid:15) ) , (77)where F syn . ( (cid:15) ) = (cid:80) mi =1 (cid:82) ∞ t i q i I i ( (cid:15), t ) d t , proving the claim.In Figure 1 (a)-(c), we show a parameter study for the total synchrotron and SSC fluence SEDs (given inarbitrary units), always considering three injections of ultrarelativistic electron populations into the emissionregion, each modeling one flare according to the estimate (77). We separately vary the nondimensional magneticfield strength b , the Doppler factor of the plasmoid D , and the first reciprocal normalized initial electron energy x , leaving the remaining free parameters fixed. Details on the values of these parameters and the specificvariations can be found in the figure caption. Note that in all three subfigures, the black, solid SEDs correspondto the same set of parameters, thus, serving as a reference curve in the parameter study. Varying only themagnetic field strength (Figure 1 (a)), toward higher values, we can see an increase of both the maximumsynchrotron 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 synchrotronemission. However, the maximum SSC energy increases inasmuch as the maximum synchrotron energy can nowreach larger values that add to the SSC scattering process. Further, shifting the Doppler factor toward highervalues results in the typical increase in the apparent luminosity (Figure 1 (b)), which is well-known from studiesof the physical effects of relativistic beaming, i.e., aberration, the Doppler effect, time dilation, and retardation.In Figure 1 (c), as we vary the reciprocal normalized initial electron energy of the first injection, towards largervalues, both the maximum synchrotron energy as well as the maximum SSC energy increase, hence, broadeningthe SED curves on their right slopes according to the specific nonlinear cooling behavior and the strength of theparameter variation. Note that we could also create a plot in which we vary the strengths q i , i ∈ { , , } , ofthe different injections. However, we refrain from showing such a plot because increasing the injection strengthstowards larger values obliges us to raise the number of grid points (i.e., the number of time steps, which iscurrently 4 × amounting to a computation time of approximately four days) in direct proportion, as theinjection strengths and the time variable are coupled by multiplication. This would result in computation timesranging from several months to years. We further remark in passing that in order to detect patterns due tovariations in the number of injections or due to their nonlinear coupling, a larger parameter study than the onepresented here is necessary. Finally, but most importantly, the resulting functional shapes of the various plotsindicate that our model can reproduce the typical fluence SEDs of blazars visible in observational data. VI. SUMMARY AND OUTLOOK We introduced a fully analytical, time-dependent leptonic one-zone model for the flaring of blazars that employscombined synchrotron and SSC radiative losses of multiple interacting, ultrarelativistic electron populations. Ourmodel assumes several injections of electrons into the emission region as the cause of the flaring, which differsfrom 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 of the order of tens of kiloparsecs and, thus, it is mostlikely that there is a pick up of more than just one particle population from interstellar and intergalactic clouds.At the same time, it further assumes both radiative cooling processes to occur simultaneously, as would be thecase in any physical scenario. In more detail, applying Laplace transformations and the method of matched3asymptotic expansions, we derived an approximate analytical solution of the relativistic kinetic equation of thevolume-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 Thomsonlimit. Using this solution, we computed the optically thin synchrotron intensity, the SSC intensity in theThomson limit, as well as the corresponding total fluences. Moreover, we mimicked finite injection durations andradiative transport by modeling flares in terms of sequences of instantaneous injections. Ultimately, we presenteda parameter study for the total synchrotron and SSC fluence SEDs for a generic three-injection scenario withvariations of the magnetic field strength, the Doppler factor, and the initial electron energies, showing thatour model can reproduce the characteristic broad-band SED shapes seen in observational data. We point outthat the SSC radiative loss term considered here is strictly valid only in the Thomson regime, limiting theapplicability of the model to at most GeV blazars. Nonetheless, it can be generalized to describe TeV blazarsby using the full Klein-Nishina cross section in the SSC energy loss rate. This leads to a model for whichsimilar yet technically more involved methods apply. Further, in order to make our simple analytical model morerealistic, terms accounting for spatial diffusion and for electron escape could be added to the kinetic equation.Also, more elaborate source functions, e.g., with a power law energy dependence, a time dependence in form ofrectangular 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 directnumerical evaluation of the associated kinetic equation. Acknowledgments The authors are grateful to Horst Fichtner, Reinhard Schlickeiser, and Michael Zacharias for useful discussionsand comments. Appendix A: Laplace Transformation Method We derive the solution R ( x, G ) of the PDE (6) by using a composition of two Laplace transformations. First,applying a Laplace transformation with respect to the reciprocal electron energy x L w ( x ) [ · ] := (cid:90) ∞ [ · ] exp ( − w x ) d x to Eq. (6) gives (cid:90) ∞ ∂R∂G exp ( − w x ) d x + (cid:90) ∞ ∂R∂x exp ( − w x ) d x = m (cid:88) i =1 q i δ ( G − G i ) (cid:90) ∞ δ ( x − x i ) exp ( − w x ) d x . Evaluating the second integral on the left-hand side via integration by parts and employing the Laplace transformof R K ( w, G ) := (cid:90) ∞ R ( x, G ) exp ( − w x ) d x , (A1)we obtain ∂K∂G + R exp ( − w x ) | x →∞ − R (0 , G ) + w K = m (cid:88) i =1 q i δ ( G − G i ) exp ( − w x i ) (cid:2) H ( x − x i ) | x →∞ − H ( x − x i ) | x =0 (cid:3) and, hence, ∂K∂G − R (0 , G ) + w K = m (cid:88) i =1 q i δ ( G − G i ) exp ( − w x i ) . (A2)4As the normalized initial electron energies are finite and bounded from above by γ i < . × b − / 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 supportsupp n ( γ, t ) = (cid:8) ( γ, t ) ∈ R ≥ | γ ≤ γ max. with γ max. := max { γ i | ≤ i ≤ m } (cid:9) . Thus, we can write n ( γ, t ) = H ( γ max. − γ ) n ( γ, t ), yielding for the second term on the left-hand side of Eq. (A2) R (0 , G ) = lim x (cid:38) (cid:18) H ( x − x max. ) n ( x, G ) x (cid:19) = 0 , where x max. := 1 /γ max. . Accordingly, we find ∂K∂G + w K = m (cid:88) i =1 q i δ ( G − G i ) exp ( − w x i ) . (A3)Secondly, applying a Laplace transformation with respect to the function G L s ( G ) [ · ] := (cid:90) ∞ [ · ] exp ( − s G ) d G to Eq. (A3) results in (cid:90) ∞ ∂K∂G exp ( − s G ) d G + w (cid:90) ∞ K exp ( − s G ) d G = m (cid:88) i =1 q i exp ( − w x i ) (cid:90) ∞ δ ( G − G i ) exp ( − s G ) d G . (A4)This Laplace transformation implies that G is non-negative. With the Laplace transform of KM ( w, s ) := (cid:90) ∞ K ( w, G ) exp ( − s G ) d G and integration by parts as before, Eq. (A4) reads K exp ( − s G ) | G →∞ − K ( w, 0) + ( w + s ) M = m (cid:88) i =1 q i exp (cid:0) − ( w x i + s G i ) (cid:1)(cid:2) H ( G − G i ) | G →∞ − H ( G − G i ) | G =0 (cid:3) . Since the Heaviside functions are not well-defined at the jump discontinuities at G = G i for i : 1 ≤ i ≤ m , thisequation reduces to − K ( w, 0) + ( w + s ) M = m (cid:88) i =1 q i exp (cid:0) − ( w x i + s G i ) (cid:1) − q exp ( − w x ) H ( G ) | G =0 , (A5)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 n = q x . Substituting this density into (A1) by employing the relation R ( x, G ) = n ( x, G ) /x , we get K ( w, 0) = n (cid:90) ∞ exp ( − w x ) x δ ( x − x ) d x = q exp ( − w x ) . Choosing H ( G ) | G =0 = 1, we can use the last term on the right-hand side of Eq. (A5) in order to compensatethis function, yielding M ( w, s ) = 1 w + s m (cid:88) i =1 q i exp (cid:0) − ( w x i + s G i ) (cid:1) . We point out that both M and the sum are positive. Hence, s > − w , which allows us to apply the inverseLaplace transformations with respect to the variables s and w to M . This gives the solution of Eq. (6) R ( x, G ) = L − w ( x ) L − s ( G ) M ( w, s ) = L − w ( x ) m (cid:88) i =1 q i H ( G − G i ) exp (cid:0) − w ( G − G i + x i ) (cid:1) = m (cid:88) i =1 q i H ( G − G i ) δ ( x − x i − G + G i ) . Appendix B: Approximation Errors and Optimal Domains The errors of the NID, IID, and FID approximations (13), (14), and (16) yield E NID u ( G ) = ( G − G u ) x u ( G − G u + x u ) (cid:20) x u + 1 G − G u + x u (cid:21) E IID v ( G ) = G − G v − x v x v ( G − G v + x v ) (cid:20) G − G v x v − G − G v + x v (cid:21) E FID w ( G ) = ( G w − x w ) G ( G − G w + x w ) (cid:20) G + 1 G − G w + x w (cid:21) . Evaluating the NID and IID errors at the NID-IID transition point G = G (N → I)T and the IID and FID errors atthe IID-FID transition point G = G (I → F)T , we obtain E NID u (cid:0) G (N → I)T ( u ) (cid:1) = 2 + 3 ξ (N → I) u x u ξ (N → I) u (1 + ξ (N → I) u ) E IID v (cid:0) G (N → I)T ( v ) (cid:1) = ξ (N → I) v − x v ( ξ (N → I) v + 1) (cid:34) ξ (N → I) v ξ (N → I) v + 1 − ξ (N → I) v (cid:35) and E IID v (cid:0) G (I → F)T ( v ) (cid:1) = ξ (I → F) v − x v ( ξ (I → F) v + 1) (cid:34) ξ (I → F) v − ξ (I → F) v + 1 (cid:35) E FID w (cid:0) G (I → F)T ( w ) (cid:1) = ( G w − x w ) x w ( ξ (I → F) w + 1) ( G w + ξ (I → F) w x w ) (cid:34) x w ( ξ (I → F) w + 1) + 2 G w + ξ (I → F) w x w (cid:35) . We derive optimal values for the interval parameters ξ (N → I) i and ξ (I → F) i by requiring that the total error at eachtransition point becomes minimal. To this end, we first impose the sufficient conditions ∂∂ξ (N → I) i (cid:16)(cid:13)(cid:13) E NID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E IID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13)(cid:17) = 0 ∂∂ξ (I → F) i (cid:16)(cid:13)(cid:13) E IID i (cid:0) G (I → F)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E FID i (cid:0) G (I → F)T ( i ) (cid:1)(cid:13)(cid:13)(cid:17) = 0 , (B1)which result in the solution ξ (N → I) i ≈ . 73 and an equation for the zeros of a seventh-order polynomial in ξ (I → F) i given by( ξ (I → F) i ) +3 ξ (I → F) i ( ξ (I → F) i +1) − − G i − x i ) ( G i + ξ (I → F) i x i ) (cid:32) x i ( ξ (I → F) i + 1) (cid:0) G i + x i [3 + 5 ξ (I → F) i ] (cid:1) ( G i + ξ (I → F) i x i ) (cid:33) = 0 , (B2)respectively. We verify that the total error (cid:13)(cid:13) E NID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E IID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) has a minimum at ξ (N → I) i ≈ . 73 by means of the second derivative test ∂ ( ∂ξ (N → I) i ) (cid:16)(cid:13)(cid:13) E NID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E IID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13)(cid:17)(cid:12)(cid:12) ξ (N → I) i =4 . > . From the Abel-Ruffini theorem, we know that there exists no general analytical solution to Eq. (B2). Hence,one has to employ a root-finding algorithm to compute numerical approximations.6One may obtain a more accurate overall approximation of Eq. (12) by using the more general IID approximation1( G − G v + x v ) = 1 x v ∞ (cid:88) n =0 n ! d n d G n (cid:32)(cid:20) G − G v x v (cid:21) − (cid:33)(cid:12)(cid:12) G − G v = a v x v × ( G − G v − a v x v ) n ≈ x v (1 + a v ) (cid:34) − 21 + a v (cid:18) G − G v x v − a v (cid:19)(cid:35) with errors E IID v (cid:0) G (N → I)T ( v ) (cid:1) = 1 x v (cid:34) ( ξ (N → I) v ) (1 + ξ (N → I) v ) − ξ (N → I) v [1 + 3 a v ] − ξ (N → I) v (1 + a v ) (cid:35) E IID v (cid:0) G (I → F)T ( v ) (cid:1) = 1 x v (cid:34) ξ (I → F) v ) − a v − ξ (I → F) v (1 + a v ) (cid:35) , where the unspecified expansion point a v constitutes another degree of freedom. Note that our original IIDapproximation (14) corresponds to expansion points which were set to a v = 1 for all v ∈ S . Then, one candetermine the optimal values for all a v such that the total errors E tot i := (cid:13)(cid:13) E NID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E IID i (cid:0) G (N → I)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E IID i (cid:0) G (I → F)T ( i ) (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) E FID i (cid:0) G (I → F)T ( i ) (cid:1)(cid:13)(cid:13) , i ∈ { , ..., j } , become minimal by making the first and second derivate tests ∂E tot i ∂a i = 0 and ∂ E tot i ∂a i > Appendix C: NID-IID and IID-FID Transition Times In the following, we present two different methods for the determination of the transition times t (N → I)T ( j ) and t (I → F)T ( j ). To this end, it is advantageous to start from the integral equation representation of the ODE (11)evaluated at the respective transition point G (N → I)T ( j ) = G j + x j /ξ (N → I) j or G (I → F)T ( j ) = G j + ξ (I → F) j x j . In theformer NID-IID case, we obtain t (N → I)T ( j ) = (cid:90) G j + x j /ξ (N → I) j G j J − ( (cid:101) G ) d (cid:101) G with J − ( (cid:101) G ) := 1 J ( (cid:101) G ) = (cid:32) D + A j (cid:88) i =1 q i (cid:0) (cid:101) G − G i + x i (cid:1) (cid:33) − . (C1)The first method applies the first mean value theorem for integration. Thus, as J − ∈ C ∞ ( R ≥ ), there exists apoint p j ∈ (cid:2) G j , G j + x j /ξ (N → I) j (cid:3) such that the transition time (C1) can be written as t (N → I)T ( j ) = x j ξ (N → I) j J − ( p j ) . Because the mean value theorem is merely an existence theorem, we have to approximate the value of p j by meansof an additional input. Since J − is strictly increasing, a possible choice for p j is the midpoint G j + x j / (2 ξ (N → I) j )of the integration interval. The second method employs a trapezoid approximation. Again due to the strictly7increasing functional shape of J − , the integral in (C1) 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 G j + x j /ξ (N → I) j of the integration interval and the height J − ( G j ) A r = x j ξ (N → I) j J − ( G j )and the area A t ( j ) of the right-angled triangle with one cathetus given by the distance between the endpointsand the other one by the height J − ( G j + x j /ξ (N → I) j ) − J − ( G j ) A t = x j ξ (N → I) j (cid:0) J − ( G j + x j /ξ (N → I) j ) − J − ( G j ) (cid:1) . This leads to the formula t (N → I)T ( j ) = x j ξ (N → I) j (cid:0) J − ( G j ) + J − ( G j + x j /ξ (N → I) j ) (cid:1) for the transition time. To obtain a more accurate approximation, one could use additional supporting pointsin the interval (cid:2) G j , G j + x j /ξ (N → I) j (cid:3) , giving rise to a finer trapezoid decomposition of the integral. Note thatin the present study, the trapezoid method yields a better approximation of the transition time. The IID-FIDtransition time t (I → F)T ( j ) = (cid:90) G j + ξ (I → F) j x j G j J − ( (cid:101) G ) d (cid:101) G can be computed similarly, resulting, on the one hand, in the expression t (I → F)T ( j ) = ξ (I → F) j x j J − (cid:0) G j + ξ (I → F) j x j / (cid:1) for the mean value theorem method and, on the other hand, in the expression t (I → F)T ( j ) = ξ (I → F) j x j (cid:0) J − ( G j ) + J − ( G j + ξ (I → F) j x j ) (cid:1) 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)) aredetermined. 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 G (cid:0) t = t (N)T ( j ) | t j ≤ t < t (N)T ( j ) (cid:1) = G (cid:0) t = t (N)T ( j ) | t (N)T ( j ) ≤ t < t (N → I)T ( j ) (cid:1) . The initial condition provides continuity of G at the transition from the ( j − j th injection domain at t = t j and the transition condition between the nonlinear and linearNID solution branches at t = t (N)T ( j ). Similarly, the IID integration constants d , d , d , and d are determinedvia the respective initial condition G (cid:0) t = t (N → I)T ( j ) (cid:1) = G (N → I)T ( j ) and d via the transition condition G (cid:0) t = t (I)T ( j ) | t (N → I)T ( j ) ≤ t < t (I)T ( j ) (cid:1) = G (cid:0) t = t (I)T ( j ) | t (I)T ( j ) ≤ t < t (I → F)T ( j ) (cid:1) , where the former condition yieldscontinuity between the NID and IID solution branches at t = t (N → I)T ( j ) and the latter between the nonlinearand linear IID solution branches at t = t (I)T ( j ). Last, the FID integration constants e , e , e , e , e , and e are specified by the initial condition G (cid:0) t = t (I → F)T ( j ) (cid:1) = G (I → F)T ( j ) and the integration constants e and e by the transition conditions G (cid:0) t = t (F)T , ( j ) | t (I → F)T ( j ) ≤ t < t (F)T , ( j ) (cid:1) = G (cid:0) t = t (F)T , ( j ) | t (F)T , ( j ) ≤ t < t j +1 (cid:1) and G (cid:0) t = t (F)T , ( j ) | t (I → F)T ( j ) ≤ t < t (F)T , ( j ) (cid:1) = G (cid:0) t = t (F)T , ( j ) | t (F)T , ( j ) ≤ t < t j +1 (cid:1) , respectively. These conditionsguarantee continuity between the IID and FID solution branches at t = t (I → F)T ( j ) and between the nonlinearand linear FID solution branches at t = t (F)T , ( j ) and t = t (F)T , ( j ). Further, the integration constants have to be8updated if (cid:8) t (N → I)T ( i ) , t (I → F)T ( i ) (cid:9) ∩ [ t j , t j +1 ) (cid:54) = ∅ for i ∈ { , ..., j − } , as elements of S switch over to S and/orelements of S switch over to S . Because these updates cause shifts in the NID, IID, and FID transition times t (N)T , t (I)T , t (F)T , , and t (F)T , towards larger values, we have to account for updating conditions that cover transitionsto other solution branches. Below, the determination of the NID integration constant c for the case S (cid:54) = ∅ (cid:54) = S is shown in detail. The remaining integration constants are computed accordingly. Applying the initial condition G (cid:0) t = t j | t j ≤ t < min (cid:0) t j +1 , t (N)T ( j ) (cid:1)(cid:1) = G j to the first branch of Sol. (31), we obtain c ( j ; S ) = G j − C ( j ; S ) t j . As c depends just on S , we have to consider only IID-FID updates. Hence, the conditions – and the updatedconstants – for all possible IID-FID transitions of the i th injection at t j < t = t (I → F)T ( i ) < min (cid:0) t j +1 , t (N)T ( j ) (cid:1) are: • First branch Sol. (31) → first branch Sol. (31) G (cid:0) t = t (I → F)T ( i ) | t j ≤ t < t (N)T ( j ; S \{ i } ) (cid:1) = G (cid:0) t = t (I → F)T ( i ) | t j ≤ t < t (N)T ( j ; S ∪ { i } ) (cid:1) c (cid:0) j ; S ∪ { i } (cid:1) = 3 (cid:16) C (cid:0) j ; S \{ i } (cid:1) − C (cid:0) j ; S ∪ { i } (cid:1)(cid:17) t (I → F)T ( i ) + c (cid:0) j ; S \{ i } (cid:1) • First branch Sol. (31) → Sol. (32) G (cid:0) t = t (I → F)T ( i ) | t j ≤ t < t (N)T ( j ; S \{ i } ) (cid:1) = G (cid:0) t = t (I → F)T ( i ) | t j ≤ t < t (N → I)T ( j ) (cid:1) c (cid:0) j ; S ∪ { i } (cid:1) = 3 (cid:16) C (cid:0) j ; S \{ i } (cid:1) − C (cid:0) j ; S ∪ { i } (cid:1)(cid:17) t (I → F)T ( i ) + c (cid:0) j ; S \{ i } (cid:1) . Since S ∪ { i } (cid:54) = ∅ and t (N)T ( j ; S \{ i } ) < t (N)T ( j ; S ∪ { i } ), 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 | ≤ 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 G SID ( t ) := H (cid:0) t (S)T (cid:1) H (cid:0) t − t (S)T (cid:1) × (cid:104) H (cid:0) t (S)T − t (cid:1) G (cid:0) t | ≤ t < t (S)T ; Sol . (57) (cid:1) + H (cid:0) t − t (S)T (cid:1) G (cid:0) t | t (S)T ≤ t < t ; Sol . (57) (cid:1)(cid:105) + H (cid:0) t (S)T − t (cid:1) G (cid:0) t | ≤ t < t ; Sol . (58) (cid:1) + H (cid:0) − t (S)T (cid:1) G (cid:0) t | ≤ t < t ; Sol . (59) (cid:1) . The more elaborate NID component is given by G NID ( t ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) G (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; Sol . (26); S = ∅ = S (cid:1) + χ ( S ) (cid:2) − χ ( S ) (cid:3) G (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; Sol . (26); S (cid:54) = ∅ , S = ∅ (cid:1) + (cid:2) − χ ( S ) (cid:3) χ ( S ) G NID (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; S = ∅ , S (cid:54) = ∅ (cid:1) + χ ( S ) χ ( S ) G NID (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; S (cid:54) = ∅ (cid:54) = S (cid:1) , G NID (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) = H (cid:0) t (N)T ( i ) − t i (cid:1) H (cid:0) t (N → I)T ( i ) − t (N)T ( i ) (cid:1) × (cid:104) H (cid:0) t (N)T ( i ) − t (cid:1) G (cid:0) t | t i ≤ t < t (N)T ( i ) ; Sol.(31) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t − t (N)T ( i ) (cid:1) G (cid:0) t | t (N)T ( i ) ≤ t < t (N → I)T ( i ) ; Sol.(31) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1)(cid:105) + H (cid:0) t (N)T ( i ) − t (N → I)T ( i ) (cid:1) G (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; Sol.(32) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t i − t (N)T ( i ) (cid:1) G (cid:0) t | t i ≤ t < t (N → I)T ( i ) ; Sol.(33) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) and χ ( S k ) := (cid:40) S k (cid:54) = ∅ S k = ∅ , with k ∈ { , , } , is the characteristic function. For the IID and FID components, we obtain similar expressions.On the one hand, we have G IID ( t ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) G (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; Sol . (38) ; S = ∅ = S (cid:1) + χ ( S ) (cid:2) − χ ( S ) (cid:3) G (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; Sol . (38) ; S (cid:54) = ∅ , S = ∅ (cid:1) + (cid:2) − χ ( S ) (cid:3) χ ( S ) G IID (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; S = ∅ , S (cid:54) = ∅ (cid:1) + χ ( S ) χ ( S ) G IID (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; S (cid:54) = ∅ (cid:54) = S (cid:1) , where G IID (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) = H (cid:0) t (I)T ( i ) − t (N → I)T ( i ) (cid:1) H (cid:0) t (I → F)T ( i ) − t (I)T ( i ) (cid:1) × (cid:104) H (cid:0) t (I)T ( i ) − t (cid:1) G (cid:0) t | t (N → I)T ( i ) ≤ t < t (I)T ( i ) ; Sol.(39) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t − t (I)T ( i ) (cid:1) G (cid:0) t | t (I)T ( i ) ≤ t < t (I → F)T ( i ) ; Sol.(39) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1)(cid:105) + H (cid:0) t (I)T ( i ) − t (I → F)T ( i ) (cid:1) G (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; Sol.(40) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t (N → I)T ( i ) − t (I)T ( i ) (cid:1) G (cid:0) t | t (N → I)T ( i ) ≤ t < t (I → F)T ( i ) ; Sol.(41) ; S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) , and on the other hand, we have G FID ( t ; i ) := (cid:2) − χ ( S ) (cid:3) (cid:2) − χ ( S ) (cid:3) G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S = ∅ = S (cid:1) + χ ( S ) (cid:2) − χ ( S ) (cid:3) G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S (cid:54) = ∅ , S = ∅ (cid:1) + (cid:2) − χ ( S ) (cid:3) χ ( S ) G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S = ∅ , S (cid:54) = ∅ (cid:1) + χ ( S ) χ ( S ) G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S (cid:54) = ∅ (cid:54) = S (cid:1) , G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S = ∅ = S (cid:1) = H (cid:0) t (F)T , ( i ) − t (I → F)T ( i ) (cid:1) H (cid:0) t i +1 − t (F)T , ( i ) (cid:1) × (cid:104) H (cid:0) t (F)T , ( i ) − t (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t (F)T , ( i ) ; Sol.(54) ; S = ∅ = S (cid:1) + H (cid:0) t − t (F)T , ( i ) (cid:1) G (cid:0) t | t (F)T , ( i ) ≤ t < t i +1 ; Sol.(54) ; S = ∅ = S (cid:1)(cid:105) + H (cid:0) t i +1 − t (I → F)T ( i ) (cid:1)(cid:104) H (cid:0) t (F)T , ( i ) − t i +1 (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; Sol.(55) ; S = ∅ = S (cid:1) + H (cid:0) t (I → F)T ( i ) − t (F)T , ( i ) (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; Sol.(56) ; S = ∅ = S (cid:1)(cid:105) as well as G FID (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) = H (cid:0) t (F)T , ( i ) − t (I → F)T ( i ) (cid:1) H (cid:0) t i +1 − t (F)T , ( i ) (cid:1) × (cid:104) H (cid:0) t (F)T , ( i ) − t (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t (F)T , ( i ) ; Sol.(46) ; S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t − t (F)T , ( i ) (cid:1) G (cid:0) t | t (F)T , ( i ) ≤ t < t i +1 ; Sol.(46) ; S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1)(cid:105) + H (cid:0) t i +1 − t (I → F)T ( i ) (cid:1)(cid:104) H (cid:0) t (F)T , ( i ) − t i +1 (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; Sol.(47) ; S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1) + H (cid:0) t (I → F)T ( i ) − t (F)T , ( i ) (cid:1) G (cid:0) t | t (I → F)T ( i ) ≤ t < t i +1 ; Sol.(48) ; S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ (cid:1)(cid:105) . Further, we determine the initial values G i for all i : 3 ≤ i ≤ m by requiring continuity between the solutionbranches of the ( i − i th injection domains. If the i th injection enters the system while the ( i − S (cid:54) = ∅ , S = ∅ G i = G (cid:0) t = t i | t i − ≤ t < t (N → I)T ( i − 1) ; Sol.(26) (cid:1) for t i < t (N → I)T ( i − S (cid:54) = ∅ , S (cid:54) = ∅ G i = G (cid:0) t = t i | t i − ≤ t < t (N)T ( i − 1) ; Sol.(31) (cid:1) for t i < t (N)T ( i − < t (N → I)T ( i − G (cid:0) t = t i | t (N)T ( i − ≤ t < t (N → I)T ( i − 1) ; Sol.(31) (cid:1) for t (N)T ( i − ≤ t i < t (N → I)T ( i − G (cid:0) t = t i | t i − ≤ t < t (N → I)T ( i − 1) ; Sol.(32) (cid:1) for t i < t (N → I)T ( i − ≤ t (N)T ( i − G (cid:0) t = t i | t i − ≤ t < t (N → I)T ( i − 1) ; Sol.(33) (cid:1) for t (N)T ( i − ≤ t i − < t i < t (N → I)T ( i − . If, however, the ( i − S (cid:54) = ∅ , S = ∅ in G i = G (cid:0) t = t i | t (N → I)T ( i − ≤ t < t (I → F)T ( i − 1) ; Sol.(38) (cid:1) for t (N → I)T ( i − ≤ t i < t (I → F)T ( i − , whereas for S (cid:54) = ∅ , S (cid:54) = ∅ , they yield1 G i = G (cid:0) t = t i | t (N → I)T ( i − ≤ t < t (I)T ( i − 1) ; Sol.(39) (cid:1) for t i < t (I)T ( i − < t (I → F)T ( i − G (cid:0) t = t i | t (I)T ( i − ≤ t < t (I → F)T ( i − 1) ; Sol.(39) (cid:1) for t (I)T ( i − ≤ t i < t (I → F)T ( i − G (cid:0) t = t i | t (N → I)T ( i − ≤ t < t (I → F)T ( i − 1) ; Sol.(40) (cid:1) for t i < t (I → F)T ( i − ≤ t (I)T ( i − G (cid:0) t = t i | t (N → I)T ( i − ≤ t < t (I → F)T ( i − 1) ; Sol.(41) (cid:1) for t (I)T ( i − ≤ t i − < t i < t (I → F)T ( i − . Finally, when the ( i − S (cid:54) = ∅ , S = ∅ or S (cid:54) = ∅ , S (cid:54) = ∅ G i = G (cid:0) t = t i | t (F)T , ( i − ≤ t < t i ; Sol.(46) (cid:1) for t (I → F)T ( i − < t (F)T , ( i − ≤ t i G (cid:0) t = t i | t (I → F)T ( i − ≤ t < t i ; Sol.(47) (cid:1) for t (I → F)T ( i − ≤ t i < t (F)T , ( i − G (cid:0) t = t i | t (I → F)T ( i − ≤ t < t i ; Sol.(48) (cid:1) for t (F)T , ( i − ≤ t (I → F)T ( i − ≤ t i and in case S = ∅ = S G i = G (cid:0) t = t i | t (F)T , ( i − ≤ t < t i ; Sol.(54) (cid:1) for t (I → F)T ( i − < t (F)T , ( i − ≤ t i G (cid:0) t = t i | t (I → F)T ( i − ≤ t < t i ; Sol.(55) (cid:1) for t (I → F)T ( i − ≤ t i < t (F)T , ( i − G (cid:0) t = t i | t (I → F)T ( i − ≤ t < t i ; Sol.(56) (cid:1) for t (F)T , ( i − ≤ t (I → F)T ( i − ≤ t i . Appendix F: The CS Function For z ∈ R ≥ , the CS function is defined by [15] CS ( z ) := 1 π (cid:90) π sin ( θ ) (cid:90) ∞ z/ sin ( θ ) K / ( y ) d y d θ = W , / ( z ) W , / ( z ) − W / , / ( z ) W − / , / ( z ) , (F1)where K a is the modified Bessel function and W a, b denotes the Whittaker function [1]. On account of the degreeof complexity of (F1), one usually employs an approximate function that is adapted to its asymptotics CS ( z ) (cid:39) (cid:40) a z − / for z (cid:28) z − exp ( − z ) for z (cid:29) , where a := 1 . 15. Standard approximations are, therefore, given by CS ( z ) := a exp ( − z ) z / , CS ( z ) := a exp ( − z ) z , and CS ( z ) := a z / (cid:0) z / exp ( z ) (cid:1) . Figure 2 shows the absolute values of the relative deviations of these approximations with respect to (F1) inpercent, that is, Dev( z ) := 100 (cid:12)(cid:12)(cid:12)(cid:12) CS ( z ) − CS n ( z ) CS ( z ) (cid:12)(cid:12)(cid:12)(cid:12) for n ∈ { , , } . Since the CS function is primarily used for the computation of the synchrotron intensity, where the spectrumcovers the energy range from radio waves to X-rays, i.e., with a lower energy limit of the order neV and anupper limit of the order keV, and z = 2 (cid:15)/ (3 (cid:15) γ ) with (cid:15) ∼ − b , initial electron energies γ i ∼ b − / , andnondimensional magnetic field strength b ∼ − up to b ∼ 10, we have to consider values z (cid:28) z (cid:29) CS ( z ), which coincides with both asymptotic endsof (F1) and has the smallest overall deviation.2 -2 -1 z D e v ( z ) [ % ] CS CS CS FIG. 2: Absolute values of the relative deviations of the approximate functions CS n for all n ∈ { , , } with respect tothe exact CS function in percent. Appendix G: Relativistic Beaming Accounting for relativistic beaming, the energy ε ∈ { (cid:15), (cid:15) s } and the time t , which are defined in the plasmoid restframe, are related to the corresponding quantities in an observer frame (denoted with an asterisk) by ε (cid:63) = D ε and t (cid:63) = t D , where D := 1Γ (cid:0) − β cos ( θ (cid:63) ) (cid:1) is the boost factor, Γ := 1 / (cid:112) − β the Lorentz factor of the plasmoid, β := v/c , and θ (cid:63) the angle between thejet axis and the line of sight of the observer. For blazars, one can assume that θ (cid:63) (cid:38) D (cid:39) (cid:115) β − β . Furthermore, because the ratio I/ε is Lorentz-invariant, i.e., I ( ε, t ) ε = I (cid:63) ( ε (cid:63) , t (cid:63) )( ε (cid:63) ) , one directly finds that the intensity and the fluence transform as I (cid:63) ( ε (cid:63) , t (cid:63) ) = D I ( ε, t ) and F (cid:63) ( ε (cid:63) ) = D F ( ε ) . Appendix H: Plotting Algorithm for the Fluence Spectral Energy Distributions The numerical implementation of G , the synchrotron and SSC intensities, as well as the corresponding totalfluence SEDs is carried out with Python. Here, we describe the functionality of the algorithm and the specificincorporation of the analytical formulas. The algorithm makes heavy use of the decimal package , which is The plotting algorithm is available on request via e-mail to [email protected]. For more detailed information on Python’s decimal package, we refer to https://docs.python.org/2/library/decimal.html. D and A , the injection strengths q i , and the reciprocalinitial electron energies x i spans several orders of magnitude, which may lead to a loss of accuracy in expressionswhere these parameters come up. This problem occurs in its most severe form in the evaluation of the formulasfor the transition times t (N)T , t (I)T , t (F)T , , and t (F)T , , yielding deviations from the expected values by more than 50%with standard floating point precision. But even with higher precision, it is preferable to avoid the evaluation ofthe transition times altogether. Thus, we compute G on a grid, using fixed time steps. This makes it possibleto read out the values of G at the grid points and directly compare them to the values of G (N)T , G (I)T , G (F)T , , and G (F)T , 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 injectionstrengths q i , and the reciprocal initial electron energies x i for i : 1 ≤ i ≤ m , the nondimensional magnetic fieldstrength b , the synchrotron and SSC cooling rate prefactors D and A , the Lorentz boost D of the plasmoid, andthe upper time boundary of the grid t end . The value of t end is chosen as one and a half times the injection timeof the final injection for a multiple-injection scenario or given by a sufficiently large value for a single-injectionscenario. In a realistic setting, however, t end corresponds to the end of the observation time. The time grid isset up homogeneously and linearly, and the number of grid points can be chosen arbitrarily. All computationsare performed in the plasmoid rest frame. For the later evaluation of the fluence SEDs, the relevant quantitiesare transformed into the observer frame (see Appendix G). Ordered lists of the initial and transition values G i , G (S)T = (cid:112) A q /D , G (N)T , G (I)T , G (F)T , , G (F)T , , G (N → I)T , and G (I → F)T , as well as a list keeping track of the elementsof the sets S , S , and S , and a list of the various solution branches are implemented. These lists are constantlyupdated during runtime. The algorithm constructs G incrementally in two loops. The first loop covers G fromthe time of the first injection t = 0 to the time of the second injection t = t or – in a scenario with only asingle 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 valueand G cur. := G ( t cur. ) is determined according to the proper solution branch, which is automatically selectedvia the above-mentioned lists. Moreover, at each grid point, the analytical expressions for G are glued togethercontinuously. In more detail, after initializing the values of C , C , C , C , as well as G (N → I)T and G (I → F)T at t = 0, the first loop starts its iteration, if 0 < G (S)T < G , in the first SID solution branch of (57). At eachstep, it evaluates G at the current time grid point and checks whether G cur. exceeds or equals G (S)T , i.e., whether G needs to be expressed by the second SID solution branch of (57). If necessary, G is altered accordingly andconnected to the previous solution branch continuously. (In case 0 < G ≤ G (S)T or G (S)T ≤ < G , the first loopmakes use of the solutions (58) or (59), respectively.) Also, if G cur. becomes larger than or equal to G (N → I)T (1)or G (I → F)T (1), 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 construction of G is completed. The second loop constructs G in a similar wayas the first loop, but now more cases, which arise from the more elaborate structure of the analytical multipleinjection solution, have to be taken into account. Once the second loop ends, the values of G are known at everypoint of the grid. This allows us to evaluate the synchrotron and SSC intensities at each grid point by simplysubstituting these values into the corresponding analytical formulas (64) and (68). The associated total fluencesare approximated by sums of the areas of rectangles, each of which is defined by the intensity at the left gridpoint and the size of the time step. Alternatively, one could implement the approximate analytical formulasderived in Section IV. The total synchrotron and SSC fluence SEDs are then computed as the product of therespective energy and fluence. Since the number of grid points can be increased arbitrarily, the precision of thesecomputations is limited only by the machine accuracy and the available CPU time. [1] Abramowitz M., Stegun I. A., “Handbook of mathematical functions with formulas, graphs, and mathematicaltables,” National Bureau of Standards (1972).[2] Aharonian F. A. et al., “An exceptional very high energy gamma-ray flare of PKS 2155-304,” The AstrophysicalJournal , L71 (2007).[3] Aharonian F. A. et al., “Simultaneous multiwavelength observations of the second exceptional γ -ray flare of PKS2155-304 in July 2006,” Astronomy & Astrophysics , 749 (2009). [4] Albert J. et al., “Variable very high energy gamma-ray emission from Markarian 501,” The Astrophysical Journal , 862 (2007).[5] Biteau J., Giebels B., “The minijets-in-a-jet statistical model and the rms-flux correlation,” Astronomy & Astro-physics , A123 (2012).[6] Blandford R. D., Payne D. G., “Hydromagnetic flows from accretion discs and the production of radio jets,” MonthlyNotices of the Royal Astronomical Society , 883 (1982).[7] Blandford R. D., Znajek R. L., “Electromagnetic extraction of energy from Kerr black holes,” Monthly Notices ofthe Royal Astronomical Society , 433 (1977).[8] Blazejowski M., Sikora M., Moderski R., Madejski G. M., “Comptonization of infrared radiation from hot dust byrelativistic jets in quasars,” The Astrophysical Journal , 107 (2000).[9] Blumenthal G. R., Gould R. J., “Bremsstrahlung, synchrotron radiation and Compton scattering of high-energyelectrons traversing dilute gases,” Reviews of Modern Physics , 237 (1970).[10] B¨ottcher M., “Modeling the emission processes in blazars,” Astrophysics and Space Science , 95 (2007).[11] B¨ottcher M., “Models for the spectral energy distributions and variability of blazars,” Proceedings “Fermi MeetsJansky - AGN at Radio and Gamma-Rays,” 41 (2010).[12] B¨ottcher M., Reimer A., Sweeney K., Prakash A., “Leptonic and hadronic modeling of Fermi-detected blazars,” TheAstrophysical Journal , 54 (2013).[13] Cerruti M., Zech A., Boisson C., Inoue S., “Lepto-hadronic modelling of blazar emission,” Proceedings of the Annualmeeting of the French Society of Astronomy and Astrophysics, 555 (2011).[14] Cerruti M., Zech A., Boisson C., Inoue S., “A hadronic origin for ultra-high-frequency-peaked BL Lac objects,”Monthly Notices of the Royal Astronomical Society , 910 (2015).[15] Crusius A., Schlickeiser R.,“Synchrotron radiation in a thermal plasma with large-scale random magnetic fields,”Astronomy & Astrophysics , 327 (1988).[16] Cui W., “X-ray flaring activity of Markarian 421,” The Astrophysical Journal , 662 (2004).[17] Dermer, C. D., Schlickeiser, R., “Model for the high-energy emission from blazars,” The Astrophysical Journal ,458 (1993).[18] Dermer, C. D., Sturner, S. J., Schlickeiser, R., “Nonthermal Compton and synchrotron processes in the jets of activegalactic nuclei,” The Astrophysical Journal Supplement Series , 103 (1997).[19] Eichmann B., Schlickeiser R., Rhode W., “On the duration of blazar synchrotron flares,” The Astrophysical Journal , 153 (2012).[20] Felten J. E., Morrison P., “Omnidirectional inverse Compton and synchrotron radiation from cosmic distributions offast electrons and thermal photons,” The Astrophysical Journal , 686 (1966).[21] Ghisellini G., Tavecchio F., Bodo G., Celotti A., “TeV variability in blazars: how fast can it be?,” Monthly Noticesof the Royal Astronomical Society , L16 (2009).[22] Giannios D., Uzdensky D. A., Begelman M. C., “Fast TeV variability in blazars: jets in a jet,” Monthly Notices ofthe Royal Astronomical Society , L29 (2009).[23] Graff P. B., Georganopoulos M., Perlman E. S., Kazanas D., “A multizone model for simulating the high-energyvariability of TeV blazars,” The Astrophysical Journal , 68 (2008).[24] Jones F. C., “Calculated spectrum of inverse-Compton-scattered photons,” Physical Review , 1159 (1968).[25] Kardashev N. S., “Nonstationariness of spectra of young sources of nonthermal radio emission,” Soviet Astronomy , 317 (1962).[26] Li H., Kusunose M., “Temporal and spectral variabilities of high-energy emission from blazars using synchrotronself-Compton models,” The Astrophysical Journal , 729 (2000).[27] Marscher A. P., “Turbulent, extreme multi-zone model for simulating flux and polarization variability in blazars,”The Astrophysical Journal , 87 (2014).[28] Pohl M., Schlickeiser R., “On the conversion of blast wave energy into radiation in active galactic nuclei and gamma-ray bursts,” Astronomy & Astrophysics , 395 (2000).[29] Reynolds S. P., “Theoretical studies of compact radio sources. II. Inverse-Compton radiation from anisotropic photonand electron distributions: General results and spectra from relativistic flows,” The Astrophysical Journal , 38(1982).[30] R¨oken C., Schlickeiser R., “Synchrotron self-Compton flaring of TeV blazars II. Linear and nonlinear electron cooling,”Astronomy & Astrophysics , 309 (2009).[31] R¨oken C., Schlickeiser R., “Linear and nonlinear radiative cooling of multiple instantaneously injected monoenergeticrelativistic particle populations in flaring blazars,” The Astrophysical Journal , 1 (2009).[32] Schlickeiser R., R¨oken C., “Synchrotron self-Compton flaring of TeV blazars I. Linear electron cooling,” Astronomy& Astrophysics , 701 (2008).[33] Schlickeiser R., “Non-linear synchrotron self-Compton cooling of relativistic electrons,” Monthly Notices of the RoyalAstronomical Society , 1483 (2009).[34] Schlickeiser R., Lerche I., “Nonlinear radiative cooling of relativistic particles under equipartition conditions I. In-stantaneous monoenergetic injection,” Astronomy & Astrophysics , 1 (2007).[35] Schlickeiser R., B¨ottcher M., Menzler U.,“Combined synchrotron and nonlinear synchrotron-self-Compton cooling of relativistic electrons,” Astronomy & Astrophysics , A9 (2010).[36] Sikora M., Begerlman M. C., Rees M. J ., “Comptonization of diffuse ambient radiation by a relativistic jet: Thesource of gamma rays from blazars?,” The Astrophysical Journal , 153 (1994).[37] Sikora M., Blazejowski M., Madejski G. M., Moderski R., “Modeling flares produced in blazars,” Proceedings of thefourth INTEGRAL workshop “Exploring the gamma-ray universe,” 259 (2001).[38] Sokolov A., Marscher A. P., McHardy I. M., “Synchrotron self-Compton model for rapid nonthermal flares in blazarswith frequency-dependent time lags,” The Astrophysical Journal , 725 (2004).[39] Urry, C. M., Padovani, P., “Unified schemes for radio-loud active galactic nuclei,” Publications of the AstronomicalSociety of the Pacific , 803 (1995).[40] Weidinger M., Spanier F., “A self-consistent and time-dependent hybrid blazar emission model,” Astronomy &Astrophysics , A7 (2015).[41] Yan D., Zhang L., “Understanding the TeV emission from a distant blazar PKS 1424 + 240 in a lepto-hadronic jetmodel,” Monthly Notices of the Royal Astronomical Society , 2810 (2015).[42] Zacharias M., Schlickeiser R., “Synchrotron lightcurves of blazars in a time-dependent synchrotron-self Comptoncooling scenario,” The Astrophysical Journal777