*Article* **Studying the Influence of External Photon Fields on Blazar Spectra Using a One-Zone Hadro-Leptonic Time-Dependent Model**

**Michael Zacharias 1,2**


<sup>2</sup> Centre for Space Science, North-West University, Potchefstroom 2520, South Africa

**Abstract:** The recent associations of neutrinos with blazars require the efficient interaction of relativistic protons with ambient soft photon fields. However, along side the neutrinos, *γ*-ray photons are produced, which interact with the same soft photon fields producing electron-positron pairs. The strength of this cascade has significant consequences on the photon spectrum in various energy bands and puts severe constraints on the pion and neutrino production. In this study, we discuss the influence of the external thermal photon fields (accretion disk, broad-line region, and dusty torus) on the proton-photon interactions, employing a newly developed time-dependent one-zone hadro-leptonic code *OneHaLe*. We present steady-state cases, as well as a time-dependent case, where the emission region moves through the jet. Within the limits of this toy study, the external fields can disrupt the "usual" double-humped blazar spectrum. Similarly, a moving region would cross significant portions of the jet without reaching the previously-found steady states.

**Keywords:** non-thermal radiation mechanisms; relativistic jets; relativistic processes; BL Lacertae objects

## **1. Introduction**

The theory of the blazar emission was transformed in the early 1990s by the introduction of the so-called *external-Compton* scenario. The scenario explains the high-energy component of the spectral energy distribution (SED) through relativistic electrons inverse-Compton (IC) scattering soft, thermal photon fields that originate outside the jet. This transformation of blazar research was significantly driven by the works of Reinhard Schlickeiser and collaborators employing the accretion disk (AD) as a source for soft external photons [1–4].

Blazars, a sub-class of active galaxies, are indeed peculiar objects with—in the words of Reinhard Schlickeiser [5]—

"properties [that] include high optical polarization, extreme optical variability, flat-spectrum radio emission associated with a compact core, and apparent superluminal motion. Such properties are thought to be produced by those few, rare extragalactic radio galaxies and quasars that are favorably aligned to permit us to look almost directly down a relativistically outflowing jet of matter expelled from a supermassive black hole."

Despite the decades of research that have passed since the advent of the external-Compton model, a clear consensus on the source of the *γ*-ray emission in blazars has not yet been reached.

The low-energy component of the double-humped SED is the least controversial part, as synchrotron emission of relativistic electrons fits all the required properties (including the aforementioned polarization). However, the nature of the high-energy component is subject of intensive discussions. Within leptonic models, it is explained through IC emission–either with the self-made synchrotron photons (synchrotron-self Compton, SSC)

**Citation:** Zacharias, M. Studying the Influence of External Photon Fields on Blazar Spectra Using a One-Zone Hadro-Leptonic Time-Dependent Model. *Physics* **2021**, *3*, 1098–1111. https://doi.org/10.3390/ physics3040069

Received: 22 September 2021 Accepted: 8 November 2021 Published: 18 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

or with external photons, such as the AD, photons from the broad-line region (BLR) [6], the dusty torus (DT) [7,8], or even the cosmic microwave background (CMB) [9,10]. However, if relativistic jets are also capable of accelerating protons, the *γ*-rays could also originate from such interactions, through direct proton-synchrotron emission, or via proton-photon interactions, causing a cascade of pairs [11–15]. Especially the production of pions would have the capability to discriminate between the leptonic and the hadronic scenario, as it would also produce neutrinos. While neutrinos have been associated with blazars recently [16,17], the significances are not yet sufficient to claim a real detection. Nonetheless, the discussion is ongoing [18], and the upcoming neutrino observatories KM3NET and IceCube-Gen2 may provide the definitive answer.

In this study, the hadro-leptonic model is combined with the external soft photons, to study their influence on the resulting pair cascade and the jet emission. A newly developed time-dependent, one-zone hadro-leptonic code—*OneHaLe*—is introduced in Section 2. It is used in Section 3 to study the influence of the external photon fields by first calculating steady-state spectra at various locations within the jet, as the region of influence of the soft photon fields on the jet is strongly distance-dependent. Subsequently, we present the case of an emission region moving outward passing through the various external photon fields. We note that the study conducted is a toy model: In order to properly identify the influence of the external fields, all other parameters of the emission region remain the same, irrespective of the location. This may have significant consequences for the emerging spectra. Section 4 provides the discussion of the results and the conclusions.

#### **2. Code Description**

The code is based on the recently developed extended hadro-leptonic steady-state code *ExHaLe-jet* [19]. In fact, the fundamental equations governing the particle and radiation processes are the same, and we only provide a brief overview here describing the free parameters. In the following, quantities in the host galaxy frame are marked with a hat, while quantities in the observer's frame are marked by the superscript "obs". Unmarked quanitites are either in the co-moving frame of the emission region or invariant.

A spherical emission region is assumed with radius *R* located a distance *z*<sup>0</sup> from the black hole within the jet, pervaded by a tangled magnetic field of strength *B*. The emission region moves with bulk Lorentz factor Γ under a viewing angle *θ*obs with respect to the observer's line-of-sight implying a Doppler factor, *<sup>δ</sup>* = [Γ(<sup>1</sup> <sup>−</sup> *<sup>β</sup>*<sup>Γ</sup> cos *<sup>θ</sup>*obs)]−1, where *<sup>β</sup>*<sup>Γ</sup> <sup>=</sup> <sup>√</sup> 1 − Γ−2.

The Fokker-Planck equation governing the time-dependent evolution of a given particle species *i* (protons, charged pions, muons, or electrons) with spectral density *ni*(*χ*) is given as

$$\begin{split} \frac{\partial n\_i(\boldsymbol{\chi},t)}{\partial t} &= \frac{\partial}{\partial \boldsymbol{\chi}} \left[ \frac{\boldsymbol{\chi}^2}{(a+2)t\_{\text{acc}}} \frac{\partial n\_i(\boldsymbol{\chi},t)}{\partial \boldsymbol{\chi}} \right] \\ &- \frac{\partial}{\partial \boldsymbol{\chi}} (\boldsymbol{\chi}\_i n\_i(\boldsymbol{\chi},t)) + Q\_i(\boldsymbol{\chi},t) - \frac{n\_i(\boldsymbol{\chi},t)}{t\_{\text{csc}}} - \frac{n\_i(\boldsymbol{\chi},t)}{\boldsymbol{\chi} t\_{i,\text{decay}}^\*}. \end{split} \tag{1}$$

For numerical reasons, we use the normalized particle momentum, *χ* = *pi*/(*mic*) = *γβ*, where *pi* = *γmiβc* is the particle momentum, *mi* is the particle mass, *c* the speed of light, *<sup>γ</sup>* the particle's Lorentz factor, and *<sup>β</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>γ</sup>*−2. The first term on the right-hand side of Equation (1) describes Fermi-II acceleration through scattering of particles on magnetohydrodynamic waves. The parametrization of [20] is used with *a* = 9*v*<sup>2</sup> *s*/4*v*<sup>2</sup> *A*, *vs* and *vA* the shock speed and Alfvèn speed, respectively, and the energy-independent acceleration time scale, *t*acc. This parametrization approximates the momentum diffusion through hard-sphere scattering.

The second term on the right-hand side of Equation (1) provides momentum changes *χ*˙ *<sup>i</sup>* through gains (Fermi-I acceleration *χ*˙FI = *χ*/*t*acc) and continuous losses. All charged particles lose energy through synchrotron radiation and adiabatic expansion of the emission region. Protons also lose energy through pion production and Bethe-Heitler pair production, while electrons suffer additional losses through IC scattering of ambient photon fields. These ambient fields consist of all intrinsically produced radiation fields—such as synchrotron—as well as the external photon fields, namely the AD, the BLR, the DT, and the CMB. Naturally, the pion production can turn a proton into a neutron. As we do not explicitly consider neutrons at this point, this effect is approximated here by a continuous loss process instead of a catastrophic loss. This channel is marked as "neutron" losses in Figures 2 and 6, while the nominal pion production cooling term is marked as "pion".

The remaining three terms on the right-hand side of Equation (1) mark the injection of particles, the escape of particles from the emission region, and the decay of unstable particles, respectively. *t* ∗ *<sup>i</sup>*,decay is the proper decay time scale, which is 2.6 <sup>×</sup> <sup>10</sup>−<sup>8</sup> s for charged pions, and 2.2 <sup>×</sup> <sup>10</sup>−<sup>6</sup> s for muons, respectively. As neutral pions decay after 2.8 <sup>×</sup> <sup>10</sup>−<sup>17</sup> s into *<sup>γ</sup>* rays, Equation (1) is not solved for neutral pions, while their radiation output is directly calculated from their injection spectrum.

While Fermi-I and II acceleration terms are considered here, we treat them merely as re-acceleration processes characterized by the acceleration time scale *t*acc = *η*acc*t*esc; namely, a multiple *η*acc of the escape time scale [21]. We do not consider the primary acceleration of protons and electrons, which may take place in small sub-regions of the larger emission region [20,22], but approximate it through the injection term *Q*(*χ*, *t*). Here, a simple powerlaw injection is used with spectral index *si* between a minimum and maximum Lorentz factor, *γ*min,*<sup>i</sup>* and *γ*max,*i*, respectively. The injection normalization, *Q*0,*i*(*t*), is given by

$$Q\_{0,i}(t) = \frac{L\_{\text{irr}\text{j}}(t)}{Vm\_iC^2} \begin{cases} \frac{2 - s\_i(t)}{\gamma\_{\text{max},i}^{2 - s\_i(t)} - \gamma\_{\text{min},i}^{2 - s\_i(t)}} & \text{if } s\_i(t) \neq 2\\ \left(\ln \frac{\gamma\_{\text{max},i}}{\gamma\_{\text{min},i}}\right)^{-1} & \text{if } s\_i(t) = 2 \end{cases} \tag{2}$$

with the injection luminosity, *L*inj,*i*, and the volume, *V*, of the spherical emission region. The injection functions for pions and muons are calculated directly from the photo-hadron interactions [23] and decays. Let us again emphasize that Equation (1) is explicitly solved for (charged) pions and muons.

The escape of particles is described by *t*esc = *η*esc*R*/*c*, a multiple *η*esc of the light travel time. As *η*esc > 1, this mimics the advective flow of particles through the emission region.

Equation (1) is solved with a Chang and Cooper routine [24]; for a detailed description, see [19,25].

The interaction of protons with photons can result in the creation of pions. Charged pions decay into muons, which in turn decay into electrons. During both decay processes, neutrinos are produced. The neutrino spectra are calculated following [19,26,27]. The secondary electrons produced in this decay chain are injected into the electron-Fokker-Planck equation along with the primary electrons. Additionally, secondary electrons are also produced from Bethe-Heitler pair production and *γ*-*γ* pair production.

We do not consider explicitly neutrons in this code. Their number density is low compared to the proton density [13]; so their effect is small. Nonetheless, we plan to rectify this issue in a future update of the code.

The photon density *n*ph within the emission region is governed by the radiation transport equation:

$$\frac{\partial n\_{\rm ph}(\nu, t)}{\partial t} = \frac{4\pi}{h\nu} j\_{\rm \nu}(t) - n\_{\rm ph}(\nu, t) \left( \frac{1}{t\_{\rm ec, ph}} + \frac{1}{t\_{\rm abs}} \right). \tag{3}$$

with the frequency *ν*, the Planck constant *h*, the emissivity *jν*, the photon escape time scale *t*esc,ph = 4*R*/3*c*, and the absorption time scale *t*abs due to synchrotron-self absorption and *γ*-*γ* pair production. For the latter, all internal and external photon fields are considered.

From the photon distribution, *n*ph, one can calculate the spectral luminosity in the observer's frame:

$$
\omega^{\rm obs} L\_{\nu^{\rm obs}}^{\rm obs} = \delta^4 \frac{h\nu^2 V}{t\_{\rm esc, ph}} n\_{\rm ph}(\nu, t) \,. \tag{4}
$$

Equations (3) and (4) hold for all radiation processes within the emission region.

On their way from the source to the observer, *γ*-ray photons are subject to further absorption processes. A routine to calculate the important cases of absorption in the BLR and DT following the prescription in [28] is implemented in the code.

The AD is described with a standard Shakura-Sunyaev disk [29] implying that the disk is fully described through the mass of the supermassive black hole *M*BH and its accretion efficiency *η*SS (or Eddington ratio). The proper transformation of the angles into the comoving frame is considered. The BLR and DT are approximated as isotropic photon fields in the host galaxy frame within an distance *R*ˆ BLR and *R*ˆ DT from the black hole, and their energy distribution is given through a grey-body spectrum of temperature *T*ˆ BLR and *T*ˆ DT normalized to a luminosity of *L*ˆ BLR and *L*ˆ DT, respectively.

The above description holds for both steady-state and time-dependent cases. The steady state is achieved if the proton and electron densities derived from Equation (1), do not vary by more than 10−<sup>4</sup> compared to the respective values of the previous two time steps. Time-dependency can be achieved by varying any of the free parameters, in which case steady states may not be achieved from time step to time step.

#### **3. Influence of the External Fields**

Table 1 provides an overview of the free parameters that have been described in the previous section. The given parameter values are a toy model, which we use to perform a small parameter study. The parameters are based upon the flat spectrum radio quasar 3C 279 [30], however a direct data comparison is beyond the scope of this paper.

**Table 1.** Free parameters of the code along with symbols, units, and toy model values. Quantities in the host galaxy frame are marked with a hat.


Instead, we wish to analyze the influence of the external fields on the SED and the particle distributions. We chose four locations: close to the AD (*z*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>16</sup> cm), within the BLR (*z*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>17</sup> cm), within the DT (*z*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>18</sup> cm), and outside the external fields (referred to as "jet", *<sup>z</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> 1019 cm). All other parameters remain unchanged, including the radius and the magnetic field of the emission region. This highlights that these are indeed toy models meant to study the influence of the external fields without any degeneracies introduced by varying other parameters.

The result is shown in Figures 1–3. While the SEDs are transformed into the observer's frame, the photon spectra are shown as they leave the emission region in the jet. The internal *γ*-*γ* absorption processes are fully considered (the corresponding optical depth *τγγ* is shown in Figure 4 left), however no additional absorption of *γ* rays while traveling through the photon field of the host galaxy (namely, BLR and DT) or through the cosmological photon fields (extragalactic background light and CMB) are shown. Any of these photon fields could additionally (and severely) attenuate the photon flux above 10 GeV. These absorption processes are, however, not important for the conclusions of this study.

**Figure 1.** SEDs in the observer's frame for the four locations: close to the AD (**top left**), within the BLR (**top right**), within the DT (**bottom left**), and outside the external fields (jet, **bottom right**). The black solid line marks the total photon spectrum, while the colored lines mark individual photon components as labeled. Only those processes are labeled, which are visible in at least one panel. No external absorption is applied, implying that photon spectra are shown as they leave the jet. The black dashed line marks the neutrino spectrum.

**Figure 2.** Steady-state particle distributions (times Lorentz factor squared), and relevant time scales as labeled as a function of Lorentz factor *γ* for the same locations as in Figure 1. For proton losses, the total, photo-pion, "neutron", and Bethe-Heitler loss time scales are shown. Adiabatic losses dominate at lower proton energies, where the loss time scale is constant, while synchrotron losses may contribute at the highest proton energies. For electron losses, the total, and IC loss time scales are shown. Synchrotron losses dominate, where IC losses are negligible, while adiabatic losses are irrelevant.

**Figure 3.** Steady-state electron injection rates *Q* (times Lorentz factor squared) as a function of the Lorentz factor *γ* as labeled for the same locations as in Figure 1.

**Figure 4.** Optical depth *τγγ* due to *γ*-*γ* pair production as a function of frequency (observer's frame) for the steady-state cases (**left**) and the moving-blob case (**right**) at the different positions within the jet, as labeled. The thin horizontal line marks *τγγ* = 1.

Close to the AD, the external fields are very intense, and are further enhanced through the large chosen bulk Lorentz factor of 50. In turn, the cooling of protons through protonphoton interactions is very strong (Figure 2), as indicated by the cooling time scales being dominated by pion production (indicated by the "pion" and "neutron" loss channels) at Lorentz factors *γ* > 105. This severely influences the proton distribution function and results in negligible proton synchrotron emission. The strong pion production, which can also be seen in the SED (Figure 1) through the neutral pion bump at PeV energies, results in a significant production of muons and highly relativistic electrons (Figure 3) with Lorentz factors *γ* > 1010. Similarly, highly energetic electrons are also injected through Bethe-Heitler pair production. These electrons produce *γ* rays through synchrotron emission, as well as through IC emission for lower-energetic electrons. The *γ* rays are absorbed through *γ*-*γ* pair production with all photon fields that permeate the emission region. The strength of the *γ*-*γ* absorption is shown in the left panel of Figure 4, and manifests itself in Figure 1 by the significant flux suppression at energies above 10 GeV. In turn, a strong electron-positron cascade is initiated. This results in an electron distribution, which is dominated by secondaries (Figure 3). The resulting electron synchrotron flux (Figure 1) extends through almost the entire frequency range, destroying the familiar double-hump shape in the SED. The peak of the flux at *γ* rays stems from IC scattering of AD photons.

Within the BLR, the proton cooling is drastically reduced at high Lorentz factors with cooling time scales being longer than the escape time scale of particles at all (relevant) energies (Figure 2). Unlike in the AD case, where the proton distribution cuts off sharply at *γ*max,p, in this case (and the following cases) the proton distribution extends beyond the injection cut-off because of the (re-)acceleration terms present in Equation (1). The change in the spectral shape between the AD and BLR cases allows for an enhanced proton synchrotron emission in the BLR case, influencing the SED at GeV energies (Figure 1). While pion and Bethe-Heitler pair production are reduced compared to the AD case, the pair cascade is still very significant (Figure 3) because of *γ*-*γ* pair production (Figure 4 left). While the process is less severe than in the AD case, the secondaries still dominate the electron distribution (Figure 2), and produce synchrotron emission beyond PeV energies. In the BLR case, IC emission is negligible.

This trend continues in the DT case, as the cascade weakens (Figures 3 and 4 left) and the more familiar double-humped SED emerges (Figure 1). At ultra-violet (UV) energies in the SED, a minor contribution from the AD itself is visible. The *γ*-ray peak is dominated by proton synchrotron emission, even though the secondary electron synchrotron emission still dominates at X-ray and TeV energies. The neutral pion bump is below the shown flux scale, indicating the reduced interaction of protons with photons. In fact, the protons are completely in a slow-cooling regime (Figure 2).

Lastly, the emission region is located outside the external photon fields in the "jet" case. While secondary pairs are still being produced (Figure 3), their number is low (Figure 2) due to low absorption (Figure 4 left), and their flux contribution only shows around photon energies on the TeV-scale, but at relatively low flux values (Figure 1). Apart from that, the SED is dominated by synchrotron emission of protons at X- and *γ* rays, and primary electrons in the optical domain. Both peaks are cleanly separated. The AD itself is clearly visible in the UV range as a big blue bump.

The changes in the cooling strength can also be seen in the energy densities of the particles, which are given in Table 2. The particle energy densities are always dominated by protons (by several orders of magnitude compared to the electrons in most cases). The strong cooling in the AD case results in a low particle energy density, while the reduced cooling in the other cases results in increased and comparable energy densities. Given the constant value of the magnetic field in all cases, the ratio of magnetic to particle energy density decreases from case to case but is always larger than unity.

**Table 2.** Energy densities in particles *u*par (in erg/cm3) and the ratio *uB*/*u*par of magnetic to particle energy density. The magnetic energy density in all cases is *uB* = 100 erg/cm3. The horizontal line separates the steady-state (top) from the moving (bottom) cases.


The different cases are also manifested in the emerging neutrino spectra. With the weakening production of pions and muons from case to case, the flux of neutrinos also decreases and drops below the scale of the plots in the "jet" case. The AD case produces not just the highest neutrino flux, but also a different neutrino spectral shape than the other cases with a flat maximum (or mildly double-humped structure) over almost three orders of magnitude in energy. In the BLR and DT case, the neutrino spectra show a single peak at about 100 PeV. Interestingly, all three cases would be detectable with the future IceCube-Gen2 instrument [31]. However, the unrealistic SEDs—especially in the AD and BLR cases—make it seem unlikely that neutrinos could be observed from a blazar–at least, under this simple set-up.

For the examples discussed above, we have used a bulk Lorentz factor of 50. Hence, if the emission region were moving, it would cover a lot of space in a relatively short amount of time because of the Lorentz contraction: *z*ˆ = *z*<sup>0</sup> + Γ*β*Γ*ct*, where *t* is the time since launch in the comoving frame, and *z*ˆ is the location of the emission region in the host galaxy frame. In turn, the external fields, and thus the conditions within the emission region may change quickly. We try to analyze this, by letting the emission region flow from the base (placed at six times the Schwarzschild radius (innermost stable circular orbit) of the black hole) downstream through the jet.

As before, none of the other parameters change, implying that also the primary injection of protons and electrons continues with the same rate *Q* and spectral shape throughout the simulation. This assumes a quasi-instantaneous acceleration of particles [32], as well as a continuous supply. This is not realistic, as the acceleration of particles also takes time [20]. Additionally, neither the magnetic field *B* nor the radius *R* vary. While the radius of the emission region may not expand as rapidly as the larger jet structure that surrounds it, it expands nonetheless while it travels through the jet [33] given the high energy densities in the emission region. While recent observational results [30,34] indicate compact emission regions beyond the BLR, and maybe even at tens of parsecs from the black hole, it is not clear whether these are indeed moving emission regions originating close to the black hole or turbulent cells within a larger flaring region. Similarly, while a high magnetic field can be expected close to the black hole, the expansion of the emission region causes a drop of the magnetic field with increasing distance. These considerations highlight once more the

toy character of this study. Applying such parameter changes are interesting avenues for future studies beyond the scope of this paper.

Having obtained the full journey of the emission region through the jet, we extract the SEDs and particle distributions at roughly the same distances as in the steady-state cases. In fact, given the finite time resolution in the simulation, we extract the SEDs and particle distribution at the time step closest to the respective distances of the steady-state models (AD: 9.83 <sup>×</sup> <sup>10</sup><sup>15</sup> cm, BLR: 9.75 <sup>×</sup> <sup>10</sup><sup>16</sup> cm, DT: 9.64 <sup>×</sup> 1017 cm, jet: 9.75 <sup>×</sup> <sup>10</sup><sup>18</sup> cm). In order to save computation time, while also properly resolving the initial steps within the BLR, an adaptive time step of <sup>Δ</sup>*tj* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>3+*j*/20 is used, where *<sup>j</sup>* is the step number. This ensures reasonable accuracy and resolution, and also explains why the time values given in Figures 5 and 6 are not simple increases by a factor 10, as one would expect. This is a reasonable trade-off. The results are shown in Figures 5–7, while the optical depth due to *γ*-*γ* pair production is shown in the right panel of Figure 4. Any times and time scales discussed below are in the comoving frame.

The changes to the SEDs and the particle spectra are profound. The emission region has passed the AD position after merely 5 ks. The bright external photon fields cause proton-photon interactions, producing a significant amount of pions (Figure 6), which decay into photons or muons and pairs. In fact, Figure 7 shows that the injection of pairs from muon decay is almost at the level as in the steady state (Figure 3), but Bethe-Heitler produced pairs are about two orders of magnitude below. Similarly, *γ*-*γ* pair production is below the steady-state level, because the internal photon fields (Figure 5) have not yet been fully developed. In turn, the optical depth due to *γ*-*γ* pair production (Figure 4 right) is not at the steady-state level—merely the absorption caused by external fields is fully present. One consequence is the reduced absorption of PeV photon energies, allowing for a very strong flux in the neutral pion bump (Figure 5). The "under-development" of the internal photon fields is a consequence of the low electron and proton densities (Figure 6) compared to the steady-state values. The consequence is the absence of the "nominal" electron synchrotron bump in the infrared domain. The *γ* rays are dominated by IC scattering of AD photons—although orders of magnitude below the steady-state case.

The situation only changes mildly until the BLR position is reached after 5.7 <sup>×</sup> <sup>10</sup><sup>4</sup> s. This is still less than the escape times of photons (2 <sup>×</sup> <sup>10</sup><sup>5</sup> s) and particles (7.5 <sup>×</sup> <sup>10</sup><sup>5</sup> s). Therefore, particle and photon densities continue to increase. The spectral shape of the SED shown in Figure 5 is somewhat similar to the steady-state case (Figure 1), but at a factor of a few reduced in flux. There are a few more details where SEDs differ. In the *γ*-ray domain, Figure 5 shows contributions from IC scattering of both the AD and the BLR. Comparing the IC/AD spectra of the top panels in Figure 5, one notices the similarity between them. Given that not even one light-crossing time scale has passed since the launch, the photons produced below have not yet vanished from the emission region, and therefore continue to contribute to the SED even though the IC/AD production has reduced a lot at this distance. The IC/BLR spectrum shows a different spectral shape and a higher flux than in the steady-state case, which can be attributed to the slightly different shapes in the electron distributions (Figure 6 vs. Figure 2). These are a consequence of the reduced *γ*-*γ* pair production at lower energies (Figure 7). As the protons have also not reached the steady-state density, their synchrotron flux is reduced compared to the steady state, while pion and muon production are similarly reduced. Most notably, the neutral pion decay flux is barely visible at PeV to EeV energies in Figure 5—a reduction of about an order of magnitude compared to the steady state.

After 5.7 <sup>×</sup> 105 s, the emission region has reached the DT position. The time the emission region has traveled is now comparable to the escape time scales of light and particles. In turn, the SED in Figure 5 is almost equal to the steady-state case (Figure 1), except for a reduced peak *γ*-ray flux by a factor of a few. This can be attributed to the still lower number of protons compared to the steady state, resulting in an equally reduced proton synchrotron flux. This is also coupled to the low efficiency of protonsynchrotron emission, implying that the flux needs more time to build compared to the

electron synchrotron flux, which is basically instantaneous–cf., the cooling time scales in Figure 6, where electron synchrotron cooling is faster than basically any other time scale (including the travel time), while the proton synchrotron cooling only dominates at energies beyond the cut-off of the proton distribution. The IC/AD and IC/BLR components visible in the SED (Figure 5) exhibit a flux about an order of magnitude below the flux at the BLR position. This corresponds very well to an exponential decay, as the photons leave the emission region without being replenished.

The following, relatively long cruise towards the "jet" position (reached after 5.8 <sup>×</sup> <sup>10</sup><sup>6</sup> s) allows for the near-complete relaxation of the emission region towards the steady state that was obtained above. At this position, SED and particle distributions are practically equal to the steady-state case.

The particle energy densities change considerably from position to position because of the accumulation of relativistic particles in the emission region. This is the reason why the particle energy density at the AD position is about a factor 35 lower than in the steady state case. This accumulation of particles continues through the other position, increasing the particle energy density along the way until the jet position, where the previous steady-state value is obtained. Similarly, the ratio of energy densities is initially very large and decreases on the way out.

The neutrino spectra shown in Figure 5 indicate as well that the interactions and distributions require time to unfold. While at the AD position lots of neutrinos are produced, their flux is a factor of a few below the steady-state flux. At the BLR position the flux reduction is almost an order of magnitude (similar to the pion flux), while it is closer to the steady-state flux at the DT position. At the "jet" position, the neutrino flux is reduced a lot, as in the steady-state case.

**Figure 5.** Same as Figure 1 but for a moving blob. In each panel, the displayed time is the time that has passed in the comoving frame since the launch.

**Figure 6.** Same as Figure 2 but for a moving blob as in Figure 5. In each panel, the displayed time is the time that has passed in the comoving frame since the launch.

**Figure 7.** Same as Figure 3 but for a moving blob as in Figure 5.

#### **4. Discussion and Conclusions**

The results of the toy study presented in this paper clearly show the importance of the external fields in case of the presence of relativistic protons in the jet. Their influence on the particle evolution is significant resulting in very different steady-state SEDs at different positions in the jet. Especially at locations within the BLR, the familiar double-humped SED structure is destroyed. At the DT position, the spectrum is already comparable to "standard" blazar SEDs, while the "jet" position outside the external fields provides the cleanest separation between the low-energy and the high-energy bump.

The situation changes entirely when the motion of the emission region is taken into account. The relatively long source time scales (particle and photon accumulation, interactions, escape) compared to the fast speed imply that the external conditions change too

fast for the emission region to adapt even until the edge of the DT. Only on "jet" scales is the previous steady state fully recovered. This, of course, is a consequence of the choice of Γ = 50, which is a rather extreme value. Lower values on the order of Γ ∼ 10 could change the situation—especially as it would also significantly reduce the energy density of the external fields within the emission region. Steady-state solutions might be achieved at positions much closer to the black hole. Testing this, and the other potential changes to the model parameters as described above, is however beyond the scope of this paper.

Within the model parameters used in this toy study, the production of neutrinos depends strongly on the external fields with practically none produced at the "jet" position. While different parameter sets of the emission region might produce better SED shapes at positions within the external photon fields, it corroborates the results obtained by other authors [18,35,36], which makes it difficult to reconcile the neutrino and photon observations within a one-zone model.

To conclude, the production of neutrinos in a blazar jet in reasonable quantities remains a challenge, as the requirement for a reasonably dense soft photon field—in order to produce the required pions—also supports the pair cascade through *γ*-*γ* absorption and Bethe-Heitler pair production. The intrusion of a gas cloud or a star into the jet [37,38] might provide sufficient numbers of cold protons for direct proton-proton interactions [39], but the consequences (efficiency of the process, developing pair cascade, etc.) would also need further studies.

**Funding:** The author acknowledges postdoctoral financial support from LUTH, Observatoire de Paris.

**Data Availability Statement:** The *OneHaLe* code is still under development and is therefore not yet meant for public use. However, the code can be shared upon reasonable request to the author.

**Acknowledgments:** I am eternally grateful to Reinhard Schlickeiser for his supervision and guidance from the Bachelor thesis up to the Ph.D. thesis and beyond. Thank you very much for all the guidance along the way—and for all the travel destinations that you have made possible. Because of you, I realized quickly that this is indeed what I want to do—being a researcher, I mean. I am also grateful for stimulating discussions with Anita Reimer, Andreas Zech, Catherine Boisson, Markus Böttcher, Anton Dmytriiev and Chris Diltz, which helped in creating and improving the code. I would like to thank the referees for constructive reports that helped a lot to improve the paper. Simulations for this paper have been performed on the TAU-cluster of the Centre for Space Research at North-West University, Potchesftroom, South Africa.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

