© ESO 2018

1 Introduction

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

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

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

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

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

2 The relativistic kinetic equation

The kinetic equation for the volume-averaged differential electron number density n = n(γ, t) (where t is the time, γ:= E_e∕(m_e c²) the normalized electron energy, and [n] = cm⁻³) of m ultrarelativistic, mono-energetic, instantaneously injected and spatially isotropically distributed electron populations in the rest frame of a nonthermal radiation source with dominant magnetic field self-generation and radiative loss rate L = L(γ, t) (with [L] =s⁻¹) reads (Kardashev 1962) $$\frac{\partial n}{\partial t}-\frac{\partial }{\partial \gamma }\left(Ln\right)={\sum }_{i=1}^m{q}_i\delta (\gamma -{\gamma }_i)\delta (t-t_i),$$(1)

where δ(⋅) is the Dirac distribution,q_i (for which [q_i] = cm⁻³) the ith injection strength, γ_i:= E_e,i∕(m_e c²) ≫ 1 the ith normalized initial electron energy, and t_i the ith injection time for i: 1 ≤ i ≤ m. In this work, radiative losses in form of both a linear, time-independent synchrotron cooling process (with a constant magnetic field B) and a nonlinear, time-dependent SSC cooling process in the Thomson limit $$L=L_{\text{syn}.}+L_{\text{SSC}}={\gamma }^2\left(D_0+A_0{\displaystyle {\int }_0^{\infty }{\gamma }^2}n\text{d}\gamma \right)$$(2)

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

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 ${({\gamma }_i,t_i)}_{i\in \{1,\mathrm{...},m\}}$. They are to be understood in the distributional sense, that is, each is a linear functional on the space of smooth test functions φ on $ℝ$ with compact support $$C_0^{\infty }(ℝ):=\left\{\phi |\phi \in C^{\infty }(ℝ),\text{supp}\phi \text{compact}\right\}.$$

More precisely, for $k\in ℝ$, one can rigorously define them as the mapping $$\delta :C_0^{\infty }(ℝ)\to ℝ,\phi \mapsto \phi (0)$$

with the integral of the Dirac distribution against a test function given by $${\displaystyle {\int }_{-\infty }^{\infty }\phi }(k)\delta (k)\text{d}k=\phi (0)\text{forall}\phi \in C_0^{\infty }(ℝ).$$

2.1 Formal solution of the relativistic kinetic equation

In terms of the function R(γ, t):= γ² n(γ, t) and the variable x:= 1∕γ, we can rewrite Eq. (1) in the form $$\frac{\partial R}{\partial t}+J\frac{\partial R}{\partial x}={\sum }_{i=1}^m{q}_i\delta (x-x_i)\delta (t-t_i),$$(3)

where $$J=J(t):=D_0+A_0{\displaystyle {\int }_0^{\infty }\frac{R(x,t)}{x^2}}\text{d}x.$$(4)

Defining the strictly increasing, continuous function G = G(t) via $$0<\frac{\text{d}G}{\text{d}t}:=J,$$(5)

Eq. (3) becomes $$\frac{\partial R}{\partial G}+\frac{\partial R}{\partial x}={\sum }_{i=1}^m{q}_i\delta (x-x_i)\delta (G-G_i)$$(6)

with G_i:= G(t_i). Although the nonlinear kinetic equation given in (1) is now transformed into a linear PDE, its solution can obviously serve as a Green’s function only for the single-injection scenario with m = 1. Consequently, in order to solve the generalized kinetic equation $$\frac{\partial n}{\partial t}-\frac{\partial }{\partial \gamma }\left(Ln\right)=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)={\sum }_{i=1}^m{q}_{i}H(G-G_i)\delta (x-x_i-G+G_i),$$(7)

where H(⋅) is the Heaviside step function $$H(k-k_0):=\{\begin{array}{ll}0\hfill & \text{for}k<k_0\hfill \\ 1\hfill & \text{for}k>k_0\hfill \end{array}$$

with $k,k_0\in ℝ$, which has a jump discontinuity at k = k₀. A detailed derivation of this solution can be found in Appendix A. In order to determine the function G, we substitute Sol. (7) into (4) and, by using (5), obtain the ordinary differential equation (ODE) $$\frac{\text{d}G}{\text{d}t}=D_0+A_0{\sum }_{i=1}^m\frac{q_{i}H(G-G_i)}{{\left(G-G_i+x_i\right)}^2}\cdot$$(8)

This equation can be regarded as a compact notation for the set of m piecewise-defined ODEs $$\begin{array}{l}\{\begin{array}{ll}\frac{\text{d}G}{\text{d}t}=D_0+A_0\frac{q_1}{{\left(G-G_1+x_1\right)}^2}\hfill & \text{for}{G}_1\le G<G_2\hfill \\ \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}=D_0+A_0\left(\frac{q_1}{{\left(G-G_1+x_1\right)}^2}+\frac{q_2}{{\left(G-G_2+x_2\right)}^2}\right)\hfill & \text{for}{G}_2\le G<G_3\hfill \\ \hfill & \hfill \\ ⋮⋮⋮\hfill & ⋮\hfill \\ \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}=D_0+A_0\left(\frac{q_1}{{\left(G-G_1+x_1\right)}^2}+\frac{q_2}{{\left(G-G_2+x_2\right)}^2}+\dots +\frac{q_m}{{\left(G-G_m+x_m\right)}^2}\right)\hfill & \text{for}{G}_m\le G<\infty .\hfill \end{array}\hfill \end{array}$$(9)

In Sect. 2.2, we present an approximate analytical solution for the general case of j injections, that is, for the interval G_j ≤ G < G_j+1, where j ∈{1,..., m} and G₁ = 0, G_m+1 = ∞, employing the method of matched asymptotic expansions. Having a separate analytical solution for each ODE of (9), in Sect. 2.3, these are connected successively requiring continuity at the transition points. Finally, by substituting (7), the electron number density results in $$n(\gamma ,t)={\gamma }^{-2}R(\gamma ,t)={\gamma }^{-2}{\sum }_{i=1}^m{q}_{i}H(G(t)-G_i)\delta \left(\frac{1}{\gamma }-\frac{1}{{\gamma }_i}-G(t)+G_i\right)={\sum }_{i=1}^m{q}_{i}H\left(t-t_i\right)\delta \left(\gamma -\frac{{\gamma }_i}{1+{\gamma }_i(G(t)-G_i)}\right)\cdot$$(10)

2.2 Solution of the relativistic kinetic equation for { t $\in {ℝ}_{\ge 0}|$ t_j ≤t < t_j+1 with j ∈{1,..., m}}

In the interval G_j ≤ G < G_j+1, Eq. (8) gives $$\frac{\text{d}G}{\text{d}t}=D_0+A_0{\sum }_{i=1}^j\frac{q_i}{{\left(G-G_i+x_i\right)}^2}\cdot$$(11)

In order to solve this equation, we introduce the time-dependent sets S₁, S₂, and S₃ that contain indices i: 1 ≤ i ≤ j belonging to either injections in the near-injection domain (NID) 0 ≤ G − G_i ≪ x_i, the intermediate-injection domain (IID) G − G_i ≈ x_i, or the far-injection domain (FID) G − G_i ≫ x_i, respectively. We may now rewrite Eq. (11) in terms of these sets by continuing their domains of validity to $G_i\le G<\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{N}\to \text{I})}(i)\right)$ for the NID, $G_{\text{T}}^{(\text{N}\to \text{I})}(i)\le G<\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{I}\to \text{F})}(i)\right)$ for the IID, and $G_{\text{T}}^{(\text{I}\to \text{F})}(i)\le G<G_{j+1}$ for the FID with the NID-IID and IID-FID transition values $G_{\text{T}}^{(\text{N}\to \text{I})}(i):=G_i+x_i/{\xi }_i^{(\text{N}\to \text{I})}$ and $G_{\text{T}}^{(\text{I}\to \text{F})}(i):=G_i+{\xi }_i^{(\text{I}\to \text{F})}{x}_i$, where ${\xi }_i^{(\text{N}\to \text{I})},{\xi }_i^{(\text{I}\to \text{F})}>1$ are positive constants, taking into account that the (j + 1)th injection can occur in each domain. However, for the computations below, we still use the former much less than, approximately equal, and much greater than relations when we carry out approximations. Then, Eq. (11) can be represented formally by $$\frac{\text{d}G}{\text{d}t}=D_0+A_0\left({\sum }_{u\in S_{\text{}1}}\frac{q_u}{{\left(G-G_u+x_u\right)}^2}+{\sum }_{v\in S_{\text{}2}}\frac{q_v}{{\left(G-G_v+x_v\right)}^2}+{\sum }_{w\in S_{\text{}3}}\frac{q_w}{{\left(G-G_w+x_w\right)}^2}\right)\cdot$$(12)

Next, we express the NID and the FID summands by the Taylor series of the function ${(1+z)}^{-2}$ evaluated at the point z = 0 (generalized geometric series) $$\frac{1}{{(1+z)}^2}={\sum }_{n=0}^{\infty }{(-1)}^n(n+1)z^n\text{for}\left|z\right|<1,$$

and theIID summands by the Taylor series of the same function evaluated at the point z = 1 $$\frac{1}{{(1+z)}^2}={\sum }_{n=0}^{\infty }{(-1)}^n\frac{n+1}{2^{n+2}}{(z-1)}^n\text{for}|z-1|<2,$$

and suitably approximate these series by considering only the leading- and next-to-leading-order terms. This results in $$\begin{array}{ll}\frac{1}{{(G-G_u+x_u)}^2}\hfill & =\frac{1}{x_u^2}{\sum }_{n=0}^{\infty }(n+1){\left(-\frac{G-G_u}{x_u}\right)}^n\approx \frac{1}{x_u^2}\left(1-\frac{2(G-G_u)}{x_u}\right)\hfill \\ \hfill & \hfill \\ \frac{1}{{(G-G_v+x_v)}^2}\hfill & =\frac{1}{x_v^2}{\sum }_{n=0}^{\infty }\frac{n+1}{2^{n+2}}{\left(1-\frac{G-G_v}{x_v}\right)}^n\approx \frac{1}{4x_v^2}\left(2-\frac{G-G_v}{x_v}\right)\hfill \\ \hfill & \hfill \\ \frac{1}{{(G-G_w+x_w)}^2}\hfill & =\frac{1}{{(G-G_w)}^2}{\sum }_{n=0}^{\infty }(n+1){\left(-\frac{x_w}{G-G_w}\right)}^n\approx \frac{1}{{(G-G_w)}^2}\left(1-\frac{2x_w}{G-G_w}\right)\hfill \end{array}$$

for all u ∈ S₁, v ∈ S₂, and w ∈ S₃. To find an approximate analytical solution to Eq. (12), we require the function G to be extricable from the respective sums. This is already achieved with the approximations (13) and (14). The latter approximation (15), however, needs further modifications as G still couples to G_w in the denominator. To this end, we employ both geometric and generalized geometric series and again consider only the leading- and next-to-leading-order terms $$\frac{1}{{(G-G_w+x_w)}^2}\approx \frac{1}{G^2}{\sum }_{k=0}^{\infty }(k+1){\left(\frac{G_w}{G}\right)}^k\left[1-\frac{2x_w}{G}{\sum }_{l=0}^{\infty }{\left(\frac{G_w}{G}\right)}^l\right]\approx \frac{1}{G^2}\left(1+\frac{2(G_w-x_w)}{G}\right)\cdot$$(16)

The approximation errors and optimal values for the interval parameters ${\xi }_i^{(\text{N}\to \text{I})}$ and ${\xi }_i^{(\text{I}\to \text{F})}$ are given in Appendix B. Substituting (13), (14), and (16) into the ODE (12), we obtain $$\frac{\text{d}G}{\text{d}t}\approx {\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})+{\mathcal{C}}_2(j;S_{\text{}1},S_{\text{}2})G+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}+\frac{{\mathcal{C}}_4(j;S_{\text{}3})}{G^3}$$(17)

with the “constants” $$\begin{array}{ll}\hfill & {\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2}):=D_0+A_0\left({\sum }_{u\in S_{\text{}1}}\frac{q_u}{x_u^2}\left[1+\frac{2G_u}{x_u}\right]+\frac{1}{2}{\sum }_{v\in S_{\text{}2}}\frac{q_v}{x_v^2}\left[1+\frac{G_v}{2x_v}\right]\right),\hfill \\ \hfill & \hfill \\ \hfill & {\mathcal{C}}_2(j;S_{\text{}1},S_{\text{}2}):=-A_0\left(2{\sum }_{u\in S_{\text{}1}}\frac{q_u}{x_u^3}+\frac{1}{4}{\sum }_{v\in S_{\text{}2}}\frac{q_v}{x_v^3}\right),\hfill \\ \hfill & \hfill \\ \hfill & {\mathcal{C}}_3(j;S_{\text{}3}):=A_0{\sum }_{w\in S_{\text{}3}}{q}_w,\text{and}{\mathcal{C}}_4(j;S_{\text{}3}):=2A_0{\sum }_{w\in S_{\text{}3}}{q}_w(G_w-x_w),\hfill \end{array}$$(18)

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

2.2.1 Near-injection domain solution

In the NID $G_j\le G<\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{N}\to \text{I})}(j)\right)$, at least the jth injection is an element of S₁. Further, one may – but does not necessarily need to – have injections in S₂ and S₃. Therefore, we distinguish the following four cases $$\begin{array}{llll}\frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1})+{\mathcal{C}}_2(j;S_{\text{}1})G\hfill & \hfill & \text{for}{S}_{\text{}2}=\varnothing =S_{\text{}3}\hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})+{\mathcal{C}}_2(j;S_{\text{}1},S_{\text{}2})G\hfill & \hfill & \text{for}{S}_{\text{}2}\cancel{=}\varnothing ,S_{\text{}3}=\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}2}=\varnothing ,S_{\text{}3}\cancel{=}\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}2}\cancel{=}\varnothing \cancel{=}{S}_{\text{}3}.\hfill \end{array}$$

Since Eqs. (19) and (20) as well as Eqs. (21) and (22) are identical except for the present constants, only two different kinds of ODEs have to be solved. Starting with Eqs. (19) and (20), we use separation of variables in order to obtain $${\displaystyle \int \frac{\text{d}G}{{\mathcal{C}}_1+{\mathcal{C}}_{2}G}}=\frac{1}{{\mathcal{C}}_2}\ln ({\mathcal{C}}_1+{\mathcal{C}}_{2}G)=t+c_1,c_1=\text{const}.\in ℝ,$$

which isequivalent to $$G(t)=\frac{1}{|{\mathcal{C}}_2(j)|}[{\mathcal{C}}_1(j)-\exp (-|{\mathcal{C}}_2(j)|(t+c_1))]\text{for}{t}_j\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N}\to \text{I})}(j)\right).$$

As this solution converges strictly against the value ${\mathcal{C}}_1(j)/|{\mathcal{C}}_2(j)|$, which can be lower than the NID transition value $G_{\text{T}}^{(\text{N}\to \text{I})}(j)$, it is not suitable for covering this domain. Including a second-order term in the ODE under consideration, that is, working with the equation $$\frac{\text{d}G}{\text{d}t}={\mathcal{C}}_1+{\mathcal{C}}_{2}G+\mathcal{C}{G}^2,$$(23)

where $\mathcal{C}=\text{const}.\in ℝ$, leads to a solution in form of a tangent function, which may have several poles in the NID, rendering it useless as well. Adding further contributions to (23) yields analytically non-solvable ODEs. Hence, we have to employ the simpler ODEs $$\begin{array}{llll}\frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1})\hfill & \hfill & \text{for}{S}_{\text{}2}=\varnothing =S_{\text{}3}\hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})\hfill & \hfill & \text{for}{S}_{\text{}2}\cancel{=}\varnothing ,S_{\text{}3}=\varnothing \hfill \end{array}$$

with the linear solution $$G(t)={\mathcal{C}}_1(j)t+c_2(j)\text{for}{t}_j\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N}\to \text{I})}(j)\right)$$(26)

and $c_2=\text{const}.\in ℝ$. For Eqs. (21) and (22), we also apply separation of variables and extend the integrand with the multiplicative identity ${\mathcal{C}}_1/{\mathcal{C}}_1$ and the neutral element ${\mathcal{C}}_3-{\mathcal{C}}_3$, leading to $${\displaystyle \int \frac{G^2\text{d}G}{{\mathcal{C}}_3+{\mathcal{C}}_1{G}^2}}=\frac{1}{{\mathcal{C}}_1}\left(G-{\mathcal{C}}_3{\displaystyle \int \frac{\text{d}G}{{\mathcal{C}}_3+{\mathcal{C}}_1{G}^2}}\right)=\frac{1}{{\mathcal{C}}_1}\left[G-\sqrt{\frac{{\mathcal{C}}_3}{{\mathcal{C}}_1}}\text{arctan}\left(\sqrt{\frac{{\mathcal{C}}_1}{{\mathcal{C}}_3}}G\right)\right]=t+c_3,$$(27)

where $c_3=\text{const}.\in ℝ$. One way of deriving an approximate analytical solution of this transcendental equation is to determine G asymptotically for small and large arguments of the arctan function $($in case both asymptotic ends exist for $G:G_j\le G<\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{N}\to \text{I})}(j)\right))$ and extrapolate these solutions up to an intermediate transition point, where they are connected requiring continuity, such that the entire domain of definition is covered. For small arguments $\sqrt{{\mathcal{C}}_1/{\mathcal{C}}_3}G\ll 1$, we use the third-order approximation $\text{arctan}{(z)}_{|z\ll 1}\approx z-z^3/3$, as the linear term in Eq. (27) and the linear contribution in the approximation of the arctan function cancel each other out. This leads to the asymptotic solution $$G(t)\simeq {(3{\mathcal{C}}_{3}t+{\mathfrak{c}}_0)}^{1/3},{\mathfrak{c}}_0=\text{const}.\in ℝ.$$(28)

For large arguments $\sqrt{{\mathcal{C}}_1/{\mathcal{C}}_3}G\gg 1$, the arctan function can be well approximated by π∕2. We directly obtain an asymptotic solution of the form $$G(t)\simeq {\mathcal{C}}_{1}t+{{\mathfrak{c}}^{\prime }}_0,{{\mathfrak{c}}^{\prime }}_0=\text{const}.\in ℝ.$$(29)

By means of the condition $\sqrt{{\mathcal{C}}_1/{\mathcal{C}}_3}G=1$, we derive the NID transition value $$G_{\text{T}}^{(\text{N})}(j)=\sqrt{\frac{{\mathcal{C}}_3(j)}{{\mathcal{C}}_1(j)}}\cdot$$(30)

This value specifies the upper bound of the domain of validity of (28) and the lower bound of the domain of validity of (29). The corresponding transition time can be directly computed by substituting (30) into (27), in which the constant c₃ had to be fixed via the initial condition G(t = t_j) = G_j resulting in $$c_3=\frac{1}{{\mathcal{C}}_1}\left[G_j-\sqrt{\frac{{\mathcal{C}}_3}{{\mathcal{C}}_1}}\text{arctan}\left(\sqrt{\frac{{\mathcal{C}}_1}{{\mathcal{C}}_3}}{G}_j\right)\right]-t_j.$$

This yields $$t_{\text{T}}^{(\text{N})}(j)=t_j+\frac{1}{{\mathcal{C}}_1(j)}\left(\sqrt{\frac{{\mathcal{C}}_3(j)}{{\mathcal{C}}_1(j)}}\left[1-\frac{\pi }{4}+\text{arctan}\left(\sqrt{\frac{{\mathcal{C}}_1(j)}{{\mathcal{C}}_3(j)}}{G}_j\right)\right]-G_j\right).$$

We point out that in general, one does not have the ordered sequence $G_j<G_{\text{T}}^{(\text{N})}(j)<G_{\text{T}}^{(\text{N}\to \text{I})}(j)$ because $G_{\text{T}}^{(\text{N})}(j)$ can in principle be larger than or equal to $G_{\text{T}}^{(\text{N}\to \text{I})}(j)$ as well as smaller thanor equal to G_j. These cases arise when only one asymptotic end of the arctan function exists. Hence, for $G_j<G_{\text{T}}^{(\text{N})}(j)<G_{\text{T}}^{(\text{N}\to \text{I})}(j)$, the solution of Eq. (27) is approximately given by $$G(t)=\{\begin{array}{ll}{(3{\mathcal{C}}_3(j)t+c_4(j))}^{1/3}\hfill & \text{for}{t}_j\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N})}(j)\right)\hfill \\ {\mathcal{C}}_1(j)t+c_5(j)\hfill & \text{for}\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N})}(j)\right)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N}\to \text{I})}(j)\right),c_4,c_5=\text{const}.\in ℝ,\hfill \end{array}$$(31)

for $G_j<G_{\text{T}}^{(\text{N}\to \text{I})}(j)\le G_{\text{T}}^{(\text{N})}(j)$ by $$G(t)={(3{\mathcal{C}}_3(j)t+c_6(j))}^{1/3}\text{for}{t}_j\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N}\to \text{I})}(j)\right),c_6=\text{const}.\in ℝ,$$(32)

and for $G_{\text{T}}^{(\text{N})}(j)\le G_j$ by $$G(t)={\mathcal{C}}_1(j)t+c_7(j)\text{for}{t}_j\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{N}\to \text{I})}(j)\right),c_7=\text{const}.\in ℝ.$$(33)

The integration constants c₂ and c₄,..., c₇ are determined via the proper initial and transition conditions in Appendix D. In order to compute the constants ${\mathcal{C}}_1(j;S_{\text{}1})$, ${\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})$, ${\mathcal{C}}_2(j;S_{\text{}1})$, ${\mathcal{C}}_2(j;S_{\text{}1},S_{\text{}2})$, and ${\mathcal{C}}_3(j;S_{\text{}3})$, we have to continually check duringthe evolution of $G\in \left[G_j,\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{N}\to \text{I})}(j)\right)\right)$ whether an injection i: 1 ≤ i ≤ j belongs to S₁, S₂, or S₃. Therefore, we specify two different kinds of updates. The first kind occurs when a new injection enters the system, while the second kind is due to either NID-IID or IID-FID transitions. We start with the initial update at the time of the jth injection, verifying $$\begin{array}{ll}\hfill & G_j-G_i<x_i/{\xi }_i^{(\text{N}\to \text{I})}\text{for the}i\text{th injection being in}{S}_{\text{}1}\hfill \\ \hfill & x_i/{\xi }_i^{(\text{N}\to \text{I})}\le G_j-G_i<{\xi }_i^{(\text{I}\to \text{F})}{x}_i\text{for the}i\text{th injection being in}{S}_{\text{}2}\hfill \\ \hfill & {\xi }_i^{(\text{I}\to \text{F})}{x}_i\le G_j-G_i\text{for the}i\text{th injection being in}{S}_{\text{}3}.\hfill \end{array}$$

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

2.2.2 Intermediate-injection domain solution

In the IID $G_{\text{T}}^{(\text{N}\to \text{I})}(j)\le G<\text{min}\left(G_{j+1},G_{\text{T}}^{(\text{I}\to \text{F})}(j)\right)$, where at least the jth injection is in S₂, the previousj − 1 injections can reside in their respective NIDs, IIDs, or FIDs. Thus, we have to discuss the four cases $$\begin{array}{llll}\frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}2})\hfill & \hfill & \text{for}{S}_{\text{}1}=\varnothing =S_{\text{}3}\hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})\hfill & \hfill & \text{for}{S}_{\text{}1}\cancel{=}\varnothing ,S_{\text{}3}=\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}2})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}1}=\varnothing ,S_{\text{}3}\cancel{=}\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}1}\cancel{=}\varnothing \cancel{=}{S}_{\text{}3}.\hfill \end{array}$$

As these coincide structurally with the ODEs (24) and (25) as well as (21) and (22) of the NID, we can directly write down their solutions. For Eqs. (34) and (35), we get $$G(t)={\mathcal{C}}_1(j)t+d_1(j)\text{for}{t}_{\text{T}}^{(\text{N}\to \text{I})}(j)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I}\to \text{F})}(j)\right),d_1=\text{const}.\in ℝ,$$(38)

whereas for Eqs. (36) and (37), we obtain in case $G_{\text{T}}^{(\text{N}\to \text{I})}(j)<G_{\text{T}}^{(\text{I})}(j)<G_{\text{T}}^{(\text{I}\to \text{F})}(j)$ $$G(t)=\{\begin{array}{ll}{(3{\mathcal{C}}_3(j)t+d_2(j))}^{1/3}\hfill & \text{for}{t}_{\text{T}}^{(\text{N}\to \text{I})}(j)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I})}(j)\right)\hfill \\ {\mathcal{C}}_1(j)t+d_3(j)\hfill & \text{for}\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I})}(j)\right)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I}\to \text{F})}(j)\right),d_2,d_3=\text{const}.\in ℝ,\hfill \end{array}$$(39)

if $G_{\text{T}}^{(\text{N}\to \text{I})}(j)<G_{\text{T}}^{(\text{I}\to \text{F})}(j)\le G_{\text{T}}^{(\text{I})}(j)$ $$G(t)={(3{\mathcal{C}}_3(j)t+d_4(j))}^{1/3}\text{for}{t}_{\text{T}}^{(\text{N}\to \text{I})}(j)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I}\to \text{F})}(j)\right),d_5=\text{const}.\in ℝ,$$(40)

and if $G_{\text{T}}^{(\text{I})}(j)\le G_{\text{T}}^{(\text{N}\to \text{I})}(j)$ $$G(t)={\mathcal{C}}_1(j)t+d_5(j)\text{for}{t}_{\text{T}}^{(\text{N}\to \text{I})}(j)\le t<\text{min}\left(t_{j+1},t_{\text{T}}^{(\text{I}\to \text{F})}(j)\right),d_5=\text{const}.\in ℝ,$$(41)

where $$t_{\text{T}}^{(\text{I})}(j)=t_{\text{T}}^{(\text{N}\to \text{I})}(j)+\frac{1}{{\mathcal{C}}_1(j)}\left(\sqrt{\frac{{\mathcal{C}}_3(j)}{{\mathcal{C}}_1(j)}}\left[1-\frac{\pi }{4}+\text{arctan}\left(\sqrt{\frac{{\mathcal{C}}_1(j)}{{\mathcal{C}}_3(j)}}{G}_{\text{T}}^{(\text{N}\to \text{I})}(j)\right)\right]-G_{\text{T}}^{(\text{N}\to \text{I})}(j)\right)$$

is the IID transition time associated with the transition value $$G_{\text{T}}^{(\text{I})}(j)=\sqrt{\frac{{\mathcal{C}}_3(j)}{{\mathcal{C}}_1(j)}}\cdot$$

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

2.2.3 Far-injection domain solution

In the FID $G_{\text{T}}^{(\text{I}\to \text{F})}(j)\le G<G_{j+1}$, the jth injection is, per definition, an element of S₃. As before, the previous j − 1 injections can be in their respective NIDs, IIDs, or FIDs depending on the initial injection parameters, and therefore we have to evaluate thefour ODEs $$\begin{array}{llll}\frac{\text{d}G}{\text{d}t}\hfill & =D_0+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}+\frac{{\mathcal{C}}_4(j;S_{\text{}3})}{G^3}\hfill & \hfill & \text{for}{S}_{\text{}1}=\varnothing =S_{\text{}2}\hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}1}\cancel{=}\varnothing ,S_{\text{}2}=\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}2})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}1}=\varnothing ,S_{\text{}2}\cancel{=}\varnothing \hfill \\ \hfill & \hfill & \hfill & \hfill \\ \frac{\text{d}G}{\text{d}t}\hfill & ={\mathcal{C}}_1(j;S_{\text{}1},S_{\text{}2})+\frac{{\mathcal{C}}_3(j;S_{\text{}3})}{G^2}\hfill & \hfill & \text{for}{S}_{\text{}1}\cancel{=}\varnothing \cancel{=}{S}_{\text{}2}.\hfill \end{array}$$

Since Eqs. (43)–(45) are also structurally identical to the ODEs (21) and (22), we can once again directly write down their solutions, yielding for $G_{\text{T}}^{(\text{I}\to \text{F})}(j)<G_{\text{T},1}^{(\text{F})}(j)<G_{j+1}$ $$G(t)=\{\begin{array}{ll}{(3{\mathcal{C}}_3(j)t+e_1(j))}^{1/3}\hfill & \text{for}{t}_{\text{T}}^{(\text{I}\to \text{F})}(j)\le t<t_{\text{T},1}^{(\text{F})}(j)\hfill \\ {\mathcal{C}}_1(j)t+e_2(j)\hfill & \text{for}{t}_{\text{T},1}^{(\text{F})}(j)\le t<t_{j+1},e_1,e_2=\text{const}.\in ℝ,\hfill \end{array}$$(46)

for $G_{\text{T}}^{(\text{I}\to \text{F})}(j)<G_{j+1}\le G_{\text{T},1}^{(\text{F})}(j)$ $$G(t)={(3{\mathcal{C}}_3(j)t+e_3(j))}^{1/3}\text{for}{t}_{\text{T}}^{(\text{I}\to \text{F})}(j)\le t<t_{j+1},e_3=\text{const}.\in ℝ,$$(47)

and for $G_{\text{T},1}^{(\text{F})}(j)\le G_{\text{T}}^{(\text{I}\to \text{F})}(j)<G_{j+1}$ $$G(t)={\mathcal{C}}_1(j)t+e_4(j)\text{for}{t}_{\text{T}}^{(\text{I}\to \text{F})}(j)\le t<t_{j+1},e_4=\text{const}.\in ℝ,$$(48)

where $$t_{\text{T},1}^{(\text{F})}(j)=t_{\text{T}}^{(\text{I}\to \text{F})}(j)+\frac{1}{{\mathcal{C}}_1(j)}\left(\sqrt{\frac{{\mathcal{C}}_3(j)}{{\mathcal{C}}_1(j)}}\left[1-\frac{\pi }{4}+\text{arctan}\left(\sqrt{\frac{{\mathcal{C}}_1(j)}{{\mathcal{C}}_3(j)}}{G}_{\text{T}}^{(\text{I}\to \text{F})}(j)\right)\right]-G_{\text{T}}^{(\text{I}\to \text{F})}(j)\right)$$