UltraHighEnergy Cosmic Rays from LowLuminosity Active Galactic Nuclei
Abstract
We investigate the production of ultrahighenergy cosmic ray (UHECR) in relativistic jets from lowluminosity active galactic nuclei (LLAGN). We start by proposing a model for the UHECR contribution from the black holes (BHs) in LLAGN, which present a jet power erg s. This is in contrast to the opinion that only highluminosity AGN can accelerate particles to energies EeV. We rewrite the equations which describe the synchrotron selfabsorbed emission of a nonthermal particle distribution to obtain the observed radio flux density from sources with a flatspectrum core and its relationship to the jet power. We find that the UHECR flux is dependent on the observed radio flux density, the distance to the AGN, and the BH mass, where the particle acceleration regions can be sustained by the magnetic energy extraction from the BH at the center of the AGN. We use a complete sample of 29 radio sources with a total flux density at 5 GHz greater than 0.5 Jy to make predictions for the maximum particle energy, luminosity, and flux of the UHECRs from nearby AGN. These predictions are then used in a semianalytical code developed in Mathematica (SAM code) as inputs for the MonteCarlo simulations to obtain the distribution of the arrival direction at the Earth and the energy spectrum of the UHECRs, taking into account their deflection in the intergalactic magnetic fields. For comparison, we also use the CRPropa code with the same initial conditions as for the SAM code. Importantly, to calculate the energy spectrum we also include the weighting of the UHECR flux per each UHECR source. Next, we compare the energy spectrum of the UHECRs with that obtained by the Pierre Auger Observatory.
keywords:
cosmic ray, UHECR, AGN, jet, Auger1 Introduction
Cosmic rays (CRs) are a direct sample of matter from outside the solar system, and their study can, for instance, provide important information on the chemical evolution of the universe or improve constraints on Galactic and extragalactic magnetic fields. They can be measured indirectly through the study of extensive air showers that are induced as the CRs hit the top of the atmosphere (known as CR events). The extensive air showers are currently observed using air fluorescence [e.g., High Resolution Fly’s Eye (HiRes) experiment^{3}^{3}3http://www.cosmicray.org] or large array, groundbased detectors [e.g., Akeno Giant Air Shower Array (AGASA)^{4}^{4}4http://wwwakeno.icrr.utokyo.ac.jp/AGASA], or both [e.g., Pierre Auger Observatory (Auger)^{5}^{5}5http://www.auger.org]. In the future, spacebased detectors might be another option. UHECR particles are mostly protons or fully ionized nuclei with energy above 50 EeV (1 EeV = 10 eV). At such high energies, the flux of UHECRs is very low and only a few dozen particles per square kilometer per century are expected. This is one of the main reasons for the difficulty posed in understanding the origin and nature of the UHECRs. Therefore, very large detector arrays are required. The Pierre Auger Observatory, by far the biggest cosmic ray detection instrument, uses air fluorescence and water detection in a hybrid instrument with an aperture of 7000 km sr.
Joint efforts have been made during the past decade by worldwide, cosmic ray experiments to help us understand from where the UHECRs come and what is their nature. It is believed that the UHECRs originate in extragalactic sources, as the gyroradius of a proton with an energy of 100 EeV is of the order of the dimension of our Galaxy, whereas most of the CR particles with energy below 50 EeV originate within our Galaxy (e.g., Berezinsky et al. 2006; Stanev 2010a, b). If the UHECR particles are protons, they are subject to energy loss by creating pions through their occasional collisions with the cosmic microwave background (CMB) photons. This process produces a suppression of the cosmic ray energy spectrum beyond 50 EeV, which is known as the GreisenZatsepinGuzmin (GZK) cutoff (Greisen 1966; Zatsepin & Kuzmin 1966). Therefore, the UHECRs would not be able to survive the propagation from their acceleration sites to us unless their sources are located within Mpc. The presence of the GZK cutoff at the expected energy in the data released by the HiRes collaboration was taken as strong evidence that the UHECR flux is dominated by protons (Abbasi et al. 2010).
A suppression of the CR flux has also been observed in the data released by the Pierre Auger collaboration (Abraham et al. 2008b, 2010b). With respect to primary composition, this collaboration has exploited the observation of the longitudinal shower development with fluorescence detectors to measure the depth of the maximum of the shower evolution, , which is sensitive to the primary mass. A gradual increase of the average mass of cosmic rays with energy up to 59 EeV is deduced when comparing the absolute values of and RMS() to air shower simulations (Abreu et al. 2010c).
The data collected by the Pierre Auger Observatory provide evidence for a correlation between the arrival directions of CR events above 55 EeV and the positions of AGN with (Abraham et al. 2007, 2008a; Abreu et al. 2010c, 2011; Aab et al. 2013b), where the region around the position of the radiogalaxy Cen A has the largest excess of arrival directions relative to the isotropic expectations. The correlation is shown for a selection of AGN from the catalog of VéronCetty & Véron (2006), which do not necessarily follow the same structure as the gammaray bursts (GRBs). We emphasize that there is no clear detection of the UHECR sources, just a strong evidence for the anisotropy in the arrival directions of UHECRs.
At highest energies, heavy nuclei may be deflected by Galactic magnetic fields, whereas proton propagation is affected by the CMB, as well as by magnetic deflection, though to a less degree compared to that of particles with higher mass number (e.g., Medina Tanco et al. 1998).
UHECRs are most probably accelerated at astrophysical shocks, for instance, through a firstorder Fermi mechanism (e.g., Gallant & Achterberg 1999), in very powerful systems that can be associated with jets and hot spots in AGN and GRBs (Waxman 1995; Vietri 1995), in largescale shocks in clusters (e.g., Farrar et al. 2006), or as iron nuclei in pulsar winds from rapidlyspinning, young neutron stars (Blasi et al. 2000; Fang et al. 2012). Numerical simulations of particle acceleration in shocks have been widely performed using different values for the shock Lorentz factor and background conditions at the shock front (e.g., Bednarz & Ostrowski 1998; Achterberg et al. 2001; Kirk et al. 2000; Keshet & Waxman 2005; Niemiec & Ostrowski 2006; Meli et al. 2008), which lead to a slope of the particle distribution of . Such shocks can also be associated with Poynting flux models for the origin of jets from forcefree magnetosphere above accretion disks (e.g., Lovelace 1976; Blandford 1976; Boldt & Ghosh 1999; Blandford 2000; Biermann et al. 2008). Magnetic reconnection in relativistic jets represents another option for UHECR acceleration cite[e.g.,][]giannios10. As an alternative, Farrar & Gruzinov (2009) showed that very intense, shortduration AGN flares that result from the tidal disruption of a star or from a disk instability can accelerate UHECRs. (See also Waxman & Loeb (2009).)
In this paper, we propose a model for the UHECR contribution from relativistic jets in LLAGN and calculate the expected energy spectrum of UHECRs using the SAM code developed by Biermann et al. (2008), Caramete et al. (2011a). For comparison, we also employ the CRPropa code, which is set up with the same initial conditions as the SAM code. The particles in the jet are powered by the BH accretion disk and then accelerated at relativistic shocks with energies up to the ultrahigh energy (UHE) domain. We limit the launching area of the jets to the innermost part of the disk located inside the BH ergosphere, where the rotational effects of the spacetime are very strong. There are several general relativistic magnetohydrodynamic (GRMHD) codes, the result of which show that the jets can be magnetically driven from a thin disk located inside the BH ergosphere via a Penroselike process (Koide et al. 1999; Nishikawa et al. 2005; Mizuno et al. 2007) or via the BlandfordZnajek mechanism (BZ, Blandford & Znajek (1977)) when a thick accretion disk is considered (McKinney 2006; McKinney & Blandford 2009). (But see the simulations by Fragile et al. (2012), where the BZ driven jet does not depend on the disk thickness.) In contrast to the BZ mechanism, where the power of the jet is proportional to the square of the BH spin (), in the model presented here the dependence of on comes through the launching area of the jets (see the equations in Appendices A and C). In the jets, the electrons lose their energy through synchrotron emission, whereas the protons, as well as heavy nuclei (here, iron nuclei), are capable of surviving the radiative cooling and, perhaps, of propagating through the intergalactic and Galactic medium towards us. Since the particle species undergo the same acceleration process, there must be a correlation between the electron synchrotron emission and the energy of the UHECR particles (protons and iron nuclei). We seek this correlation to make predictions for maximum energy, luminosity, and flux of the UHECRs from nearby LLAGN. This is in contrast to the opinion that only highluminosity AGN can accelerate particles to UHE domain (e.g., Zaw et al. 2009). Taking into account the deflection of the trajectories of the UHECRs in the intergalactic and Galactic magnetic fields, we calculate the distribution of the arrival direction at the Earth and the energy spectrum of the UHECRs. The latter is then compared with the energy spectrum obtained by the Pierre Auger Observatory. We point out that LLAGN as sources of UHECRs were also proposed by Moskalenko et al. (2009), where discussions about the implication of AGN jet power and intergalactic, magnetic field configurations for the observed statistical correlation between AGN and UHECR events are presented. Our work is a step further to that of Moskalenko et al. (2009), as we include quantitative estimations of UHECR flux using its correlation to the AGN jet power, as well as simulations of UHECR particle propagation in the intergalactic and Galactic magnetic fields. To obtain the energy spectrum of UHECRs, we use a complete sample of 29 LLAGN with a total radio flux density larger than 0.5 Jy (Biermann et al. 2008; Caramete et al. 2011a). About 80% of our sample is contained in the allsky catalog of local radio galaxies of van Velzen et al. (2012), which is used to seek for the correlation between the UHECRs and LLAGN (van Velzen & Falcke 2013). The fact that some sources of our sample are not included in the catalog by van Velzen et al. (2012) might be attributed to the difference in the data; i.e., the frequency at which the radio flux density was measured: 5 GHz in our case and 1.5 GHz and 843 MHz in the case of van Velzen et al. (2012).
In Section 2, we provide a description of the model for the UHECR contribution from relativistic jets in LLAGN. We derive the luminosity and flux of the UHECRs based on the relation between the jet power and the observed radio flux density for a flatspectrum core source (see C) and calculate the particle maximum energy taking into account the spatial limit and synchrotron emission losses. In Section 3, we provide the predictions for nearby galaxies as possible sources of UHECRs by employing the SAM code. For comparison, we also use the CRPropa code with the same input setup as for the SAM code. We then compare the results of the two codes with those obtained by the Pierre Auger Observatory. In Section 4, we present a summary of the key points and discuss the implication of this model for further studies of UHECRs.
2 Description of the model for UHECR source
2.1 Model conditions

We assume that the UHECRs are accelerated by shocks in AGN jets, which are launched from the inner accretion disk which is located inside the BH ergosphere, where the rotational effects of the spacetime become much stronger (Duţan & Biermann 2005). (The inner disk extends from the stationary limit inward to the innermost stable orbit .) In the model by Duţan & Biermann (2005); Duţan (2010), the rotation of the spacetime channels a fraction of the disk energy (i.e., the gravitational energy of the disk plus the rotational energy of the BH that is deposited into the disk by closed magnetic field lines, which connect the BH to the accretion disk) into a population of particles that escape from the disk surfaces, carrying away mass, energy, and angular momentum in the form of jets, allowing the remaining disk gas to accrete. The power of the jets can have two components; the one which comes from the accretion power and the other one which comes from the BH rotational energy. The latter component dominates when the mass accretion rate^{6}^{6}6Here, the term mass accretion rate refers to the mass flow rate through the disk up to the ergosphere (), where is the Eddington accretion rate. is defined from the Eddington luminosity as , where is the efficiency of converting the gravitational energy of the accretion disk into radiation and denotes the Thomson opacity. depends on the BH spin parameter as (Thorne 1974), so that for a Schwarzschild BH and for an extremely spinning Kerr BH. If we scale the BH mass to , then , where g s. relative to the Eddington accretion rate is below .

Being driven from the accretion disk, the jet propagates along a cylinder of length (see Fig. 1), where the poloidal and toroidal components of the comoving magnetic field become approximately equal, and then extends into a conical shape with a constant opening angle , as a consequence of the free adiabatic expansion of the jet plasma (similar to Markoff et al. (2001)). The cylindrical surface represents the envelope of the magnetic field lines, which close to the BH are mainly parallel. (We do not include here a monopolelike configuration of the magnetic field.) We suppose that a shock is produced at the jet height . (Production of a shock associated with the twisting of the magnetic field lines, where the toroidal component of the magnetic field dominates, was observed in GRMHD simulations. See, e.g., fig 3 in Mizuno et al. (2004).)

As a result of the shock produced at , a powerlaw energy distribution of the particles is established. The number density of the electrons in the energy interval , [or , )] has, in terms of the Lorentz factor, the following powerlaw form: . In the case of a conical jet, the normalization of the electron number density is (Blandford & Königl 1979):
(1) We set the slope of the particle energy distribution to , and .

The calculations are performed for the case when the UHECRs are composed of 90% protons and 10% iron nuclei (e.g., Allard et al. 2008).

We consider the strength of the BH magnetic field to have its maximum value. This condition provides, in turn, the minimum values of the particle maximum energy, luminosity, and flux of the UHECRs.
2.2 Magnetic field scaling along a steady jet
To describe the jet physics, we use the following reference frames: (i) the frame comoving with the jet and (ii) the (rest) frame of the observer, in which the relativistic jet moves with the bulk Lorentz factor. In a frame comoving with the jet, the poloidal component of the magnetic field is considered to vary as . This variation follows from the conservation of magnetic flux along the axis . To keep the field divergencefree, the toroidal component must vary as . This topology of was first derived by Parker (1958) for the magnetohydrodynamics solution of a sphericalsymmetric flow (so that, a jet can be considered a conical cut along the flow surfaces). [See also Blandford & Königl (1979).] At the distance , the poloidal and toroidal components of the comoving magnetic field become approximately equal . The distance might be a few tens of gravitational radii,^{7}^{7}7The gravitational radius is defined as cm, where is the Newtonian gravitational constant, is the BH mass, and is the speed of light. based on the fact that the VLBI observation, for instance, of the jet in M87 at 43 GHz gives evidence on the jet collimation (by the toroidal magnetic field) on scales of 60200 (Biretta et al. 2002). More recently, Hada et al. (2011) argued that 43 GHz VLBI core is located at from the central engine, where they performed the core shift measurement by using multifrequency, phasereferencing Very Long Baseline Array observations, and that the measured frequency dependence of the core shift is in good agreement with a synchrotron selfabsorbed jets. The structure of the M 87 jet, from submiliarsec to arcsec scales, was also investigated by Nakamura & Asada (2013), where the (bulk) acceleration and collimation of the jet are correlated in a parabolic streamline. The question is whether the streamline can be extrapolated continuously towards the jet launching region. To answer such question future observations with high sensitivity, e.g., a space VLBI and/or VLBI at higher frequency than 86 GHz are needed. Next, if we consider the ray emission from M 87 and Cen A, the absorption of the TeV photons produced at about and , respectively, from the BH is small, which means that the observed radiation could come from the BH vicinity (Acciari et al. 2009; Brodatzki et al. 2011; Saba et al. 2013). Nevertheless, in calculating the luminosity and flux of the UHECRs we do not make use of the exact location of , we rather extract the ratio of from the expression of the optical depth (Eq. 23). A largescale and predominantly toroidal magnetic field can exert an inward force (hoop stress), confining and collimating the jet (e.g., BisnovatyiKogan & Ruzmaikin 1976; Blandford & Payne 1982). The magnetic hoop stress is balanced either by the gas pressure of the jet or by centrifugal force if the jet is spinning. From upward, the poloidal component of the magnetic field becomes weaker, so that the field lines are soon wound up in the azimuthal direction by the jet rotation. Thus, above , the magnetic field along the jet is nearly azimuthal (for a steady jet) and varies inversely proportional to the distance along the jet:
(2) 
where is the strength of the magnetic field at the height above the equatorial plane of the BH. This dependence of the magnetic field appears to be contradicted by radiopolarization observations (Bridle & Perley 1984). These observations strongly suggest that the magnetic field is predominantly parallel to the jet axis initially and only later becomes perpendicular to the jet axis, with some parallel magnetic field left over. Becker & Biermann (2009) argued that the basic pattern of the magnetic field is indeed and that the observational evidence for a parallel magnetic field is due to highly oblique shocks. Their argument is based on the observations of the jet structure which might be explained through the occurrence of the moving shocks between and (Marscher et al. 2008), while the first stationary, strong shock can be produced in the approximate range of (Markoff et al. 2001, 2005) or as seen in blazars (Marscher et al. 2008).
The strength of the magnetic field in the comoving frame can be related to the poloidal magnetic field in the BH frame (e.g., Drenkhahn 2002) as
(3) 
where the maximum value of the BH magnetic field is given by
(4) 
which is obtained in a similar manner as the calculation performed by Znajek (1978), with the difference that we set the BH potential drop to the specific energy of the particles at the innermost stable orbit, whereas Znajek (1978) makes use of the fact that the Eddington luminosity sets an upper bound on the radiation pressure (as the disk is radiatively efficient). The maximum value of the BH magnetic field corresponds to the time when the accretion rate was as high as the Eddington accretion rate. In this case, the BH spin parameter^{8}^{8}8The BH spin parameter is defined as , where is the angular momentum of the BH about the spinning axis per unit mass and per speed of light and is the maximal angular momentum of the BH. Furthermore, the BH spin parameter obeys the condition: . is limited to (Thorne 1974). Although this limit might be even closer to the maximal value of the spin parameter , this will introduce just a small variation of the maximum value of the BH magnetic field.
2.3 Luminosity and flux of the UHECR
In this section, we seek for the UHECR luminosity flux () specified as a function of the observed radio flux density. First, we consider the UHECR luminosity defined as
(5) 
where it is assumed that the UHECR luminosity is a fraction () of the jet power, with . If we were to adopt the point of view that the jet power is shared equally in a comoving frame between the baryonic matter, magnetic field, and cosmic rays extending to the highest energy, . In the jetdisk model of Falcke & Biermann (1995), the energy equipartition in the comoving frame appears to be a good approximation. It would also suggest that AGN driven by the BH rotational energy, and suffering from a low mass accretion rate, may attain a higher Lorentz factor, consistent with some observations.
Now, we can obtain the cosmic ray luminosity by dividing the expression of the jet power in Eq. 35 by 3. In section 2.5, we write down the UHECR luminosity for the case of . Given the UHECR luminosity, we can easily obtain the UHECR flux:
(6) 
where we do not include the cosmological distance as we refer to nearby radio sources with a flatspectrum core and a redshift down to .
2.4 Maximum particle energy of the UHECR
Now, we look for the maximum energy of the UHECR in the case of the spatial (geometrical) limit (Falcke & Biermann 1995); i.e., the jet particle orbits must fit into the Larmor radius. Conform to Gallant & Achterberg (1999), the maximum particle energy in the downstream rest frame can be written as
(7) 
where is the Lorentz factor of the shock and is the particle mass number. Equation 7 corresponds to relativistic shocks, being larger by a factor than that resulting from a conventional geometrical comparison of the gyration radius (Hillas 1984).
Using the expression for the magnetic field along the jet (Eqs. 2 and 3) and the fact that , the maximum energy of the UHECR particles (in the observer frame) becomes:
(8) 
For protons,
(9) 
where was used.
Next, we look for the maximum energy of the UHECR in the case of the synchrotron loss limit (Biermann & Strittmatter 1987). Setting synchrotron losses equal to diffusive shock acceleration gains, Biermann & Strittmatter (1987) showed that a ubiquitous cutoff in the nonthermal emission spectra of AGN can be explained. This requires that the protons to be accelerated near eV. The frequency cutoff () might be produced at about . Rewriting the expression for the maximal proton energy derived by Biermann & Strittmatter (1987),
(10) 
and using the expression for the magnetic field along the jet (Eqs. 2 and 3), the maximal proton energy in the loss limit reads:
(11) 
2.5 Model set of equations
Taking and , the equations for the maximum particle energy in the spatial (Eq. 9) and loss (Eq. 11) limits, as well as for the UHECR luminosity (Eq. 5 with 35) for , become:
(12) 
For the expression of in Eq. (12), we can choose , based on the results obtained by Becker & Biermann (2009), which show that a first large steady shock can be produced at about [following the work by Markoff et al. (2001)]. Nevertheless, for the calculation of the UHECR energy spectrum, we use the expression of , which provides higher values for the particle maximum energy. We note that if , we should make use of .
3 Predictions for nearby galaxies as ultrahighenergy cosmic ray sources
In order to formulate predictions for the sources of UHECRs, we use a twostep analysis, as follows:

we elaborate a physical model of the UHECR sources (previous section) and then

we consider the deflection of the trajectories of the UHECRs by the intergalactic magnetic fields and calculate the distribution of the arriving directions and the energy spectrum of the UHECRs using both the SAM code and the CRPropa code (see below).
We apply the model proposed in this paper to a complete sample of extended, steep spectrum radio sources^{9}^{9}9For radio galaxies, a flatspectrum core source has a spectral index which is typically at the very center and at a larger radius (the compact radio core can be extremely weak compared with the very bright and luminous radio lobes, suggesting that the radio core might suffer absorption). The latter spectral index corresponds to a powerlaw index of the accelerated particle of . The extended emission, which includes the emission from the radio lobes, can have a steep spectrum. (Biermann et al. 2008; Caramete et al. 2011a), at redshift (about 100 Mpc), with a total radio flux density larger than 0.5 Jy. The numbers for the estimated flux and maximal energy exclude the GZK effect, but includes the distance effect. The selection criteria used by the authors are presented in more detail in their papers. Table 1 lists the predictions for the maximum energy and flux of the UHECRs. To calculate the errors of these quantities, we use the method of propagation of uncertainty. We emphasize that there could be a common scaling limit, such as a condition that the Larmor radius has to fit three times or five times into the jet. The scaling limit is not critical to our predictions as long as we refer the quantities to, say, those of M 87; therefore, the jet parameters (, , and ) for all the sources in Table 1 are assumed to be the same as for M 87. In the simulations described below, we use columns 5 and 7 of Table 1.
For the SAM code, we start with a sufficient initial number of particles at the level of the detected UHECR flux in the selected energy interval (above eV) and then distribute these particles to each source in the list (Table 1) taking into account the calculated flux of the UHECRs at the source. This is performed by using the MonteCarlo method in which we consider the calculated UHECR fluxes as weights to randomly distribute the particles to each source in the list. This is equivalent to drawing numbers from a distribution according to a list of weights which is attached to the distribution. In this way, the most active source produces more UHECRs and, therefore, we can relate such source with more particles from the total number of the injected particles (that propagate through intergalactic and Galactic magnetic fields). Next, to associate a value of the energy to each UHECR, we randomly generate a range for the energy of the particles that come from each source using three power laws: , and , where we consider the estimated maximal energy of the particles as an upper limit of the energy range. Here, we assume the initial distribution of the UHECRs at the source to be constituted of 90% protons and 10% iron nuclei, as considered by Allard et al. (2008). Future upgrade of the SAM code will include the interactions of the UHECR particles with different backgrounds (e.g., CMB, IR, Optical, and UV radiation).
To obtain the energy spectrum and the distribution of the arrival directions at the Earth of the UHECRs, we estimate the deflections in the intergalactic magnetic fields using the method applied in the numerical simulations performed by (Das et al. 2008). This method provides a distribution of deflection angles of protons in the intergalactic magnetic fields corresponding to the following intervals of the energy of the particles: (i) from EeV to EeV, (ii) from EeV to EeV, and (iii) beyond EeV. When taking into account the deflection of UHECR trajectories by the magnetic fields, we have to consider separately the effects produced by the intergalactic magnetic fields (Ryu et al. 1998; Dolag et al. 2005; Ryu et al. 2008) from those produced by the magnetic field of our Galaxy (Breitschwerdt et al. 1991; Everett et al. 2008, 2010). For the SAM code, we use the scattering laws from the Fig.7 of Das et al. (2008) with the following fitting formulas:
where denotes the deflection angle and represents the fraction of events per deflection angle bin, which can be easily translated into a probability to have a certain number of events with a specific deviation angle. Moreover, the core of the distribution of angle deflection is chosen such that to simulate deflections in the Galactic disk, while the rest of the distribution of the deflection angle reflects either the scattering in the cosmological magnetic fields or the scattering in a Galactic magnetic wind which extends into the halo of our Galaxy.
Again, using a MonteCarlo weighted choice for each energy range, we associate to each UHECR a deflection angle from the direction of the source on the Earth sky, and by uniformly random selection we chose a position of the UHECR event in a ring around the source. This leads us to obtaining a sky map of the distribution of the arrival directions of the UHECRs (Fig. 2). We can easily notice that almost all particles come from Cen A, which is the most active source in the list in terms of the UHECR flux. Near the position of this source, a clustering of events is observed, a picture that resembles the one obtained by the Pierre Auger Observatory (Abraham et al. 2007). Such clear cluster of arrival directions of UHECR events around the Cen A point to their source.
To consistently check the results of the SAM code, we compare them with those obtained from simulations that we perform with the CRPropa code (Kampert et al. 2013). The CRPropa code is freely available and designed to study the propagation of UHECRs in the intergalactic space by including the effects of background interactions and magnetic fields. For the simulations run with the CRPropa code, we use the same initial conditions as for those performed with the SAM code; that is, the same list of sources, initial composition, and initial flux of UHECRs. However, there are two main differences between the two simulation codes that concern (i) the interaction with the backgrounds (which is present in the CRPropa code, but not in the SAM one) and (ii) the models for the intervening magnetic fields between the sources of UHECRs and the Earth’s atmosphere (for the SAM code, we consider both the extragalactic and the Galactic magnetic fields, as described above, whereas for the CRPropa code, we use the provided standard intergalactic magnetic field. In the present public version of the CRPropa code, the Galactic magnetic field is not included).
Figure 3 shows the energy spectrum resulted from the simulations performed with the SAM code (the grey area), which is compared with that obtained with the CRPropa code (red crosses) and with the observed by the Pierre Auger Collaboration (blue dots). Here, we show only the plot corresponding to the injected energy spectrum of particles of , which provides the best approximation to the observed spectrum. The other two slopes ( and ) generate much higher energy spectra for UHECRs than . In order to estimate the simulation errors, a set of 1000 simulations are performed and the global error is computed according to D. Here, we use from the all sky distribution of arriving particles only the ones which fall on the general field of view of the Pierre Auger detector. We note that the experiment makes a certain cutoff in the events that it detects which depends on many parameters like the arrival angle of the particles with respect of the Pierre Auger detector or other factors which are not considered here (for a discussion see Abreu et al. (2011)). These factors could explain the large errors that appear at the highest energies, as well as the offset of the end points. However, a slight overestimate of the maximal energy of UHECRs at the source might also account for the large errors at the highest energies in the spectrum. Within errors, the simulation results obtained with the SAM code are in good agreement with the observed energy spectrum. The CRPropa simulation provides an energy spectrum which, for the energy range eV, is slightly above both the SAM energy spectrum and the observed one. We correlate this discrepancy with the fact the CRPropa code does not include deflection of particles in the Galactic magnetic field.
Using a special developed representation software in Mathematica, we also include here a 3D representation of the trajectories of UHECRs (Fig. 4), generated from the 29 LLAGN in Table 1. We obtain this representation from simulations with the CRPropa code where we use the same initial setup as for the SAM code. We marked in red the trajectories of the particles that are generated from Cen A. The simulation box size is Mpc, where our Galaxy is located at the center of the box.
4 Summary and conclusions
In this paper, we address the possibility that nearby LLAGN can be sites of production of UHECRs. The analysis that we perform in the work presented here is twofold:

First, we elaborate a physical model of the source of UHECRs, where the particles are accelerated in LLAGN jets, for which the acceleration regions can be sustained by the magnetic energy extraction from the spinning BH at the heart of the AGN. We relate the observed radio flux density to the luminosity and flux of the UHECRs and calculate the maximum particle energy in both spatial and loss limits. Next, we calculate the energy and flux of UHECR that can be generated from a complete sample of 29 LLAGN (Table 1). The results of these calculations are then used as initial conditions for the simulations performed with both the SAM code and the CRPropa one.

Second, using the SAM code we estimate the deflection of the trajectories of the UHECRs (90% protons and 10% iron nuclei) by the intergalactic and Galactic magnetic fields, using appropriate distributions of the angle deflections (Section 3) and calculate the distribution of the arriving directions and the energy spectrum of the UHECR particles. We also compare the energy spectrum obtained from the SAM simulations with that resulted from the simulations with the CRPropa code and from the measurements performed by the Pierre Auger Observatory.
From the simulations performed with the SAM code and CRPropa, the best approximation to the observed spectrum measured by the Pierre Auger Observatory is obtained for the injected energy spectrum of particles corresponding to . Although the two simulated spectra are in good agreement with the observed spectrum at the highest energies, the spectrum obtained with the CRPropa code overestimates the one measured by the Pierre Auger Observatory in the energy range eV. The different behavior of the simulated energy spectra care be attributed to the differences between the two codes. The main difference between the codes regards the deflection of the particles in the Galactic magnetic field, which is not included in the current public version of the CRPropa code, whereas such deflection is already implemented in the SAM code. More importantly, to calculate the energy spectrum, we also include the weighting of the UHECR flux per each source; i.e., we attach to each UHECR source a weigh in flux based on the value of its UHECR flux relative to a canonical value of the flux, which we choose it to be the flux corresponding to M 87. (The use of the relative flux to a canonical value can eliminate the systematic errors.) This weighting of the UHECR flux is not taken into account in the CRPropa code, where all sources can have the same UHECR flux. Therefore, as in the CRPropa code there is not possible to run a single simulation with different number of particles per each source, we have to run separate simulations per each UHECR source (or sublist of sources) from the sample in Table 1. Thus, we run (i) one simulation for Cen A using of the total injected particle number (), (ii) one simulation for M 87 and the sources with an UHECR flux at the same level of that of M 87 using of , and (iii) one simulation for the rest of the sources in Table 1 using the remaining of . Next, we merge these three simulations together to obtain the energy spectrum of UHECRs. The SAM code has the weighting of the UHECR flux already implemented, so that we just use the UHECR flux provided by the model to give the initial number of particle per source.
Although we presented a simple model for relativistic jets in LLAGN as sites of production of UHECRs, without including complex phenomena (e.g., magnetic reconnection and instabilities in the jet plasma) that can change the spectrum of the emitted radiation from the jets, the uncertainty on several parameters of the model are still present. For example, we do not know exactly at which distance from the BH a strong shock can be produced. We might overcome this problem by studying each source individually and using observations of high resolution, which may not be available. However, in this paper we have to make assumptions based on the performed measurements. We also do not have observational data for the jet parameters (the Lorentz factor, the angle to the line of site, and the semiopening angle) of the all sources in Table 1. We supposed that these jet parameters have similar values as those of M 87. (We note that the difference between the value of the UHECR flux using the observed parameters of the jet for Cen A and that for Cen A but using the parameters of the jet for M 87 is within one order of magnitude.) Thus, the maximum particle energy and flux of UHECRs will be scale only by the BH mass, the radio power, and the distance to the source:
where the errors are given by the uncertainties in the radio core flux density, the distance to the AGN, and the BH mass measurements. Therefore high resolution observations in different frequency domains are needed. Moreover, increasing the number of BHs with measured mass, the distribution function of the BH mass of Caramete & Biermann (2010, 2011b) can be improved and then use to determine the masses of supermassive BHs more precisely.
We plan to extend the analysis presented in this paper. We will implement the calculations for neutrino and gammaray energy spectra. We also intend to use a different sample of LLAGN, e.g., the catalog by van Velzen et al. (2012) and to modify the chemical composition of the UHECRs. These calculations can provide a further test on the exact composition of the UHECRs as different particles behave differently when they propagates through the intergalactic and Galactic magnetic fields.
Acknowledgments
We would like to thank Peter L. Biermann. We also thank the anonymous referees for their thorough review and highly appreciate their comments and suggestions, which significantly contributed to improving the quality of the paper. I.D. was supported through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. She appreciates the support from MPIfR during the last phase of this work. I.D. was also supported by STAR, CONTR. 12/19.11.2012, and L.I.C. was partially supported by PNIIRUTE201130184 and STAR, CONTR. 12/19.11.2012. We make use of the NASA/IPAC Extragalactic Database (ned.ipac.caltech.edu) and the CRPropa simulation code (https://crpropa.desy.de/Main_Page).
References
 Aab et al. (2013b) Aab, A. et al. (Pierre Auger Collaboration) 2013b, arXiv: astroph/1307.5059
 Abbasi et al. (2010) Abbasi, R. et al. (HiRes Collaboration) 2010, Phys. Rev. Lett., 104, id. 161101
 Acciari et al. (2009) Acciari V. A. et al. (VERITAS, MAGIC, VLBA M87 and H.E.S.S. Collaborations) 2009, Science, 325: 444
 Achterberg et al. (2001) Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W. 2001, MNRAS, 328, 393
 Allard et al. (2008) Allard, D., Busca, N. G., Decerprit, G., Olinto, A. V., & Parizot, E. 2008, JCAP, 10, 33
 Abraham et al. (2007) Abraham, J. et al. (Pierre Auger Collaboration) 2007, Science, 318, 938
 Abraham et al. (2008a) Abraham, J. et al. (Pierre Auger Collaboration) 2008a, Astropart. Phys., 29, 188
 Abraham et al. (2008b) Abraham, J. et al. (Pierre Auger Collaboration) 2008b, Phys. Rev. Lett., 101, id. 061101
 Abraham et al. (2010a) Abraham, J. et al. (Pierre Auger Collaboration) 2010a, Phys. Lett. B, 685, 239
 Abraham et al. (2010b) Abraham, J. et al. (Pierre Auger Collaboration) 2010b, Phys. Rev. Lett., 104, id. 091101
 Abreu et al. (2010c) Abreu, P. et al. (Pierre Auger Collaboration) 2010c, Astropart. Phys., 34, 314
 Abreu et al. (2011) Abreu, P. et al. (Pierre Auger Collaboration) 2011, arXiv: astroph/1107.4809
 Abreu et al. (2013a) Abreu, P. et al. (Pierre Auger Collaboration) 2013a, ApJL, 762, L13
 Bardeen (1970) Bardeen, J. M. 1970, Nature, 226, 64
 Becker & Biermann (2009) Becker, J. K. & Biermann, P. L. 2009, APh, 31, 138
 Bednarz & Ostrowski (1998) Bednarz, J. & Ostrowski, M. 1998, Phys. Rev. Lett., 80, 3911
 Berezinsky et al. (2006) Berezinsky, V., Gazizov, A., & Grigorieva, S. 2006, Phys. Rev. D, 74, id. 043005
 Biermann et al. (2008) Biermann, P. L. et al. 2008, Proc. for Origin, Mass, Composition and Acceleration Mechanisms of UHECRs (CRIS 2008), Malfa, Italy (arXiv: astroph/0811.1848)
 Biermann & Strittmatter (1987) Biermann, P. L. & Strittmatter, P. A. 1987, ApJ, 322, 643
 Biretta et al. (2002) Biretta, J. A., Junor, W., & Livio, M. 2002, NewAR, 46, 239
 Biretta et al. (1999) Biretta, J. A., Sparks, W. B., & Macchetto, F. 1999, ApJ, 123, 351
 BisnovatyiKogan & Ruzmaikin (1976) BisnovatyiKogan, G. S. & Ruzmaikin, A. A. 1976, Ap&SS, 42, 401
 Blandford (1976) Blandford, R. D. 1976, MNRAS, 176, 465
 Blandford (2000) Blandford, R. D. 2000, Physica Scripta, T85, 191
 Blandford & Königl (1979) Blandford, R. D. & Königl, A. 1979, ApJ, 232, 34
 Blandford & Payne (1982) Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883
 Blandford & Znajek (1977) Blandford R. D. & Znajek R. L. 1977, MNRAS, 179, 433
 Blasi et al. (2000) Blasi, P., Epstein, R., & Olinto, A. V. 2000, ApJ, 533, L123
 Boldt & Ghosh (1999) Boldt, E. & Ghosh, P. 1999, MNRAS, 307, 491
 Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Völk, H. J. 1991, A&A, 245, 79
 Bridle & Perley (1984) Bridle, A. & Perley, A. 1984, ARAA, 22, 319
 Brodatzki et al. (2011) Brodatzki, K. A., Pardy, D. J. S., Becker, J. K., & Schlickeiser, R.
 Cappellari et al. (2009) Cappellari, M. et al. 2009, MNRAS, 394, 660
 Caramete et al. (2011a) Caramete, L. I. et al. 2011, arXiv: astroph/1106.5109
 Caramete & Biermann (2011b) Caramete, L. I. & Biermann, P. L. 2011, arXiv: astroph/1107.2244
 Caramete & Biermann (2010) Caramete, L. I. & Biermann, P. L. 2010, A&A, 521, 8
 Das et al. (2008) Das, S., Kang, H., Ryu, D., & Cho, J. 2008, ApJ, 682, 29
 Dolag et al. (2005) Dolag, K., Grasso, D., Springel, V., & Tkachev, I. 2005, JCAP, 1, 9
 Drenkhahn (2002) Drenkhahn, G. 2002, A&A, 387, 714
 Duţan (2010) Duţan, I., arXiv: astroph/1001:5434
 Duţan (2011) Duţan, I., arXiv: astroph/1104.0825
 Duţan & Biermann (2005) Duţan, I. & Biermann, P. L. 2005, in Proc. of the 14th Course of the ISCRA, Erice, Italy, 2004. Eds. M. M. Shapiro, S. Todor, and J. P. Wefel. NATO science series II: Mathematics, physics and chemistry, vol. 209, p.175
 Everett et al. (2010) Everett, J. E., Schiller, Q. G., & Zweibel, E. G. 2010, ApJ, 711, 13
 Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, ApJ, 674, 258
 Falcke & Biermann (1995) Falcke, H. & Biermann, P. L. 1995, A&A, 293, 665
 Falcke et al. (1995) Falcke, H., Malkan, M. A., & Biermann, P. L. 1995, A&A, 298, 375
 Fang et al. (2012) Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 16
 Farrar et al. (2006) Farrar, G. R., Berlind, A. A., & Hogg, D. V. 2006, ApJ, 642, L89
 Farrar & Gruzinov (2009) Farrar, G. R. & Gruzinov, A. 2009, ApJ, 693, 329
 Fragile et al. (2012) Fragile, P. C., Wilson, G., & Rodriguez, R. 2012, MNRAS, 424, 524
 Gallant & Achterberg (1999) Gallant, Y. A. & Achterberg, A. 1999, MNRAS, 305, L6
 Giannios (2010) Giannios, D. 2010, MNRAS, 408, L46
 Giovannini et al. (2001) Giovannini, G., Cotton, W. D., Feretti, L., Lara, L., & Venturi, T. 2001, ApJ, 552, 508
 Greisen (1966) Greisen, K. 1966, Phys. Rev. Lett., 16, 748
 Hada et al. (2011) Hada, K. et al. 2011, Nature, 477 , 185
 Heinz & Sunyaev (2003) Heinz, S. & Sunyaev, R. A. 2003, MNRAS, 343, L59
 Hillas (1984) Hillas, A. M. 1984, ARA&A, 22, 425
 Kaiser (2006) Kaiser, C. R. 2006, MNRAS, 367, 1083
 Kampert et al. (2013) Kampert, K.H., Kulbartz, J., Maccione, L., et al. 2013, Astroparticle Physics, 42, 41
 Keshet & Waxman (2005) Keshet, U. & Waxman, E. 2005, Phys. Rev. Lett., 94, id. 111102
 Kirk et al. (2000) Kirk, J. G., Guthmann, A. W., Gallant, Y. A., & Achterberg, A. 2000, ApJ, 542, 235
 Koide et al. (1999) Koide S., Shibata K., & Kudoh T. 1999, ApJ, 522, 727
 Lee et al. (2008) Lee, S.S. et al. 2008, ApJ, 136, 159
 Longair (1994) Longair, M. 1994, High Energy Astrophysics, Vol. 2 (Cambridge Uni. Press), 259
 Lovelace (1976) Lovelace, R. V. E. 1976, Nature, 262, 649
 Macchetto et al. (1997) Macchetto, F. et al. 1997, ApJ, 489, 579
 Macri et al. (1999) Macri, L. M. et al. 1999, ApJ, 521, 155
 Marconi & Hunt (2003) Marconi, A. & Hunt, L. K. 2003, ApJ, 589, 21
 Markoff et al. (2001) Markoff, S., Falcke, H., & Fender, R. 2001, A&A, 372, L25
 Markoff et al. (2005) Markoff, S., Nowak, M. A., & Wilms, J. 2005, ApJ, 635, 1203
 Marscher et al. (2008) Marscher, A. P. et al. 2008, Nature, 452, 966
 McKinney (2006) McKinney J. C. 2006, MNRAS, 368, 1561
 McKinney & Blandford (2009) McKinney J. C. & Blandford R. D. 2009, MNRAS, 394, L126
 Medina Tanco et al. (1998) Medina Tanco, G. A., de Gouveia dal Pino, E. M., & Horvath, J. E. 1998, ApJ, 492, 200
 Meisenheimer et al. (2007) Meisenheimer, K. et al. 2007, A&A, 471, 453
 Meli et al. (2008) Meli, A., Becker, J. K., & Quenby, J. J. 2008, A&A, 492, 323
 Mizuno et al. (2007) Mizuno Y., Nishikawa N., Koide S., Hardee P., & Fishman G. J. 2007, in Proc. of the VI Microquasar Workshop: Microquasars and Beyond, Como, Italy. Ed. B. Tomasso (Trieste: SISSA), p.45
 Mizuno et al. (2004) Mizuno Y., Yamada, S., Koide, S., & Shibata, K. 2004, 615, 389
 Moskalenko et al. (2009) Moskalenko, I. V., Stawarz, L., Porter, T. A., & Cheung, C. C. 2009, ApJ, 693, 1261
 Nagar et al. (2001a) Nagar, N. M., Wilson, A. S., & Falcke, H. 2001, ApJ, 559, L87
 Nagar et al. (2001b) Nagar, N. M., Falcke, H., & Wilson, A. S. 2005, A&A 435, 521
 Nakamura & Asada (2013) Nakamura, M. & Asada, K. 2013, ApJ, 775, 11
 Niemiec & Ostrowski (2006) Niemiec, J. & Ostrowski, M. 2006, ApJ, 641, 984
 Nishikawa et al. (2005) Nishikawa K.I., Richardson G., Koide S., Shibata K., Kudoh T., Hardee P., & Fishman G. J. 2005, ApJ, 625, 60
 Novikov & Thorne (1973) Novikov, I. D. & Thorne, K. S. 1973, in Black Holes Les Astres Occlus, Ed. DeWitt, C. and DeWitt, B. S. (NewYork: Gordon & Breach), 359
 Parker (1958) Parker, R. D. 1958, ApJ, 128, 664
 Protheroe & Biermann (1996) Protheroe, R. J. & Biermann, P. L. 1996, Astropart. Phys., 6, 45
 Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, Radiative Processes in Astrophysics (John Wiley & Sons)
 Ryu et al. (1998) Ryu, D., Kang, H., & Biermann, P. L. 1998, A&A, 335, 19
 Ryu et al. (2008) Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Science, 320, 909
 Saba et al. (2013) Saba, I., Becker Tjus, J. K., & Halzen, F. 2013, Astropart. Phys., 48, 30
 Sanders (1983) Sanders, R. H. 1983, ApJ, 266, 73
 Sigl et al. (2003) Sigl, G., Miniati, F., & Ensslin, T. A. 2003, Phys. Rev. D, 68, 043002
 Slee et al. (1994) Slee, O. B. et al. 1994, MNRAS, 269, 928
 Stanev (2010a) Stanev, T. 2010a, Review presented at the 2010 Vulcano workshop, Vulcano, Italy (arXiv: astroph/1011.1872)
 Stanev (2010b) Stanev, T. 2010b, High Energy Cosmic Rays (Praxis Publishing Ltd, Chichester, UK, 2nd edition (First edition published in 2004.))
 Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507
 Tingay et al. (1998) Tingay, S. J. et al. 1998, AJ, 115, 960
 van Velzen & Falcke (2013) van Velzen, S. & Falcke, H., The Intriguing Life of Massive Galaxies, Proceedings of the International Astronomical Union, IAU Symposium, Volume 295, pp. 271271
 van Velzen et al. (2012) van Velzen, S., Falcke, H., Schellart, P., Schellart, N., & Kampert, K.H. 2012, A&A, 544, 18
 VéronCetty & Véron (2006) VéronCetty, M. P. & Véron, P. 2006, A&A 455, 773
 Vietri (1995) Vietri, M. 1995, ApJ, 453, 883
 Vila & Romero (2010) Vila, G. S. & Romero, G. E. 2010, MNRAS, 403, 1457
 Waxman (1995) Waxman, E 1995, Phys. Rev. Lett., 75, 386
 Waxman & Loeb (2009) Waxman, E. & Loeb, A. 2009, JCAP, issue 8, id. 026
 Whysong & Antonucci (2003) Whysong, D. & Antonucci, R. 2003, NewAR, 47, 219
 Zatsepin & Kuzmin (1966) Zatsepin, G. T. & Kuzmin, V. A. 1966, Pis’ma Zh. Eksp. Teor. Fiz., 4, 114, [English translation in JETP Lett., 4, 78 (1966)]
 Zaw et al. (2009) Zaw, I., Farrar, G. R., & Greene, J. E. 2009, ApJ, 696, 1218
 Znajek (1978) Znajek, R. L. 1978, MNRAS, 185, 833
Appendix A Electron and proton number densities
First, we consider the case that the jet is composed mainly of electrons, positrons, and protons. Then, we chose the case where the mass composition of the UHECRs consists of 90% protons and 10% iron nuclei instead of 100% protrons. We denote by the ratio of the electron to proton number densities, where the number densities are measured in a frame comoving with the jet plasma. Unless otherwise noted, should be assumed to include the positron number density as well. It is straightforward to generalize to a mixed chemical composition, including many heavy nuclei. Furthermore, both electrons and protons can have thermal and nonthermal populations before being accelerated at the shock. There may also be a substantial number of positrons from pion production and decay processes (also called secondaries).
Now, we look for the expression of the proton and electron number densities injected into the accelerating region. First, we consider the mass flow rate into the jets, which in the comoving frame is given by
(13) 
where is the restmass density of the jet, is the comoving volume of the jet, is the launching area of the jet, is the length of the cylinder along which the jet propagates before expanding freely in a conical geometry, and is the bulk velocity of the jet.
The surface area between two equatorial surfaces of a Kerr BH can be calculated as
(14) 
where the Kerr metric functions are:
(15) 
where is the coordinate radius. Next, we use normalizations to the gravitational radius, so that is the dimensionless radius. The surface area is then:
(16) 
where the factor increases from 2 to 80 as the BH spin parameter increases from 0.95 to 1. For the first equality, we use the fact that the inner disk, from where the jet is launched, has its inner and outer radii at the innermost stable orbit^{10}^{10}10Once the accretion flow reaches the innermost stable orbit, it drops out of the disk and falls directly into the BH. The expression for the radius of the innermost stable orbit is given by eq. (2.21) in Bardeen (1970). and stationary limit , respectively.
The comoving density of the jet can be expressed in terms of the ratio of the electron to proton number densities:
(17) 
For protons dominating over the electrons, , where electrons and positrons can partially occur as secondaries. This fraction of CR is from data of CRs at 1 GeV. One can get a ratio of unity assuming that the spectra go down to rest mass, which is implausible [see, e.g., Protheroe & Biermann (1996) and references therein].
Substituting Eqs. (16) and (17) for (13), we obtain the mass flow rate into the jet in the observer frame by including :
(18) 
This expression provides the proton number density, which we use to derive the electron number density:
(19) 
Now, if we consider that the UHECR particles are composed of 90% protons and 10% iron nuclei, the denominator of Eq. 19 is modified through a multiplying factor of , as we replace with . We shall use this result later for evaluating the selfabsorbed synchrotron emission of the jets (Section B).
Appendix B Selfabsorbed synchrotron emission of the jets
The spectra from compact radio sources can be explained by selfabsorbed synchrotron emission of the jets produced by electrons with a powerlaw energy distribution. In this section, we rewrite the quantities which describe the selfabsorbed synchrotron emission, derived in Rybicki & Lightman (1979), and express them under the considerations of the model presented here. Using Eqs. 1 and 2, the absorption coefficient [eq. 6.26 in Rybicki & Lightman (1979)] becomes:
(20) 
where
(21) 
To calculate the observed distance along the jet where the jet becomes selfabsorbed, we first determine the optical depth of the jet material. The averaged path of a photon through the jet has the length , which is a reasonable approximation for a jet observed at large inclination angle (e.g., Kaiser 2006). We introduce a factor in the expression of the path length to account for a small inclination angle. Thus, we can write the optical depth as
(22) 
For conical jets, the intrinsic halfopening angle is given by . With the absorption coefficient specified through Eq. (20), the optical depth can be written as
(23) 
One can define the distance along the jet where the jet becomes selfabsorbed as the distance for which . Using Eq. (23), one obtains:
(24) 
The total power radiated per unit volume per unit frequency by a nonthermal particle distribution equals:
(25) 
where [Eq. 6.36 in Rybicki & Lightman (1979)]. Using Eqs. (2) and (1), as well as the method to calculate the averaged pitch angle employed in Longair (1994), the expression of the total power becomes:
(26) 
where
(27) 
Next, we derive from Eq. 23 when . With this, the expression for the total power takes the form:
(28) 
and the emission coefficient is simply .
At low frequencies, the emitting region is opaque to synchrotron radiation and the observed intensity of radiation is proportional to the source function, while at high frequencies, the region is transparent and the observed intensity is proportional to the emission coefficient. This transition corresponds to an optical depth . Using Eqs. (20) and (26), the source function () when becomes:
(29) 
where .
To obtain the emission spectrum, one needs to solve the equation for the radiative transfer through a homogeneous medium. Because the angular sizes of the jets are small, instead of the specific intensity of the radiation, one usually measures the flux density (energy per unit time, per unit frequency interval, that passes through a surface of unit area). Thus,
(30) 
where , with the distance from the observer to the jet source and . If we integrate Eq. 30 from to , we obtain the flux density of the synchrotron emission in the case of as
(31) 
where the second term in the last squared bracket can be neglected with respect to the first term for (where the jet emission becomes selfabsorbed). Using Eqs. (24) and (29), the flux density is then:
(32) 
where . The radio flux density in Eq. (32) does not depend on the emitted frequency of the radiation since we already adopted the case of flatspectrum core sources when .
For a powerlaw synchrotron spectrum of the form from a continuous jet, the observed flux density is related to the intrinsic flux density as
(33) 
where is the Doppler factor of the jet and is the inclination angle of the jet axis with respect to the line of sight .
Appendix C Relation between the jet power and the observed radio flux density for a flatspectrum core source
In the previous section, we established the expression for the radio flux density from sources with a flatspectrum core (Eq. 32). This quantity reflects the radiative property of the jet, as the radiated energy is replaced by dissipation of the jet power (e.g., Blandford & Königl 1979). In this section, we seek the relation between the jet power and the observed radio flux density. First, we consider the jet power in the observer frame defined as
(34) 
which follows, e.g., from Falcke & Biermann (1995) [see also Vila & Romero (2010)], and for which we need to evaluate using Eq. 19. An upper limit for the electron density is specified by . So, we can substitute Eq. (19) for the expression of the observed radio flux density (Eq. 33) and find the mass flow rate into the jet . The strength of the magnetic field follows from Eqs. (3) and (4). This procedure yields the power of the jet:
(35) 
where
(36) 
where for a continuous jet emission (Eq. 33) and . We use a normalization value for the Lorentz factor of the jet, say 5, although this factor can range from to , as observational data suggest. We adopt (and then ), which means that there is, in average, one hundred electrons/positrons for at least one proton (or one heavy nucleus) in a jet that is powered by a very rapidly spinning BH () and observed at a large angle (). Flat spectrum cores are predicted for any angle to the line of sight (Blandford & Königl 1979). They are pointing close to the line of sight only if the cores dominate over the extended emission.
Using Eq. (35) for the observed radio flux density of a conical jet, we obtain:
(37) 
where the electron number density in the jet scales as , , and . Since the observed radio flux density in Eq. (37) is not dependent on the distance along the jet, the expression can be applied to microquasars as well.
Appendix D Error bars in Monte Carlo computations
Monte Carlo computations are likely to have large errors and it is important to estimate the order of magnitude of the error in an easy way. Therefore, all Monte Carlo computations should report error estimates.
If we take as a random variable for times and then we consider the value of by
(38) 
Next, if we use the central limit theorem from the probability theory, which states that the mean of a sufficiently large number of independent random variables, each with a welldefined mean and welldefined variance, will be approximately normally distributed, then:
(39) 
where is the standard deviation of and . We calculate further that is , where . Then:
(40) 
and finally,
(41) 
The Monte Carlo data is usually described like and graphically we represent the error using a bar of length with the estimated value of , in the center and we evaluate the value of in the interval
Source  D  M  /  

(Mpc)  ()  (mJy)  ( eV)  
(1)  (2)  (3)  (4)  (5)  (6)  (7) 
ARP 308  4.081.02  0.170.06  0.280.04  
CGCG 114025  5.631.33  0.230.09  3.120.04  
ESO 137G006  12.43.1  0.520.20  2.130.26  
IC 4296  12.923.23  0.540.21  0.180.02  
IC 5063  5.781.44  0.240.09  0.650.08  
NGC 0193  5.781.44  0.240.96  0.070.004  
NGC 0383  9.582.35  0.400.15  0.100.004  
NGC 1128  5.781.44  0.240.09  0.090.01  
NGC 1167  8.762.19  0.360.14  0.050.007  
NGC 1316  12.43.1  0.520.20  0.0080.0001  
NGC 1399  7.081.77  0.290.11  0.0050.0007  
NGC 2663  10.092.48  0.420.16  0.130.02  
NGC 3801  6.061.51  0.250.10  2.270.49  
NGC 3862  8.572.14  0.350.14  6.521.41  
NGC 4261  9.320.89  0.390.09  0.360.05  
NGC 4374  12.921.29  0.540.13  0.070.001  
NGC 4486  23.843.5  1  1  
NGC 4651  2.580.64  0.100.04  0.030.005  
NGC 4696  7.081.77  0.290.11  0.070.01  
NGC 5090  11.122.78  0.460.18  0.290.04  
NGC 5128  3.030.82  0.120.05  32.605.25  
NGC 5532  13.433.35  0.560.22  0.060.01  
NGC 5793  4.831.2  0.200.08  0.930.16  
NGC 7075  6.461.55  0.270.10  0.080.01  
UGC 01841  4.081.02  0.170.06  1.140.17  
UGC 02783  8.372.09  0.350.13  0.650.09  
UGC 11294  6.961.68  0.290.11  0.010.002  
VV 201  4.081.02  0.170.06  0.390.05  
WEIN 45  6.711.61  0.280.10  1.850.23 
NOTE: Col. (1) Source name; Col. (2) Distance to the source; Col. (3) BH mass; Col. (4) Core flux density at 5 GHz; Col (5) Maximum particle energy (spatial limit) for ; Col. (6) Maximum particle energy relative to that of M 87 (spatial limit); (7) UHECR flux relative to that of M 87. To calculate the errors for the quantities in columns (5), (6), and (7), we use the method of propagation of uncertainty. REFERENCES for Table 1: For Col. (2): The NASA/IPAC Extragalactic Database (NED), ned.ipac.caltech.edu (exception: NGC 5128 has unspecified error, we take it as 10%); For Col. (3): we use the method from Caramete & Biermann (2010) (exceptions: NGC 4261, NGC 4374, and NGC 4486 from Marconi & Hunt (2003) and NGC 5128 from Cappellari et al. (2009)); For Col. (4) The NED (exceptions: ARP 308, NGC 5532, and UGC 01841 have unspecified error, we take it as 10%; NGC 1167 and NGC 4261 have unspecified error, we take it as 10%; NGC 3801 and NGC 3862 have unspecified error, we take it as 10%; NGC 1167 and NGC 4261 from Nagar et al. (2001b) have unspecified error, we take it as 10%; NGC 4651 have unspecified error, we take it as 10%; UGC 02783 has unspecified error, we take it as 15%. For CGCG 114025, ESO 137G006, IC 5063, NGC 5793, NGC 7075, UGC 11294, and WEIN 45 we estimate the core radio flux from the total one using the fitting formula from Giovannini et al. (2001)).