Studying the influence of external photon fields on blazar spectra using a one-zone hadro-leptonic time-dependent model

The recent associations of neutrinos with blazars require the efficient interaction of relativistic protons with ambient soft photon fields. However, along side the neutrinos gamma-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.


Introduction
The theory of blazar emission was transformed in the early 1990's 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 Prof. Reinhard Schlickeiser and collaborators employing the accretion disk (AD) as a source for soft external photons [1][2][3][4].
Blazars, a sub-class of active galaxies, are indeed peculiar objects with -in the words of Prof. Schlickeiser [5] properties [that] include high optical polarization, extreme optical variability, flatspectrum 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.
And 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 is not yet reached.
The low-energy component of the famous 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) 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][12][13][14][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 -will be introduced (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. 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 comoving frame of the emission region or invariant.

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.
We assume a spherical emission region with radius R located a distance z 0 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 δ = [Γ(1 − β Γ cos θ obs )] −1 . Here, β Γ = The Fokker-Planck equation governing the time-dependent evolution of a given particle species i (protons, charged pions, muons, or electrons) with spectral density n i (χ) is given as For numerical reasons, we use the normalized particle momentum χ = p i /(m i c) = γβ, where m i is the particle mass, c the speed of light, γ the particle's Lorentz factor, and β = 1 − γ −2 . The first term on the right-hand side of Eq. (1) describes Fermi-II acceleration through scattering of particles on magnetohydrodynamic waves. We use the parametrization of [20] with a = 9v 2 s /4v 2 A , v s/A the shock and Alfvèn speed, respectively, and the energyindependent acceleration time scale t acc . This parametrization approximates the momentum diffusion through hard-sphere scattering.
The second term on the right-hand side of Eq. (1) provides momentum changesχ i 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 1 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.
The remaining three terms on the right-hand side of Eq. (1) mark the injection of particles, the escape of particles from the emission region, and the decay of unstable particles, respectively. t * i,decay is the proper decay time scale, which is 2.6 × 10 −8 s for charged pions, and 2.2 × 10 −6 s for muons, respectively. As neutral pions decay after 2.8 × 10 −17 s into γ rays, we do not solve Eq. (1) for neutral pions, but calculate their radiation output directly from their injection spectrum.
While we consider Fermi-I and II acceleration terms, 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, we use a simple power-law injection with spectral index s i between a minimum and maximum Lorentz factor, γ min,i and γ max,i , respectively. The injection normalization Q 0,i (t) is given by 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. We emphasize again that Eq. (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.
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: with the emissivity j ν , the photon escape time scale t esc,ph = 4R/3c, 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 , we can calculate the spectral luminosity in the observer's frame 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. We consider the important cases of absorption in the BLR and DT following the prescription in [28].
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 distanceR BLR andR DT from the black hole, and their energy distribution is given through a grey-body spectrum of temperatureT BLR andT DT normalized to a luminosity ofL BLR andL 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 Eq. (1), do not vary by more than 10 −4 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. 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. 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 0 = 1 × 10 16 cm), within the BLR (z 0 = 1 × 10 17 cm), within the DT (z 0 = 1 × 10 18 cm), and outside the external fields (referred to as "jet", z 0 = 1 × 10 19 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.

Influence of the external fields
The result is shown in Figs. 1, 2 and 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 Fig. 4 left), however we do not show the 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). 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. 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 proton-photon interactions is very strong (Fig. 2), as indicated by the cooling time scales being dominated by pion production (indicated by the "pion" and "neutron" loss channels) at Lorentz factors γ > 10 5 . 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 ( Fig. 1) through the neutral pion bump at PeV energies, results in a significant production of muons and highly relativistic electrons ( Fig. 3) with Lorentz factors γ > 10 10 . 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 Fig. 4, and manifests itself in Fig. 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 (Fig. 3). The resulting electron synchrotron flux (Fig. 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 (Fig. 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 Eq. (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 (Fig. 1). While pion and Bethe-Heitler pair production are reduced compared to the AD case, the pair cascade is still very significant (Fig.  3) because of γ-γ pair production (Fig. 4 left). While the process is less severe than in the AD case, the secondaries still dominate the electron distribution (Fig. 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 ( Fig. 3 and Fig. 4 left) and the more familiar double-humped SED emerges (Fig. 1). At 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 (Fig. 2).
Lastly, the emission region is located outside the external photon fields in the "jet" case. While secondary pairs are still being produced (Fig. 3), their number is low (Fig. 2) due to low absorption (Fig. 4 left), and their flux contribution only shows around photon energies on the TeV-scale, but at relatively low flux values (Fig. 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 Tab. 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.
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 3 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 Lorentz contraction:ẑ = z 0 + Γβ Γ ct, where t is the time since launch in    the comoving frame, andẑ 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  Figure 5. Same as Fig. 1, but for a moving blob. In each panel, the time in the comoving frame is given that has passed since the launch.
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 the same distances as in the steady-state cases 2 . The results are shown in Figs. 5, 6, and 7, while the optical depth due to γ-γpair production is shown in the right panel of Fig. 4. We note that 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 protonphoton interactions producing a significant amount of pions (Fig. 6), which decay into photons or muons and pairs. In fact, Fig. 7 shows that the injection of pairs from muon decay is almost at the level as in the steady state (Fig. 3), but Bethe-Heitler produced pairs are about 2 orders of magnitude below. Similarly, γ-γ pair production is below the steady-state level, because the internal photon fields (Fig. 5) have not yet been fully developed. In turn, the optical depth due to γ-γ pair production (Fig. 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 at PeV photon energies allowing for a very strong flux in the neutral pion bump (Fig. 5). The p loss (tot) p loss (pion) p loss (neu) p loss (bethe-heitler) e loss (tot) e loss (ic) escape acceleration Figure 6. Same as Fig. 2, but for a moving blob as in Fig. 5. In each panel, the time in the comoving frame is given that has passed since the launch.
"under-development" of the internal photon fields is a consequence of the low electron and proton densities (Fig. 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 -though orders of magnitude below the steady-state case. The situation only changes mildly until the BLR position is reached after 5.7 × 10 4 s. This is still less than the escape times of photons (2 × 10 5 s) and particles (7.5 × 10 5 s). Therefore, particle and photon densities continue to increase. The spectral shape of the SED shown in Fig. 5 is somewhat similar to the steady-state case (Fig. 1), but at a factor of a few reduced in flux. There are a few more details where SEDs differ. In the γ-ray domain, Fig. 5 shows contributions from IC scattering of both the AD and the BLR. Comparing the IC/AD spectra of the top panels in Fig. 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 much reduced 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 (Fig. 6 vs. Fig. 2). These are a consequence of the reduced γ-γ pair production at lower energies (Fig. 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 Fig. 5 -a reduction of about an order of magnitude compared to the steady state.
After 5.7 × 10 5 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 Fig. 5 is almost equal to the steady-state case (Fig. 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 proton-synchrotron 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 Fig. 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 (Fig. 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 × 10 6 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 Fig. 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 much reduced as in the steady-state case.

Discussion and conclusion
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, the previous steady state is 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 situationespecially 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] that it is 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 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 Prof. Dr. 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 destination 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, 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.