© The Authors 2023

1. Introduction

Blazar emission (with blazars defined as active galaxies whose jets point toward the Earth) is typically explained on the basis of the so-called one-zone model, where a single zone is responsible for most of the source’s radiative output. The spectral energy distribution (SED) is characterized by two broad humps. The low-energy one is explained by electron-synchrotron emission, while the high-energy hump can be explained by inverse-Compton emission or hadronically induced processes, such as proton synchrotron or synchrotron emission from the leptonic cascade. An important role is played by photon fields external to the jet, such as the accretion disk (AD), the broad-line region (BLR), and the dusty torus (DT), as these fields may provide ample seed photons for particle-photon and photon-photon interactions (see, e.g., Böttcher 2019; Cerruti 2020, for detailed reviews). The one-zone model is well justified in flares, where the variability timescale restricts the size of the emission region. However, while the particle flow in this region is relativistic, the emission region itself is typically assumed to remain stationary with respect to the black hole (e.g., H.E.S.S. Collaboration 2019), which allows us to ignore certain complications, such as varying (external) photon fields, adiabatic expansion and cooling, and so on. On the other hand, the observation of stationary features in radio very-long-baseline interferometry (VLBI) observations (and other wavelengths where jets have been resolved) suggests standing recollimation shocks (e.g., Weaver et al. 2022), where particles may be accelerated and emit radiation. In such a situation, the emission region would indeed be stationary with respect to the black hole.

The same radio VLBI observations have further revealed moving components launched somewhere upstream with respect to the radio core. The interaction of such moving features with the standing features has been connected with multiwavelength flaring events (e.g., Ahnen et al. 2017; H.E.S.S. Collaboration 2021). Indeed, numerical simulations by Fichet de Clairfontaine et al. (2021, 2022) have shown that the interaction of moving and standing shocks can induce rapid flaring events. Thus, it is important to study the flaring characteristics of a moving, expanding emission region (or “blob”).

Recently, Boula & Mastichiadis (2022) and Tramacere et al. (2022) discussed this model in detail and derived the expected time delays between the γ-ray and the radio band expected from the progressive optical thinning at radio frequencies of the expanding source (see also Saito et al. 2015). In the present work, we specifically consider the effect of γ − γ pair production on the evolution of the particle distribution and the photon fluxes. In particular, bright external photon fields can have a major influence, as they provide a significant amount of absorption for the γ rays of the emission region (e.g., Zacharias 2021). However, as these photon fields are only present in the direct vicinity of the black hole, the optical thickness for γ rays changes over time, as the blob moves down the jet. Additionally, the opening angle will strongly influence the escape time of both photons and particles owing to the increase in the radius of the emission region. This could potentially prolong the pair production as the (internal) photons have more time to interact before escape. This naturally competes with the thinning of the photon density due to expansion. In any case, it is possible to obtain variability simply via the motion of an expanding blob through the various radiation fields due to the varying injection of pairs.

Firstly, we derive a simple semi-analytical model that describes the time-dependent evolution of the particle distribution, given in Sect. 2. In this section, we also introduce the basic assumptions about the escape timescale and the implications of the expansion. We then derive the particle distribution without pair injection, followed by a very simple linear cascade model for a few exemplary cases. In Sect. 3, we use the numerical code OneHaLe (Zacharias 2021; Zacharias et al. 2022) to derive the light curves of a realistic simulation based on the parameters of PKS 1510−089, which is known for its bright external photon fields. We present our conclusions in Sect. 4.

2. A simple electron evolution model

We derived simple semi-analytical models for the time-dependent evolution of the electron distribution and while ignoring all energy dependencies. This is indeed a major simplification, however, it allows us to study three distinct cases with a linear cascade (at most). Using the time-dependent, one-zone radiation code OneHaLe, we reproduced two of the three cases well. We discuss the failure of the third case, as well as the influence of the choice of parameters on the solutions.

2.1. Escape time

The blob travels with constant speed, β_Γc, corresponding to the bulk Lorentz factor $\mathrm{\Gamma }={(1-{\beta }_{\mathrm{\Gamma }}^2)}^{-1/2}$ along the z-axis. Within a conical jet, the opening angle is constant. The VLBI observations suggest an opening angle of ∝α/Γ with α ∼ 0.26 (Pushkarev et al. 2017), but we used α as a free parameter. The radius, R, of the blob thus evolves as a function of comoving time¹, t, and jet coordinate, z′:

$$\begin{array}{cc}\hfill R(t)& =R_0+(z\prime (t)-z_0^{\prime })tan(\alpha /\mathrm{\Gamma })\hfill \\ \hfill & =R_0+\mathrm{\Gamma }{\beta }_{\mathrm{\Gamma }}cttan(\alpha /\mathrm{\Gamma })\hfill \\ \hfill & \approx R_0+\alpha ct,\hfill \end{array}$$(1)

where, in the last line, we present approximations for small angles and β_Γ ≈ 1. Boula & Mastichiadis (2022) expressed this equation by an “expansion speed”, u_exp, which relates to our equation as u_exp = αc. In their paper, they used α between 0.01 and 0.2. It should be noted that one of the common estimates for jet expansion, α = 1 (e.g., Ackermann et al. 2016), implies a radial expansion with the speed of light, while α > 1 implies a superluminal expansion and thus a causal disconnection of regions within the blob or jet.

On average, particles escape the emission region on an energy-independent timescale,

$$\begin{array}{c}\hfill t_{\mathrm{esc}}(t)={\eta }_{\mathrm{esc}}R(t)/c,\end{array}$$(2)

with η_esc > 1. This resembles an advective motion of the plasma below the speed of light mimicking the trapping of charged particles in the magnetized blob. Calculating the average escape timescale for photons, we obtain η_esc, ph = 3/4 (Böttcher & Chiang 2002).

In an expanding blob, the escape time increases. In turn, particles and photons take longer to escape and can only efficiently do so once the intrinsic time, t, since launch surpasses ∼t_esc(t); namely,

$$\begin{array}{cc}\hfill t& >\frac{{\eta }_{\mathrm{esc}}}{c}(R_0+\alpha ct)\hfill \\ \hfill \iff t& >\frac{t_{\mathrm{esc}}(0)}{1-{\eta }_{\mathrm{esc}}\alpha },\hfill \end{array}$$(3)

with t_esc(0) = η_escR₀/c. In the case of η_escα → 1, particles are effectively trapped in the blob without a chance for a meaningful escape. In this case, the escape of photons is also significantly slowed down. Hence, a significant cascade could still materialize at far distances from the black hole. It also suggests that a cascade that has begun developing closer to the black hole (say, within the external photon fields) can continue to grow even at far distances.

2.2. Particle density evolution without secondary injection

We are interested in the time-dependent evolution of the total particle density, while neglecting its detailed energy-dependent evolution, as the latter does not change the overall density. Considering only time-dependent primary injection and escape, but no secondary injections, the kinetic equation for the particle distribution n(t) becomes:

$$\begin{array}{c}\hfill \frac{\partial n(t)}{\partial t}+\frac{n(t)}{t_{\mathrm{esc}}(t)}=Q(t),\end{array}$$(4)

with the analytical solution:

$$\begin{array}{c}\hfill n(t)=exp{\displaystyle (-{\int }^t\frac{\mathrm{d}t\prime }{t_{\mathrm{esc}}(t\prime )})\underset{0}{\overset{t}{\int }}Q(t″)exp\left({\int }^{t″}\frac{\mathrm{d}t\prime }{t_{\mathrm{esc}}(t\prime )}\right)\mathrm{d}t″·}\end{array}$$(5)

The particle injection rate Q(t) is coupled to the particle injection luminosity L_inj(t), taking the form:

$$\begin{array}{c}\hfill Q(t)=\frac{L_{\mathrm{inj}}(t)}{V(t)E(\gamma )}=q_0{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-3-p},\end{array}$$(6)

where E(γ) is a function of the injection energy spectrum, which is of no concern to us, while V(t) is the volume of the spherical blob. With the assumption $L_{\mathrm{inj}}(t)=L_0{[R_0/R(t)]}^{\mathrm{p}}=L_0{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\mathrm{p}}$ and the definition q₀ := L₀/[V(0)E(γ)], the second equality of Eq. (6) can readily be achieved. The power-law index, p, describes the decrease of the injection luminosity as a function of R(t). As the jet power is proportional to R(t)², p = 2 implies a constant jet particle injection power. While p is a free parameter, we will mostly use p = 2 below. For constant injection and escape – that is a straight jet, α = 0 – the solution of Eq. (5) is n = q₀t_esc(0), as expected.

With Eq. (2), the integrals in the exponentials of Eq. (5) are easily solved:

$$\begin{array}{cc}\hfill n(t)& =q_0{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}{\displaystyle \underset{0}{\overset{t}{\int }}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t″)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(3+p)}\mathrm{d}t″}\hfill \\ \hfill & =\frac{q_0{t}_{\mathrm{esc}}(0)}{1-{\eta }_{\mathrm{esc}}\alpha (2+p)}[{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-(2+p)}-{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}]·\hfill \end{array}$$(7)

This equation is positive for all times and p. For p = (η_escα)⁻¹ − 2, the density becomes

$$\begin{array}{c}\hfill n(t)=\frac{q_0{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}ln(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)·\end{array}$$(8)

Equation (7) holds for all particle species as long as there is no secondary injection or particle destruction. Below, we refer to Eq. (7) as the “standard solution”.

Figure 1 exemplifies Eq. (7) for two values of p and various values of η_escα. The parameters q₀ and R₀ are chosen to resemble typical one-zone parameters, while η_esc = 3 implies $\alpha \in [0.0\overline{3},0.3]$, resembling the typical range of this parameter (Pushkarev et al. 2017; Boula & Mastichiadis 2022). Interestingly, both p and η_escα have a similar effect. A larger p implies a faster decrease in the injection luminosity and thus a faster decrease of the density. More importantly though, larger opening angles imply a slower escape of particles from the emission region, and thus a slower decline of the density (see Eq. (3)). The peak in each curve is attained at roughly t ∼ t_esc(0) (red vertical line), however with a slightly earlier peak time for larger opening angles because of the quicker injection density decrease. For α → 0, the lines would approach a constant as the blob approaches the steady-state. The solid lines in Fig. 1 are derived from simulations using OneHaLe, confirming the analytical result.

Fig. 1. Particle density as a function of time for various cases of ηescα and p as labeled. The dashed lines are the analytical result, Eq. (7). Vertical red line is tesc(0). Dash-double-dotted line marks t−3. Further parameters are q0 = 1 cm−3 s−1, R0 = 5 × 1015 cm, and ηesc = 3. Solid lines are from simulations using OneHaLe.

The dash-double-dotted line in Fig. 1 indicates n ∼ t⁻³, which is the expected evolution of the density if it were solely due to the increase in volume of the blob. Most model lines are harder than this line implying a continuous increase of particle number in the emission region, even for p = 2. It is only for small opening angles (black and blue cases) for p = 2 that the model lines are softer and the total number of particles decreases at large times. This can be easily understood from Eq. (7), as the evolution at late times is governed by the power-law index: (2 + p) or 1/η_escα, whichever is smaller. If the injection switches off entirely, the density drops ∝t^{−1/η_escα}, which for large opening angles is the same behavior as with continuous injection. This shows the aforementioned trapping of particles in the (rapidly) expanding blob.

2.3. Linear cascade evolution

In this section, we treat the additional injection of electrons through γ − γ pair production of a γ ray colliding with a soft photon. We do so by first deriving three simple (semi-)analytical scenarios. These are compared to simulations in Sec. 2.3.4, where we discuss both the success and failure of the approach. We continue to neglect the energy dependency of the process n this work, as we are only interested in the rough time dependency of the total particle density. In other words, we assume that pair production takes place. This clearly is a very strong simplification, which may not hold in many cases. In this scenario, the injection rate of pairs can be written as:

$$\begin{array}{c}\hfill Q_{\gamma \gamma }(t)=\xi n_{\mathrm{soft}}(t)n_{\gamma }(t),\end{array}$$(9)

with the soft photon density, n_soft(t), the γ-ray photon density, n_γ(t), and the correlation factor, ξ ∼ σ_Tc, absorbing all constants and energy dependencies.

Both photon distributions can have various underlying production processes. The soft photons are most likely electron-synchrotron photons, or thermal photons from external sources, such as the AD, the BLR, and the DT. The γ rays can be produced from electron-inverse-Compton emission or proton-dependent processes, such as proton-synchrotron or neutral-pion decay. The internal photon densities thus depend on the underlying particle distribution, implying n_phot(t)∼n_e/p(t), with the exception of SSC radiation, which depends quadratically on the electron distribution. As we only want to treat a linear cascade evolution – that is, Q_γγ shall depend at most linearly on the electron distribution – we will ignore from now on the following combinations: electron-synchrotron and electron-inverse-Compton, as well as external photon fields and electron-SSC.

In fact, the simplest case is hadronically induced γ rays absorbed by external photon fields, as this does not depend on the electron distribution at all. Cases with a linear dependency on the electron distribution are electron-synchrotron absorbing hadronically induced γ rays, as well as external photons absorbing electron-inverse-Compton radiation. This leaves us with three distinct cases:

$$\begin{array}{cc}\hfill Q_{\gamma \gamma }(t)& \sim \xi n_{\mathrm{ext}}(t)n_{\mathrm{p}}(t),\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill Q_{\gamma \gamma }(t)& \sim \xi n_{\mathrm{e}}(t)n_{\mathrm{p}}(t),\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill Q_{\gamma \gamma }(t)& \sim \xi n_{\mathrm{ext}}(t)n_{\mathrm{e}}(t).\hfill \end{array}$$(12)

The external photon fields comprising AD and isotropic sources, such as the BLR, are by themselves assumed to be time-independent. However, the motion of the blob changes the distance to these external sources, implying:

$$\begin{array}{c}\hfill n_{\mathrm{ext}}(t)=\frac{n_{\mathrm{AD}}}{{(1+\frac{\mathrm{\Gamma }c}{R_{\mathrm{AD}}}t)}^2}+n_{i}H\left[\frac{z_i}{\mathrm{\Gamma }c}-t\right]·\end{array}$$(13)

The first summand describes the AD field, which is roughly constant close to the accretion disk and falls off with distance-squared once the blob has traveled a distance that roughly corresponds to the AD radius, R_AD (Dermer & Schlickeiser 2002). The second summand represents isotropic photon fields (in the black hole frame), such as the BLR and the DT, within a given distance z_i from the black hole.

2.3.1. Case 1: External and proton-induced photons

In this case, Eq. (10) is simply added to Eq. (4) and in Eq. (5). The solution is thus separated into the primary injection following Eq. (7) and the secondary injection leads to:

$$\begin{array}{cc}\hfill n_{\mathrm{e}}^{\gamma \gamma }(t)& =\frac{\xi q_{0,\mathrm{p}}{t}_{\mathrm{esc}}(0)}{1-{\eta }_{\mathrm{esc}}\alpha (2+p)}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{t}{\int }}}[{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t″)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p)}-1]\hfill \\ \hfill & \times [\frac{n_{\mathrm{AD}}}{{(1+\frac{\mathrm{\Gamma }c}{R_{\mathrm{AD}}}t″)}^2}+n_{i}H\left[\frac{z_i}{\mathrm{\Gamma }c}-t″\right]]\mathrm{d}t″.\hfill \end{array}$$(14)

While the isotropic contribution can be easily integrated, the accretion disk contribution leads to a hypergeometric integral without a simple analytical solution.

For isotropic external photons, that is, n_AD = 0, the analytical result for Eq. (14) becomes

$$\begin{array}{cc}\hfill n_{\mathrm{e}}^{\gamma \gamma }(t)& =\frac{\xi q_{0,\mathrm{p}}{[t_{\mathrm{esc}}(0)]}^2{n}_i{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}}{[1-{\eta }_{\mathrm{esc}}\alpha (2+p)][1-{\eta }_{\mathrm{esc}}\alpha (1+p)]}\hfill \\ \hfill & \times \{\begin{array}{cc}[({(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(1+p)}-1)-t]\hfill & \text{for}t\le \frac{z_i}{\mathrm{\Gamma }c}\hfill \\ [({(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}\frac{z_i}{\mathrm{\Gamma }c})}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(1+p)}-1)-\frac{z_i}{\mathrm{\Gamma }c}]\hfill & \mathrm{for}t>\frac{z_i}{\mathrm{\Gamma }c}\hfill \end{array}.\hfill \end{array}$$(15)

The sum of Eqs. (7) and (15) is shown in Fig. 2. Compared to the standard solution, Eq. (7), there are significantly more particles injected within the boundaries of the external field, z_i. The boundary of z_i is noticeable by the break in the dash-dotted and solid lines. The peak of the particle density is attained later for larger opening angles compared to the standard solution. Beyond z_i, the decay of the density again depends strongly on the opening angle. The standard solution is only reached for small opening angles. The semi-analytical solutions of Eq. (14) for the AD, that is n_i = 0, are shown in Fig. 3. This behavior is comparable to the isotropic case except for an earlier peak owing to R_AD ≪ z_i.

Fig. 2. Electron density as a function of time for various cases of ηescα as labeled. Secondary electrons are produced from proton-induced γ rays absorbed by isotropic external photons. The thick dash-dotted line is the analytical solution of Eq. (14), assuming only one isotropic photon field with ni = 1 × 1010 cm−3, and zi = 1 pc. The dashed lines are the standard solution, Eq. (7). The vertical red line is tesc(0), while the horizontal dotted lines mark the maximum level of the standard solution and when they are reached by the dash-dotted lines. Further parameters are p = 2, q0, e = q0, p = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10. The solid lines mark simulations with OneHaLe.

Fig. 3. Same as Fig. 2, except that the absorber is AD photons with nAD = 1 × 1011 cm−3, and RAD = 1 × 1016 cm. Further parameters are p = 2, q0, e = q0, p = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.

In both figures (Figs. 2 and 3) the horizontal dotted lines mark the maximum level of the standard solution. They indicate the time delay for the semi-analytical solution to drop below this density. Naturally, this happens later for the isotropic field than for the AD field and depends strongly on the opening angle. In fact, the time delay can be orders of magnitude implying a significantly more pronounced and longer high state.

2.3.2. Case 2: Electron-synchrotron and proton-induced photons

In adding Eqs. (11) to (4), we obtain with a slight rearrangement:

$$\begin{array}{c}\hfill \frac{\partial n_{\mathrm{e}}(t)}{\partial t}+(\frac{1}{t_{\mathrm{esc}}(t)}-\xi n_{\mathrm{p}}(t))n_{\mathrm{e}}(t)=Q_{\mathrm{e}}(t).\end{array}$$(16)

The linear cascade acts as a “catastrophic” injection with time scale t_γγ(t) = [ξn_p(t)]⁻¹. We directly obtain the generalization of Eqs. (5) and (7) as:

$$\begin{array}{cc}\hfill n_{\mathrm{e}}(t)& =exp{\displaystyle (-{\int }^t\frac{\mathrm{d}t\prime }{t_{\mathrm{esc}}(t\prime )}+\xi {\int }^t{n}_{\mathrm{p}}(t\prime )\mathrm{d}t\prime )}\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{t}{\int }}{Q}_{\mathrm{e}}(t″)exp({\int }^{t″}\frac{\mathrm{d}t\prime }{t_{\mathrm{esc}}(t\prime )}-\xi {\int }^{t″}{n}_{\mathrm{p}}(t\prime )\mathrm{d}t\prime )\mathrm{d}t″}\hfill \\ \hfill & =q_{0,\mathrm{e}}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}exp{\displaystyle \left(\xi {\int }^t{n}_{\mathrm{p}}(t\prime )\mathrm{d}t\prime \right)}\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{t}{\int }}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t″)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(3+p)}exp(-\xi {\int }^{t″}{n}_{\mathrm{p}}(t\prime )\mathrm{d}t\prime )\mathrm{d}t″.}\hfill \end{array}$$(17)

As the protons follow Eq. (7), the integrals in the exponentials can be performed:

$$\begin{array}{cc}\hfill {\displaystyle {\int }^t}{n}_{\mathrm{p}}(t\prime )\mathrm{d}t\prime & =\frac{q_{0,\mathrm{p}}{t}_{\mathrm{esc}}{(0)}^2}{[1-{\eta }_{\mathrm{esc}}\alpha (2+p)](1+p)(1-{\eta }_{\mathrm{esc}}\alpha )}\hfill \\ \hfill & \times [(1-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }){(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-(1+p)}\hfill \\ \hfill & +(1+p){(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{1-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}]·\hfill \end{array}$$(18)

With this solution, the remaining integral in Eq. (17) can be solved numerically. The result is shown in Fig. 4.

Fig. 4. Same as Fig. 2, but for secondary electrons produced from proton-induced γ rays absorbed by electron-synchrotron photons, Eq. (17). Additional parameters are: p = 2, q0, e = 1 cm−3 s−1, q0, p = 1 × 103 cm−3 s−1, R0 = 5 × 1015 cm, and ηesc = 3.

The cascade develops only after t_esc(0), while the opening angle dictates the further evolution. Small opening angles result in a high number density and a quick decay, while large opening angles seemingly keep the cascade going for very long times without any significant decay. This is also indicated by the dotted lines.

2.3.3. Case 3: External and electron-IC photons

This case is similar to the case 2, except that we have to replace n_p(t) with n_ext(t) and the γ rays stem from IC scattering of the external photon fields. Then the electron density becomes

$$\begin{array}{cc}\hfill n_{\mathrm{e}}(t)& =q_{0,\mathrm{e}}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}\hfill \\ \hfill & \times exp(-\frac{\xi n_{\mathrm{AD}}{R}_{\mathrm{AD}}}{\mathrm{\Gamma }c(1+\frac{\mathrm{\Gamma }c}{R_{\mathrm{AD}}}t)}+\xi n_{i}t\stackrel{\sim }{H}\left[t\frac{z_i}{\mathrm{\Gamma }c}\right])\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{t}{\int }}}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t″)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(3+p)}\hfill \\ \hfill & \times exp(\frac{\xi n_{\mathrm{AD}}{R}_{\mathrm{AD}}}{\mathrm{\Gamma }c(1+\frac{\mathrm{\Gamma }c}{R_{\mathrm{AD}}}t″)}-\xi n_{i}t″H\left[\frac{z_i}{\mathrm{\Gamma }c}-t″\right])\mathrm{d}t″,\hfill \end{array}$$(19)

with the modified Heaviside function $\stackrel{\sim }{H}[xa]=1$ for x ≤ a and $\stackrel{\sim }{H}[xa]=a/x$ for x > a.

For isotropic external photon fields, that is, n_AD = 0, Eq. (19) can be analytically integrated:

$$\begin{array}{cc}& n_{\mathrm{e}}(t\le \frac{z_i}{\mathrm{\Gamma }c})=\frac{q_{0,\mathrm{e}}}{\xi n_i}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}{\left(\frac{{\eta }_{\mathrm{esc}}\alpha }{\xi n_i{t}_{\mathrm{esc}}(0)}\right)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(3+p)}\hfill \\ \hfill & \times exp\left[\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)\right]\hfill \\ \hfill & \times [\mathrm{\Gamma }(\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p),\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha })\hfill \\ \hfill & -\mathrm{\Gamma }(\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p),\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t))]\hfill \end{array}$$(20)

$$\begin{array}{cc}& n_{\mathrm{e}}(t>\frac{z_i}{\mathrm{\Gamma }c})=q_{0,\mathrm{e}}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{-\frac{1}{{\eta }_{\mathrm{esc}}\alpha }}exp\left(\xi n_i\frac{z_i}{\mathrm{\Gamma }c}\right)\hfill \\ \hfill & \times \{exp\left(\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }\right)\frac{{\left(\frac{{\eta }_{\mathrm{esc}}\alpha }{\xi n_i{t}_{\mathrm{esc}}(0)}\right)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(3+p)}}{\xi n_i}\hfill \\ \hfill & \times [\mathrm{\Gamma }(\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p),\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha })\hfill \\ \hfill & -\mathrm{\Gamma }(\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p),\frac{\xi n_i{t}_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}\frac{z_i}{\mathrm{\Gamma }c}))]\hfill \\ \hfill & +\frac{t_{\mathrm{esc}}(0)}{1-{\eta }_{\mathrm{esc}}\alpha (2+p)}[{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p)}\hfill \\ \hfill & -{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}\frac{z_i}{\mathrm{\Gamma }c})}^{\frac{1}{{\eta }_{\mathrm{esc}}\alpha }-(2+p)}]\},\hfill \end{array}$$(21)

with the incomplete Gamma-function Γ(q, x). The result is shown in Fig. 5, where compared to Fig. 2, a strong pile-up is visible, which results from ξn_it > 1 close to z_i. This may be due to an extreme choice of parameters, however, it does display the importance of the immediate feedback of the cascade on the γ rays in this case compared to case 1. We note that the pile-up for η_escα = 0.1 in Fig. 5 might be influenced by some approximation inaccuracies for small values of η_escα and so, a smaller pile-up is to be expected in this case.

Fig. 5. Same as Fig. 2, but for secondary electrons produced from IC-induced γ rays absorbed by isotropic external photons, Eqs. (20) and (21), with ni = 4 × 107 cm−3, and zi = 1 pc. Further parameters are p = 2, q0, e = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.

For AD photons, we need to numerically integrate Eq. (19) with the result shown in Fig. 6. The result is fairly similar to case 1, which may, again, be due to parameter choices, but also due to the argument of the exponential is a decreasing function with time. The latter effect reduces the amount of absorbing photons before the γ rays have fully developed, resulting in a reduced cascade compared to the isotropic case.

Fig. 6. Same as Fig. 2, but for secondary electrons produced from IC-induced γ rays absorbed by AD photons with nAD = 1 × 1010 cm−3, and RAD = 1 × 1016 cm. Further parameters are p = 2, q0, e = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.

In both cases, the dotted lines are roughly as long as in case 1. In the isotropic case, this is significantly influenced by the pile-up close to z_i. With a slightly smaller size for the isotropic field and a resulting smaller pile-up for case 3, the duration of the high-state in case 3 would be shorter than in case 1. However, no such strong influence is expected in the AD case.

2.3.4. Simulations

We use the OneHaLe code (Zacharias 2021; Zacharias et al. 2022, see the appendix for a brief description) to model the time-dependent evolution of the blob. In order to comply with some of the assumptions above, we have to slightly tweak the code. The γ − γ pair production injection rate is given by Aharonian et al. (1983), Cerruti et al. (2021)

$$\begin{array}{cc}\hfill Q_{\gamma \gamma }(\gamma )=& 2\frac{3{\sigma }_{\mathrm{T}}c}{32}{\displaystyle \underset{\gamma }{\overset{\infty }{\int }}\mathrm{d}ϵ\frac{n_{\mathrm{ph}}(ϵ)}{{ϵ}^3}\underset{\frac{ϵ}{4\gamma (ϵ-\gamma )}}{\overset{1}{\int }}\mathrm{d}\stackrel{\sim }{ϵ}\frac{n_{\mathrm{ph}}(\stackrel{\sim }{ϵ})}{{\stackrel{\sim }{ϵ}}^2}}\hfill \\ \hfill & \times [\frac{4{ϵ}^2}{\gamma (ϵ-\gamma )}ln\left(\frac{4\gamma \stackrel{\sim }{ϵ}(ϵ-\gamma )}{ϵ}\right)-8ϵ\stackrel{\sim }{ϵ}\hfill \\ \hfill & +\frac{2{ϵ}^2(2ϵ\stackrel{\sim }{ϵ}-1)}{\gamma (ϵ-\gamma )}-(1-\frac{1}{ϵ\stackrel{\sim }{ϵ}}){\left(\frac{{ϵ}^2}{\gamma (ϵ-\gamma )}\right)}^2],\hfill \end{array}$$(22)

with the photon distribution, n_ph, photon energies, ϵ and $\stackrel{\sim }{ϵ}$, normalized to the electron rest energy, and the electron Lorentz factor, γ. Equation (22) requires $ϵ\gg 1/\stackrel{\sim }{ϵ}$, ensuring, together with $\stackrel{\sim }{ϵ}<1$, an absorber with low photon energies for the γ rays.

In order to suppress severe cooling (due to the extreme parameter settings; see Table 1), we switch off electron IC cooling. While this is a massive change, it is necessary to create sufficient synchrotron or IC photons. Furthermore, the magnetic field is kept constant throughout these simulations. Lastly, as we are only interested in the additional injection of pairs from γ − γ processes, the secondary injections from muon decay and Bethe-Heitler pair production have also been disabled. The input parameters for the simulations are given in Table 1, while the results are shown as solid lines in Figs. 2–6. We emphasize that these parameters are not chosen for realism but merely to reproduce the semi-analytical models as closely as possible.

In this regard, the external fields labeled “BLR” and “DT” should just be considered as two separate external fields without any intended resemblance to actual BLR and DT photon fields. In case 1, the external fields act as absorbers. In case 3, they are needed for the IC process, while the absorption is done by electron-synchrotron emission.

Cases 1 and 3 can be reproduced reasonably well. In case 1, the isotropic photon field is made of two thermal fields with different temperature and luminosity, while in the corresponding case in case 3, only one thermal field is needed. In the latter case, the pile-up close to the edge of the external field is present, but not as strong as in the analytical model². The accretion disk field is represented by a Shakura-Sunyaev disk (Shakura & Sunyaev 1973), with luminosity L_AD = η_EddL_Edd. In case 3, the AD field absorbs IC emission scattering both isotropic fields (but not the AD field).

We have found no parameter set that could reproduce the semi-analytical result in case 2. The simulations with OneHaLe with the parameters given in Table 1 show a pronounced maximum shortly after t_esc(0) followed by a decrease depending on the opening angle. This seems reasonable, as the dashed lines roughly correspond to the evolution of the γ rays (p-synchrotron). In turn, the number of “absorbees” (i.e., γ rays) is quickly reduced. The long evolution suggested by the semi-analytical model would only be achievable if the absorber were to increase in such a way that a higher fraction of the γ rays would be progressively absorbed. Then, at late times, practically all γ rays would be absorbed to keep the cascade going. This is implausible and marks the limit of the simple analytical model that ignores all the energy-dependencies of the cross-section.

The influence of the parameters on the simulations is shown in Fig. 7, where three variations for each case have been derived. It is clear that reproducing the analytical model requires a specific set of parameters. This also shows the limits of the approach, however, it is nonetheless reassuring that we have managed to reproduce two out of the three cases overall.

Fig. 7. Small parameter study using OneHaLe. The black lines correspond in style to those in Figs. 2–6, but for the case ηescα = 0.5. The solid colored lines show simulations varying the input parameter as labeled. In the bottom panels, the variation in the electron parameters (magenta lines) preserves the initial injection rate of ∼1 cm−3 s−1.

2.4. Interlude

The sections above have demonstrated the influence of the expansion of the moving blob on the development of the pair cascade. The effects are two-fold. First, a significantly higher density can be achieved through the pair production, as compared to no secondaries. Second, a larger opening angle can drag out the peak of the density evolution. While the densities for the standard solution peak around t_esc(0), the curves with pair injection typically peak later than this timescale. The exception to this rule are secondary injections involving AD photons, because of the timescale, namely, R_AD/Γ/c < t_esc(0). For small values of z_i, this could become the case as well for the isotropic photon field. For internal photon fields (case 2), however, the peak is always attained after t_esc(0), as the internal fields only reach their maximum at this time and the developing cascade reinforces the interaction. Non-linear cascades, which we have not treated here, might increase this effect, as both absorber and absorbee are increased at later times.

The dotted lines in Figs. 2–6 mark the peak density of the standard solution and indicate the time at which the models with pair injection drop below this level. This typically occurs at least an order of magnitude after t_esc(0). Of course, this also strongly depends on the amount of pairs produced, the parameters, and the opening angle – highlighting the importance of pair production and the geometry of the jet for the development and duration of a flux high state.

Naturally, the inclusion of the energy dependency on the particle distribution, as well as the pair production process, can have a significant influence. In the absence of efficient reacceleration of particles, the particle cooling will lead to a drop in the pair production process, as less energetic γ-ray photons will be produced. The adiabatic expansion of the blob will also cause a reduction in the energy density of the particles and the photon fields, similarly reducing the pair production. Hence, the time evolution might be shorter than envisaged in the present semi-analytical model.

3. Modeling based on PKS 1510−089

The flat spectrum radio quasar (FSRQ) named PKS 1510−089, located at a redshift of 0.361, is one of the most famous blazars known to date. It was detected at very-high-energy γ rays with H.E.S.S. in 2009 (H.E.S.S. Collaboration 2013) and has since received significant attention (e.g., Ahnen et al. 2017; H.E.S.S. Collaboration 2021). Its modeling typically requires both the BLR and DT photon fields (Barnacka et al. 2014) or multiple zones (Nalewajko et al. 2012; Prince et al. 2019).

While no modeling of the source or a specific light curve is approached in this work, we use this source to study the impact of the pair cascade in a moving, expanding blob in terms of enhancing flaring events in a realistic setting. The modeling parameters are given in Table 2 and are based upon the works by Nalewajko et al. (2012) and Barnacka et al. (2014). We derived multiple modelings: (1) a steady-state model similar to the aforementioned works to provide a baseline flux that may describe PKS 1510−089 on average; (2) a moving blob using only leptonic processes similar to Saito et al. (2015); and (3) a moving blob including hadronic interactions. For the moving blobs we again derive models for various opening angles as before. We also derive specifically models without γ − γ pair production to show the difference. For the moving-blob models 2 and 3, the launching position is chosen relatively close to the black hole, where the jet has not expanded much (Zacharias et al. 2022). At this position the blob is still small in order to fit into the jet, and its magnetic field is high compared to the steady-state one-zone model, which is located further down in the jet.

One important addition to the simulations of the previous section is the consideration of magnetic flux conservation, that is, B(z) = B₀ R₀/R(z), assuming a dominating toroidal guide field (Kaiser 2006). This also implies that the larger the opening angle, the faster the strength of the magnetic field drops. This may have a significant influence on the synchrotron emission, compared to the expectation from the simple model above.

3.1. Light curves

The multi-wavelength light curves are shown in Fig. 8. The wavebands are the very-high-energy γ-ray (VHE) band (E > 30 GeV³, corresponding to the energy band of the forthcoming CTA observatory), the high-energy γ-ray (HE) band (100 MeV < E < 100 GeV, corresponding to the energy band of the Fermi-LAT), the X-ray band (2 − 10 keV, corresponding to the range of Swift-XRT), and the optical R-band. It is immediately obvious that the flare evolution depends strongly on the wavelength and the presence of relativistic protons.

Fig. 8. Model light curves based upon parameters describing PKS 1510−089 for VHE γ rays (top row), HE γ rays (2nd row), X-rays (3rd row), and the optical R-band (bottom row) for purely leptonic processes (left column) and including hadronic processes (right column). Displayed models are either with (thick dashed) or without (solid) γ − γ pair production. Colors are for various opening angles as labeled. The vertical lines mark tesc(0) (red), the crossing of the BLR edge (blue) and of the DT edge (magenta); all times in the observer’s frame. For clarity, the launch of the blob takes place at tobs = 0.1 d.

Models without pair production (solid lines) do not show major differences between the leptonic and hadronic simulations. Except for the VHE band, which is fully absorbed at early times, all light curves rise quickly to a maximum at the initial escape timescale, t_esc(0), followed by a subsequent decay. The switch from IC-BLR cooling to IC-DT or synchrotron cooling during the crossing of the BLR edge⁴ results in a minor bump in the X-rays and the R-band, which is absent in the HE band. The VHE light curve shows a minor bump at this point, as the environment is now optically thin to VHE photons. Except for the HE band, the decay of the light curve is faster for larger opening angles, which is a consequence of the diminishing magnetic field, and hence a faster decay of the synchrotron emission. We should note here that the X-ray emission from the blob is also influenced by electron-synchrotron emission, as the (initial) high magnetic field shifts the synchrotron peak to higher energies compared to the steady-state model. As the HE band is governed by inverse-Compton emission of external photon fields, the magnetic field has basically no effect, and the decay of the light curve is fully governed by the evolution of the electron density. In the hadronic model and for large opening angles, the VHE light curve shows an additional bump around the time of the DT-edge crossing. This radiation is electron-synchrotron emission from secondary particles injected from muon decay.

The inclusion of the pair production has almost no consequence for the VHE band, except for hadronic simulations of large opening angles, where the DT-edge-crossing bump is much more pronounced and shows a plateau phase lasting for a few hours. At that time, a peak also shows up in the X-ray domain for large opening angles. It is induced by pair production and is also visible (although less pronounced) in the leptonic simulations. The peak flux in the main flare is significantly higher compared to the case without pair production. Similarly, the HE band exhibits higher peak fluxes compared to the solid lines. In both bands, HE and X-rays, the peak fluxes are much more pronounced in the hadronic simulations indicating the significant impact of relativistic protons on the pair production. The R-band shows no change between pair or no pair production, implying that pairs are produced predominantly at higher energies. While cooling will reduce the energy of the secondary pairs, their density is insufficient to add to the synchrotron emission of the primaries in the R-band.

The most important effects in the light curves shown in Fig. 8 take place within a few days in the observer’s frame. Figure 9 displays the first 100 days in the evolution of the radio band at 15 GHz. The rise of the light curve begins after a few days, when the flares in the other bands have already ceased (see also Boula & Mastichiadis 2022; Tramacere et al. 2022). For small opening angles, the radio flare diminishes within 50 days, while for large opening angles, the light curve keeps rising. While it is not shown, for η_escα = 0.9 the peak is reached only after more than 1000 days. There is barely a difference between the cases with and without pair production for the same reason as in the R-band.

Fig. 9. Same as Fig. 8, but for the radio band at 15 GHz.

3.2. Luminosity evolution

An important measure for any blazar model is the luminosity of the emission region. Power is contained in particles, magnetic field, and radiation. The evolution of the luminosities as a function of time is shown in Fig. 10. By design, the magnetic field luminosity is constant with a value of $L_B^{\mathrm{obs}}=1.8\times {10}^{44}$ erg s⁻¹, which is why we do not display it in Fig. 10. The total luminosity in the bottom row of that figure is the sum of the electron, proton, magnetic field, and radiation luminosities. The individual luminosities are calculated from the respective energy densities in the comoving frame u_i(t) according to $L_i^{\text{obs}}(t)=u_i(t){\mathrm{\Gamma }}^2\pi R{(t)}^{2}c$ deriving the bulk Lorentz factor Γ from the assumption Γ = δ.

Fig. 10. Electron luminosity (top row), proton luminosity (second row), radiation luminosity (third row), and total luminosity (bottom row) as a function of time for pure electron (left column) and electron-proton simulations (right column) for different opening angles (color code as labeled) using the parameters of PKS 1510−089 as in Fig. 8. Simulations without pair production are marked by solid lines, while dashed lines refer to simulations with pair production. The vertical lines have the same meaning as in Fig. 8. The black dot-dashed and black dash-double-dotted lines mark the Eddington and AD luminosities, respectively. We note the logarithmic time axis.

For small opening angles (black and blue lines), the respective luminosities quickly rise to the maximum as the source is being filled with particles and radiation, and then decrease. In turn, the total luminosity reaches a bottom value that is given by the constant magnetic field luminosity. It also implies that the initial particle and radiative powers are much larger than the magnetic field power. For η_escα = 0.5, the particle powers reach a constant level implying that the injection and escape of particles is fully countered by the expansion. In the presence of protons (right column in Fig. 10), the electrons slowly continue to rise in the case where pair production is enabled (magenta dashed line). This implies continuous pair creation even at very late times. For large opening angles (orange and yellow lines), the particle luminosities keep rising, as the rapid expansion does not allow for a meaningful escape of particles. The radiative luminosity in this case continues to decrease initially, but reverses this trend for late time, as the continuing accumulation of particles also increases the amount of photons, which also remain longer and longer in the emission region.

Interestingly, the total luminosity is initially dominated by radiation and only after about 1 day in the observer’s frame does it change to particle or magnetic field dominance. It implies that transfer of power into radiation is initially very efficient and that the particles are efficiently cooled. With the decrease in the magnetic field strength, and the decrease of density of internal and external photon fields due to expansion and motion, the radiation production becomes less efficient over time and its luminosity drops more rapidly than that of the particles.

In all cases, the comparison between particles and magnetic fields suggests that the emission region is initially particle dominated. For large opening angles, this is maintained throughout, while for small opening angles, the magnetic field dominates the particles at late times. This is contrary to most models, where an initial magnetic field dominance is expected which changes to particle dominance later (e.g., Zacharias et al. 2022, and references therein). The reason is that we chose a relatively modest initial magnetic field of only 20 G along with a linear decrease with radius. If the magnetic field were initially mostly poloidal, one would expect a quadratic decrease with radius (Kaiser 2006), and a much faster decrease. In this case, the magnetic power would also drop as a function of time. While our initial value is in line with the measurement in M 87 (Event Horizon Telescope Collaboration 2021), theoretical works suggest much larger values along with a rapid decrease (Zdziarski et al. 2022). Testing such scenarios is left for a future work.

The black horizontal lines in Fig. 10 mark the Eddington luminosity and the AD luminosity, respectively. As the accretion process fuels the jet, the jet luminosity must be compared with these values. Apparently, all models remain initially close to or below the limit of the Eddington luminosity. In case of protons and pair production (right column, dashed lines), the Eddington luminosity is initially surpassed by a factor 5. This is, however, due to the large radiative power, which includes also the IC scattering of the external fields. While the BLR and the DT are less luminous than the AD, they are strongly beamed in the comoving frame. Thus, the short excess of the Eddington luminosity can be expected and is within reasonable bounds. The continued increase of jet power for large opening angles has already been discussed and naturally surpasses the Eddington limit by far. Clearly, the model is not realistic in these cases.

4. Conclusions

In this paper, we have discussed the particle evolution within a moving, expanding emission region within the jet of a blazar (Boula & Mastichiadis 2022; Tramacere et al. 2022). With a simple semi-analytical model, we have shown that the opening angle has an important influence on the particle distribution and the electron-positron cascade. With an increasing opening angle the high-state (including the cascade) may last much longer owing to the increased escape timescale for both the particles and photons.

With the help of the OneHaLe code (Zacharias 2021; Zacharias et al. 2022), we have shown that the simple model can be reproduced for linear cascades induced by external photon fields, while the linear cascade resulting from internal processes was not reproduced successfully. Naturally, due to the neglect of the energy dependencies of the particle evolution and the cascade, it required quite specific, and extreme parameters in the simulation. Nonetheless, the exercise has demonstrated the importance of the opening angle on the evolution of the particle distribution.

The simple model may be improved by considering delta-function approximations to the energy dependency of the particle distribution and the cascade. This may provide a more reasonable view on the evolution providing a broader range of application. These calculations are left to the interested reader.

In order to obtain realistic simulations, we have made use of parameters from PKS 1510−089 (Nalewajko et al. 2012; Barnacka et al. 2014). The leptonic and hadro-leptonic simulations were conducted for various opening angles. From the resulting light curves in Fig. 8, we can derive the following conclusions. A moving blob induces a rapid flare. The evolution depends critically on the blob’s speed and on the initial size at launch. Observations in the X-ray and VHE γ-ray are crucial to determine the opening angle and the particle composition, as for large opening angles these bands may (or may not) show a delayed flare when the blob moves out of the DT. These secondary flares are significantly influenced by the pair production process and are more pronounced if a hard and highly energetic proton distribution is present. Interestingly, the optical band is also useful to obtain information on the opening angle, as a small opening angle induces a longer decay phase due to the slower decay of the magnetic field. In the HE γ-ray band, the flare is bright and short with barely a difference between the various opening angles. Importantly, the HE band peaks slightly before the other bands (cf. the red vertical line in Fig. 8), whereas the X-ray and R-band peak at the time corresponding to t_esc(0). In all bands, the (primary) flare is asymmetric with a fast rise and a slow decay (Saito et al. 2015; Boula & Mastichiadis 2022; Tramacere et al. 2022).

Truly simultaneous and high-cadence observations are paramount to obtain all the potential information. The main flare lasts just for a few hours, and the subtle differences between peak times and decay profiles might be hidden in the statistical errors or a delay in observation. It should also be noted that it is worthwhile to continue observing for a few days in the X-ray and VHE γ-ray band to search for the secondary flare.

Radio observations may also provide indications for the opening angle. Large values of η_escα result in years-long variations, while small opening angles result in flares lasting a few tens of days. Single-dish radio observations every few days for up to 100 days after the event at higher energies should provide sufficient information on the opening angle. The downside could be the higher chance of multiple events superposing each other.

The jet luminosities are mostly within reasonable bounds compared to the Eddington luminosity. However, for large opening angles, the particle luminosities keep increasing due to the lack of escape, which results in unreasonable jet powers. Additionally, the magnetic field evolution is such that its power remains constant. In turn, the initially partical-dominated emission region becomes magnetic-field-dominated at later times and small opening angles. This is contrary to standard expectations. A scenario, which is more in line with the usual jet evolution, will be discussed elsewhere.

The above mentioned evolutionary details of the flare depend on the type of object. The given results are based on parameters describing the FSRQ PKS 1510−089. In BL Lac objects with much weaker or even absent external fields, the effects induced from the pair cascade are probably much less pronounced (Zacharias et al. 2022). In this case, the solid lines in Fig. 8 might already be an adequate description and the secondary flares in the X-ray and VHE γ-ray band are absent. Similarly, the details of the magnetic field decay may have a significant influence on the evolution of the flux, especially in the R-band and needs to be discerned from the expansion profile.

Lastly, using the standard one-zone description implies that at all times an instantaneous particle spectrum was used. While shock acceleration can quickly accelerate particles (Böttcher & Baring 2019), it would be useful to properly include the particle acceleration in this model as it might have an important influence on the light curve evolution at early, but also at late times. At early times, the delayed injection of high-energetic particles would slow down the flare evolution, and might also reduce the strength of the cascade development. At late times, the weakened magnetic field might not be able to accelerate particles to the highest energies, which might reduce or even inhibit the development of the secondary flare. It would also substantially reduce the particle luminosities. In order to adequately account for these effects, a two-zone model is needed (Weidinger & Spanier 2015; Chen et al. 2015; Dmytriiev et al. 2021), which is beyond the scope of this paper.

Acknowledgments

I would like to thank Andreas Zech, Markus Böttcher, Anton Dmytriiev, and Stella Boula for fruit- and beerful discussions on the content of the manuscript. I extend my gratitude to the anonymous referee for valuable comments that helped to improve the manuscript. M.Z. acknowledges postdoctoral financial support from LUTH, Observatoire de Paris. Simulations for this work have been performed on the TAU-cluster of the Centre for Space Research at North-West University, Potchesftroom, South Africa. The OneHaLe code is available upon request to the author.

References

Appendix A: The OneHaLe code

The OneHaLe code is a time-dependent, one-zone, hadro-leptonic radiation code that solves the Fokker-Planck equation of the particle distribution for protons, charged pions, muons, and electrons (including positrons). In each time step, the radiation transport equation is solved allowing for the direct feedback of the particle and photon interactions. The Fokker-Planck equation for the particle distribution of species, i, is:

$$\begin{array}{cc}\hfill \frac{\partial n_i(\chi ,t)}{\partial t}& =\frac{\partial }{\partial \chi }\left[\frac{{\chi }^2}{(a+2)t_{\mathrm{acc}}}\frac{\partial n_i(\chi ,t)}{\partial \chi }\right]\hfill \\ \hfill & -\frac{\partial }{\partial \chi }\left({\dot{\chi }}_i{n}_i(\chi ,t)\right)+Q_i(\chi ,t)-\frac{n_i(\chi ,t)}{t_{\mathrm{esc}}}-\frac{n_i(\chi ,t)}{\gamma t_{i,\mathrm{decay}}^{\ast }}·\hfill \end{array}$$(A.1)

The distributions are given as a function of normalized momentum χ = γβ, with the particle Lorentz factor γ and its corresponding speed β normalized to the speed of light. The first term on the right-hand-side is momentum diffusion representing Fermi-II acceleration using hard-sphere scattering with the ratio a of shock to Alfvèn speed. The second term marks continuous momentum gains and losses. Gains are achieved through Fermi-I acceleration, while losses depend on the particle species and include synchrotron, adiabatic, Bethe-Heitler, pion production, and inverse-Compton processes. The third term marks the injection term, while the forth term marks the catastrophic escape of particles. The last term is the decay term for unstable particles.

The acceleration timescale is given as a multiple of the escape timescale: t_acc = η_acct_esc. It merely marks the reacceleration of particles in the acceleration zone and does not provide “first-principle” acceleration, which typically requires a much smaller zone (e.g., Weidinger & Spanier 2015; Chen et al. 2015; Dmytriiev et al. 2021). The initial acceleration is mimicked through the primary injection term Q_i, which for protons and electrons takes the form of a power-law between minimum and maximum Lorentz factors, γ_i, min and γ_i, max, respectively, with spectral index, s_i. The injection of pions and muons is directly calculated from the respective interactions (using the template approach of Hümmer et al. 2010) and decays. Secondary electrons are injected from muon decay, Bethe-Heitler pair production, and γ − γ pair production. Further details are given in Zacharias (2021) and Zacharias et al. (2022). These references provide the gory details about the particle and photon integrals solved in the code.

A.1. Acceleration and cooling timescales

In this section, we briefly discuss the influence of the expanding blob scenario on the acceleration and cooling timescales for illustrative purposes. These are not necessarily considered in the discussion and simulations of the main part. For the (re-)acceleration of particles, OneHaLe includes acceleration terms based on hard-sphere scattering. This is parameterized by an energy-independent timescale, namely, t_acc(t) = η_acct_esc(t), where η_acc is a free parameter. Hence,

$$\begin{array}{c}\hfill t_{\mathrm{acc}}(t)={\eta }_{\mathrm{acc}}{t}_{\mathrm{esc}}(0)(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)·\end{array}$$(A.2)

The cooling of particles also becomes time-dependent, as several important variables are a function of radius. For instance, with the conservation of magnetic flux, the magnetic field B(t) can be written as:

$$\begin{array}{c}\hfill B(t)=\frac{B_0}{{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^b},\end{array}$$(A.3)

with the initial magnetic field B₀, and the power b ∈ [1, 2] of the radial dependence (b = 1 for a toroidal guide field, b = 2 for poloidal guide field, Kaiser 2006). The synchrotron cooling timescale then becomes:

$$\begin{array}{c}\hfill t_{\mathrm{syn}}(\gamma ,t)=\frac{6\pi m_e{c}^2}{c{\sigma }_T{B}_0^2}{\left(\frac{m}{m_e}\right)}^3{\gamma }^{-1}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{2b},\end{array}$$(A.4)

with the particle Lorentz factor, γ, the rest mass of the electron, m_e, the rest mass of the particle, m, (m_p for protons), and the Thomson cross section, σ_T.

Following the prescription of Schlickeiser (2009), the SSC cooling timescale in the Thomson regime can be written as:

$$\begin{array}{c}\hfill t_{\mathrm{ssc}}(\gamma ,t)=\frac{m_e{c}^2}{3c_1{c}_2^2{\sigma }_T{P}_0{R}_0{B}_0^2}{\gamma }^{-1}F{(t)}^{-1}{(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t)}^{2b-1},\end{array}$$(A.5)

with c₁ = 0.684, c₂ = 1.856 × 10⁻²⁰ erg^1/2 cm^3/2, and P₀ = 2 × 10²⁴ erg⁻¹ s⁻¹. The function

$$\begin{array}{c}\hfill F(t)={\displaystyle \underset{0}{\overset{\infty }{\int }}}{\gamma }^2{n}_e(\gamma ,t)\mathrm{d}\gamma ,\end{array}$$(A.6)

depends on the instantaneous electron distribution function n_e(γ, t), which may have a complicated time-dependency resulting in a non-linear evolution of the particle distribution (Zacharias & Schlickeiser 2012). If we assume that the energy part of the particle energy distribution is time-independent, the time-dependency derived above for the particle distribution, Eq. (7), can be included here (cf. Fig. A.1).

Fig. A.1. Timescales in arbitrary units as labeled versus time in the comoving frame. The lines are indicative of the time-dependency. The vertical red line marks tesc(0). Free parameters are ηescα = 0.5, ηesc = 3, R0 = 5 × 1015 cm, b = 1, p = 2, Γ = 10, and $z_{\text{ext}}^{\prime }=3\times {10}^{18}$ cm.

The expansion of the blob results in adiabatic cooling with timescale:

$$\begin{array}{c}\hfill t_{\mathrm{adi}}(\gamma ,t)=\frac{\mathrm{\Gamma }}{3(1-{\gamma }^{-2})}\frac{t_{\mathrm{esc}}(0)}{{\eta }_{\mathrm{esc}}\alpha }(1+\frac{{\eta }_{\mathrm{esc}}\alpha }{t_{\mathrm{esc}}(0)}t),\end{array}$$(A.7)

where we followed the prescription of Zdziarski et al. (2014).

As we assume Γ = const., the cooling due to external isotropic photon fields such as the BLR is constant in time, as long as the blob is within the distance of the external field:

$$\begin{array}{c}\hfill t_{\mathrm{ext}}(\gamma ,t)=\frac{m_e{c}^2}{c{\sigma }_T{\mathrm{\Gamma }}^2{u}_{\mathrm{ext}}^{\prime }}{\gamma }^{-1}H\left[\frac{z_{\mathrm{ext}}^{\prime }}{\mathrm{\Gamma }c}-t\right],\end{array}$$(A.8)

with the energy density of the external photon field in the black hole frame, $u_{\text{ext}}^{\prime }$, and the radius of the external photon field, $z_{\text{ext}}^{\prime }$. Again, we approximated β_Γ ≈ 1.

Figure A.1 shows the evolution with time of the cooling and escape timescale. The adiabatic and acceleration (not shown) timescales behave exactly like the escape timescale, while the external-Compton cooling simply stops once the blob leaves the respective region. The synchrotron cooling timescale (with b = 1) increases rapidly for t > t_esc(0). The SSC cooling timescale (including the time-dependency of Eq. (7)) initially decreases due to the increase in particle density. This reverses at t ∼ t_esc(0), as the maximum density has been passed and particles overall leave the emission region.

All Tables

All Figures

	Fig. 1. Particle density as a function of time for various cases of ηescα and p as labeled. The dashed lines are the analytical result, Eq. (7). Vertical red line is tesc(0). Dash-double-dotted line marks t−3. Further parameters are q0 = 1 cm−3 s−1, R0 = 5 × 1015 cm, and ηesc = 3. Solid lines are from simulations using OneHaLe.
In the text

Fig. 2. Electron density as a function of time for various cases of ηescα as labeled. Secondary electrons are produced from proton-induced γ rays absorbed by isotropic external photons. The thick dash-dotted line is the analytical solution of Eq. (14), assuming only one isotropic photon field with ni = 1 × 1010 cm−3, and zi = 1 pc. The dashed lines are the standard solution, Eq. (7). The vertical red line is tesc(0), while the horizontal dotted lines mark the maximum level of the standard solution and when they are reached by the dash-dotted lines. Further parameters are p = 2, q0, e = q0, p = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10. The solid lines mark simulations with OneHaLe.

In the text

	Fig. 3. Same as Fig. 2, except that the absorber is AD photons with nAD = 1 × 1011 cm−3, and RAD = 1 × 1016 cm. Further parameters are p = 2, q0, e = q0, p = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.
In the text

	Fig. 4. Same as Fig. 2, but for secondary electrons produced from proton-induced γ rays absorbed by electron-synchrotron photons, Eq. (17). Additional parameters are: p = 2, q0, e = 1 cm−3 s−1, q0, p = 1 × 103 cm−3 s−1, R0 = 5 × 1015 cm, and ηesc = 3.
In the text

	Fig. 5. Same as Fig. 2, but for secondary electrons produced from IC-induced γ rays absorbed by isotropic external photons, Eqs. (20) and (21), with ni = 4 × 107 cm−3, and zi = 1 pc. Further parameters are p = 2, q0, e = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.
In the text

	Fig. 6. Same as Fig. 2, but for secondary electrons produced from IC-induced γ rays absorbed by AD photons with nAD = 1 × 1010 cm−3, and RAD = 1 × 1016 cm. Further parameters are p = 2, q0, e = 1 cm−3 s−1, R0 = 5 × 1015 cm, ηesc = 3, and Γ = 10.
In the text

	Fig. 7. Small parameter study using OneHaLe. The black lines correspond in style to those in Figs. 2–6, but for the case ηescα = 0.5. The solid colored lines show simulations varying the input parameter as labeled. In the bottom panels, the variation in the electron parameters (magenta lines) preserves the initial injection rate of ∼1 cm−3 s−1.
In the text

Fig. 8. Model light curves based upon parameters describing PKS 1510−089 for VHE γ rays (top row), HE γ rays (2nd row), X-rays (3rd row), and the optical R-band (bottom row) for purely leptonic processes (left column) and including hadronic processes (right column). Displayed models are either with (thick dashed) or without (solid) γ − γ pair production. Colors are for various opening angles as labeled. The vertical lines mark tesc(0) (red), the crossing of the BLR edge (blue) and of the DT edge (magenta); all times in the observer’s frame. For clarity, the launch of the blob takes place at tobs = 0.1 d.

In the text

	Fig. 9. Same as Fig. 8, but for the radio band at 15 GHz.
In the text

Fig. 10. Electron luminosity (top row), proton luminosity (second row), radiation luminosity (third row), and total luminosity (bottom row) as a function of time for pure electron (left column) and electron-proton simulations (right column) for different opening angles (color code as labeled) using the parameters of PKS 1510−089 as in Fig. 8. Simulations without pair production are marked by solid lines, while dashed lines refer to simulations with pair production. The vertical lines have the same meaning as in Fig. 8. The black dot-dashed and black dash-double-dotted lines mark the Eddington and AD luminosities, respectively. We note the logarithmic time axis.

In the text

	Fig. A.1. Timescales in arbitrary units as labeled versus time in the comoving frame. The lines are indicative of the time-dependency. The vertical red line marks tesc(0). Free parameters are ηescα = 0.5, ηesc = 3, R0 = 5 × 1015 cm, b = 1, p = 2, Γ = 10, and $z_{\text{ext}}^{\prime }=3\times {10}^{18}$ cm.
In the text