Scrutinizing FR 0 Radio Galaxies as Ultra-High-Energy Cosmic Ray Source Candidates

Fanaroff-Riley (FR) 0 radio galaxies compose a new class of radio galaxies, which are usually weaker but much more numerous than the well-established class of FR 1 and FR 2 galaxies. The latter classes have been proposed as sources of the ultra-high-energy cosmic rays (UHECRs) with energies reaching up to ${\sim}10^{20}$ eV. Based on this conjecture, the possibility of UHECR acceleration and survival in an FR 0 source environment is examined in this work. In doing so, an average spectral energy distribution (SED) based on data from the FR 0 catalog (FR0CAT) is compiled. The resulting photon fields are used as targets for UHECRs, which suffer from electromagnetic pair production, photo-disintegration, photo-meson production losses, and synchrotron radiation. Multiple mechanisms are discussed to assess the UHECR acceleration probability, including Fermi-I order and gradual shear accelerations, and particle escape from the source region. This work shows that in a hybrid scenario, combining Fermi and shear accelerations, FR 0 galaxies can contribute to the observed UHECR flux, as long as $\Gamma_\mathrm{j}\gtrsim 1.6$, where shear acceleration starts to dominate over escape. Even in less optimistic scenarios, FR 0s can be expected to contribute to the cosmic-ray flux between the knee and the ankle. Our results are relatively robust with respect to the realized magnetic turbulence model and the speed of the accelerating shocks.


Introduction
The idea of radio galaxies (RGs) contributing to the overall flux of ultra-high-energy cosmic rays (UHE-CRs), as measured on Earth, is several decades old and has been revisited many times in the past (e.g., [1,2,3,4,5]). Indeed, radio galaxies are among those source populations that are considered capable of accelerating charged nuclei to energies up to ∼10 20 eV, and hence, they fulfill the famous Hillas criterion [6] as well as energetics criteria [7,8] -as do their beamed counterparts, blazars. In this latter case, a few powerful sources are considered as the origin of detected UHECRs, and a correspondingly anisotropic cosmic ray (CR) arrival distribution is expected (e.g., [9]). Indeed, recently the Pierre Auger Observatory (PAO) collaboration reported a ∼5σ dipole signal above 8 EeV [10].
Another search strategy for UHECR sources considers instead a sufficiently numerous class of (comparatively) low-luminosity sources more isotropically distributed in the nearby universe, such as starburst galaxies (see, e.g., [11,12]) or radio galaxies. The crucial influence of diffuse CR transport in the intergalactic magnetic field on the expected arrival distribution has been examined, e.g., in [5,13,14]. Studies of either scenario have typically focused on CR transport between the source distribution and Earth, thereby providing constraints on the required spectra and composition of UHECRs that have escaped the sources (e.g., [15,16]). Only a few provide a study of acceleration and energy losses of CRs in the putative CR source environments (e.g., [17,18]), which is required to answer questions regarding the energization and survival probability of charged particles in a particular source class.
Radio galaxies belong to the radio-loud jetted class of Active Galactic Nuclei (AGN) and have been historically divided into Fanaroff-Riley (FR) 1 and FR 2 radio galaxies depending on their morphology in the radio band [19]. While FR 2s display prominent hot spots at the end of their powerful radio lobes and bright outer edges, FR 1 jets are more diffuse and weaker. Low power FR 1-type jets are usually linked with radiatively inefficient accretion flows, while radiatively efficient accreting objects produce FR 2-type jets.
Thanks to the ever-increasing sensitivity of large-area radio and optical surveys, a new population of very weak radio galaxies emerged by cross-matching the SDSS with the NVSS sample [20] and has been named FR 0 RGs [21].
So far FR 0s have never been proposed as contributing to the observed UHECR flux 1 . Their core radio properties turned out to be similar to typical FR 1 sources (see, e.g., [23,24,25,26,27] for radio properties of FR 1/2), as well as their optical classification as low-excitation RGs (LERGs); however, their extended radio emission shows a pronounced deficit (by a factor on the order of 100) of radio power (e.g., Baldi and Capetti [20]). Indeed, in FIRST and VLA radio images, they appear rather compact down to a scale of ≤ 1 kpc [28]; only on the pc-scale VLBI images could radio jets be resolved in most cases [29], and Doppler boostings in the range 1.7 to 6 were inferred. Their Eddington-ratios, L bol /L edd , of ∼10 −3 . . . 10 −5 indicate radiatively very inefficient accretion taking place in FR 0s, compared to FR 1/2 (see, [30,27]). Their spectral information in the X-ray band suggests a circum-nuclear environment depleted of a dusty torus [31]. Their host galaxies are classified as (slightly less) massive ellipticals ( 10 8 . . . 10 9 M ), i.e., contain a rather old stellar population with a higher than solar metallicity. FR 0s seem to prefer a lower galaxy density environment than that of FR 1s [32]. Additionally, their intrinsic space density is about a factor of ∼5 times higher than that of FR 1s in the local universe, n FR0 ∼ 5n FR1 [20,33], which seems to feed into the expectations of the cosmic down-sizing scenario (e.g., [34]).
A sample of FR 0 RGs ("FR0CAT "), limited to a small fraction of the sky, has been published in [35], which is used in the present study. 2 The population of FR 0s qualifies as a contributor to the observed UHECR flux from an energetics perspective. The required per source CR-power of U obs,UHECR /n FR0 ∼10 40.5 erg s −1 (with U obs,UHECR being the measured energy density of the UHECRs at Earth; see, e.g., [37]) lies well below their typically available jet power of ∼10 42...43.5 erg s −1 (as estimated from their radio power following [30]). Hence, FR 0s seem sufficiently numerous within the cosmic-ray horizon to account for an appreciable part of the observed UHECR flux measured on Earth.
In this work, we scrutinize this class of RGs as a contributor to the observed UHECR flux with a particular focus on their environment. For this purpose, we first (Section 2) study various acceleration mechanisms for CR nuclei -Fermi-I/II, gradual shear, and relativistic blast wave acceleration -in tandem with all relevant particle losses in the radiative environment of AGN jets -Bethe-Heitler pair production, photo-disintegration, photo-meson production synchrotron losses, and diffusive particle escape. This is then applied to the average radiative environments of FR 0 RGs in Section 3. Here, we assess the potential of these sources to accelerate nuclei to the highest energies observed so far; while at the same time surviving intact in this environment. We conclude this work by summarizing our results in Section 4.

Particle Energy Losses and Gains in Jets
To assess whether a class of astrophysical objects can be considered a cosmic-ray accelerator -the total energetics (as discussed in Section 1) and the source environment must be evaluated. Models of the relevant processes that influence the maximal acceleration capability of a source are summarized here. These processes can be divided into three different competing groups: 1) Particle acceleration, 2) Energy losses of CRs, and 3) Particle transport out of the region of interest. Whichever process happens on the shortest time scale will dominate the cosmic-rays' energy gains or losses. Comparing the three associated time scales (τ acc , τ loss , τ esc ) will allow the assessment of the maximum energies that can be reached in the sources. Accordingly, the first part of Section 2 deals with established particle acceleration models for jetted AGN. In the second part, the loss time scales are discussed, where the focus is on losses due to electromagnetic pair production (Bethe-Heitler), photo-disintegration, and photo-meson production. Finally, a parameter scan of the loss length is provided for reasonable target field energies, , and nuclei energies, E N . The corresponding look-up plots will allow quick evaluation of approximate values for the acceleration chances by comparison with the derived acceleration time scales.
In this paper, we refrain from detailed modeling of cosmic-ray energy spectra and secondary spectral energy distributions (SEDs) in the source region but rather attempt to answer whether FR 0 can potentially reach the required energies to contribute to the UHECR flux.

Acceleration Time Scales
In this section, we will briefly discuss possible acceleration scenarios for the core region of FR 0 radio galaxies. We focus on the time scales and refer interested readers to the original publications for the proposed acceleration mechanisms' technical details. In general, the acceleration time scale is given by where η ≥ 1 is a model-dependent scaling factor, and λ is the scattering length. In general, this scaling factor's dependencies can complex, e.g., the scattering center speed (see Section 2.1.2) or the shock, which is neglected in this work. The shortest acceleration time scale is given by η = 1 and λ = r Larmor (see, e.g., [38]).

Fermi 2nd Order
More than 70 years ago, Fermi proposed a model for cosmic ray acceleration in turbulent magnetized clouds [39]. The stochastic nature of this process allows for energy gains and losses due to the randomly moving scattering centers. This fact makes the second-order process relatively slow. On the other hand, it is easily realized since no large-scale shock structure is needed, only a turbulent magnetic field.
Second-order Fermi acceleration can be described by a corresponding (momentum diffusion) transport equation (see, e.g., Rieger et al. [40]). From this equation, the acceleration time scale can be derived to be where the momentum diffusion coefficient is given by κ p , and p is the absolute value of the momentum. This time scale now depends on the specific diffusion model that is applicable in the source region.
In the turbulent cascading regime, where a power law describes the power spectrum, the following relation between the spatial and momentum diffusion coefficient holds (see, e.g., Drury and Strong [41], Sigl [42]): where the Alfvén speed is v A = B/ √ 4π ρ, with ρ being the plasma density, and B the magnetic field strength, while α is the power-law index of the spatial diffusion coefficient D x (p) ∝ p α . The relevant time scale can then be calculated by inserting Eq. (3) into Eq. (2) and assuming a proportionality constant of 4/3 [42], which yields:

Fermi 1st Order
The concept of Fermi acceleration was later extended to include a more efficient acceleration mechanism, which can be seen as the predecessor of diffusive shock acceleration [43], and is often the preferred acceleration model. An updated description of diffusive shock acceleration can be found, e.g., in [44].
The first-order Fermi acceleration time-scale is determined by the conditions on both sides of the shock. When the shock is expanding with velocity u s in the laboratory frame, the time scale is given by (see, e.g., [45,46]): Here, u 1 = −u s and u 1 = s u 2 are the upstream and downstream velocity components parallel to the shock normal measured in the shock frame, respectively, and s is the compression ratio. The spatial diffusion coefficients describing the transport on either side of the shock are denoted with κ i . Given that the FR 0 source class environment is not well known, we choose to consider the most optimistic scenario. This includes the assumptions of a strong non-relativistic shock with s = 4 and Bohm diffusion, which lead to the shortest acceleration time scales. This ansatz then yields: where q = Ze is the particle charge, and Z is the charge number. Here, the ultra-relativistic approximation for particle speeds (γ N 1, v ≈ c) is applied throughout this paper, which implies E ≈ pc. Comparing the time scales of first-and second-order Fermi acceleration, it becomes clear that the Fermi-II process is slower by approximately a factor of ∼(v A /u s ) 2 since the Alfvén speed is usually much lower than typical shock speeds.
In general, the acceleration efficiency can hence be expressed by η Fermi−I = 20/(3α) β −2 s , where α is the diffusion coefficient power law index and β s = u s /c.

Gradual Shear Acceleration
Gradual shear acceleration can be interpreted as a stochastic acceleration with scattering centers moving in an ordered direction. Therefore, the acceleration time scale can be calculated based on Eq. (2), when the corresponding effective diffusion coefficient is known. In Rieger and Duffy [47], this is quantified for flows with a cylindrical profile as which yields an acceleration time scale Here, ∂u z /∂r is the gradient of the flow profile, τ 0 is the energy-independent scattering time, and Γ j is the jet's Lorentz factor. The scattering time is related to the mean free path λ of the transport via λ = τ 0 p α c. This time scale corresponds to an acceleration efficiency of η shear = (c/λ) 2 /(Γ(4 + α)). More details, and the fully relativistic derivation, including other flow profiles, can be found in Rieger and Duffy [48]. Rieger and Duffy assume that the particle's Larmor-radius is smaller than the width of the shear layer r Larmor < ∆r, introducing an upper limit on the applicable energy range. When the gyroradius of the CR is larger, the shear layer can be approximated by a discontinuity; numerical results for the acceleration rate, in that case, can be found, e.g., in [49]. Assuming a linear flow profile (∂u z /∂r = ∆u/∆r), where u z is the flow speed along the jet axis, and Bohm diffusion for the scattering time, Eq. (8) becomes One interesting fact to notice is the energy dependency of shear acceleration: The acceleration time scale decreases with increasing CR energy -the inverse of Fermi-type acceleration. This change in energy dependence behavior opens the possibility that the dominant acceleration mechanism changes with energy if the environment allows. Depending on the exact jet properties (∆u, ∆r, and Γ j ), this could lead to unphysically fast acceleration, corresponding to η(E) < 1 at high energies. Therefore, this description of gradual shear acceleration can only be applied up to a maximum energy E shear max , making sure that r Larmor (E shear max ) < ∆r, and η(E shear max ) > 1. Comparing results from Rieger and Duffy to [49,50] suggests a transition of the acceleration time-scale energy-dependence, from ∝ E −1 to ∝ E +1 , around r Larmor ≈ ∆r. Since the acceleration process in the transition is poorly known, we will apply the most optimistic scenario there: η = 1 for energies above the threshold E > E shear max . As shown in Section 3.1, larger values for η ≈ 10, as estimated in [50], would not significantly change the resulting maximum energy.
In the following calculations it will be assumed that the central jet speed is given by u z (r = 0) = c 1 − Γ −2 j and vanishes at the edge of the jet, u z (r = ∆r) = 0.

Relativistic Blast Wave
The relativistic blast wave scenario [51], sometimes referred to as Espresso shot acceleration [52,53], differs from the other scenarios as it incorporates a drastic change of the energy gain per encounter over time. While the CRs' first crossing of the ultra-relativistic shock front can give an energy boost on the order of Γ 2 b (Γ b is the Lorentz factor of the blast wave), all other crossings will give only an increase of a factor ∼2. The downstream residence time t d approximately gives the time needed for the first boosting. In all subsequent cycles, the time scale is given by the sum of the up-and downstream residence time, which is mainly dominated by the upstream time scale t u [51]: Here, l c is the correlation length of the turbulent magnetic field, and r g is the gyro-radius. Therefore, this acceleration scenario is especially interesting when an already pre-accelerated particle distribution is available at the source site. This population can then be significantly boosted in energy on a very short time scale. However, it is not much more efficient as a sole acceleration mechanism than ordinary Fermi-I-type acceleration. Furthermore, the boosting factors are relatively low for FR 0 jets. Therefore, this mechanism is not further investigated in this work.

Loss Time Scales
In this section, the principle loss processes relevant to the source region are discussed. Since abundant literature already exists on this topic (see, e.g., [54,55] in application to AGN models or [56,57,58] for a general description), we refrain from detailed discussions and only present the necessary information to reproduce our results.
The most relevant loss processes in the radiatively dominated source region are: 1) Electromagnetic Bethe-Heitler pair production, 2) Photo-disintegration of the primary nucleus, 3) Photo-meson production, and 4) Synchrotron losses when the magnetic field is relatively high (∼1 G). In this section, the losses are discussed for mono-energetic target fields to estimate the resulting time scales independently of the target photon field shape. In Section 3, the loss time scales are evaluated for a concrete target photon field based on the average SEDs of FR 0 radio galaxies.
The energy loss length for photo-meson production (Section 2.2.2) and photo-disintegration (Section 2.2.1) is calculated based on the interaction rate τ −1 Nγ , weighted with the mean energy loss per interaction. For Bethe-Heitler pair production [59] and synchrotron losses, a direct (semi-) analytical description is possible and given in Section 2.2.3 and Section 2.2.4, respectively.
In general, the interaction rate for particle-photon processes is described by Here, √ s is the center-of-mass energy for the interaction between a photon with energy and a nucleus with mass m N and energy E: where r is the photon energy in the nucleus rest frame. The interaction angle is denoted with θ, and β N = v N /c is the normalized nucleus velocity. Quantities observed in the rest frame of the nucleus are denoted with a subscripted 'r'.
Assuming an isotropic 3 photon distribution n γ ( , Ω γ ) = n γ ( )/4π allows integration over the azimuthal component, which results in where r,thr is the process threshold, which will be discussed in the following sections. For a monoenergetic photon field, n γ ( ) = n 0 δ( − 0 ), one gets The energy loss length for photo-disintegration and photo-pion production is calculated as: where the interaction rate τ −1 is defined according to Eq. (12), and ∆E is the mean energy loss per interaction. For Bethe-Heitler pair production, only Eq. (16) is used.

Photo-disintegration
In this work, we refer to the photo-disintegration regime when the target photon's nucleus rest-frame energy is in the range of 10 −2 ≤ r /GeV ≤ 0.14, where the upper boundary is approximately the pion production threshold. The left plot in Fig. 1 shows the total cross-sections used for the four relevant tracer elements 4 , helium, nitrogen, silicon, and iron.  The cross-sections are taken from tabulated values, as published in the CRPropa data repository 5 . The elements with A > 12 are based on the TALYS 1.8 Monte-Carlo simulations of the photo-hadronic cross-section [61]. The cross-section for lower mass numbers A ≤ 12 are taken from various references, including [62,63,64,65]. 6 The energy loss length is calculated as a weighted average over all possible disintegration channels Here, b i is the branching ratio, which is the normalized exclusive interaction rate of the specific process: b i = τ i (E)/ j τ j (E). The relative energy loss is proportional to the mass ratio of the fragment(s) and primary particle ∆E/E = m sec /m N , where ∆E is the energy loss and m sec is the mass of the fragment(s). As an example, we show the exclusive cross-sections for iron in the right plot of Fig. 1.

Photo-Meson Production
For interaction energies above r,thr ≈ 140 MeV photo-meson production sets in. The cross-section, see Fig. 2, is first dominated by several hadronic resonances before multi-pion-production dominates at higher center-of-mass energies. The neutron cross-section deviates from that of protons, mainly in the secondary resonance region. In this work, the cross-section for more massive particles is approximated using the superposition model: σ N = A α (N p σ pγ + N n σ nγ ). The variable A α describes a shielding factor, determining which portion of the nucleons are shielded and not available for interactions. The exponent α can vary from unity to 2/3, depending on the mass number A of the nucleus.
In this simple superposition model, exactly one nucleon takes part in the interaction, making the energy loss calculation a simple scaling with the mass number: ∆E = E/A. The energy loss for protons and neutrons is approximated as occurring in the ∆-resonance regime, leading to a decrease by a factor of ∆E p/n = E(1 − m p/n /m ∆ ) ≈ 0.24E.

Bethe-Heitler Pair Production
The calculation of the electron-positron pair production loss length is based on the semi-analytical approach described in Chodorowski et al. [69]. The assumptions made there -γ N 1 and r m e c 2 5 https://github.com/CRPropa/CRPropa3-data 6 See [66,67], for an explanation and the data repository for the exact listings. The photo-pion production cross-section for protons, neutrons, and iron is shown. In this work, we used a simple superposition model to construct the cross-section of heavier elements. The cross-section data comes from the CRPropa [68] implementation of the SOPHIA code [56].
-are usually very well justified for UHECR. Assuming isotropy of the target field, the energy loss is then described by: Here, α S ≈ 1/137 is Sommerfeld's fine-structure constant, r 0 is the classical electron radius, and Z is the nucleus charge number. The nuclear rest frame photon energy is parameterized in units of the electron momentum as κ = 2γ N /(m e c). The function φ(κ) is a semi-analytical description of the cross-section integral (see Eq. 3.12 of the original paper [69]).

Synchrotron Radiation
The average energy loss of an isotropic nuclei ensemble due to synchrotron radiation can be written as: Here, σ T is the Thomson cross-section, and u mag is the magnetic energy density. Furthermore, v is the nucleus speed, γ its Lorentz factor, Z the charge number, m e is the electron mass, and m N the nucleus mass.

Combined Loss Length
The total energy loss length is the combination of all four processes given by L −1 loss,tot = i L −1 loss,i . Figure 3 shows the individual loss lengths for the four loss processes and the total loss length depending on the nucleus energy E N for iron; all other elements can be found in Fig. A.9. Here, an exemplary monoenergetic target photon field with an energy of 0 = 1 meV and a number density of n 0 = 400 cm −3 , close to the CMB monochromatic approximation, was chosen. The magnetic field strength (B = 4 µG) was chosen such that equipartition between the magnetic and the target photon field energies is reached. Depending on the nucleus energy, different processes dominate, with photo-disintegration leading to the highest losses at a given energy above its threshold. The position of the minimum changes, of course, with changing energy of the target photon field.  which is combined with photo-disintegration (PDI, orange-dashed), photo-pion-production (PPP, green-dash-dotted), and synchrotron radiation (synch., red-dash-dot-dotted) to give the total loss length scale L loss (violet-solid).
As an initial step, a parameter scan for reasonable target field energies was performed in the range 10 −3 ≤ 0 /eV ≤ 10 2 . These scans are shown in Fig. 4 and Fig. A.10 for a monochromatic target field number density of n 0 = 1 cm −3 as an example. Note, that synchrotron radiation has not been included in these plots. However, the simple inverse scaling with energy for any given nucleus allows adding this loss to the combined loss rate of the other processes.
The gradient is found along the major diagonal, which corresponds to an increasing nuclear rest frame photon energy r . The target field energies in the optical regime (∼eV) have the largest influence on nuclei with energies around the ankle (∼10 18 eV). For higher nuclei energies, the most important target field energies are in the range of the CMB. Furthermore, at least for a monochromatic target field, only a relatively narrow range in nuclei energies will suffer from the highest losses. A broader target photon field, such as the SEDs discussed in Section 3, will smooth the relatively peaked loss rate.

Escape Time Scale
Cosmic-ray escape from the acceleration region is not yet completely understood. Ideally, a transport equation (including advection, streaming, and diffusion) should be solved -though this is beyond this paper's scope. With a velocity profile along the jet axis, perpendicular escape is only possible by diffusion (τ ⊥ = R 2 /(2κ)). In the parallel direction, escape could be dominated by advection (τ = L/u).
However, for the modeled jet structures (see Section 3), perpendicular diffusive escape is the more relevant loss process, when the jet is longer than L 10 pc. 7 Additionally, at the highest energies, where gradual shear acceleration dominates the energy gains (see Section 3.2), perpendicular diffusion is the most likely escape process (see, e.g., [70,71]). Therefore, we focus here only on diffusive escape with a time scale of where R is the size of the acceleration region, and κ is the diffusion coefficient. Assuming again Bohm diffusion, which gives a lower limit on the escape time scale, Eq. (18) yields

Maximum Energy
The maximum energy for any multiplicative acceleration process is given by the Hillas criterion, which yields for the maximum energy: Here, R is the size and B the magnetic field strength of the acceleration region. The speed of the scattering centers is βc. However, in most cases, the maximum energy will be smaller than defined by the Hillas criterion. Neglecting losses -which is a valid assumption in some cases (see Section 3) -the maximum energy can be calculated by comparing the escape time scale with the acceleration time scale.
For Fermi-I acceleration -assuming a strong shock and the same diffusion coefficient κ on both sides of the shock (see Eq. (5)) -and an escape process driven by the same diffusion process (see Eq. (18)), the maximum energy becomes which yields for the maximal particle energy. Here, κ 0 is the diffusion coefficient normalization with κ(E 0 ) = κ 0 . In the case of Bohm diffusion Eq. (22) can be approximated by Since shear acceleration has the same energy dependence as the escape process, a general condition for efficient acceleration is given by Γ j 1.6, where τ shear acc < τ esc . The gradual shear acceleration breaking point is given by η(E shear max ) = 1. For Bohm diffusion and a linear shear profile (see Section 2.1.3) the maximal shear acceleration energy is E shear max = 3/(Γ 2 − 1) (qBc)/Γ ∆r. Above this energy the acceleration process cannot be described by Eq. (8), therefore we assume the most optimistic case (τ acc = λ/c) for higher energies.
Comparing the acceleration length scale with the scales of losses and escape leads to a convenient definition for the maximum energy E Acc max by Here, i and j are wildcards for all possible acceleration and loss processes, respectively. Based on Eq. (24) we define the acceleration probability: which describes the chances for acceleration, with P acc E Acc max = 0.5.

FR 0 source environment
In this section we aim to find a reasonable description of the target-photon field expected for CRinteractions in the environment of FR 0s. For this purpose we derive the average source SED in the relevant target field energy range. This includes internal (jet) and external photon fields, where both particle acceleration and energy losses take place. We note that a complete and self-consistent description of the source physics is beyond the scope of this work which focuses only on the nuclei and their acceleration chances, which are examined in Section 3.2.

Average source SEDs
The FR 0s sample used here consists of the FR0CAT [35] (in its revised version, containing 104 sources, more details can be found in the appendix of [28]) to which were added 10 of the FR 0s studied in [31] 8 for their X-rays properties. In total, there are then 114 FR 0s for which we collected SED data from the NASA Extra-galactic Database 9 and SSDC Sky Explorer 10 . The selection criteria for the 104 sources included in FR0CAT are the following [35]: • FIRST flux ≥ 5 mJy and no extended emission.
• Observed major-axis in FIRST images ≤ 6.7", at redshift z = 0.05, this corresponds to a radius of ∼2.5 kpc. • Radio sources with a maximum offset of 2" from the optical center.
Similar criteria (less constrained to enlarge the sample) were applied to select the other ten sources [31]: • FIRST flux ≥ 30 mJy (to ensure X-ray emission) and no extended emission.
While collecting data for the 114 sources, cuts were made on some of the available data points. First, spectral data points without corrections from other bright sources in the instruments' beam were excluded. Further, a comparison between the sources' position ellipses and the IRAS observations' targeted positions led us to not consider these IRAS data points, as the ellipses did not overlap for the sources concerned.
Due to the low spatial resolution of the UV images provided by Galaxy Evolution Explorer (Galex ) [72], we follow [73] and treat these data points as upper limits.
As mentioned before, the case of Tol 1326-379 (which was claimed to be the first FR 0 associated with a γ-ray counterpart) is uncertain since it is no longer associated with a γ-ray emission in the last version of the Fermi catalog 4LAC [22]. For this reason, and because they are not relevant as target photon fields for this work, the γ-ray data obtained for Tol 1326-379 are not included in the global SED.
Based on all the source sample data described above, we construct an average SED, from which we derive the average target photon field in the jet frame. We note that the FIR to UV band is likely dominated by the host galaxy, whereas the rest of the multi-wavelength data is associated with the jet. Hence, two different photon fields, one external and one internal to the jet, need to be considered. From now on, the unprimed quantities refer to the ones in the galaxy frame. In this frame, each source's luminosity is L ν = (4πd 2 L )/(1 + z)F ν,obs = 4πd 2 L F ν , where d L is the luminosity distance, and z is the redshift. The mean value for the luminosity distance of all the sources in our sample is d L ≈ 5.6 × 10 26 cm and z ≈ 0.04.

Jet component
The Synchrotron-Self-Compton (SSC) model can be considered a minimal model to describe the jet component as synchrotron radiation and Inverse Compton scattering inevitably occurs in any magnetized region that hosts relativistic particles. This motivates us to use a simple steady-state SSC model, in the following, to describe the typical jet component of FR 0s at low energies. The apparent positive slope in the X-ray SED of FR 0s may indicate the high-energy component's onset. Different parameters were investigated, focusing on those that impact the maximum energy particles may reach (see Section 2.4). Specifically, the product of the emission region size and the magnetic field strength r em B was maximized in Table 1: Parameters of the SSC models. The following parameters are the same for all models: γ min = 100, p 1 = 2, p 2 = 3, and θ = 20 • . The luminosity is given in logarithmic units log 10 (L i /(erg s −1 )), where a luminosity ratio between protons and electrons of ξ = 10 is assumed. order to optimize the maximal energy limited by the Hillas criterion. All other parameters have been tuned in order not to overestimate the observed fluxes. In doing so, it was considered that some of the parameters, e.g., the maximal electron Lorentz factor γ max , might be limited; the theoretical upper limit for γ max can be estimated by comparing acceleration and synchrotron losses for electrons.
Here we present six different models, computed using the SSC/EC Simulator by ASI 11 . These models are able to represent the synchrotron component of the jet and do not violate the X-ray data. The models use a broken power-law for the electron differential density n (γ ), described with the parameters γ min , γ max , γ cut , p 1 and p 2 , such that: where γ' is the electrons' Lorentz factor. The parameters for each model can be found in Table 1.
In Fig. 6, the model SEDs are shown together with the data. All models agree well on the description of the low energy data ν 10 16 Hz. Given the different sizes of the emission region -and therefore different photon-energy densities -the second bump is quite different for each of the models. Especially, models A, B, E, and F need an additional (possibly hadronic) component to explain the observed flux above ∼10 16 Hz.
A Lorentz transformation is used to evaluate the jet radiation field in the comoving frame. For an internal photon field, if the emission is isotropic in the comoving frame (denoted with primed quantities), the jet frame specific energy density is [74]: where F ν is the spectral flux density, δ D = (Γ j (1 − βµ)) −1 is the Doppler factor, µ is the cosine of the angle between the line of sight and jet direction, Γ j is the jet's Lorentz factor, β = (1 − Γ −2 j ) 1/2 is the normalized velocity, and r em is the emission region comoving size.
For all models, the energy densities of the three major constituents are listed in Table 1: 1) The magnetic field energy density (where a homogeneous turbulent field strength is assumed), 2) The total radiation energy density u γ , and lastly, 3) The CR energy densities u CR = u e + u p , with contributions from electrons u e and protons u p .
We estimate the protons' contribution in two ways: 1) Assuming equipartition of the magnetic field energy density and the cosmic-ray energy density u B = u CR (which makes sense only if the former is  larger than the electron energy density u B > u e ) and, 2) Assuming a constant energy density ratio of protons and electrons u p = ξu e . Additionally, assuming charge neutrality for the non-thermal proton and electron distribution and the same power-law index α ≈ 2.2 for their spectra, the luminosity ratio becomes ξ ≈ 10 − 100. The ratio ξ can change significantly when both populations follow different spectral slopes (see, e.g., [75,76]). However, in this work, we fix the ratio to ξ = 10 for our calculations.
The total jet luminosity is given by where u i are all relevant energy densities, and is calculated for all models (see Table 1).
Comparing the derived values with the jet's average mechanical power, L jet 10 43.5 erg s −1 (cf. Section 1), shows that only models E and F comply with that boundary. However, the underlying correlation is not well restricted [30]; therefore, it is plausible that some of the FR 0s will have enough jet power to support models B, C, and D, and only model A will be neglected for the rest of this work.

Host galaxy component
Following the treatment of [77], the spectral and radial dependencies are separated to express the energy density provided by the host galaxy inside a given region of radius R: u(E, R) = u(E) u(R). A Sersic profile can adequately reproduce the radial dependence of the surface brightness from elliptical galaxies [78]: where η = R/R e and R e is the half-light radius of the elliptical galaxy that encloses half of the total luminosity. For elliptical galaxies, R e ∼ 1 kpc, and with m = 4, we recover the De Vaucouleur's profile [79]. The advantage of this formulation is that it can be solved analytically [80] to give the projected monochromatic luminosity inside a given radius R: with b ≈ 7.669. Here γ(α, x) is the incomplete gamma function. The total luminosity can then be calculated as which gives its normalization with the complete gamma function Γ(α) = γ(α, ∞).
For the spectral part, we used a diluted black-body spectrum to reproduce the FIR to UV component, for which the normalized energy density is: with R = r em δ D being the size of the emitting region in the host galaxy frame. The luminosity is related to the energy density by We then determine the dilution factor dil by adjusting ∞ 0 L(E, R)dR = L tot u(E) to match the data with L(E, R) = u(E) L(R).
Finally the energy density for a region of size R is then given by: We use the invariance of u(E, Ω)/E 3 for the transformation of this external photon field energy density into the jet's frame: u (E , µ , r em ) = u(E, R)/(2Γ 3 j (1 + βµ ) 3 ) for an isotropic photon field in the galaxy frame. Replacing u(E, R) by Eq. (33) and integrating over µ one gets:

CR Acceleration in Loss Environment
Based on our estimates above for the typical target photon fields expected in FR 0 environments, it is now possible to evaluate their role as potential cosmic-ray sources. The target photon field's modeling fixes many parameters of the overall source region model, as explained above. It should be noted that some of the parameters are degenerate in the sense that other combinations of the fitting parameters might have led to equally good approximations of the observed photon fluxes. However, the FR 0 SEDs spread is much larger than the changes introduced by varying some of the degenerated parameters.
In the following, two scenarios will be tested based on Bohm and Kolmogorov diffusion. Bohm diffusion is the limiting case with the shortest time scale of all Fermi-I acceleration processes and does not need the introduction of parameters in addition to those given by the SSC model. However, at least on Galactic scales, observational evidence points towards a Kolmogorov diffusion κ(E) = κ 0 (E/(ZE 0 )) 1/3 (see, e.g., [81]). We fix the normalization by requiring that a particle with an energy according to the Hillas criterion E Hillas max has a mean free path smaller than the emission region size. 12 This will lead by construction to comparable time scales for the Bohm and Kolmogorov case and any difference will be mainly due to the different spectral indices. If another turbulence model or normalization were realized, the time scales will change significantly. However, the currently available observations are not able to discriminate such speculations from each other. Furthermore, we are assuming a wide shear layer which is given by the size of the emission region r em = ∆r and constant gradient ∂u z /∂r = β jet c/∆r (see Eq. (9)).
In addition, we will focus on Fermi-I and gradual shear acceleration for the rest of this work. These two cases are among the most promising processes considering FR 0s as sources of CRs. Classical Fermi-II acceleration can be neglected since it is not fast enough compared with the former two. Blast wave acceleration only provides a fast energy increase for the acceleration cycle's first encounter. For this reason, it provides a better mechanism for re-acceleration than for the full acceleration cycle. One of the neglected scenarios could, however, still play some role in the acceleration. Therefore, two acceleration scenarios are considered. For the first acceleration scenario, we consider pure diffusive shock acceleration, whereas the second scenario assumes a hybrid approach that models the energy gain by a combination of Fermi-I and gradual shear acceleration. Figure 7 shows all relevant length scales for iron in a single plot 13 . Exemplarily, the Bohm diffusion scenario is presented here using model E. The complexity of the processes involved makes it hard to draw quantitative conclusions about the maximum energy that might be reached in the underlying model E from a plot like this alone. Due to the rather broad-banded target photon field, all four loss processes are no longer sharply distinguished but smoothed into one another over several orders of magnitude in energy. As expected from the results presented in Section 2.2, the most relevant process for heavy nuclei is photo-disintegration. In contrast to the monochromatic target fields discussed above, photo-disintegration has a non negligiblesometimes dominant -influence on the total loss length even at the highest nuclei energies E > 1000 EeV.
Gradual shear acceleration can become the dominant acceleration mechanism at high energies, as reported in, e.g., [70]. The transition energy can be derived as Figure 8 shows ten models for the acceleration probability of iron (B-1 to F-2) for comparison. For the models based on pure Fermi-I acceleration, the probability will drop below P acc < 0.5 at some point. This statement holds independently of the chosen target field mode, as was already shown in Section 2.1. The maximum energy might be even smaller when the losses become significant.
In contrast, the hybrid scenario could lead to a saturation of the acceleration probability at larger energies lim E→∞ P acc (E) > 0.5, when E shear max ≥ E Hillas max . This saturation is because the loss length depends only weakly on energy at the highest energies, and the escape and gradual shear acceleration rates have the same energy dependence ∝ E −α (and are stronger compared with the losses). When shear acceleration is L acc -Fermi-I -shear L min acc L esc Figure 7: A summary of all relevant length scales for iron nuclei is shown. The parameters for the target photon field, magnetic field strength, and size of the emission region used here are that of model E. It is clear that the interplay between losses, escape, and acceleration can become complex. The total loss length L loss (blue line) combines pair production (BH), photopion production (PPP), photo-disintegration (PDI), and synchrotron losses. The same is done for gradual shear and Fermi-I acceleration (orange lines). The escape length scale is shown in green, and the minimum acceleration length, corresponding to η = 1, is depicted in red. The gray shaded area shows the region above the Hillas energy calculated for β = 1.
cut off by the maximum acceleration rate (E shear max < E Hillas max ) a steep decrease of the acceleration probability can be observed for E > E shear max . This is due to the fact that the sign of the power-law index α of energy dependence of the acceleration time scale flips at that energy.
The maximum energy in the jet frame for all the models has been derived by calculating the upper limit of the Hillas energy E Hillas max assuming β = 1, and the maximum acceleration energy given by P acc (E acc max ) = 0.5; the smaller of the two defines the maximum energy. The results for all models, and the five tracer elements, are given in Table 3 for the Bohm diffusion ansatz and in Table C.4 for Kolmogorov diffusion. From these five tracer elements, the average maximal rigidity ζ was also calculated. The shock speed u s , has of course, some influence on the time scale of the Fermi-I acceleration. However, the effect on the maximum energy is almost negligible for 0.03 < u s /c < 0.1. Therefore, all numbers shown in this work are calculated for u s = 0.1 c.
Models based on the hybrid ansatz lead, on average, to an increase of the maximum rigidity ζ max (which can be as large as a factor of ∼20). Furthermore, a higher charge number Z generally leads to a higher maximum energy as expected from the acceleration mechanism. However, when the target photon field density is high, and the FR 0 SED shows a larger second bump (e.g., in model C ), a clear proportionality with Z cannot be found. This non-proportionality effect is more pronounced in the hybrid scenario than Fermi-type acceleration alone. In general, it can be observed that only in a hybrid scenario, energies above the ankle can be reached. Fermi-type acceleration alone is only efficient enough to reach energies below the ankle.
Looking at the magnetic turbulence model's influence, it is interesting that pure Fermi-1 order acceleration in a Kolmogorov turbulence spectrum does not generate UHECR at all. In this case, the maximum energy is determined by balancing escape and acceleration and is, therefore, independent of the SED. However, gradual shear acceleration leads to the highest maximum energies that we have found within our model; see model B-2 and model E-2 in Table C    ζ = E/q is the rigidity of the particle in units of volt.

Summary & Conclusions
In this work, the relatively recently discovered class of FR 0 radio galaxies has been discussed as a potential source class for UHECRs. FR 0s are much more numerous than more energetic radio galaxies such as FR 1s or FR 2s (see Section 1). In this case, each source does not need to provide a high CR flux, which can be advantageous when modeling the rather isotropic CR arrival directions. In doing so, a brief summary of the most relevant processes has been given in Section 2, where different acceleration, energy loss, and escape models have been explained.
To answer the question of whether UHECRs can be accelerated, and survive, in the source region of FR 0s the available multi-wavelength data was compiled into an average SED to describe the expected target field for CR-photon interactions. This SED was decomposed into a jet component (based on six different SSC models) and a host galaxy component. The host galaxy was assumed to be a typical elliptical galaxy and represented by a de-Vaucouleurs'-radial-luminosity-profile combined with a diluted blackbody energy spectrum. Based on the Doppler and Lorentz factor constraints of the jetted regions the jet frame's expected photon target field has been derived. This modeling allowed a reasonable estimation of the time scales, of all relevant particle-photon interactions, to be expected in the environment of FR 0s, in tandem with acceleration and escape.
We define the probability P acc as a measure for the acceleration potential to derive the expected maximal energy for all different target field models. Since the nature of the magnetic field is not known, two turbulence models (Bohm or Kolmogorov) have been tested, as well as different shock speeds u s . From our calculations, we can draw the following main conclusions: 1. Acceleration of UHECR, up to the highest measured energies, in FR 0s, is possible under certain conditions: a) Acceleration due to the hybrid scenario, where HECR are predominantly accelerated by Fermi-I acceleration and the most energetic particles are driven by gradual shear acceleration. b) A jet Lorentz factor Γ j 1.6, so gradual shear acceleration happens on faster time scales than the escape (given a linear flow profile). c) According to our models an emission and acceleration region size of ∼5 × 10 17 cm (as described in, e.g., model E ) is the most promising size scale.
2. Less optimistic scenarios, based on pure Fermi-I acceleration with Bohm diffusion, still reach energies large enough to contribute to the CR flux, between the knee and the ankle. Models based on pure Fermi-I acceleration within Kolmogorov turbulence are not able to reach energies significantly above the knee.
3. The influence of the realized turbulence model -Bohm or Kolmogorov -is found to be relatively weak. This is partly due to the similar normalization of the two turbulence models. In reality, the diffusion coefficient might be different and may be incapable of keeping the most energetic particles confined in the acceleration region.
4. Within the hybrid scenario, the shock speed's u s influence on the maximum energy is weak as well. However, it has an impact on the transition energy between the two acceleration modes.

5.
A scenario based on gradual shear acceleration alone is not realistic, as the acceleration time scales at the lowest energies are too large due to the small boosting factors of FR 0 jets. In the hybrid approach, this problem is overcome by Fermi-I acceleration but could also be realized by other mechanisms such as, e.g., magnetic reconnection [83].
6. Despite the fact that FR 0s are radiatively weak sources, radiative losses not only play an important role for the maximum CR energy of FR 0s but they will have a significant impact on the expected cosmic-ray spectrum at the source (see Fig. A.11, Fig. A.12, Fig. A.13, and Fig. A.14). However, we find that synchrotron losses of nuclei are largely negligible. This suggests that a simple unbroken power-law with a cut off is probably insufficient to describe the cosmic-ray energy distribution escaping the source region.
We can summarize that FR 0s make an interesting candidate source class for UHECRs that deserves further investigations.
The code developed for this work is available upon request in a repository for reference.
Appendix A. Loss, Escape, and Acceleration Lengths  Loss lengths for the tracer elements hydrogen, helium, nitrogen, and silicon is shown. It is clearly visible that pair production (blue dotted line) can almost be neglected for a CMB-like monoenergetic target field. For the chosen magnetic field strength, yielding equipartition between the magnetic and the target photon field umag = utarget, synchrotron losses (red dash-dash-dotted) can be neglected as well. The short loss length is due to photo-disintegration (orange dashed line), while at the very highest energies, the losses are dominated by photo-pion production (green dash-dotted line).      Figure B.15: Acceleration probability for hydrogen, helium, nitrogen, and silicon based on Bohm diffusion is shown. Apart from a shift approximately corresponding to the nuclei's charge, the curves are similar to each other. The higher the influence of the photo-disintegration resonances the more pronounced the local minimum is for the acceleration probability; see, e.g., the model C-2 for nitrogen.  Figure B.16: Acceleration probability for hydrogen, helium, nitrogen, silicon, and iron for Kolmogorov diffusion is shown. Compared with the Bohm diffusion approach the acceleration probability is much smoother.