© ESO 2021

1. Introduction

The launching of astrophysical jets from accreting black holes (BHs) involves fundamental gravitomagnetic processes. The frame-dragging effect on the plasma of the inner accretion flow and its advected magnetic field creates a rotating magnetosphere that can lead to a vacuum voltage drop of 10²⁰ eV for a typical spinning supermassive BH in the centres of galaxies (Levinson & Rieger 2011). It is commonly thought that thermal pair production (PP) in photon-photon collisions provides charge carriers in sufficient numbers to supply the Goldreich-Julian density, which allows electric currents to remove the space charges resulting from the electric potential and to reshape the magnetosphere into a force-free configuration. Any imbalance in the supply of charged particles, however, results in the immediate creation of vacuum (electrostatic) gaps. In vacuum gaps the actual charge density, provided by photon-photon PP (see, however, Petropoulou et al. 2019, for an assessment of the importance of triplet PP in limiting the gap growth), is different from the Goldreich-Julian charge density. Because of this deviation, the force-free condition is not satisfied and the electric field component parallel to the magnetic field lines is non-vanishing and able to accelerate intruding charge carriers to relativistic velocities, extracting energy from the central rotating BH (Blandford & Znajek 1977; Beskin et al. 1992; Hirotani & Okamoto 1998; Levinson 2000).

The interaction of a particle beam with a low-energy radiation field results in electromagnetic cascades. Inverse-Compton (IC) scattering of the primary electrons off low-energy photons produces gamma rays that then undergo PP by interacting with the same photon field leading to a secondary generation of electrons. This sequence repeats until the gamma rays can escape freely (Lovelace et al. 1979; Akharonian et al. 1985; Svensson 1987; Zdziarski 1988; Aharonian & Plyasheshnikov 2003).

Cascades have been studied in a variety of astrophysical environments. In outer magnetospheric vacuum gaps of pulsars, curvature radiation is thought to initiate the pair cascade, limiting the gap growth by supplying charged particles (cf. Cheng et al. 1986; Hirotani 2005; Tang et al. 2008). The gamma rays produced by the cascade and leaving the pulsar magnetosphere in a confined beam are made responsible for the observed high-energy (HE) and very high energy (VHE) pulsed emission. In the galactic centre, vacuum gaps are hypothesised to drive cascades that are responsible for PeV cosmic rays and for the VHE gamma-ray point source coincident with Sagittarius A* (cf. e.g., Levinson & Rieger 2011; Katsoulakos et al. 2020).

In active galactic nuclei (AGN) with high accretion rates such as Seyfert galaxies, IC pair cascades in the coronal plasma of the accretion disk are thought to shape their X-ray spectrum up to MeV energies (Svensson 1987; Zdziarski 1988). In AGN with jets such as blazars, unsaturated IC pair cascades can develop from which gamma rays at energies well above MeV can escape. The cascades can be initiated by relativistic particles accelerated in vacuum gaps in the BH magnetosphere when the accretion rate is too low to supply the Goldreich-Julian density but their magnetic field is strong enough to launch a jet (Neronov & Aharonian 2007; Neronov et al. 2009; Levinson & Rieger 2011; Hirotani et al. 2017; Katsoulakos & Rieger 2020; Hirotani 2018, for a review). The gaps are thought to be located at the poles of the spinning BH magnetosphere near to the null surface where the Goldreich-Julian charge density is equal to zero (Beskin et al. 1992; Hirotani & Okamoto 1998; Ptitsyna & Neronov 2016; Hirotani et al. 2016; Levinson & Segev 2017; Ford et al. 2018; Chen & Yuan 2020; Kisaka et al. 2020). It was also found that gaps form around the divide of the magnetohydrodynamic (MHD) flow, called the stagnation surface (Vincent & Lebohec 2010; Broderick & Tchekhovskoy 2015; Hirotani & Pu 2016; Aharonian et al. 2017), or near to the inner light surface of the magnetosphere (Crinquand et al. 2020). Time-dependent one-dimensional, general relativistic MHD simulations by Levinson & Segev (2017), Levinson & Cerutti (2018), Chen et al. (2018), and Chen & Yuan (2020) and two-dimensional simulations by Crinquand et al. (2020, for computational reasons, the latter authors had to restrict themselves to unrealistic parameters, however) have recently found inherently non-stationary gap solutions. Levinson & Cerutti (2018) and Kisaka et al. (2020) indicate that the gap activity relaxes to low-amplitude quasi-steady oscillations after an initial spark event, while Chen et al. (2018) and Chen & Yuan (2020) find cyclically enduring gap activity.

Gap-accelerated particles then interact with ambient photon fields, that is, photons from the accretion flow or from the broad line region (BLR), initiating the onset of a cascade (cf. Blandford & Levinson 1995; Broderick & Tchekhovskoy 2015), which gives rise to an additional component of IC up-scattered photons that are emitted along a beam in the same direction as the original electrons were moving. When these cascaded photons can escape from the AGN gamma sphere without being absorbed by PP, and when the line of sight coincides with the beam, the emission could cause an observable imprint on the spectral energy distribution (SED) of the AGN. Intermittency of the particle flux injected from the gap or of the target photon field results in time variability of the emerging radiation and in transient features of the total SED. This could serve as an explanation of the observed minute- and hour-scale variability of HE and VHE emission in blazars (cf. Aharonian et al. 2007; Albert et al. 2007; Foschini et al. 2013; Pittori et al. 2018; Liao 2018). Similar short-term variability has been found in AGN with jets at inclination angles larger than a few degrees, which supports a similar interpretation (cf. Aleksić et al. 2014; MAGIC Collaboration 2018). Alternative scenarios to explain the short-time variability are the interaction of the jet with the radiation field from a stellar envelope, for instance, or a relativistically propagating emission zone in a relativistically moving jet (cf. Aharonian et al. 2017; Rieger 2019).

In this paper, we examine IC pair cascades in blazars emphasising the interaction of beam electrons from (short-lived) vacuum gaps with recombination-line photons. Blazars are subdivided into two subclasses, flat-spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects (for a recent review, see e.g., Foschini 2017). These two types were originally distinguished by the presence of broad optical emission lines in the case of FSRQs, with equivalent width |W_λ| > 0.5 nm, and the absence of broad optical emission lines for BL Lac objects (Urry & Padovani 1995). However, even in BL Lac objects, in which by definition |W_λ| < 0.5 nm, optical and ultraviolet recombination lines from the BLR are often present (cf. Stocke et al. 2011; Becerra González et al. 2020, for Markarian 501 and S4 0954+65, respectively), and might imprint a diagnostically important signature on the gamma-ray spectra emerging from the innermost regions.

This paper is organised as follows: in Sect. 2 we describe the physical model for the beam-photon interaction scenario and its underlying assumptions, including the set of coupled kinetic equations that were employed to numerically evaluate the model. We outline the numerical procedure for solving them to obtain the steady-state particle and photon distributions in Sect. 3. In Sect. 4 we apply this numerical scheme to fit the observational SED of the BL Lac object Markarian 501 (Mrk 501), which exhibited a peculiar spectral feature at ≈ 3 TeV during a strong flaring activity in 2014 July (MAGIC Collaboration 2020), and discuss the implications of the proposed scenario for the jet formation in Mrk 501. Finally, Sect. 5 summarises our findings.

In the following m_e, c, e, σ_Th, k_B, and h denote the electron rest mass, the velocity of light, the elementary charge, the Thomson cross section, the Boltzmann constant, and Planck’s constant, respectively. SI units are employed unless noted otherwise.

2. Unsaturated inverse-Compton pair cascade model

The electron beam from an active spark gap region close to the accreting BH can initiate IC pair cascades (Broderick & Tchekhovskoy 2015; Ptitsyna & Neronov 2016; Ford et al. 2018; Levinson & Cerutti 2018; Chen & Yuan 2020). When the beam electrons propagate further along the jet axis, they can encounter the intense photon field of an emission-line cloud as depicted in Fig. 1, and initiate IC pair cascades. Because the threshold energy for PP in photon-photon collisions with the optical and ultraviolet photons characteristic for recombination-line spectra is high, the cascade spectra do not saturate at the electron rest mass energy delivering X-ray spectra, but instead produce spectra peaking at VHE gamma rays. Such cascades are also linear because the reservoir of low-energy target photons is not modified by the cascade in a major way. For their SEDs, it is essential to consider an escape term accounting for the energy-dependent gamma sphere resulting from the pair creation optical depth, which decreases towards the edge of the emission region.

Fig. 1. Sketch of the geometry (not to scale) assumed for the model. Seed electrons are accelerated in the gap and multiplied in and behind the gap through ADAF photons. For simplicity of drawing, photons are not shown in the in-gap and post-gap multiplication. After propagating away from the gap, the electron beam meets an emission-line photon field, driving a cascade, whose steady-state gamma-ray photon density is nγ and whose escaping flux is Fcasc. To convert nγ into Fcasc, cf. Eq. (3).

2.1. Unsaturated inverse-Compton pair cascade with escape

For the model of the emission region, we treated electrons and positrons identically and assumed three reservoirs of particle species to be present: firstly, we have relativistic electrons with spectral number density N(γ), where γ denotes the Lorentz factor of the electron. Secondly, we have highly energetic photons with spectral number density n_γ(x_γ), where x_γ denotes the photon energy in units of the electron rest mass energy m_ec². In what follows, the term “highly energetic photon” does not only denote photons with energies in the established HE range from 100 MeV to 100 GeV, but more generally photons with energies well above m_ec², in contrast to low-energy photons with energies well below m_ec². Thirdly, we have external low-energy photons with spectral number density n₀(x) with x₁ ≤ x ≤ x₀ ≪ 1.

We assumed that all number densities are time-independent, isotropic, and homogeneous throughout the interaction region, and that interactions take place only inside of it. The interactions were assumed to be dominated by IC scattering and photon-photon PP. The cascade is then driven by PP from collisions in which the low-energy target photons turn the highly energetic photons into electron-positron pairs that then IC scatter off the same low-energy target photons to become the next generation of highly energetic photons, and so on (cf. Fig. 2). We considered four additional mechanisms acting as sinks or sources for the highly energetic photons or the electrons: first, there is the injection of highly energetic photons with spectral rate ṅ_γ,i(x_γ) and x_γ ≫ 1. Additionally, we demanded that the threshold condition for PP in collisions with the isotropic target radiation field with x_γ x >  1 was fulfilled. Second, there is the injection of electrons with spectral rate Ṅ_i(γ). Henceforth, ṅ_γ,i and Ṅ_i are called spectral injection rate, and are defined as the number of particles, injected per unit space volume interval, per unit time interval, and per unit energy interval. We required that Ṅ_i is non-vanishing only for γ ≫ 1, hence the electrons have to be highly relativistic. Furthermore, we required that γ x >  1 was satisfied for the injected electrons. This is equivalent with primary IC scattering events being in the Klein-Nishina (KN) regime. According to Zdziarski (1988), this type of cascade is denoted IC-Klein-Nishina pair cascade. Third, there is escape of highly energetic photons from the interaction region with timescale T_ph esc(x_γ). Fourth, there is electron escape with timescale T_e esc(γ). This is the main difference of our treatment to that by Zdziarski (1988) and Wendel et al. (2017), who considered saturated cascades; it makes our scenario more realistic. A finite escape timescale is equivalent to a bounded interaction region and is thus closer to reality than an infinite escape timescale, which is equivalent to an infinitely large interaction volume. The cascades considered here are therefore not necessarily saturated.

Fig. 2. Sketch of the considered linear IC pair cascade with escape terms.

We assumed that the reservoir of target photons was not modified by the cascade. In other words, n₀(x) is kept constant. This assumption is equivalent to assuming that a cascade is linear (Zdziarski 1988). Linearity of a cascade means that the low-energy photons alone serve as targets. Self-interactions of the highly energetic photons and of the electrons are neglected. Zdziarski (1988) showed by comparing cross sections that when the energy density of low-energy photons is at least similar to that of the highly energetic photons, the condition for a cascade to be linear is met whose interactions occur in the KN regime. Repetition of IC scattering events successively decreases the electron energy, until IC scattering events take place in the Thomson regime (γ x <  1). Here, cascades are always non-linear because x_γ x <  1 implies that PP is not possible and thus PP is only possible through highly energetic photon self-interactions.

2.2. Description through kinetic equations

To determine the change rates of the spectral number densities of photons and electrons, we used the spectral interaction rates of IC scattering and of PP given by Zdziarski (1988).

For IC events with incident electrons of energy γ ≫ 1 scattering off photons with energy x ≪ 1 and spectral number density n₀(x), which result in down-scattered electrons with energy γ′ and highly energetic photons with energy x_γ, the spectral IC scattering interaction rate is denoted by C(γ, γ′) and determined with Eq. (A.1) by Zdziarski (1988), cf. Appendix A. For such events, the threshold γ_IC,  th(x_γ, x) as defined in Appendix A is the minimum required value of γ. The maximum possible γ is γ_IC,  max(γ′,x), as outlined in Appendix A. Furthermore, we define ${\gamma }_{\text{IC}, \text{min}}^{\prime }(\gamma ,x)$ as the minimum possible γ′ and x_γ,  max(γ, x) as the maximum energy of the up-scattered photon.

Because of the interaction of highly energetic photons of energy x_γ with low-energy photons of energy x and spectral number density n₀(x), resulting in the PP of electrons of energy γ, the spectral PP interaction rate is called P(x_γ, γ) and is computed with Eq. (B.1) of Zdziarski (1988), cf. Appendix A. The PP threshold, that is, the minimum required x_γ, is defined as x_{γ,  PP,  th}(γ, x) in Appendix A. The values of γ that can be reached are limited by γ_PP,  min(x_γ, x) and γ_PP,  max(x_γ, x), as described in Appendix A.

The kinetic equation of the electrons and of the photons ensues after all the relevant sinks and sources are specified. The change rate of the spectral number density of the electrons (electron energy distribution) is affected by the electron spectral injection rate Ṅ_i, the spectral loss rate (as the sum of the escape rate with timescale T_e esc and of the IC down-scattering rate to lower energies), and by the production rate (as the sum of the IC down-scattering rate from higher energies and of the PP rate). Similarly, the change rate of the spectral number density of the photons is determined from the photon spectral injection rate ṅ_γ,i, the loss rate (which is the sum of the escape rate with timescale T_ph esc and of the pair absorption rate), and the production rate ṅ_γ,IC.

Assuming steady state, we obtain, as shown in Appendix A,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill n_{\gamma }(x_{\gamma })=& \underset{:=n_{\gamma ,\mathrm{i}}(x_{\gamma })}{\underbrace{\frac{{\dot{n}}_{\gamma ,\mathrm{i}}(x_{\gamma })}{\frac{1}{T_{\mathrm{ph}\mathrm{esc}}(x_{\gamma })}+{\int }_{{\gamma }_{\mathrm{PP},\min }(x_{\gamma },x_0)}^{{\gamma }_{\mathrm{PP},\max }(x_{\gamma },x_0)}P(x_{\gamma },\gamma )\mathrm{d}\gamma }}}\hfill \\ \hfill & +\underset{:=n_{\gamma ,\mathrm{IC}}(x_{\gamma })}{\underbrace{\frac{{\int }_{{\gamma }_{\mathrm{IC},\mathrm{th}}(x_{\gamma },x_0)}^{\infty }N(\gamma )C(\gamma ,\gamma -x_{\gamma })\mathrm{d}\gamma }{\frac{1}{T_{\mathrm{ph}\mathrm{esc}}(x_{\gamma })}+{\int }_{{\gamma }_{\mathrm{PP},\min }(x_{\gamma },x_0)}^{{\gamma }_{\mathrm{PP},\max }(x_{\gamma },x_0)}P(x_{\gamma },\gamma )\mathrm{d}\gamma }}},\hfill \end{array}\end{array}$$(1)

where we have defined two contributions to the photon spectral number density. n_γ,  i describes the contribution that arises from the interplay of injection and escape and absorption losses, while n_γ,  IC describes the contribution that is achieved from the interplay of production through IC up-scattering and both losses. Exchanging γ and γ′ in Eq. (1) and plugging this into Eq. (A.8), we can eliminate n_γ,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill N(\gamma )=& \frac{{\dot{N}}_{\mathrm{i}}(\gamma )+{\int }_{\gamma }^{\gamma {\prime }_{\mathrm{IC},\max }(\gamma ,x_0)}N(\gamma \prime )C(\gamma \prime ,\gamma )\mathrm{d}\gamma \prime +{\int }_{x_{\gamma ,\mathrm{PP},\mathrm{th}}(\gamma ,x_0)}^{\infty }({\dot{n}}_{\gamma ,\mathrm{i}}(x_{\gamma })}{\frac{1}{T_{\mathrm{e}\mathrm{esc}}(\gamma )}}\hfill \\ \hfill & \frac{+{\int }_{\gamma {\prime }_{\mathrm{IC},\mathrm{th}}(x_{\gamma },x_0)}^{\infty }N(\gamma \prime )C(\gamma \prime ,\gamma \prime -x_{\gamma })\mathrm{d}\gamma \prime )\frac{2P(x_{\gamma },\gamma )}{\frac{1}{T_{\mathrm{ph}\mathrm{esc}}(x_{\gamma })}+{\int }_{\gamma {\prime }_{\mathrm{PP},\min }(x_{\gamma },x_0)}^{\gamma {\prime }_{\mathrm{PP},\max }(x_{\gamma },x_0)}P(x_{\gamma },\gamma \prime )\mathrm{d}\gamma \prime }\mathrm{d}{x}_{\gamma }}{+{\int }_{\gamma {\prime }_{\mathrm{IC},\min }(\gamma ,x_0)}^{\gamma }C(\gamma ,\gamma \prime )\mathrm{d}\gamma \prime }\hfill \end{array}\end{array}$$(2)

3. Numerical solution procedure

In this section, we outline explicit solutions for N(γ) and n_γ(x_γ).

The main task is to find solutions for N. This is done by solving Eq. (2) iteratively. The right-hand side of Eq. (2) directly incorporates γ, T_ph esc, T_e esc, N, Ṅ_i, ṅ_γ,i, C, P, and the functions in the integration borders, which have been defined in the previous section. C and P are known as soon as n₀ has been specified. Hence, Eq. (2) can briefly be written as N(γ) = ℱ′(n₀, Ṅ_i, ṅ_γ,i, T_{ph esc}, T_{e esc}, N, γ). As described in Sect. 2.1, the first five arguments of ℱ′ are quantities that describe the physical setting. Hence, they should be understood as input functions and have to be prescribed to specify the physical problem. When this is done, we obtain the equation N(γ) = ℱ(N, γ), where ℱ was defined accordingly. Now, we chose a sequence (γ_k)_{k = 1, …, κ} of κ discrete values γ_k along which the kinetic equation is to be solved, in other words, along which N(γ) is to be determined. We write N(γ_k) = N_k for convenience. As a starting point of the iteration, we guessed an initial course N_init(γ) of the function N(γ) and determined the corresponding sequence of values (N_k)_{k = 1, …, κ;  jinit}. In Appendix B, we discuss how the converged sequence is determined numerically.

We performed the iteration using the python3 language. To numerically compute the integrals with an infinite integration range, we restricted the integration range we used to a finite range. This is possible when n₀(x), ṅ_γ,i(x_γ), and Ṅ_i(γ) are equal to zero everywhere above an upper cut-off, which we denote by x₀, x_γ,  0, and by γ_i,  0, respectively.

To summarise, the physical setting is specified as soon as the following functions are defined: The low-energy photon spectral number density n₀, the highly energetic photon spectral injection rate ṅ_γ,i, the electron spectral injection rate Ṅ_i, the highly energetic photon escape time T_ph esc, and the electron escape time T_e esc. By iterating, our code then determined the steady state N. When N was known, we used Eq. (1) to determine n_γ, and from this, we obtained the spectral flux density F_casc. This approach differs from that of Wendel et al. (2017). Therein, although T_ph esc = T_e esc = ∞ (no escape) was used in the kinetic equation, it was assumed that the IC up-scattered photons leave the interaction region from its shell-like boundary. Accordingly, Wendel et al. (2017) directly converted the photon spectral production rate (last term in Eq. (A.7)) into the observed spectral flux density. This appears to be a very crude approach because first it neglects that the injected highly energetic photons can also leave the interaction region (which we take into account by n_γ,  i in Eq. (1)), and second it assumes that all IC up-scattered photons escape from the boundary of the interaction region and that no IC up-scattered photons escape from its interior. On the one hand, this second assumption contradicts the T_ph esc = T_e esc = ∞ assumption and on the other hand it neglects the energy-dependent behaviour of photon attenuation (which we took into account through the division of the photon spectral injection and production rate (numerator of Eq. (1)) by the spectral loss rate (denominator of Eq. (1))).

Based on n_γ, we determined the spectral flux density F_casc of highly energetic photons (i.e. the number of photons per unit energy, per unit area, and per unit time) through

$$\begin{array}{c}\hfill F_{\mathrm{casc}}(x_{\gamma })=n_{\gamma }(x_{\gamma })·\frac{4\pi c{R}_{\mathrm{ph}\mathrm{esc}}{(x_{\gamma })}^2}{\mathrm{\Omega }(\varphi )D^2{m}_{\mathrm{e}}{c}^2}·\end{array}$$(3)

We assumed that the photons leave the interaction region, whose approximate radial size is R_ph esc(x_γ), only through a conical beam of opening angle ϕ and solid angle Ω = 4π ⋅ sin²(ϕ/4), and that D is the distance to the observer. The size of the interaction region is linked with the escape time of the photons and electrons. A natural choice for this would be R_ph esc(x_γ) = T_ph esc(x_γ) c. The division by m_ec² arises because F_casc is measured in m⁻² s⁻¹ J⁻¹, while n_γ is in units of m⁻³.

4. Indication of gap activity in Markarian 501

The numerical procedure described above is used below to explain a peculiar feature in the VHE spectrum of the blazar Mrk 501 during a flaring activity on 2014 July 19.

The broadband emission from blazars, and in particular from the high-synchrotron-peaked BL Lac object Mrk 501, is characterised by a two-bump structure. This emission is usually described within the one-zone synchrotron-self-Compton (SSC) framework (see e.g., Maraschi et al. 1992; Dermer & Giebels 2016). In this model, a blob of relativistic electrons is assumed to move along the jet that points towards the observer. The blob is permeated by a randomly distributed magnetic field, producing synchrotron radiation from the electrons. This radiation is responsible for the low-energy SED bump in blazars. The origin of the HE bump is typically ascribed to the IC scattering of the relativistic electrons with the synchrotron photons. In addition to the SSC scenario, other theoretical models have been proposed to explain the HE and VHE emission in blazars, involving hadronic components and IC pair cascading. Minute-scale variability spectral components mean that the standard one-zone SSC scenario may be insufficient and suggest that an additional sporadic emission component is emitted from the innermost regions of the jet, see Ahnen et al. (2018).

Recently, the MAGIC Collaboration (2020) has presented a two-week flare from Mrk 501, mainly detected in X-rays with the Swift satellite and in VHE gamma rays observed with the major atmospheric gamma-ray imaging Cherenkov (MAGIC) telescopes. During this flaring state, the source displayed hard X-ray spectra as well as the highest X-ray count rate detected by the Swift X-ray Telescope during its 15 years of operation. In coincidence with the peak of the X-ray count rate, which took place on 2014 July 19 (MJD 56857.98), a narrow feature in the TeV spectrum from the MAGIC telescopes was detected with a significance of 3σ − 4σ. This feature, located around 3 TeV, is inconsistent with the standard functions to describe the VHE spectra from blazars (power law, log-parabola, and log-parabola with exponential cut-off) at more than 3σ confidence level. The VHE spectrum on 2014 July 19 is better explained assuming a narrow component (a narrow log-parabola or a Gaussian function) in addition to a broad log-parabola. Such a double function fit is preferred with respect to a single log-parabola at 3.2σ − 4.5σ, depending on the a priori assumptions (details can be found in MAGIC Collaboration 2020). This was the first time that a narrow feature (not consistent with the usual smooth spectral functions) was found with a significance higher than 3σ in the VHE spectrum of Mrk 501, and any TeV blazar in general. Indications for similar (but statistically insignificant) spectral features at ≈ 2 TeV can be found by inspection of flare spectra reported by the very energetic radiation imaging telescope array system (VERITAS; cf. Fig. 8 of Abdo et al. 2011) and the MAGIC collaboration (cf. Fig. 10 of Ahnen et al. 2017).

Under the assumption that the spectral feature described above is real, three different theoretical scenarios were proposed by the MAGIC Collaboration (2020) to reproduce the narrow spectral feature: (a) a pileup in the electron energy distribution produced by stochastic acceleration, (b) a two-zone SSC model, and (c) a narrow emission from an IC pair cascade scenario induced by electrons accelerated in a magnetospheric vacuum gap. Details for the model describing scenario (c) are provided in this work. In the following, we evaluate the model we fitted to the observational data to infer its implications and consequences.

4.1. Specification of the modelled setting

We assumed a vacuum gap near to the poles of the Mrk 501 BH magnetosphere. The potential drop in the gap accelerates seed electrons to ultra-relativistic energies. The seed electrons enter the gap both through PP of gamma rays from the hot advection-dominated accretion flow (ADAF), and possibly through direct leakage from it (Neronov et al. 2009). Seed electrons that have been accelerated in the gap, leave it in the direction of the electric field (i.e. in the direction of the magnetic field), and propagate along a beam that is well collimated to an opening angle ≈ ϕ. During the propagation away from the gap, the seed electrons still encounter the ADAF photon field, and a post-gap cascade (Broderick & Tchekhovskoy 2015; Chen & Yuan 2020) occurs and increases the number of electrons by repeated curvature-radiation and IC emission and PP events. This post-gap cascade ceases when curvature-radiation emission becomes negligible and when IC scattering and PP can no longer be sustained by ADAF photons. Several tens or hundreds of Schwarzschild radii away from the gap, the beam of electrons is assumed to enter a region in which background photons from ionised gas clouds are present. This induces a cascade in this interaction region. This model is depicted in Fig. 1. For this cascade, the model of Sect. 2 was applied. The electron beam from the gap serves as electron injection mechanism in our model. We describe the electron spectral injection rate with the following cut-off Gaussian:

$$\begin{array}{c}\hfill {\dot{N}}_{\mathrm{i}}(\gamma )=\{\begin{array}{cc}\frac{K_{\mathrm{G}}}{\varsigma \sqrt{(}2\pi )}·exp(-\frac{{(\gamma -{\gamma }_{\mathrm{mean}})}^2}{2{\varsigma }^2})\hfill & \mathrm{if}{\gamma }_{\mathrm{i},1}\le \gamma \le {\gamma }_{\mathrm{i},0},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(4)

Here, γ_mean denotes the mean Lorentz factor of the injected electrons, and ς parametrises the width of the distribution. For the cut-off values we chose γ_i,  1 = γ_mean − 3.0 ς and γ_i,  0 = γ_mean + 3.0 ς, satisfying the condition γ ⋅ x >  1 (cf. Sect. 2.1). The normalisation of the Gaussian, which describes the total number of electrons that are injected per unit time and per unit space volume, is called K_G.

Furthermore, we assumed that the external low-energy photons of the cascade are represented by emission-line photons from photo-ionised gas clouds in the surroundings of the Mrk 501 BH magnetosphere. Although Mrk 501, as a BL Lac object, does not have a prominent BLR, its BLR has clearly been confirmed observationally, possibly due to outshining by the boosted jet emission, due to a weak accretion rate, or due to the lack of ambient gas (cf. Moles et al. 1987; Stocke et al. 2011). After they are triggered by minor galaxy mergers, interstellar gas clouds, consisting mainly of hydrogen and helium (Wilms et al. 2000), can migrate from the host galaxy into its centre and thus into the central part of Mrk 501. Radiation from the accretion flow or from OB-type stars in the dense galactic centre will inevitably ionise passing clouds. As an effect of recombination, emission-line photons will be abundant in the proximity of the clouds. The external low-energetic photon field can also stem from the envelope of a post-main-sequence star with low metallicity (cf. e.g., Bednarek & Protheroe 1997; Barkov et al. 2010). In any case, we assumed that the interaction region was both penetrated by the electron beam and pervaded by emission-line photons, and consequently, that an IC pair cascade developed there. It is natural to describe the spectral number density by a sum of Dirac delta distributions,

$$\begin{array}{c}\hfill n_0(x)=K_{\mathrm{lines}}·{\displaystyle \sum _{i=1}^4\frac{K_{\mathrm{line},i}}{x_{0,i}}·{\delta }_{\mathrm{Dirac}}(x-x_{0,i})}\end{array}$$(5)

Here, x_0,  i = h/(λ_0,  i m_e c), and λ_0,  i are the dimensionless energy, and the wavelength of the ith line, respectively. K_line,  i gives the relative flux density contribution of the ith line. Precisely, K_line,  i gives the flux density of the emission line i with respect to the flux density of a hypothetical hydrogen Balmer-β line. Accordingly, K_line,  i/x_0,  i describes the relative number density contribution of the ith line. The parameter K_lines is a measure for the photon total number density. We restricted the sum to the three typically most important hydrogen lines plus the most important helium line. We compiled these lines and their respective K_line,  i based on the ultraviolet spectroscopic blazar observations by Pian et al. (2005) and on the synthetic photo-ionisation spectra by Abolmasov & Poutanen (2017). We list them in Table 1.

From an energetic point of view, the emission-line radiation is just a reprocessed fraction of the more luminous ADAF radiation. Nevertheless, we did not include the ADAF radiation in the target photon field. This is justified by the following geometrical argument. The emission of line photons in the ionised cloud is approximately isotropic. The interaction angles of injected electrons and photons are therefore equally distributed; the interaction angle is on average 90°. This also justifies the usage of the interaction rates C and P. In contrast, for interactions of electrons with ADAF photons the collision angle would be at best 90°, and on average smaller than 90°. Thus, glancing collisions occur, and for such interactions, the interaction rates are lower than for right-angled collisions. This results in increased interaction rates for the line photons in comparison to ADAF photons, and thus we can neglect the latter. Moreover, except for the reprocessed photons from the ionised cloud, reprocessed ADAF photons might be scarce in the surroundings of the beam due to low ambient gas density in the BL Lac object, which means that they drop out as target photons as well. We have to mention furthermore, that ADAF spectra usually extend to higher energies and x ≪ 1, which was assumed in the derivation of C and P, is satisfied only very critically.

We did not have injection of highly energetic photons in this setting. Therefore we set the photon spectral injection rate

$$\begin{array}{c}\hfill {\dot{n}}_{\gamma ,\mathrm{i}}(x_{\gamma })=0.\end{array}$$(6)

To choose the escape timescales, we determined

$$\begin{array}{c}\hfill T_{\mathrm{ph}\mathrm{esc}}(x_{\gamma })=T_{\mathrm{e}\mathrm{esc}}(\gamma )=\frac{R}{c}:=T_{\mathrm{esc}}.\end{array}$$(7)

In other words, we approximated the interaction volume as a spherical region with radius R. Electrons and photons have equal, and energy independent, that is, constant, escape timescales. The function R_ph esc(x_γ) in Eq. (3) is therefore also set equal to the parameter R.

4.2. Modelling results

According to the procedure outlined in Sect. 3, we determined N(γ), n_γ(x_γ), and F_casc(x_γ) of the cascade that evolves through interaction of the electrons originally stemming from the gap with the emission-line photons, such that F_casc(x_γ) was fitted to the narrow SED feature (MAGIC telescopes’ data) and a background SSC emission is fitted to the multi-wavelength SED. We understand K_G, K_lines, γ_mean, ς, R, and ϕ as fitting parameters, while we used the luminosity distance D = 149.4 Mpc (corresponding to redshift z = 0.034, dimensionless density parameter Ω_m = 0.3 of matter, dimensionless density parameter Ω_Λ = 0.7 of dark energy, Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹) and the line contributions K_line,  i given in Table 1. Because the electron beam penetrates into the interaction region from one direction, the kinematics of IC scattering and PP ensure that the highly energetic photons leave the interaction region in this same direction, hence along a beam of opening angle ϕ.

Interim results of the numerical procedure are shown in Figs. 3–5 and are discussed now.

Fig. 3. γ N(γ), which is a proxy for the electron spectral energy density, versus the product of the Lorentz factor with the dimensionless energy of the highest energetic line, the helium II Lyman-α line. The iterated points Nk,  j are drawn in red.

The results of our iterative determination of N(γ) are shown in Fig. 3. The initial function of N was chosen to be

$$\begin{array}{c}\hfill N_{\mathrm{init}}(\gamma )=\frac{2·{\dot{N}}_{\mathrm{i}}(\gamma )}{{\int }_{\gamma {\prime }_{\mathrm{IC},\min }({\gamma }_{\mathrm{i},1},x_0)}^{{\gamma }_{\mathrm{i},1}}C({\gamma }_{\mathrm{i},1},\gamma \prime )\mathrm{d}\gamma \prime +{\int }_{\gamma {\prime }_{\mathrm{IC},\min }({\gamma }_{\mathrm{i},0},x_0)}^{{\gamma }_{\mathrm{i},0}}C({\gamma }_{\mathrm{i},0},\gamma \prime )\mathrm{d}\gamma \prime },\end{array}$$(8)

cf. Eq. (A.6). This choice is a compromise between a function that can be evaluated fast and is accurate. It is a good approximation for the final N in the highest energetic regime because the electron spectral production rate is negligible here. For decreasing γ, a higher number of iteration steps are necessary because the energy transfer in one IC scattering event decreases. We chose the convergence criterion as follows: One iteration was considered completed when the relative change |N_{k,  j − 1}−N_k,  j|/N_k,  j between two successive values was smaller than 0.001 at the points ${\gamma }_k>10x_{0,1}^{-1}$ and smaller than 0.001 ⋅ (γ_kx_0,  1/10)² at the points ${\gamma }_k\le 10x_{0,1}^{-1}$. In other words, the demanded accuracy increases quadratically the deeper we iterate into the Thomson regime. There, about j_final = 100 iteration steps were necessary to meet the convergence criterion, while only two or three steps were necessary in the KN regime, cf. Fig. 3. The shallow ridge of N at $\gamma \approx 2x_{0,1}^{-1}$ is caused by PP of highly energetic photons with hydrogen Lyman photons (mainly hydrogen Lyman-α photons), while the very shallow ridge at $\gamma \approx 0.5x_{0,1}^{-1}$ is caused by PP of highly energetic photons with the helium II Lyman-α photons. In other words, these ridges in N are the effect of peaks of the spectral PP rate (fourth summand in Eq. (A.6)). In the Thomson regime, with decreasing γ the course of N deviates more and more from the ∼ γ⁻² dependence of the saturated cascade case considered by Zdziarski (1988) because particles escape in our scenario. The finite value of the escape time prevents the electrons from penetrating deep into the Thomson regime.

Figure 4 clearly shows that the photon spectral loss rate is dominated by pair absorption in the approximate range $3x_{0,1}^{-1}<x_{\gamma }<800x_{0,1}^{-1}$ corresponding to a photon energy between 19 GeV and 5.1 TeV, and by escape outside this energy range. The maximum of the photon spectral loss rate occurs at $x_{\gamma }\approx 12x_{0,1}^{-1}$, corresponding to the maximum of the hydrogen Lyman-α absorption. The optical depths, which correspond to the respective spectral loss rate contributions, are shown in Fig. 4 on the right-hand side ordinate and were determined by multiplication of the respective spectral loss rate contribution with the escape time T_esc.

Fig. 4. Various contributions (cf. Eq. (A.7)) to the photon spectral loss rate on the left-hand side ordinate in dependence on the product of the dimensionless highly energetic photon energy with the dimensionless energy of the highest energetic line. On the right-hand side ordinate, we plot the optical depth that corresponds to the respective contribution to the spectral loss rate.

We show the photon spectral production rate ṅ_γ,IC in Fig. 5. From N, it inherits the peak around $500x_{0,1}^{-1}$ as well as the deviation from a power law in the Thomson regime due to escape (cf. Zdziarski 1988). In this figure we also show the spectral number density n_γ of highly energetic photons. The contribution n_γ,  i vanishes because of the choice Eq. (6). The contribution n_γ,  IC is non-vanishing and ensues as a consequence of IC scattering, pair absorption, and escape. The pronounced trough above $x_{\gamma }\approx 4x_{0,1}^{-1}$ is due to absorption by PP with photons from the hydrogen Lyman-α line (cf. Fig. 4). Because of the smallness of K_line,  2 and K_line,  3, absorption features of the hydrogen Lyman-β line and series are almost invisible. The shallow dip of the spectral number density above $x_{\gamma }\approx x_{0,1}^{-1}$ is due to absorption on the helium II Lyman-α line. Thus, the slight bump at $x_{\gamma }\approx x_{0,1}^{-1}$ results from cascaded photons that are immediately below the threshold for PP and hence are not absorbed. The peaky feature around $500x_{0,1}^{-1}$ is more pronounced in n_γ than in ṅ_γ,IC because the spectral loss rate (i.e. the denominator of n_γ,  IC) is higher below $500x_{0,1}^{-1}$ than at $500x_{0,1}^{-1}$. If n₀ were continuous and distributed along a wider range of x instead of consisting of a few strong lines, the peak around $500x_{0,1}^{-1}$ would be wider and extend to lower energies.

Fig. 5. Spectral production rate and number density of the cascade in Mrk 501. The product of the squared dimensionless photon energy with the spectral production rate of photons is plotted in red (which is determined by the last term of Eq. (A.7)). We show the product of the squared dimensionless photon energy with the spectral number density of photons in green (which is determined by Eq. (1)). Both quantities are plotted vs. the product of the highly energetic photon energy with the energy of the helium II Lyman-α line.

The total energy density of highly energetic photons is obtained by integrating the spectral energy density x_γ n_γ(x_γ)

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{x_{\gamma ,\max }({\gamma }_{\mathrm{i},0},x_{0,1})}{x}_{\gamma }{n}_{\gamma }(x_{\gamma })\mathrm{d}{x}_{\gamma }·m_{\mathrm{e}}{c}^2\approx 4.5\mathrm{J}{\mathrm{m}}^{-3}.}\end{array}$$(9)

The total energy density of the low-energy photons is $K_{\mathrm{lines}}·\sum _{i=1}^4{K}_{\mathrm{line},i}·m_{\mathrm{e}}{c}^2\approx 6.5$ J m⁻³ and thus of similar magnitude. Our cascade is therefore indeed linear (cf. Sect. 2.1).

The cascaded spectral flux density F_casc of highly energetic photons adds to the spectral flux density F_SSC of photons that is produced farther downstream the jet in a common one-zone SSC scenario. Both these components contribute to the observed spectral flux density F. We fit both components, F_SSC and F_casc, to the multi-wavelength data of Mrk 501 of the night of MJD 56857.98 as presented by the MAGIC Collaboration (2020). As shown in Fig. 10 of MAGIC Collaboration (2020) and in Fig. 6 in the present paper, we can describe the narrow SED feature without the need of extreme SSC parameters. The parameters found by the fit are given for the cascaded emission component in Table 2 of the present paper and for the SSC emission component in Table 6 and Fig. 7 in MAGIC Collaboration (2020). The emission from a cascade initiated by electrons originating from a vacuum gap can thus explain the observed narrow peak-like feature.

Fig. 6. High-energy bump of the SED of Mrk 501 from 2014 July 19 (MJD 56857.98). The dashed grey line depicts the IC bump of a one-zone SSC model (contribution FSSC). The dot-dashed red line depicts the emission of the cascade (contribution Fcasc, cf. green line in Fig. 5). The sum F of both contributions is depicted by the solid black line. The spectral data from the MAGIC telescopes are shown as red circles, while measurements and upper limits by the Fermi Large Area Telescope are drawn by black and yellow triangles. The details on the data analysis and the SSC modelling can be found in MAGIC Collaboration (2020).

4.3. Physical inferences

In this subsection we discuss the results of modelling Mrk 501. We assess the findings in a broader context.

4.3.1. Number and size of emission-line clouds

We compared the assumed hydrogen Lyman-α photon abundance with the hydrogen Lyman-α luminosity of L_{Ly α,  obs} = 5.2 × 10³³ W that was detected from Mrk 501 by Stocke et al. (2011). The total number density of all low-energy photons is given by $K_{\mathrm{lines}}·\sum _{i=1}^4{K}_{\mathrm{line},i}/x_{0,i}\approx 3.2\times {10}^{18}$ m⁻³. This number density describes the photon field of the ionised cloud assumed to be interacting with the electron beam from the gap. The total number density of hydrogen Lyman-α photons is determined by K_lines ⋅ K_line,  4/x_0,  4 ≈ 2.6 × 10¹⁸ m⁻³. We further assumed that there are N_cl clouds of approximate radial size R_cl in the BLR. The hydrogen Lyman-α luminosity L_{Ly α,  cl} of all these clouds is responsible for the observed hydrogen Lyman-α luminosity. We can set

$$\begin{array}{c}\hfill L_{\mathrm{Ly}\alpha ,\mathrm{obs}}\stackrel{!}{=}{L}_{\mathrm{Ly}\alpha ,\mathrm{cl}}=N_{\mathrm{cl}}·4\pi R_{\mathrm{cl}}^2{K}_{\mathrm{lines}}\frac{K_{\mathrm{line},4}}{x_{0,4}}c{x}_{0,4}{m}_{\mathrm{e}}{c}^2.\end{array}$$(10)

As BL Lac objects are considered to be evolved objects running short of cold gas supplies (Liodakis 2018; Foschini 2017), we can assume that N_cl is quite small. On the other hand, to obtain a Gaussian line profile for the emission lines resulting from the Doppler broadening, the number of clouds must not be too small. Then, solving Eq. (10) for the radial size yields an estimate of R_cl. Choosing N_cl = 10 yields R_cl ≈ 1.8 × 10¹¹ m, which roughly agrees with our value used for the radial size of the interaction region R (cf. Table 2). This size of the cloud corresponds to the size of a red giant, which is suggestive of cloud formation by photospheric abrasion through beam interaction. In terms of Schwarzschild radii r_S for a central BH of mass M = 10⁹ M_⊙, it is R_cl = 0.061 r_S. From the modelling we infer that in Mrk 501 only few major gas clouds are present in the BLR. The encounter of a major cloud with electrons from the gap-beam must therefore be a rare event. Presuming that the gas clouds revolve around the BH on unperturbed Keplerian orbits, an encounter of a cloud with the beam of electrons will at best repeat with the Keplerian orbital period of the cloud. These conclusions are in line with the fact that the peak-like SED feature has not been found elsewhere in Mrk 501 data.

4.3.2. Reprocessing of the accretion flow luminosity

The total line luminosity of these N_cl clouds of radius R_cl is determined as

$$\begin{array}{c}\hfill L_{\mathrm{lines},\mathrm{tot}}=N_{\mathrm{cl}}·4\pi R_{\mathrm{cl}}^2{K}_{\mathrm{lines}}·{\displaystyle \sum _{i=1}^4{K}_{\mathrm{line},i}c{m}_{\mathrm{e}}{c}^2.}\end{array}$$(11)

For the above chosen N_cl = 10 and resulting R_cl = 1.8 × 10¹¹ m, it is L_{lines,  tot} ≈ 7.8 × 10³³ W, a factor 5 lower than the upper limit of the Mrk 501 BLR luminosity found by Celotti et al. (1997). We recall that this finding was obtained from our modelling and from the observed hydrogen Lyman-α luminosity. In canonical BLR models, the BLR luminosity is assumed to be about 10% of the accretion flow luminosity (Ghisellini et al. 2017). As Mrk 501 is a gas-depleted object, its BLR reprocesses even less of the incident accretion flow radiation. Thus, the reprocessing fraction ξ = L_{lines,  tot}/L_{ADAF,  tot} might be lower than 10%, say 1%. Here, L_{ADAF,  tot} is the total luminosity of the ADAF. To compare the emission-line luminosity with the accretion flow luminosity, we applied the synthetic ADAF spectrum computed by Mahadevan (1997). There, the spectral luminosity of an ADAF as a function of frequency was parametrised by the viscosity parameter α, the pressure ratio β, the inner boundary R_in of the hot accretion flow, its outer boundary R_out, the central BH mass M, the electron temperature T_e, and by the dimensionless mass accretion rate ṁ, which is defined as the ratio of the mass accretion rate to the Eddington accretion rate. We integrated numerically, and for the purpose of cross checking also analytically, along the photon energy ϵ over the spectral luminosity L_ADAF(ϵ) to obtain L_{ADAF,  tot}. We integrated essentially over the superposition of the cyclosynchrotron emission, its Comptonised part, and of the bremsstrahlung contribution. We used the canonical choices α = 0.3, β = 0.5, R_in = 3 r_S, and R_out = 1000 r_S as well as M = 10⁹ M_⊙. Furthermore, we chose the energy ϵ_ADAF,  0 = 300 k_B T_e as an upper integration border because the spectrum essentially vanishes above this energy. We let ṁ and T_e be free parameters. The chosen requirement L_lines,tot ≈ 0.01 · L_ADAF,tot(ṁ, T_e) can be met with the exemplary ṁ, T_e) pairs listed in row 1 of Table 3. This shows that the present scenario fits the general picture of accretion flows as illuminators of BLRs. Moreover, it gives evidence that Mrk 501 is indeed fed by a low-rate, hot ADAF. When the reprocessing fraction ξ is increased from 1% to the canonical 10%, the values of T_e decrease by about 25% (cf. Table 3).

4.3.3. Constraints on the in-gap and post-gap multiplication

The parameter K_G expresses the number density of electrons that intrude into the interaction region per time interval. In our model, these electrons enter the interaction region after their number was increased in the gap and during the post-gap cascade. The origin of the electron beam along which the multiplication takes place is the vacuum gap (cf. Fig. 1), where in turn the origin of the seed electrons is PP through self-interaction of ADAF radiation. In other words, photons that are emitted by the ADAF can collide in the BH magnetosphere and pair-produce seed electrons. We can determine the number density of ADAF photons capable of PP with the help of the ADAF model by Mahadevan (1997). The spectral number density of ADAF photons is given by n_ADAF(ϵ) = L_ADAF(ϵ)/(c ⋅ ϵ ⋅ A_ADAF), where $A_{\mathrm{ADAF}}=2\pi R_{\mathrm{out}}^2+2\pi R_{\mathrm{in}}^2+2\pi ·\left(R_{\mathrm{out}}\sqrt{(}1.25R_{\mathrm{out}}^2)-R_{\mathrm{in}}\sqrt{(}1.25R_{\mathrm{in}}^2)\right)$ is a rough estimate of the surface area of the ADAF. We determined the total surface area as a sum of the outer cylindrical lateral area and the inner cylindrical lateral area, and twice the outer conical lateral was subtracted by the inner conical lateral area. Additionally, we assumed the height of the ADAF to be r at radius r. As mentioned above, the ADAF spectrum is non-vanishing up to the energy ϵ_ADAF,  0. Because PP is a threshold process, the lowest possible ADAF photons capable of PP have the energy ${ϵ}_{\text{ADAF},\text{PP th}}=m_{\text{e}}^2{c}^4/{ϵ}_{\text{ADAF}, 0}$. Consequently, the total number density of ADAF photons capable of PP is obtained by

$$\begin{array}{c}\hfill n_{\mathrm{PP},1}\approx {\displaystyle {\int }_{{ϵ}_{\mathrm{ADAF},\mathrm{PP}\mathrm{th}}}^{{ϵ}_{\mathrm{ADAF},0}}{n}_{\mathrm{ADAF}}(ϵ)\mathrm{d}ϵ.}\end{array}$$(12)

As a comparison, we determined the total number density of pair-producing MeV photons based on the approximation

$$\begin{array}{c}\hfill n_{\mathrm{PP},2}\approx 1.4\times {10}^{17}{\dot{m}}^2\frac{{10}^9{M}_{\odot }}{M}{\mathrm{m}}^{-3}\end{array}$$(13)

from Levinson & Rieger (2011), which is motivated by findings for the bremsstrahlung luminosity of an ADAF by Narayan & Yi (1995). For the ADAF parameters we used above, the assumed ξ = 1%, and the constrained pairs of ṁ and T_e, we give the correspondingly values of n_PP,  1 and n_PP,  2 in the second and third row of Table 3. Except for the last pair, we have an order-of-magnitude agreement. In contrast, for the canonical ξ = 10% and the corresponding lower T_e, we obtain values for n_PP,  1 that are about two orders of magnitude lower (cf. Table 3) than in the 1% reprocessing fraction case.

From the number density of pair-producing ADAF photons, we estimated the materialisation rate in the gap. With 0.2 σ_Th as an approximation for the PP cross section, the number density of pairs that are produced in the gap region per unit time is $K_{\text{PP}, \text{gap}}\approx 0.2{\sigma }_{\text{Th}}{n}_{\text{PP}, 1}^{2}c$. We show the values of K_PP,  gap for the ṁ, T_e) pairs in the fourth row of Table 3. All four estimates in the ξ = 1% case are about six orders of magnitude smaller than our determined K_G and all four estimates in the ξ = 10% case are about ten orders smaller. This does not necessarily cast doubt on our model.

These pair-produced seed electrons are accelerated in the gap and multiplied in the gap and during the propagation to the interaction with the emission-line photons. The multiplication is mediated by the emission of synchro-curvature radiation in the magnetospheric magnetic field, by the emission of IC up-scattered radiation in the ADAF photon field, and by the subsequent materialisation of this radiation through PP. In net effect, this multiplication increases the number of electrons by the number 𝒩 and also lowers the energy per electron by 𝒩. Because we neither included the magnetic field and the emission of synchro-curvature radiation nor the exact spatial distribution of the ADAF photon field into our model, a detailed analysis of the in-gap and post-gap multiplication is beyond the scope of this work. We can, however, estimate from comparison of K_PP,  gap with K_G that it has to be 𝒩 ≈ 10⁶ in the 1% reprocessing fraction case and 𝒩 ≈ 10¹⁰ in the 10% case. The finding 𝒩 ≈ 10⁶ is three orders of magnitude larger than the multiplication found by Broderick & Tchekhovskoy (2015). These authors, however, determined the multiplication only in the post-gap cascade. Considering that up to 10⁷ curvature photons can be emitted by one electron in the initial spark event of the gap (Levinson & Cerutti 2018), 𝒩 ≈ 10⁶ after the post-gap cascade is plausible and does not inhibit gap stability. Furthermore, to assess the magnitude of 𝒩, we have to consider the magnetic field structure (especially the divergence of the magnetic field), which can dilute the number density of the electron bunch.

The maximum potential drop in the gap is connected with the maximum Lorentz factor of the electrons. Knowing from our modelling (cf. Eq. (4) and Table 2) that the maximum energy per electron after the multiplication is γ_i,  0 m_ec² = (γ_mean + 3 ς) m_ec² = 5.7 × 10¹² eV, we can estimate that the maximum energy that one original seed electron has attained in the gap voltage is ≈ 𝒩 ⋅ γ_i,  0 m_ec². In the ξ = 1% case, this gives ≈ 5.7 × 10¹⁸ eV, which agrees with estimated (Levinson & Rieger 2011) and simulated (Hirotani & Pu 2016) values for the potential drop in a vacuum gap. In the ξ = 10% case, this gives ≈ 5.7 × 10²² eV, which seems hard to reach even for extreme magnetospheric gap conditions. Thus, the canonical choice ξ = 10% for the reprocessing fraction does not fit a consistent picture of Mrk 501. Levinson & Segev (2017), Levinson & Cerutti (2018), Chen & Yuan (2020), Crinquand et al. (2020), and Kisaka et al. (2020) have recently cast doubt on the prior attitude of vacuum gaps being steady-state phenomena. Instead, they might be of intermittent or cyclic character. In particular, Levinson & Cerutti (2018) found that charge-starved gaps, after ignition by injection of gamma rays or electrons abruptly discharge in a gamma-ray flaring event on a timescale r_S/(2c), followed by relaxation towards a state of weak-amplitude oscillations of the electric field and plasma. In this case, rule-of-thumb estimates for the maximum potential drop would not be adequate for conclusions about the attained electron Lorentz factor. We therefore leave detailed inferences from our input parameters about conditions in the vacuum gap to future studies.

If the electron beam is emitted from an oscillating gap, we might expect that this also leads to cyclic behaviour of the cascade with the emission-line photons. However, Broderick & Tchekhovskoy (2015) found that if the gap width exceeds its height, then the gap consists of a number of causally independent tubes. The cyclic behaviour after the initial spark can therefore hardly be traced because the emitted electron beam is the transverse average of many independently fluctuating gap tubes.

4.3.4. Comparison of energy budgets

The total power that is injected into the interaction volume is determined to be $P_{\mathrm{i}}=4\pi R^3/3·{\int }_{{\gamma }_{\mathrm{i},1}}^{{\gamma }_{\mathrm{i},0}}\gamma {\dot{N}}_{\mathrm{i}}(\gamma )\mathrm{d}\gamma ·m_{\mathrm{e}}{c}^2=1.9\times {10}^{33}$ W, where R is the radial size of the interaction region. This power is due to the electrons that have inherited their energy from gap acceleration. Hence, P_i has to be smaller than the maximum power P_gap that can be extracted from the gap, which itself cannot exceed the corresponding Blandford-Znajek power P_BZ. Moreover, we determine the power, that is radiated by the electrons in the SSC emitting zone to be P_SSC = 4.4 × 10³⁶ W.

The magnetic flux density at the BH horizon upon assuming an open gap was estimated by Aharonian et al. (2017) to be B_BH = 0.19 ⋅ (β_m/M₉)^4/7 T, where β_m is the magnetisation of the ADAF. Assuming the maximum possible electric field strength, gap extension, BH spin, and charge density, Aharonian et al. (2017) obtained $P_{\mathrm{gap}}=4.8\times {10}^{37}{\beta }_{\mathrm{m}}^{8/7}\kappa h{M}_9^{6/7}{sin}^2\theta$ W, where we have substituted for the variability timescale by the height h of the gap in units of Schwarzschild radii, and where θ is the polar angle, and κ is the multiplicity in the gap as defined by Aharonian et al. (2017). For the values θ = π/4, κ = 0.1 (weakly screened gap), h = 0.25 (thickness is half a gravitational radius), and β_m = 0.5 (magnetic pressure is half the gas pressure), we find P_gap = 2.7 × 10³⁵ W. For the same parameters and for the dimensionless BH spin parameter a_* = 0.99, we find $P_{\text{BZ}}={10}^{38}{a}_{\ast }^2{M}_9{(B_{\text{BH}}/\text{T})}^2\text{W}=1.6\times {10}^{36}$ W (Hirotani et al. 2016), which is twice the accretion flow luminosity (cf. Sect. 4.3.2).

To conclude, we find P_i/P_BZ ≈ 0.001. This agrees with estimates by Chen & Yuan (2020) and Kisaka et al. (2020) and especially with the finding of Levinson & Cerutti (2018): a fraction ≈ 0.001 of P_BZ is dissipated in a time span ≈ 7 r_S/(2 c) after the initial discharge of a gap. When we again assume M₉ = 1, this time span is ≈ 9.6 h, the approximate duration of the night of observation of the narrow feature. The comparison P_gap/P_BZ ≈ 0.2 is due to the upper limit nature of the estimates by Aharonian et al. (2017). The ratio P_SSC/P_BZ ≈ 3 is acceptable within the uncertainty of M and changes of the magnetic field strength.

5. Summary

We have modelled the interaction of particle beams from vacuum gaps in the magnetospheres of spinning BHs in blazars with recombination-line photons from surrounding gas clouds. Such clouds are commonly found in the BLR of quasars, but could also arise in early-type galaxies as a result of red giants that cross the particle beam. Numerically solving the coupled kinetic equations corresponding to linear IC pair cascades with escape terms, we obtained the steady-state solution of the SED of the photons escaping from the interaction zone. Compared to previous studies of cascades in AGN, our approach considered three main ramifications: Firstly, the inclusion of an escape term corresponds to a source that is finitely extended and not necessarily optically thick. In contrast to the work of Zdziarski (1988), which dealt with saturated and hence optically thick cascades, our scenario can therefore be applied to astrophysical objects that are optically thin to the emission of HE gamma rays. This is only a first step towards the self-consistent treatment of a spherically symmetric emission region with a radially varying emission and absorption coefficients. For the SSC mechanism without pair cascades, such a treatment was developed by Gould (1979). In our framework, the photon density reacts back upon itself through the electron density. We therefore did not treat the emission coefficient as prescribed. Secondly, Wendel et al. (2017) determined the photon spectral number density n_γ essentially as the accumulation of all the IC up-scattered photons. In contrast, we determined n_γ by considering both IC up-scattering and losses through escape and pair absorption (cf. Eq. (1)). Thirdly, we used an efficient numerical iteration scheme (cf. Sect. 3) that saves computational cost by iterating from values γ_k with high k to lower k, and so using the optimal course of N(γ) in each iteration step.

Applying the model to the blazar Mrk 501, we showed that the recently observed peculiar narrow spectral feature at 3 TeV can readily be explained when rather plausible physical parameters for the beam and recombination-line clouds are adopted. In the scenario supported by this observation, photons from a weakly accreting hot ADAF materialise as pairs in the central BH magnetosphere. These seed electrons are accelerated in a vacuum gap and are multiplied by post-gap cascades while travelling away from the central region and carry about 0.1% of the Blandford-Znajek luminosity. The ADAF also irradiates and ionises ambient gas clouds. In agreement with measurements of the hydrogen Lyman-α luminosity in Mrk 501, these clouds reprocess about 1% of the original ADAF luminosity as emission-line photons. They act as a target colliding with the electron beam. The resulting cascade photons are emitted in the direction of the primary beam electrons, shaping the intermittent narrow bump at about 3 TeV in the SED. In combination with the SSC spectrum produced by particles presumably accelerated at shock waves further downstream, the model explains the broadband SED of Mrk 501 at MJD 56857.98, assuming no further pair attenuation within the host galaxy to occur. Predictions about the duty cycle of this feature are difficult considering the intermittency of gap formation and randomness of nearby cloud passages.

Acknowledgments

We thank Dorit Glawion, and Amit Shukla for fruitful discussions on the topic. We are also thankful to the anonymous referee for constructive criticism on the manuscript. Moreover, we thank Astrid Peter and Julian Sitarek for comments on this work as well as J. D. Hunter and co-workers for the development of matplotlib (Hunter 2007). C. W. gratefully acknowledges support by the project “Promotion inklusive” of the Universität zu Köln and the German Bundesministerium für Arbeit und Soziales. J. B. G. acknowledges the support of the Viera y Clavijo program funded by ACIISI and ULL.

References

Appendix A: Derivation of the kinetic equations

As an attachment to Sect. 2.2, here we give a precise derivation of the kinetic equations. To quantify the change rates of the spectral number densities of the photons and of the electrons, it is necessary to introduce the spectral interaction rates of IC scattering and of PP. We do this in line with Zdziarski (1988).

We assume an IC scattering event of an incident electron with original energy γ ≫ 1 off a photon with energy x ≪ 1 that results in a down-scattered electron with final energy γ′ and a highly energetic photon with energy x_γ. Energy conservation yields γ ≈ γ + x = γ′+x_γ. From IC scattering kinematics, it follows

$$\begin{array}{c}\hfill \gamma x>E_{\ast }\end{array}$$(A.1)

with the abbreviation E_* = (γ/γ′−1)/4 (cf. Zdziarski 1988). This ensures that the photon energy has a maximum possible value x_γ,  max(γ, x):=(4γ²x)/(1 + 4γx). This function is <  γ for all realistic values of γ, which reflects that the electron cannot transfer all its energy to the photon. Moreover, we remark that x_γ,  max(γ, x) increases with increasing γ as well as with increasing x. Hence, if N has an upper cut-off at γ_i,  0 and n₀ has an upper cut-off at x₀, then the distribution of highly energetic photons vanishes above x_γ,  max(γ_i,  0, x₀) and is non-vanishing below this.

The spectral IC scattering interaction rate C for such events on the field of low-energy photons with spectral number density n₀ is given by Zdziarski (1988) as an approximation of the exact one found by Jones (1968),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill C(\gamma ,\gamma \prime )=& {\displaystyle {\int }_{max(x_1,E_{\ast }/\gamma )}^{x_0}{n}_0(x)\frac{3{\sigma }_{\mathrm{Th}}c}{4{\gamma }^{2}x}}\hfill \\ \hfill & \times [r+(2-r)\frac{E_{\ast }}{\gamma x}-2{\left(\frac{E_{\ast }}{\gamma x}\right)}^2-2\frac{E_{\ast }}{\gamma x}ln\frac{\gamma x}{E_{\ast }}]\mathrm{d}x.\hfill \end{array}\end{array}$$(A.2)

Here it is r = (γ/γ′+γ′/γ)/2 and the second argument in the lower integration border is due to Eq. (A.1). N(γ)C(γ, γ′) is the number of IC scattering events with original energy γ and final energy γ′ per unit time, per unit space volume, per unit original energy, and per unit final energy.

Substituting γ′ with γ − x_γ in C, we realise that C(γ, γ − x_γ) describes the probability for an IC scattering event with original electron energy γ to produce a photon with energy x_γ, per unit time, and per unit x_γ-interval. N(γ)C(γ, γ − x_γ) stands for the number of IC scattering events with original electron energy γ that result in a photon with energy x_γ, per unit time, per unit space volume, per unit γ-interval, and per unit x_γ-interval.

Now, we consider a collision of a highly energetic photon of energy x_γ with a photon of energy x, resulting in the PP of an electron of energy γ and a positron with energy x_γ + x − γ ≈ x_γ − γ. Based on Agaronyan et al. (1983), Zdziarski (1988) gave the corresponding kinematic relationship

$$\begin{array}{c}\hfill x_{\gamma }x>E_{\ast }>1\end{array}$$(A.3)

where $E_{\ast }=x_{\gamma }^2/(4\gamma (x_{\gamma }-\gamma ))$ as well as the spectral PP interaction rate P, again for a photon field n₀,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(x_{\gamma },\gamma )=& {\displaystyle {\int }_{max(x_1,E_{\ast }/\gamma )}^{x_0}{n}_0(x)\frac{3{\sigma }_{\mathrm{Th}}c}{4x_{\gamma }^{2}x}}\hfill \\ \hfill & \times [r-(2+r)\frac{E_{\ast }}{x_{\gamma }x}+2{\left(\frac{E_{\ast }}{x_{\gamma }x}\right)}^2+2\frac{E_{\ast }}{x_{\gamma }x}ln\frac{x_{\gamma }x}{E_{\ast }}]\mathrm{d}x.\hfill \end{array}\end{array}$$(A.4)

Here it is r = (γ/(x_γ − γ)+(x_γ − γ)/γ)/2 and the second argument in the lower integration border is due to the first part of Eq. (A.3). n_γ(x_γ)P(x_γ, γ) is the number of PP events with photon energy x_γ that create an electron with energy γ per unit time, per unit space volume, per unit photon energy, and per unit electron energy.

To infer the kinetic equation of the electrons, we have to specify all relevant sinks and sources of the electron energy distribution in terms of change rates of the electron spectral number density N: First, the spectral injection rate Ṅ_i(γ) of electrons directly enters the kinetic equation. Second, electrons leave the interaction region on the escape timescale T_e esc. This is expressed with the term N(γ)/T_e esc(γ).

Third, an electron of original energy γ loses energy by becoming IC down-scattered to a lower energy γ′< γ, resulting in a sink at γ. We remark, however, that for given x and γ, the minimum possible value of γ′ is not equal to 1 but has a lower limit at ${\gamma }_{\text{IC}, \text{min}}^{\prime }(\gamma ,x):=\gamma /(1+4x\gamma )$ as a result of scattering kinematics. From Eq. (A.1) the limitation ${\gamma }^{\prime }>{\gamma }_{\text{IC}, \text{min}}^{\prime }$ follows. It always is ${\gamma }_{\text{IC}, \text{min}}^{\prime }>1$, reflecting the fact that the electron retains an amount of kinetic energy in every case. ${\gamma }_{\text{IC}, \text{min}}^{\prime }$ decreases with increasing x. Now, if the low-energy photons are not mono-energetic, but are distributed in energy space with an upper cut-off of n₀(x) at some x₀, then the value of γ′ that can at least be reached for an electron with original energy γ, is ${\gamma }_{\text{IC}, \text{min}}^{\prime }$(γ, x₀). Therefore we describe the number of IC scattering events with original energy γ and any allowed final energy γ′ per unit time, per unit space volume, and per unit original energy by ${\int }_{{\gamma }_{\text{IC},\text{min}}^{\prime }(\gamma ,x_0)}^{\gamma }N(\gamma )C(\gamma ,{\gamma }^{\prime })\text{d}{\gamma }^{\prime }$. The lower integration border is in contrast to the one mistakenly used by Zdziarski (1988) in the second term on the right-hand side of his Eq. (1) as well as in the first integral in A21.

Fourth, an electron loses energy by becoming IC down-scattered to a lower energy γ′< γ, resulting in a source at γ′. We can infer from Eq. (A.1) that a final electron energy γ′ cannot be reached from any original electron energy γ. For γ′≥1/(4x) (scattering into the KN regime), the original electron energy is indeed unbounded. However, for γ′< 1/(4x) (scattering into the Thomson regime), the original electron energy has to be lower than γ′/(1 − 4xγ′), which increases with increasing x. Again, if the low-energy photon spectral number density is extended with an upper cut-off at x₀, then the upper boundary of γ is expressed by γ′/(1 − 4x₀γ′). In what follows we define

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{IC},\max }(\gamma \prime ,x):=\{\begin{array}{cc}\gamma \prime /(1-4x\gamma \prime )\hfill & \mathrm{for}\gamma \prime <1/(4x),\hfill \\ +\infty \hfill & \mathrm{for}\gamma \prime \ge 1/(4x).\hfill \end{array}\end{array}$$(A.5)

With this, we can quantify the number of IC scattering events with final energy γ′ and any permitted original energy γ per unit time, per unit space volume, and per unit final energy by ${\int }_{{\gamma }^{\prime }}^{{\gamma }_{\text{IC},\text{max}}({\gamma }^{\prime },x_0)}N(\gamma )C(\gamma ,{\gamma }^{\prime })\text{d}\gamma$. The upper integration border is in contrast to the one erroneously used by Zdziarski (1988) in the third term on the right-hand side of his Eq. (1). When this term is plugged into the kinetic equation, we have to exchange γ with γ′ in this term because the kinetic equation refers to the variable γ.

Fifth, PP of an electron with energy γ is a source at this energy. However, PP on a photon with energy x that produces an electron with energy γ is possible only if the photon energy x_γ exceeds the PP threshold x_{γ,  PP,  th}(γ, x):=γ/(1 − (1/(4γx))). This threshold can be inferred from Eq. (A.3). x_{γ,  PP,  th} adopts finite and positive values only for γ >  1/(4x), expressing that electrons cannot be pair-produced below 1/(4x). In this range, x_{γ,  PP,  th} has a minimum at γ = 1/(2x). This minimum corresponds to x_{γ,  PP,  th}(1/(2x),x) = 1/x, which is the well-known set-in of PP. For γ >  1/(4x), it is of course x_{γ,  PP,  th} >  γ. Furthermore, we can show that x_{γ,  PP,  th} decreases with increasing x. If the low-energy photons are not mono-energetic, but are distributed with n₀(x), which vanishes above an upper cut-off x₀, then PP is therefore possible as soon as x_γ >  x_{γ,  PP,  th}(γ, x₀). The number of kinematically allowed PP events that create an electron with energy γ per unit time, per unit space volume, and per unit electron energy, is accordingly expressed by ${\int }_{x_{\gamma ,\text{PP},\text{th}}(\gamma ,x_0)}^{\infty }{n}_{\gamma }(x_{\gamma })P(x_{\gamma },\gamma )\text{d}{x}_{\gamma }$. As P is normalised to the production of only one electron, this term has to be multiplied by 2 in the kinetic equation.

The electron kinetic equation is obtained by setting the rate of change Ṅ(γ) of the electrons’ spectral number density equal to the sum of all source terms subtracted by all sinks,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{N}(\gamma )=& {\dot{N}}_{\mathrm{i}}(\gamma )-N(\gamma )(\frac{1}{T_{\mathrm{e}\mathrm{esc}}(\gamma )}+{\displaystyle {\int }_{\gamma {\prime }_{\mathrm{IC},\min }(\gamma ,x_0)}^{\gamma }}C(\gamma ,\gamma \prime )\mathrm{d}\gamma \prime )\hfill \\ \hfill & +{\displaystyle {\int }_{\gamma }^{\gamma {\prime }_{\mathrm{IC},\max }(\gamma ,x_0)}}N(\gamma \prime )C(\gamma \prime ,\gamma )\mathrm{d}\gamma \prime \hfill \\ \hfill & +2·{\displaystyle {\int }_{x_{\gamma ,\mathrm{PP},\mathrm{th}}(\gamma ,x_0)}^{\infty }}{n}_{\gamma }(x_{\gamma })P(x_{\gamma },\gamma )\mathrm{d}{x}_{\gamma }.\hfill \end{array}\end{array}$$(A.6)

In what follows the first term of this equation is called electron spectral injection rate, the part within round brackets is called electron spectral loss rate (which is the sum of the spectral escape rate and of the spectral IC down-scattering rate to lower energies), and the sum of the third and fourth term are called electron spectral production rate (which is the sum of the spectral IC down-scattering rate from higher energies and of the spectral PP rate). We note that the dimension of the electron spectral injection rate and of the electron spectral production rate is [γ]⁻¹[t]⁻¹[l]⁻³, respectively, and that of the electron spectral loss rate is [γ]⁻¹[t]⁻¹.

Analogously, we identify all source and sink terms that determine the rate of change ṅ_γ(x_γ) of the highly energetic photon spectral number density: First, the injection term ṅ_γ,i(x_γ) has to be included. Second, the spectral loss rate due to escape from the interaction volume is described by n_γ(x_γ)/T_ph esc(x_γ).

Third, PP is a sink of the highly energetic photon distribution at the energy x_γ of the incident photons. For given x and x_γ, the energies γ that can be reached by the created electrons, are limited due to the interaction kinematics Eq. (A.3). The electron energy γ that can be obtained obeys γ_PP,  min <  γ <  γ_PP,  max with the limits γ_PP,  min(x_γ, x):=x_γ(1 − (1 − 1/(x_γx))^1/2)/2 and γ_PP,  max(x_γ, x):=x_γ(1 + (1 − 1/(x_γx))^1/2)/2. With increasing x, the limit γ_PP,  min decreases and γ_PP,  max increases. Consequently, if n₀ is extended and non-vanishing only below and at an upper cut-off x₀, then the electron energies that are kinematically allowed satisfy γ_PP,  min(x_γ, x₀) < γ <  γ_PP,  max(x_γ, x₀). With this in mind, the number of PP events with incident photon energy x_γ and any possible electron energy per unit time, per unit space volume, and per unit x_γ-interval can be written as ${\int }_{{\gamma }_{\text{PP},\text{min}}(x_{\gamma },x_0)}^{{\gamma }_{\text{PP},\text{max}}(x_{\gamma },x_0)}{n}_{\gamma }(x_{\gamma })P(x_{\gamma },\gamma )\text{d}\gamma$.

Fourth, the IC up-scattering of photons of energy x to energy x_γ is a source of photons. The energy x_γ can, however, only be reached, when the original electron energy γ exceeds a threshold γ_IC,  th. This threshold is defined by γ_IC,  th(x_γ, x):=x_γ(1 + (1 + 1/(x_γx))^1/2)/2, as can be inferred from Eq. (A.1), and is >  x_γ for all realistic x_γ. The constraint by this threshold is again the result of the electron not being able to transfer all its energy to the photon, and hence the original electron energy has to exceed x_γ by γ′, which cannot fall below the value γ_IC,  th − x_γ. We note that γ_IC,  th decreases with increasing x. Thus, if n₀ is extended over a range of energies and has an upper cut-off at x₀, then IC scattering to the energy x_γ sets in as soon as γ >  γ_IC,  th(x_γ, x₀). Then, ${\int }_{{\gamma }_{\text{IC},\text{th}}(x_{\gamma },x_0)}^{\infty }N(\gamma )C(\gamma ,\gamma -x_{\gamma })\text{d}\gamma$ quantifies the number of IC scattering events with any kinematically allowed original electron energy that result in a photon with energy x_γ per unit time, per unit space volume, and per unit energy.

The highly energetic photon kinetic equation is obtained by setting the rate of change ṅ_γ(x_γ) of the photon spectral number density equal to the sum of all source terms subtracted by all sinks,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{n}}_{\gamma }(x_{\gamma })=& {\dot{n}}_{\gamma ,\mathrm{i}}(x_{\gamma })-n_{\gamma }(x_{\gamma })(\frac{1}{T_{\mathrm{ph}\mathrm{esc}}(x_{\gamma })}\hfill \\ \hfill & +{\displaystyle {\int }_{{\gamma }_{\mathrm{PP},\min }(x_{\gamma },x_0)}^{{\gamma }_{\mathrm{PP},\max }(x_{\gamma },x_0)}}P(x_{\gamma },\gamma )\mathrm{d}\gamma )\hfill \\ \hfill & +\underset{:={\dot{n}}_{\gamma ,\mathrm{IC}}(x_{\gamma })}{\underbrace{{\displaystyle {\int }_{{\gamma }_{\mathrm{IC},\mathrm{th}}(x_{\gamma },x_0)}^{\infty }}N(\gamma )C(\gamma ,\gamma -x_{\gamma })\mathrm{d}\gamma }}.\hfill \end{array}\end{array}$$(A.7)

In what follows the first term of this equation is called photon spectral injection rate, the part within round brackets is called photon spectral loss rate (which is the sum of the spectral escape rate and of the spectral pair absorption rate), and the last term ṅ_γ,IC is called photon spectral production rate. We remark that the dimension of the photon spectral injection rate and of the photon spectral production rate is [x_γ]⁻¹[t]⁻¹[l]⁻³, respectively, and that of the photon spectral loss rate is [x_γ]⁻¹[t]⁻¹.

We assume the spectral number densities to be in steady state. Thus, the respective rates of change vanish in Eqs. (A.6) and (A.7). Then, we can solve the respective equation for the respective spectral number density in the second term on the right-hand side of the respective equation. This yields

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill N(\gamma )=& \frac{{\dot{N}}_{\mathrm{i}}(\gamma )+{\int }_{\gamma }^{\gamma {\prime }_{\mathrm{IC},\max }(\gamma ,x_0)}N(\gamma \prime )C(\gamma \prime ,\gamma )\mathrm{d}\gamma \prime }{\frac{1}{T_{\mathrm{e}\mathrm{esc}}(\gamma )}}\hfill \\ \hfill & \frac{+2·{\int }_{x_{\gamma ,\mathrm{PP},\mathrm{th}}(\gamma ,x_0)}^{\infty }{n}_{\gamma }(x_{\gamma })P(x_{\gamma },\gamma )\mathrm{d}{x}_{\gamma }}{+{\int }_{\gamma {\prime }_{\mathrm{IC},\min }(\gamma ,x_0)}^{\gamma }C(\gamma ,\gamma \prime )\mathrm{d}\gamma \prime }\hfill \end{array}\end{array}$$(A.8)

and Eq. (1).

Appendix B: Discussion of numerical solution procedure

As an extension to Sect. 3, in this appendix we concisely explain how explicit converging solutions for N are found from iteratively solving N(γ) = ℱ(N, γ) along the sequence of points (γ_k)_{k = 1, …, κ} after prescribing the initialisation N_init(γ).

We could iterate the following way: We could determine ℱ at all the points γ_k based on the sequence (N_k)_jinit. Then, we would assign the attained values to a new sequence (N_k)₀. In short, the initialisation step of the iteration would be described by N_k,  0 = ℱ((N_k)_jinit, γ_k). The same technique would be pursued in the subsequent iteration steps. In each iteration step, we would compute a complete new sequence (N_k)_j (with integer j running along j_init,  0,  1, …,  j_final and denoting the number of the iteration step, or equivalently, the number of the sequence of values N_k,  j). In other words, in the jth iteration step, we would determine the κ values N_k,  j = ℱ((N_k)_j − 1, γ_k) based on the sequence (N_k)_j − 1 of the previous iteration. Iterating would have to proceed until all points N_k,  j converge with a value called N_{k,  jfinal} (and hence the function N(γ) converges), which is always the case according to the Banach fixed-point theorem.

Wendel et al. (2017) used this iteration scheme. However, this scheme is not favoured because it is computationally inefficient. This is shown as follows. By inspection of Eq. (2) we realise that a computation of ℱ at a given point γ_k′ evaluates (and thus needs to know) N only at points γ_k ≥ γ_k′. This observation is obvious in the case of the second term in the numerator of Eq. (2) because the integration range stretches from γ upwards. In the case of the third term, we have to consider the values of the lower integration borders of the nested integrals. First, we make use of the fact that the lower border x_{γ,  PP,  th}(γ, x₀) of the outer integral is the PP threshold and thus >  γ (for all realistic values of γ). Second, as mentioned in Sect. 2.2, the lower border γ_IC,  th(x_γ, x₀) of the inner integral is always >  x_γ. As an effect, the lowest value, at which N has to be evaluated during the computation of the third term, is γ_IC,  th(x_{γ,  PP,  th}(γ, x₀),x₀), which is always >  γ. We consider the computation of a certain value N_k′, j, in other words, the jth computation of the value of N at γ_k′. This computation is based on all the values N_{k,  j − 1} with k ≥ k′, in other words, on the sequence (N_k)_{k = k′+1, …, κ;  j − 1}. As long as the values with k >  k′ have not yet converged, the computation of N_k′, j is based on inaccurate prerequisites and thus becomes inaccurate by itself. It is therefore reasonable to implement the following iteration scheme:

In contrast to above, where the complete sequence (N_k)_j of κ values was determined in one iteration step, we iterate point-wise. In the first step, only the value N_κ,  j at the largest γ is iterated. Briefly written, we calculated N_κ,  0 = ℱ((N_k)_jinit, γ_κ) and N_κ,  j = ℱ((N_k)_j − 1, γ_κ). Not until convergence of N_κ,  j was achieved with the final value being N_{κ,  jfinal}, did we iterate at the next lower point γ_κ − 1. In short, we computed N_{κ − 1,  0} = ℱ((N_k)_jinit, γ_κ − 1) and N_{κ − 1,  j} = ℱ((N_k)_j − 1, γ_κ − 1) until N_{κ − 1,  j} converges with the final value being N_{κ − 1,  jfinal}. Then, we switch to γ_κ − 2 and iterate again until convergence. This was done successively at all points until convergence was reached at the point k = 1. This procedure is advantageous because during the iteration at a point γ_k′, the final converged sequence (N_k)_{k = k′+1, …, κ;  jfinal} of values is already known and can be used, which reduces inaccuracies at γ_k′.

In the range γ >  1/(4x₀) pairs can be produced. This is reflected in the third term in the numerator of Eq. (2) being non-zero. This makes the computation of ℱ expensive in this regime. At γ ≤ 1/(4x₀), PP does not happen and the third term vanishes, which speeds up the computations. On the other side, more iteration steps are necessary for lower values of γ due to the decreasing energy transfer in one IC scattering event. In the deep Thomson regime (i.e. for γ ≪ 1/(4x₀)) Eq. (A.6) in the steady state can be expressed as a continuity equation and could be solved by a differential equation solver (e.g., a Runge-Kutta method, cf. Zdziarski 1988). However, as we are mainly interested in the course of N in the KN regime, we used the iteration procedure over the complete range.

All Tables

All Figures

Fig. 1. Sketch of the geometry (not to scale) assumed for the model. Seed electrons are accelerated in the gap and multiplied in and behind the gap through ADAF photons. For simplicity of drawing, photons are not shown in the in-gap and post-gap multiplication. After propagating away from the gap, the electron beam meets an emission-line photon field, driving a cascade, whose steady-state gamma-ray photon density is nγ and whose escaping flux is Fcasc. To convert nγ into Fcasc, cf. Eq. (3).

In the text

	Fig. 2. Sketch of the considered linear IC pair cascade with escape terms.
In the text

	Fig. 3. γ N(γ), which is a proxy for the electron spectral energy density, versus the product of the Lorentz factor with the dimensionless energy of the highest energetic line, the helium II Lyman-α line. The iterated points Nk,  j are drawn in red.
In the text

	Fig. 4. Various contributions (cf. Eq. (A.7)) to the photon spectral loss rate on the left-hand side ordinate in dependence on the product of the dimensionless highly energetic photon energy with the dimensionless energy of the highest energetic line. On the right-hand side ordinate, we plot the optical depth that corresponds to the respective contribution to the spectral loss rate.
In the text

Fig. 5. Spectral production rate and number density of the cascade in Mrk 501. The product of the squared dimensionless photon energy with the spectral production rate of photons is plotted in red (which is determined by the last term of Eq. (A.7)). We show the product of the squared dimensionless photon energy with the spectral number density of photons in green (which is determined by Eq. (1)). Both quantities are plotted vs. the product of the highly energetic photon energy with the energy of the helium II Lyman-α line.

In the text

Fig. 6. High-energy bump of the SED of Mrk 501 from 2014 July 19 (MJD 56857.98). The dashed grey line depicts the IC bump of a one-zone SSC model (contribution FSSC). The dot-dashed red line depicts the emission of the cascade (contribution Fcasc, cf. green line in Fig. 5). The sum F of both contributions is depicted by the solid black line. The spectral data from the MAGIC telescopes are shown as red circles, while measurements and upper limits by the Fermi Large Area Telescope are drawn by black and yellow triangles. The details on the data analysis and the SSC modelling can be found in MAGIC Collaboration (2020).

In the text