**1. Introduction**

Active galactic nuclei (AGN) are among the most enigmatic objects in the universe. With luminosities in excess of 1037 W (1044 erg s−1), they represent the most luminous, continuous sources of radiation. With their central supermassive black holes (SMBHs), they provide an environment that can help us to understand black holes at work. The question of how energy is transferred from the black hole and/or the accretion disk to launch gigantic radio jets is still largely unsolved and subject to ongoing research. AGN also constitute one of the few sites in the universe that provide enough energy on total to serve as a candidate for the observed flux of ultra-high energy cosmic rays (UHECRs) and might be a key to understand particle acceleration up to the highest observed particle energies; see, e.g., [1,2] for summaries. As one of the very few source classes, active galaxies provide an astrophysical, extreme environment which might be suited to accelerate particles to an incredible amount of 1020 eV (see, e.g., [3]). These sources are therefore also considered to contribute to the diffuse astrophysical high-energy neutrino flux as measured by IceCube

**Citation:** Becker Tjus, J.; Hörbe, M.; Jaroschewski, I.; Reichherzer, P.; Rhode, W.; Schroller, M.; Schüssler, F. Propagation of Cosmic Rays in Plasmoids of AGN Jets-Implications for Multimessenger Predictions. *Physics* **2022**, *4*, 473–490. https:// doi.org/10.3390/physics4020032

Received: 31 January 2022 Accepted: 1 April 2022 Published: 28 April 2022

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

**Copyright:** © 2022 by the authors. 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/).

at Earth [4]. In particular, the sub-class of blazars is known for their strong short- and longterm variability, especially (but not exclusively) at gamma-ray energies. In the literature, blazar jets are often discussed to be dominated by leptonic particle processes, which fit observational quiescent data of blazar spectral energy densities (SEDs; see [5–7]). Yet, the understanding of the complex, time-variable structure of the SEDs is still far from being understood in full detail. Additionally, the picture of an electron-positron plasma in the jet is not the only possibility, and an electron–proton (hadron) plasma is certainly a realistic option. Thus, in the past decades, AGN jets and particularly blazar emissions have also been argued to naturally contain hadronic components, which would not only contribute to the emission of electromagnetic radiation in blazars [8], but also lead to the production of secondary high-energy neutrinos [9–11]. The general detection of up to PeV high-energy neutrinos of astrophysical origin observed by IceCube [4] has started to shed more light on the non-thermal high-energy universe. From the distribution of events in the projected sky, it is clear that the flux is not focused in the galactic plane, and that it is therefore likely to contain a significant extragalactic component, see, e.g., Becker Tjus and Merten [12] for a summary. In the past few years, several possible associations of neutrinos with astrophysical objects have been identified. Each of these is at a ∼3σ confidence level so far and serves as first evidence. A first hint for an association of a high-energy neutrino with a blazar comes from the source TXS 0506+056. In September 2017, a neutrino of ∼300 TeV could be associated with an exceptional gamma-ray flare at GeV energies from this source. This gamma-neutrino correlation can be estimated to have a significance of 3σ by combining data of the IceCube neutrino observatory and Fermi-LAT [13]. This detection initiated a dedicated search for neutrino clustering in the ∼10 years of IceCube data around the position of TXS 0506+056, with the result that there was an enhanced flux of neutrinos in a half-year period from September 2014 to March 2015, also with a statistical significance of ∼3σ of being incompatible with the background hypothesis [14].

In order to explain the neutrino signatures detected from the direction of TXS 0506+065, hadronic jet models have been applied (e.g., [15–19]). The two potential flares appear quite different in their evolution: the neutrino signal above the atmospheric background detected in 2014/2015 lasted about 100 days and consisted of about 8–18 neutrinos with energies of 10–100 TeV, and the gamma-ray light curve at GeV is in a minimum state [14]. The 2017 detection is based on one single high-energy neutrino with extreme energy (∼300 TeV). A gamma-ray flare was observed in coincidence with the arrival of the neutrino. It has been noted by Kun et al. [20], however, that at the exact time of the neutrino detection, even here the GeV gamma-ray flux was at a local minimum and only rose to a high emission state shortly after the neutrino detection. It was argued in [20] that this observation is consistent with a model for which the neutrinos are produced in a high-density medium in which gamma-rays are absorbed, which either becomes less dense with time or for which the gamma-rays take some time to cascade down to GeV energies before escaping.

By now, there are tens of possible associations of neutrinos and AGN jets [21–23]. With IceCube in continued operation, this number will further increase in the upcoming years, with the expectation to finally confirm at least some of these sources at the >5 sigma level to be neutrino emitters. One common conundrum in all of the different high-energy neutrino detections (diffuse and potential sources) is the lack of TeV, but even GeV gamma-ray emission. Hadronic interactions that are responsible for neutrino production inevitably lead to the co-production of high-energy gamma rays. The reason is that neutrinos are produced from the subsequent decay of charged pions and kaons, which in turn are produced in a fixed ratio of neutral pions and kaons, leading to the production of gamma rays with an energy threshold at the mass of the pion, *Eγ*, min = *mπ*<sup>0</sup> *c*2/2 = 70 MeV, with *c* being the speed of light. Comparing the detected diffuse neutrino flux with the measured extragalactic, diffuse component of gamma rays leads to the conjecture of a source environment that must absorb the gamma rays at energies >GeV [24,25]. The absorption of photons can be due to a strong accretion disk, as it has been discussed in, e.g., Brodatzki et al. [26] for the case of TeV emission. The potential neutrino source fluxes

from the 2014/2015 TXS signature also indicates that, in order to match the observed gamma-ray flux, it must be diminished significantly at >GeV energies to make the neutrino production model work. Such a model of gamma-ray absorption can even be used as a tracer in the searches for associations of high-energy neutrinos with blazars. That is, rather than searching for an enhanced gamma-ray flux, the neutrinos can actually arrive at times of reduced gamma-ray activity [20]. Such a scenario can be produced for regions of extreme gas or photon densities. In the first case, the photons will interact with the dense gas via Compton scattering. In the second case, gamma–gamma interactions will lead to electromagnetic cascades. It is clear that in both scenarios, the energy of the highenergy gamma rays will be deposited at much lower energies. If these environments only become transparent at MeV energies as suggested by, e.g., Halzen et al. [15], this is not observable now, as a dedicated mission for MeV detection of the gamma-ray sky is currently missing. Future missions like e-Astrogam, MeVCube, or AMEGO will shed more light on these questions. At this point, the theoretical models are being challenged by GeV-TeV measurements, which indicate that there is no significant increase in the energy output connected to the neutrinos.

In general, the modeling of the steady-state emission, but even more the modeling of the flares is complex and requires a complete consideration of the jet physics, including different scenarios for the acceleration region, gas and photon targets, as well as the magnetic field structure. The latter is highly important for the proper modeling of the diffusive cosmic-ray transport, which is relevant for the evolution of the flares, both for the leptonic and hadronic signatures [19]. Quantitative theoretical modeling is necessary to establish a physical connection of the neutrinos to the blazars. Mere directional coincidence is not enough, because the angular uncertainty of the neutrino events is larger than 1◦, making source identification difficult without theoretical input.

Models of high-energy neutrino and electromagnetic up to gamma-ray emission in the jets of AGN cover a variety of scenarios and parameter spaces. Blazars are known to be highly variable across the electromagnetic spectrum, with a variety of models put forward to explain these flares. Such models include, among others, particle acceleration via internal shocks in the jet (see, e.g., [3,9,10,27]) and reconnection driven plasmoids (see, e.g., [19,28,29]). These latter relativistic and compact structures have been discussed in the literature since the 1960s [30,31]. What here is referred to as *plasmoid* is often called *blob* in the literature, meaning compact, dense structures traveling with relativistic speeds along the jet axis. The term "plasmoid" is preferred here, as it is typically used in the context of the plasmoid creation via magnetic reconnection events. In this scenario, the injection of a relativistic plasma into the system (here the AGN jet) can lead to reconnection events that, under certain circumstances, lead to the plasmoid instability that breaks down the streaming plasma into small blobs, i.e., the plasmoids. In this scenario, charged particles can be pre-accelerated in the reconnection events. While non-relativistic reconnection is limited to below-knee energies [32], relativistic reconnection can be much more efficient [33], also by further acceleration via a Fermi second-order process when the particles scatter in between the plasmoids. Such an acceleration scenario can solve the long-standing injection-problem. They also justify the assumption that the cosmic-ray population is distributed homogeneously in the plasmoid, as the turbulent field in the plasmoid is used to isotropize the direction of the incoming particles. The assumption of a homogeneously distributed population is implicit in those models that do not resolve the blob, but work with timescales. In test-particle simulations, it is a reasonable approach to start with a homogeneous distribution, as done in this paper.

The modeled hadronic component of proton-proton interaction was discussed by Eichmann et al. [27]; the modeling of leptonic and (lepto-)hadronic emission by Christie et al. [34], Keivani et al. [35] for the case of TXS 0506+056. Other flare models based on external factors include gas clouds entering the base of jets [36–38] and jets of former binary AGN drilling through their own dust tori and/or accretion disks after being redirected by the merger of their host black holes [39]. Such scenarios of high density all tend to be more hadron-dominated due to the nature of their occurrence.

Figure 1 shows a sketch of an AGN with a jet. Also shown are the photon and/or gas targets (yellow/red). The structure of the large-scale magnetic field is shown in blue. A turbulent component of the magnetic field exists as well and is not drawn in the figure, but only indicated by purple text. The structure of the gas/photon fields is highly relevant for the particle interactions and therefore needs to be included in the models in three dimensions. The same is true for the magnetic field structure as the synchrotron radiation is sensitive to the direction of the field. For a diffusive transport description, it is also highly relevant as it is often assumed that the propagation is dominantly along the field lines with a smaller component perpendicular to the field. In fact, the perpendicular diffusion coefficient can even dominate the picture if the magnetic field turbulence level is *δB*/*B* > 1, something that is certainly possible in these extreme environments.

**Figure 1.** Structure of an active galactic nucleus (AGN) with a jet. Yellow/red components are targets for cosmic-ray interactions with gamma rays or gas. Blue colors indicate the likely structure of the magnetic field along the jet. Particle acceleration happens either at the shock fronts or in context with the relativistic plasmoids.

Current state-of-the-art of numerical codes include many of the necessary features. The codes and their most important properties are summarized in Table 1. Cerruti et al. [40] compared the codes for blazar hadronic models (with the exception of the CRPropa code). The models are typically designed to numerically solve the transport equation including loss processes from which the secondary particle radiation from electrons and protons can be calculated. A loss term for the particle transport is usually included via a timescale, i.e., a term −*n*/*τ*, where *n* is the particle density and *τ* is the characteristic timescale. In [8,41,42], the timescale is chosen to be the ballistic one, *τ* = *c*/*R* (where *R* is the radius of plasmoid), thus being energy-independent. Gao et al. [43] has argued that the propagation is of diffusive nature so that the escape time of particles is assumed to be a factor of 10 times longer than the ballistic one, but still assumed to be energy-independent. However, there are models that implement the escape of particles with a diffusive, energy dependent approach Rodrigues et al. [44] modeled the particle escape by distinguishing between the two most extreme cases, i.e., escape via diffusion and escape via advection and argue that the physical escape spectrum will emerge in between these two scenarios. All of these models are designed to model propagation and particle interaction in blazars. Due to the strong variability of the sources, it can be deduced that the signal must come from the very compact region on the order of 1014 m. This makes the plasmoid model favorable over an approach of shock acceleration and explains why all codes focus on such an approach.


**Table 1.** Basic properties of state-of-the-art blazar propagation codes. A detailed comparison of the codes in this table for blazar hadronic models was presented by Cerruti et al. [40], with the exception of CRPropa code.

A new code for propagation of particles in relativistic plasmoids has been developed recently [19]. This new framework has been derived from the public transport code CR-Propa 3.1 [45]. The advantage with this numerical approach is that ballistic propagation can be performed in a test particle approach, thus not relying on the simplifying assumption of time-scales. Further, a second transport framework is integrated in CRPropa 3.1, which is the solution of the transport equation via the stochastic differential equations (SDEs) approach. This approach enables us to solve the transport equation via pseudo-particle propagation, which makes it compatible to be used in the ballistic test-particle propagation of CRPropa. The SDE approach is designed to include a full diffusion tensor, which is also an improvement when compared to other codes. A CRPropa modification presented in Hörbe et al. [19] makes use of the modular structure of CRPropa to create a propagation environment of a plasmoid traveling along the jet axis. The photon field of a thin accretion disk is implemented at the foot of the jet for gamma–gamma and proton–gamma interactions. Technically, the propagation is done in the reference frame of the plasmoid and then transferred into the observer's frame. The plasmoid itself contains a plasma with a constant density *n*plasma, which is considered as a target for cosmic-ray interactions as well. The magnetic field, in which the particles propagate, was assumed to be of purely turbulent nature of Kolmogorov type in [19]. Due to the modular structure of the code, it can easily be changed to include a regular component as well, to change the nature of the turbulence, etc.

In this paper, the spotlight is put on the propagation regime in the plasmoids of blazars. In Section 2, the energy, at which a transition between diffusive and ballistic propagation happens and the consequences such a transition has for the description of SEDs and lightcurves of blazars are quantified. In Section 3, the influence of a first phase of ballistic propagation before the limit of diffusion is reached in time is investigated and the necessity to go from a diffusive approach to the description via the telegraph equation is discussed. In Section 4, first test simulations to investigate a possible difference in the flaring behavior in the diffusive vs. ballistic description is performed. Conclusions and outlook are given in Section 5.

#### **2. The Space-Domain: Diffusive vs. Ballistic Propagation**

As discussed above, propagation of charged particles in a turbulent (plus regular) magnetic field can be of fundamentally different nature depending on the astrophysical setting, in particular concerning the parameters of the particle energy, *E*, the ratio, *δB*/*B*, of the turbulent to regular magnetic field, and the correlation length, *lc*, of the field as the lower boundary for the deterministic description of the magnetic field lines. Reichherzer et al. [46] quantified five propagation regimes with respect to the particle's reduced rigidity *ρ* = *rg*/*lc*, with *rg* = *E*/(*cqB*) as the relativistic gyro radius and *lc* ≈ *l*max/5 as the correlation length. Here, *l*max = 2*π*/*k*min is the maximum scale of the magnetic turbulence spectrum, connected to the lowest wave number *k*min, defined by the turbulence injection scale. Diffusive propagation corresponds to the resonant-scattering regime

(RSR). This regime is only valid for particles that can scatter with the entire spectrum of wavelengths *k*min < *k* < *k*max, where *k*max = 2*π*/*l*min is the dissipation scale. The general scheme of resonant scattering then breaks down toward the lowest and highest reduced rigidities. At the lower boundary, when scattering does not happen with the entire angular spectrum anymore, mirroring occurs more and more often, altering the diffusion coefficient. It is expected that this regime is not relevant for particles in AGN plasmoids, as the gyro radius of ∼TeV-PeV particles that are considered here in magnetic fields of ∼G strength fulfills the boundary condition *ρ* > *l*min/(*π lc δB*/*B*) [46]. The situation is different toward large reduced rigidities, for which particles propagate in a *quasi-ballistic* way: for gyro radii that start to reach the correlation length of the system, meaning for plasmoids also coming closer to the actual size of the system, only a few gyrations are performed by the particles before leaving the source. This is happening close to the Hillas limit of the source. It is clear that the number of gyrations then does not suffice for a diffusive description.

The quasi-ballistic regime becomes relevant at a reduced rigidity of *ρ* = *rg*/*lc* 5/(2*π*) [46]. Inserting the relativistic gyro radius, the energy at which the ballistic regime becomes dominant is given as

$$E \gtrsim \frac{5}{2\pi} \cdot l\_c \cdot c \cdot q \cdot B \,. \tag{1}$$

For a given parameter set of magnetic field strengths, coherence lengths of the turbulence, and particle energies, Equation (1) can be applied to determine in what regimes the particles are propagating in typical astrophysical sources of cosmic rays [47]. Normalized to a standard set of parameters, the equation becomes

$$E \gtrsim Z \cdot \left(\frac{l\_c}{10^{11}\,\mathrm{m}}\right) \cdot \left(\frac{B}{0.42\,\mathrm{G}}\right) \cdot 10^{15}\,\mathrm{eV} \,. \tag{2}$$

This implies that for protons (the atomic number *Z* = 1) in a source region with a parameter set *lc* = 10<sup>11</sup> m and *B* = 0.42 G, diffusive propagation is happening below energies of 1015 eV, ballistic propagation needs to be applied above 1015 eV. For other parameter combinations, this transition energy can be calculated accordingly, always with ballistic propagation above the energy, diffusive transport below.

Figure 2 shows the energy limits for protons as a function of the product *B* · *lc*. The grey shaded area illustrates diffusive and the blue area ballistic propagation, with the transition between the resonant scattering regime and the quasi-ballistic regime indicated as the solid line in between, following Equation (1). The area of ballistic propagation in blue is bounded by the maximum possible proton energies according to the Hillas-Limit in each parameter space. The horizontal lines indicate the energies for the knee (1015 eV), ankle (1018.5 eV), and maximum observed energy (1020 eV).

The parameter space covered by the plasmoids is approximated to be in the range 1010 m < *lc* < 10<sup>14</sup> m. This range is based on the assumption that the plasmoids have a radius on the order of *<sup>R</sup>* <sup>∼</sup> <sup>10</sup>12–1016 m, using a correlation length of *lc* <sup>=</sup> 0.01 · *<sup>R</sup>* as described above. As the plasmoids are launched at the foot of the jet, magnetic fields are large, on the order of 10−<sup>3</sup> G < *B* < 10 G. What we want to understand is how the propagation of particles needs to be performed to describe the multimessenger emission from AGN in the plasmoid-model. The energy range of interest for high-energy photons reaches from GeV energies up to approximately 1016 eV, neutrino detection happens in an energy range corresponding to proton energies of approximately 2 · 1013 eV to 1017 eV. Figure 2 shows the relevant parameter space, displayed as *B* · *lc* on the x-axis, a fraction will be diffusive (grey area) at lower energies. The high-energy part before reaching the Hillas limit (colored, thick line) needs to be performed in the ballistic limit (blue area). A first extreme example is a combination *<sup>B</sup>* · *lc* <sup>=</sup> 108 G m, where diffusive propagation happens up to <sup>∼</sup>1012.5 eV, ballistic propagation up to the Hillas limit at around 1014.5 eV. These would be sources with a relatively low acceleration limit, as the combination of *R* and *B* only allows for maximum energies below the knee. More realistic parameter combinations that would allow the sources to reach the maximum observed energy would

be a combination of *<sup>B</sup>* · *lc* <sup>=</sup> <sup>10</sup><sup>14</sup> G m. In this case, diffusive propagation needs to be performed up to 1018.5 eV, ballistic propagation needs to be performed up to the Hillas limit at 1020.5 eV.

**Figure 2.** Illustration of the transition between the quasi-ballistic regime (in blue) and the resonant scattering regime (in grey) in the dependence of the magnetic field strength, *B*, correlation length of the magnetic field, *lc*, and particle energy, *E*. In a simplified approach, it is assumed here that particles in the quasi-ballistic regime propagate ballistically and in the resonant scattering regime diffusively. The position of the diagonal line represents the transition energy from diffusive (below) to ballistic (above) for protons (the atomic number *Z* = 1), determined by Equation (1). Diffusive propagation is expected for the result of a parameter combination below the line, while ballistic propagation lies above the line.

This result has immediate consequences on the observed energy spectra: a break in the energy behavior of the timescales applied in simplified transport equation approaches, where the term −*n*/*τ*esc describes the escape, is expected to be observed: For ballistic transport, the timescale needs to be chosen as a constant value, *τ*ballistic esc = *R*/*c*, while it becomes energy dependent in the case of diffusive propagation, *τ*diffusive esc <sup>=</sup> *<sup>R</sup>*<sup>2</sup> <sup>√</sup> /*κ* ∝ *R E*−*δ*, with *<sup>κ</sup>* <sup>=</sup> *<sup>κ</sup>*<sup>0</sup> · *<sup>E</sup><sup>δ</sup>* as the diffusion coefficient, for which the energy dependence can be approximated as a power-law behavior with an index *δ* that depends on the underlying magnetic turbulence, in this description being in the limit of quasi-linear theory, i.e., *δB*/*B* 1. Such a change in the timescale directly induces a change in the shape of the spectral energy distribution: applying the leaky box model, the emitted proton spectrum follows approximately *n*(*E*) ∼ *Q*(*E*) · *τ*esc(*E*). Those secondary photons and neutrinos that are induced by (photo-)hadronic interactions basically mirror that behavior in the energy range above the threshold for the process, so that even these are in first approximation proportional to the escape time and the primary injection spectrum *Q*(*E*), *nγ*,*<sup>ν</sup>* ∝ *Q*(*E*) · *τ*esc. Assuming an injection spectrum *Q*(*E*) ∝ *E*−2.3 and Kolmogorov-type turbulence, *τ*diffusive esc ∝ *E*−0.3, the SED is expected to behave as

$$m\_{\gamma,\nu}(E) \approx \begin{cases} \begin{array}{c} E^{-2.6}\_{\gamma,\nu} \\ E^{-2.3}\_{\gamma,\nu} \end{array} & E < E^{\gamma,\nu}\_{\text{transitive}} \\ \begin{array}{c} E^{-2.3}\_{\gamma,\nu} \end{array} & E > E^{\gamma,\nu}\_{\text{transitive}} \end{array} \tag{3}$$

That means for a typical parameter combination of *lc* = 10<sup>11</sup> m and *B* = 0.42 G, the transition energy for protons is *E*CR transition = 10<sup>15</sup> eV. This translates (see, e.g., [1]) into a transition energy for photons of *E<sup>γ</sup>* transition <sup>≈</sup> 1/10 · *<sup>E</sup>*CR transition <sup>≈</sup> 1014 eV. These 100 TeV photons are currently barely accessible for gamma-ray telescopes, and so the propagation in a purely diffusive regime is reasonable. Purely ballistic transport, however, leads to a spectrum that is too flat, as the steepening by the diffusive escape timescale due to the diffusion tensor is neglected.

For neutrinos, the transition energy is *E<sup>ν</sup>* transition <sup>≈</sup> 1/20 · *<sup>E</sup>*CR transition <sup>≈</sup> <sup>5</sup> · <sup>10</sup><sup>13</sup> eV. As the neutrinos detected by IceCube are in the energy range 10 TeV to a few PeV, this transition region can fall right into the relevant parameter range, so that a combination of ballistic and diffusive propagation needs to be considered. A break in the observed neutrino energy spectrum from a steeper to a flatter behavior is therefore expected in such a scenario. If such a break is observed, it can be used to estimate the parameter combination *B* · *lc*. To summarize the above result in the context of the propagation of particles in the plasmoids of blazars, it is shown that it is of high importance to evaluate the transport regime and adjust it accordingly to the problem under consideration in order to receive reliable results.

#### **3. The Time-Domain**

The result from the previous section applies to steady-state sources with *δn*/*δt* ≈ 0, where it is implicitly assumed that all particles have already had the time to reach a steadystate limit in their propagation; here, *t* denotes the time. However, blazars are highly variable objects and individual flares are often modeled by injecting a high-energy particle population on short timescales. Particle acceleration time-scales in reconnection events responsible for the blob creation in relativistic sources are on order of *<sup>τ</sup>*acc <sup>∼</sup> *<sup>E</sup>*/(*qBc*2), which is typically much shorter than the escape timescales (see, e.g., [48]), motivating a two-zone scenario where acceleration is typically performed in a first step, propagation afterwards.

While there is scientific consensus that the description of the transport process of particles in turbulent fields is ensured by the general concept, discussed in Section 2 in the limit of infinitely large times, the question arises under which conditions and on which timescale such a limit consideration is appropriate. In this section, criteria are derived for which, in a given plasmoid setup, the diffusive approximation still holds.

The general problem of assuming diffusive propagation during the initial propagation process, for which the diffusive limit is not yet granted, is expressed in the following points:


#### *3.1. Timescale for Transition to Diffusive Propagation*

Whereas the consideration of particle transport via the diffusion equation and its solution of a Gaussian particle distribution cannot distinguish between the initial, ballistic propagation and the subsequent diffusive propagation, the telegraph equation has recently been attributed this ability [49–52]:

$$\frac{\partial f}{\partial t} + \tau \frac{\partial^2 f}{\partial t^2} = \kappa \left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} \right), \tag{4}$$

where *x*, *y* and *z* are space coordinates.

The telegraph timescale *τ* describes the transition between these two propagation phases. This timescale enables us to make a statement about when the diffusive phase is established and when the description of the particle transport via the diffusion equation is sufficiently accurate. If the initial ballistic phase is neglected, *τ* disappears and the telegraph equation turns into the known diffusion equation. By neglecting adiabatic focusing, the telegraph timescale yields [49]:

$$\tau = \frac{3v}{8\lambda\_{\parallel}} \int\_{-1}^{1} \mathrm{d}\mu \left( \int\_{0}^{\mu} \mathrm{d}\mu' \frac{1-\mu'^2}{D\_{\mu\mu}(\mu')} \right)^2 \,\mathrm{}^{\prime} \tag{5}$$

with *μ* being the cosine of the pitch angle and *Dμμ* being the pitch-angle Fokker-Planck coefficient that, in a negligible background field (*δB B*0), yields [53]

$$D\_{\mu\mu} = (1 - \mu^2)D. \tag{6}$$

Here, *D* is the pitch-angle Fokker–Planck coefficient at 90◦. The mean-free path *λ* is defined as

$$
\lambda\_{\parallel} = \frac{3v}{8} \int\_{-1}^{1} \mathrm{d}\mu \frac{(1-\mu^2)^2}{D\_{\mu\mu}(\mu)} \,\mathrm{}.\tag{7}
$$

where *v* denotes the particle velocity.

#### *3.2. Quantifying the Time Needed to Achieve Certain Levels of Diffusivity*

Since the diffusion equation assumes diffusive transport at all times and, furthermore, the solution uses the final diffusion coefficient, *κ*, over all timescales, the norm *N* of the solution remains constant in time:

$$N = 4\pi \int\_0^\infty \mathrm{d}r \, r^2 f\_{\mathrm{diff}}(r) = 4\pi \int\_0^\infty \mathrm{d}r \, \frac{r^2}{(4\pi \mathrm{tr}\,\mathrm{t})^{3/2}} \exp\left(-\frac{r^2}{4\mathrm{st}}\right) = 1\,\tag{8}$$

The norm may be interpreted as the fraction of particles participating in diffusion [52]. On the other hand, the solution of the isotropic telegraph equation yields:

$$f\_{\text{telgraph}}(r,t) = \frac{e^{-t/2\tau}}{4\pi\kappa^{3/2}} \left[ \frac{\delta\left(t - r\sqrt{\pi/\kappa}\right)}{r/\sqrt{\kappa}} I\_0\left(\frac{1}{2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa\tau}}\right) \right.$$

$$+ \frac{\Theta\left(t/\sqrt{\tau} - r/\sqrt{\kappa}\right)}{2\tau^{3/2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa\tau}}} I\_1\left(\frac{1}{2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa\tau}}\right) \right]. \tag{9}$$

Here, *δ*(...) is the Dirac's delta function, Θ(...) is the Heaviside step function, and *Iν*(...) is the modified Bessel function. The norm can be computed individually on each of the two terms of the function:

$$\begin{split} N = \int f\_{\text{telogram}}(r, t) \, \mathrm{d}r &= \int \frac{e^{-t/2\tau}}{4\pi \kappa^{3/2}} \left[ \frac{\delta\left(t - r\sqrt{\tau/\kappa}\right)}{r/\sqrt{\kappa}} I\_0\left(\frac{1}{2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa \tau}}\right) \right] \mathrm{d}r \\ &+ \int \frac{e^{-t/2\tau}}{4\pi \kappa^{3/2}} \left[ \frac{\Theta\left(t/\sqrt{\tau} - r/\sqrt{\kappa}\right)}{2\tau^{3/2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa \tau}}} I\_1\left(\frac{1}{2}\sqrt{\frac{t^2}{\tau^2} - \frac{r^2}{\kappa \tau}}\right) \right] \mathrm{d}r. \end{split} \tag{10}$$

The delta function simplifies the first part to *t*/*τ* and the second part can be solved by employing Bessel function integration rules. The details of the calculation are omitted. The norm yields

$$N = 1 - \exp\left(-\frac{t}{\tau}\right). \tag{11}$$

In this scenario, no particles are diffusive at the beginning. With time, the number of diffusive particles increases exponentially until the value for large propagation times approaches the maximum value with all particles in a diffusive state.

Rearranging the equation leads to a calculation rule for the propagation time *t*diff,*<sup>N</sup>* required to establish a certain diffusion level:

$$t\_{\rm diff,N} = -\ln\left(1 - N\right)\tau\,. \tag{12}$$

The relation (11) is shown in Figure 3 in comparison with the constant number of diffusing particles in the case of modeling the transport with the diffusion equation. The initial phase of a flare of cosmic rays is therefore in a non-diffusive state until the steadystate limit is reached as shown in Figure 3. Note that the derivations made here for a purely turbulent field also apply to the generalized case of an additional directional magnetic field component, since the timescales for reaching the diffusive phase during transport parallel and perpendicular to the directional background field are identical for a large parameter space [54].

**Figure 3.** Comparison of the time evolution of the ratio of particles that are already diffusively propagating on average. The solution of the diffusion equation leads to the fact that particles always propagate diffusively. The solution of the telegraph equation shows an increase in the diffusively propagating particles with time and an approach to the maximum value.

In Section 4, this scenario is evaluated for the conditions in a plasmoid.

#### *3.3. Conditions for Plasmoid Settings*

In what follows, the critical time to reach diffusive propagation is expressed as a function of the blob properties and the particle energy. The resulting estimates give an overview of the expected type of propagation for special parameter combinations. For this purpose, let us start with the definition of the timescale connected to the mean free path, *λ*, of the particle,

$$
\tau = \frac{\lambda}{\upsilon} = \frac{3\kappa}{\upsilon^2} \,, \tag{13}
$$

where *λ* = 3*κ*/*v* is expressed as a function of the diffusion coefficient. Here, the case of isotropic turbulence without background field is considered in which case Bohm diffusion to be applied with a linear energy dependence in the resonant-scattering regime, *κ* ∝ *E*. In the quasi-ballistic regime, the diffusion coefficient becomes *κ* ∝ *E*2. The transition between the two regimes is given at a reduced rigidity of *ρ* ≈ 5/(2*π*), so that the diffusion coefficient can be written as

$$\kappa = \kappa\_0 \left(\frac{\rho}{\rho\_0}\right)^{\epsilon} \text{ with } \begin{cases} \epsilon = 1 & \text{for } \rho \lessapprox 5/(2\pi), \\ \epsilon = 2 & \text{for } \rho \gg 5/(2\pi), \end{cases} \tag{14}$$

where *κ*<sup>0</sup> and *ρ*<sup>0</sup> are normalisation constants and *ρ*<sup>0</sup> = 5/(2*π*) is chosen, so that *κ*<sup>0</sup> becomes the value of the diffusion coefficient at the chosen value of *ρ*0.

It follows for the timescale:

$$\tau = \frac{3\kappa\_0}{v^2} \left(\frac{\rho}{\rho\_0}\right)^{\varepsilon} \text{ with } \begin{cases} \varepsilon = 1 & \text{for } \rho \le 5/(2\pi), \\ \varepsilon = 2 & \text{for } \rho \gg 5/(2\pi), \end{cases} \tag{15}$$

finally leading to an expression for the propagation time required to reach a certain diffusion level *N* (from Equation (12))

$$t\_{\rm diff,N} = -\ln\left(1 - N\right) \frac{3\kappa\_0}{v^2} \left(\frac{2\pi\rho}{5}\right)^{\varepsilon} \text{ with } \begin{cases} \varepsilon = 1 & \text{for } \rho \le 5/\left(2\pi\right), \\ \varepsilon = 2 & \text{for } \rho \gg 5/\left(2\pi\right). \end{cases} \tag{16}$$

By inserting the definition of the reduced rigidity, this relation can be expressed as functions of *E*, *B* and *l*c:

$$t\_{\rm diff,N} = -\ln\left(1 - N\right) \frac{3\kappa\_0}{v^2} \left(\frac{2\pi E}{5qcBl\_\odot}\right)^\varepsilon \text{ with}\begin{cases} \varepsilon = 1 & \text{for } \rho \lessapprox 5/(2\pi),\\ \varepsilon = 2 & \text{for } \rho \gg 5/(2\pi). \end{cases} \tag{17}$$

The time *t*diff,*<sup>N</sup>* needed for the fraction *N* of the particles to be diffusive depends on the parameters *E*, *B* and *l*c. Figure 4 shows this condition for different plasmoid parameters. The figure shows the influence of the particle energy, the magnetic field properties, and the trajectory on the fraction of already diffusively propagating particles. The vertical lines indicate the timescales required for particles on ballistic trajectories to travel one plasmoid radius. This timescale has the same order as typical escape times of charged particles during initial ballistic propagation. If there is not yet a significant fraction of diffusive particles at the vertical lines for the respective plasmoid radii, the particles must be considered completely ballistic. For example, charged particles with *E* = 1 TeV, *B* = 1 G, and *l*<sup>c</sup> = 1010 m can be treated diffusively in plasmoids with *R* = 10<sup>16</sup> m, but must be treated via equation-of-motion approaches at smaller radii such as *R* = 10<sup>14</sup> m, *R* = 1012 m, and *R* = 1010 m.

**Figure 4.** Fraction of particles, *N*, that are diffusive as a function of propagation time, *t*, for different blob parameters and particle energies. The vertical lines illustrate the time required for ballistic particle propagation to traverse the respective blob radii.

#### **4. Simulations: Ballistic vs. Diffusive Simulation Results**

In this Section, the effect of ballistic vs. diffusive propagation is investigated by performing simulations of cosmic-ray transport in the plasmoid of is applied an AGN traveling along the jet axis. Here, interactions with photon and gas targets are considered to be switched off in order to focus on the effects coming from cosmic-ray propagation, but otherwise follow the procedure described in [19]. This is done as a first step to understand the influence of the change in the propagation regime, so that in future studies, one can calculate the neutrino and gamma-ray spectra, adding the effect of hadronic cascades. This allows us to disentangle the two effects of the change in the flare due to propagation vs. interaction. For now, focus is on the effect of the change in propagation regime, future work will then also quantify changes through interactions. The parameter set used in the simulation is summarized in Table 2. In addition, not only simulations with the equationof-motion [55] are performed, but diffusive propagation applied by using the module DiffusionSDE in CRPropa 3.1 [45], which solves the transport equation with a diffusion term *κ*Δ*n*. In order to have a quantitative comparison, first, one needs to calculate the diffusion coefficient *κ* as an input to the diffusive simulation from the ballistic part. This is done here for energies from 105 GeV up to 108 GeV. Here, the Taylor Green Kubo (TGK) formalism is applied as described in [46] (see also references therein). The simulation parameters from Table 2 are chosen to minimize numerical errors such as the interpolation effect of turbulence [56,57].


**Table 2.** Simulation parameters, given in the rest frame of the plasmoid. All simulations are performed in a modified version of CRPropa 3.1 as presented in Hörbe et al. in [19] with further additions made for this paper as described above.

Figure 5 shows the result of the running diffusion coefficient *κ*(*t*). For low energies, i.e., 10<sup>5</sup> GeV (brown), 105.5 GeV (purple), 10<sup>6</sup> GeV (red), and 106.5 GeV (green), particles reach a plateau after an initial ballistic phase as described in Section 3, which can be used as the steady-state diffusion coefficient. For larger energies (10<sup>7</sup> GeV, orange, and 108 GeV, blue), such a convergence is not observed. The reason is that the particles leave the plasmoid before being able to reach a steady-state diffusion limit. Therefore, the first three values of the steady-state diffusion coefficient are used to determine the energy dependence *κ*(*E*).

**Figure 5.** Running diffusion coefficient, *κ*, for energies 10<sup>5</sup> GeV to 10<sup>8</sup> GeV. A plateau is built up for energies between 10<sup>5</sup> GeV < *E* < 10<sup>6</sup> GeV. Toward higher energies, the coefficient breaks off as the particles leave the plasmoid before they can reach the steady-state diffusion coefficient. This is consistent with the calculated energy for a transition between a ballistic and diffusive behavior at 10<sup>6</sup> GeV.

In Figure 6, the values averaged from the plateau in Figure 5 are shown with the corresponding error bars. In the simulation setup, the magnetic field is a purely turbulent one. That means Bohm diffusion is at work, and the energy dependence is expected to be *κ* = *κ*<sup>0</sup> · (*E*/ GeV); see, e.g., [12]. We therefore perform a linear regression and find *κ*<sup>0</sup> = 1013.54±0.06. For the simulations considered, the calculated values are used for the three energies, where this was possible in a reliable way (105 GeV to 106 GeV). For larger values, the result from the linear regression is used to determine the diffusion coefficient. From the findings of Section 2, one expects the ballistic and diffusive flares to provide approximately the same results until the transition energy is reached according to Equation (1), for the set of parameters used (see Table 2); this happens at around 106 GeV. From the results of Section 3, a deviation between the diffusive and ballistic approach is expected to be largest at small times and the two approaches should converge toward large times, when the steady-state diffusion coefficient is reached. As can be seen from Figure 5, this effect is energy dependent and for low energies (10<sup>5</sup> GeV), the diffusive steady-state is reached at around 103 s, while it takes <sup>10</sup><sup>4</sup> s for the highest energies (*<sup>E</sup>* <sup>&</sup>gt; 107 GeV).

**Figure 6.** Steady-state diffusion coefficient as a function of energy for particles between 10<sup>5</sup> GeV and 106 GeV. A linear regression for the function *<sup>κ</sup>*(*E*) = *<sup>κ</sup>*<sup>0</sup> · (*E*/GeV) is performed, with the uncertainty band depicted as gray shaded area. Here, *κ*<sup>0</sup> is the normalisation constant. The linear behavior with energy is based on the assumption that Bohm diffusion is dominant in the purely turbulent field.

Figure 7 shows the flaring behavior for a monochromatic energy flare at *E* = 105 GeV. The diffusive description shows an especially large enhancement with respect to the equation-of-motion approach at early times below 103 s, in accordance with the findings of Section 3. The behavior of d*N*/d*t* is similar for both approaches, with a small shift that can be explained by the uncertainties in the numerical determination of the diffusion coefficient used here for the diffusive approach. In the diffusive transport regime for a constant diffusion coefficient, one expects Δ*x* ∝ *t* 1/2. This results in the differential particle number d*N*/d*t* ∝ *t* <sup>−</sup>1/2 of escaping particles. Note that, due to the steady escape, the decrease in the number of remaining particles in the plasmoid leads to a strong cut-off at large times. Since, in the diffusive approach, more particles initially leave the plasmoid, the particle density in the plasmoid is lower than in the equation-of-motion approach, so that an earlier cut-off is visible.

Figure 8 shows the flare for diffusive and equation-of-motion behavior at 108 GeV. Here, one obtains a visible difference between the two flares, with the diffusive approach yielding a dominant contribution at early times. In contrast to the diffusive regime with d*N*/d*t* ∝ *t* <sup>−</sup>1/2, one expects a constant differential particle number for the ballistic transport regime with Δ*x* ∝ *t*. The initial slight drop for the equation-of-motion may be explained by statistical deviations from the initial homogeneous particle distribution in the plasmoid, especially when slightly more particles are in the outer spheres of the plasmoid at *t* = 0. Since, in the diffusive case, significantly more particles initially leave the plasmoid, the particle density in the plasmoid is much lower than in the equation-of-motion approach, so

that a cut-off is visible much earlier in the diffusive approach. Thus, these test simulations emphasize the importance of propagating the particles in the proper transport regime. Only a thorough analysis of the transport properties will lead to a prediction that can be compared to the observation of non-thermal emission from blazars.

**Figure 7.** Cosmic-ray flare (*E* = 105 GeV) from a blazar in the diffusive propagation model (orange downward triangle) and in the ballistic propagation model (blue upward triangle) as differential particle number per unit time, d*N*/d*t*, over time. The total number of particles injected into the simulation is *N*inj = 105.

**Figure 8.** Cosmic-ray flare (*E* = 108 GeV) from a blazar in the diffusive propagation model (orange downward triangle) and in the ballistic propagation model (blue upward triangle) as differential particle number per unit time d*N*/d*t* over time. The number of injected particles is *N*inj = 105.

#### **5. Conclusions**

In this paper, the propagation regimes in plasmoids of blazars are investigated as sources of high-energy cosmic rays, which in turn become emitters of high-energy gammarays and neutrinos. To explain the spectral energy distributions and lightcurves of this high-energy emission, it is shown that it is necessary to distinguish between the different energy and time regimes of ballistic and diffusive transport. At early times and at high energies, the particles are still in the ballistic regime. At late times or in scenarios, for which the injection of high-energy cosmic rays is significantly longer than the characteristic timescale *τ R*/*c*, with *R* being rthe radius of plasmoid and *c* the speed of light, the diffusive approach needs to be applied. The details of this transport modeling have an impact on both the spectral energy distribution, and on the temporal evolution of a flare. For the energy behavior, the diffusive part of the spectrum is steepened by the diffusion timescale, which is dominated by the diffusion coefficient. The ballistic part, on the other hand, is connected to an energy-independent escape timescale, thus leading to an emission spectrum close to the acceleration spectrum. When looking at the flaring behavior, it is shown that diffusive approach and transport with the equation-of-motion approach yield about the same result at low energies around 10<sup>5</sup> GeV, where the diffusive approach is accurate at times above <sup>∼</sup>10<sup>3</sup> s. That means that if the equation-of-motion approach is performed with the correct parameter setting, the same result is expected at times larger than 103 s, which is confirmed within a factor of <sup>∼</sup>2. When approximating this behavior with an escape timescale, the diffusive timescale needs to be applied.

At large energies (108 GeV), it can be shown that the diffusive and equation-of-motion approaches lead to quite different flaring behaviors. Here, only the equation-of-motion approach yields a correct result, as it can reproduce the ballistic behavior. It can be approximated in a transport equation approach by applying an energy-independent (ballistic) escape time.

As for the observation of neutrinos, we predict that spectral energy distributions of blazars like TXS0506+056 should have a break in the spectrum at a neutrino energy around 5 · <sup>10</sup><sup>13</sup> eV. If blazars are responsible for the diffuse neutrino emission, a break at the same energy is predicted. For the gamma-ray and neutrino lightcurves of blazars like TXS0506+056, one would expect the behavior of the variability time to be affected by the effect of a ballistic behavior at early times that turns into a diffusive behavior at later times. In particular, when simulating these flares, the shape of the flare will be misinterpreted when using the diffusion approximation at the highest energies, because it breaks down as the particles escape faster than the diffusion time scale.

**Author Contributions:** Conceptualization, J.B.T.; methodology, J.B.T., M.H., I.J., P.R., W.R., M.S., F.S.; writing—original draft preparation, J.B.T., P.R., I.J.; writing—review and editing, J.B.T., M.H., I.J., P.R., W.R., M.S., F.S.; visualization, P.R., I.J., M.S.; supervision, J.B.T., F.S.; funding acquisition, J.B.T., W.R., F.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Science Foundation DFG via the Collaborative Research Center SFB1491 "Cosmic Interacting Matters—From Source to Signal". Further funding was received from the DFG via the Grant "Multi-messenger probe of Cosmic Ray Origins (MICRO)" , Grant numbers TJ 62/8-1. We would also like to thank the Research Department for Plasmas with Complex Interactions for support.

**Data Availability Statement:** Data presented in this article can be made available upon request.

**Acknowledgments:** J.B.T. and W.R. would like to use this opportunity to thank Reinhard Schlickeiser for the long, pleasant journey through physics, administration, and the Ruhr area during the past decades. We hope the journey will continue for long-even if its path might shift in its character. We would like to thank Rainer Grauer for discussions on the plasma physics of blazars, particularly concerning the launching and evolution of plasmoids and their plasma parameters. We also thank Imre Bartos, Peter Biermann, Anna Franckowiak, Francis Halzen, Emma Kun, and Walter Winter for discussions on the modeling of non-thermal blazar emission

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Physics* Editorial Office E-mail: physics@mdpi.com www.mdpi.com/journal/physics

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9141-4