Neutrinos and Ultra-high-energy Cosmic-ray Nuclei from Blazars

In spite of the experimental efforts to measure cosmic rays (CRs) at the highest energies, i.e., above 10¹⁸ eV, their origin is not yet clear. These ultra-high-energy CRs (UHECRs) are accelerated in extragalactic sources (Aab et al. 2017a) that are not yet resolved. The Pierre Auger Observatory has measured observables that are sensitive to chemical composition, favoring the interpretation of a mixed CR composition (Aab et al. 2017b) and motivating the investigation of the interactions of nuclei heavier than hydrogen in the sources and their transport through extragalactic space. An independent test can be obtained from the observations of secondary messengers produced during CR propagation: neutrinos and γ-rays (see, e.g., Heinze et al. 2016; Supanitsky 2016). Active galactic nuclei (AGNs) are considered as sites where the acceleration of protons and nuclei might take place. In this study we investigate the possibility that blazars, i.e., AGNs whose jet points in the direction of the observer, are the sources of the UHECRs and neutrinos, which are expected to be produced through interactions of CRs with photons in the source and its surroundings.

High-energy astrophysical neutrinos, recently detected by IceCube (Aartsen et al. 2013a, 2013b), are the smoking-gun signature of hadronic interactions in astrophysical environments. However, the origin of these neutrinos remains unknown. Blazars, as the dominating extragalactic gamma-ray sources (Ajello et al. 2015), have long been proposed as promising neutrino emitters if the observed gamma rays are generated in hadronic interactions (Stecker et al. 1991). Indeed, if all blazars belonged to this type, their energy budget alone would be sufficient to power the astrophysical neutrino flux observed by IceCube. However, Aartsen et al. (2017) have limited the fraction of this contribution to the neutrino flux from blazars to 7%–27%, after a study on the spatial correlation between the directions of the IceCube neutrinos and the positions of the sampled blazars from the Fermi-2LAC catalog (Fermi-LAT 2011). Correlation studies between IceCube neutrinos and GeV and X-ray flares from point sources, including blazars, have been performed as well, but no statistically significant result has been found (Peng & Wang 2017), except from some hints (Kadler et al. 2016; Padovani et al. 2016; Resconi et al. 2017). A few significant GeV blazar flares, such as 3C 279 and PKS B1424-418, have been found to be inconsistent with a purely protonic origin from theoretical modelings (Gao et al. 2017; Petropoulou et al. 2017). Therefore, the important remaining questions are the following: what fraction of the IceCube neutrino flux can be attributed to blazars, given the possibility that the blazar jet may contain a mixture of leptons, protons, and nuclei; what the characteristic neutrino energies as a function of those ingredients are; and whether blazars can be the sources of the UHECRs in spite of possible neutrino constraints. This article cannot address all of these questions, but will rather focus on the neutrino and CR production efficiencies of individual sources.

The neutrino spectrum and production efficiencies are highly dependent on the spectral energy distribution (SED) of the seed photons (Murase et al. 2014). Therefore, a meaningful prediction of the diffuse neutrino spectrum from blazars must incorporate the differences in the SEDs among the blazar family. Based on the work of Fossati et al. (1998), the description of SEDs from observed blazars can be reasonably approximated by a sequence as a function of the gamma-ray luminosity L_γ: the brighter the source is in gamma rays, the "redder" the SED becomes (i.e., with the spectral peaks at lower energies).¹ Based on this prescription, we study the neutrino and CR emissivity for a number of representative cases out of the sequence: from low-luminosity, high-frequency-peaked BL Lacs (HBL) to high-luminosity, low-frequency-peaked BL Lacs (LBL) and flat-spectrum radio quasars (FSRQs). As a function of both the luminosity and jet composition, we reveal which type of blazar is the best candidate for neutrino and UHECR sources.

Our blazar model is an extended version of Murase et al. (2014; see also Anchordoqui et al. 2008, for an earlier approach), where the ejected neutrino spectra have been calculated for a sample of blazars of different luminosities considering a pure proton injection composition. In that model, blazar emission cannot simultaneously explain the observed diffuse PeV neutrino flux and the UHECR flux, since a different CR loading would be required to explain either one; furthermore, the diffuse neutrino spectra obtained were found to be too high compared to IceCube observations at multi-PeV energies. In comparison with that work, in the present model we include the injection of nuclei heavier than protons including all relevant radiation processes, such as photonuclear disintegration. The latter process refers to all interactions exhibited by nuclei below the pion production threshold, which lead to the disintegration of the nucleus into lighter isotopes and one or multiple lighter fragments (see, e.g., Boncioli et al. 2017). At higher energies these processes are accompanied by inelastic interactions of target photons with individual nucleons, resulting in photohadronic production of mesons that can subsequently decay into neutrinos. Given the enhanced star-forming rates in the centers of the host galaxies, it is natural to expect high abundance of these heavy elements, as is observed from the emission lines from the accretion disk (see, e.g., Fabian et al. 2000). In this work we assume that these elements can be transported to the jet from either the accretion disk or other locations, such as in a star-in-jet scenario (Barkov et al. 2010). We demonstrate that for heavier injection compositions the emitted neutrino spectrum has a lower cutoff and overall normalization, which may help solve both aforementioned conflicts. We also include a more refined discussion of possible CR escape mechanisms from the jet zone (diffusion versus advection) and show their impact on the production of secondary particles in the broad-line region (BLR) and dust torus (DT) in FSRQs. For the sake of simplicity, we focus on pure injection compositions into the jet, ranging from protons to iron.

The paper is structured as follows: In Section 2, we introduce the blazar jet model and the simulation methods. In Section 3, we analyze the effect of the source parameters on the nuclear cascade and identify two qualitatively different nuclear disintegration regimes in the jet. Neutrino production is discussed in Section 4, and the CR escape mechanisms are discussed in Section 5, where the emitted CR spectra are also presented. In Section 6, we introduce a model for FSRQs that includes the BLR and DT as additional radiation zones. In Section 7 we apply our model to a sample of sources across the blazar sequence and calculate their neutrino and UHECR emissivity. Our conclusions are presented in Section 8.

For the simulation of blazars with the injection of nuclei we use similar techniques to those in Boncioli et al. (2017) and Biehl et al. (2017), where details can be found in Biehl et al. (2017). Our jet model closely follows Murase et al. (2014) as far as geometric setup and SEDs are concerned. Until Section 6 we will consider only radiation processes in the jet, which includes CR interactions with the nonthermal SED, neutrino production, and CR escape. In Section 6, the model will be extended to include external radiation typically associated with FSRQs.

2.1. Model Geometry

We assume that CRs are accelerated in a plasma blob of radius² ${r}_{\mathrm{blob}}^{{\prime} }\simeq 3\,\times \,{10}^{16}\,\mathrm{cm}$ unless noted otherwise, where the size can be estimated from a flaring time ${t}_{\mathrm{flare}}^{{\prime} }={10}^{6}\,{\rm{s}}$ such that ${r}_{\mathrm{blob}}^{{\prime} }\simeq {\rm{\Gamma }}\,c\,{t}_{\mathrm{flare}}$. The blob is spherical in its rest frame with volume $V^{\prime} =(4\pi /3){r}_{\mathrm{blob}}^{{\prime} 3}$ and travels along the jet with Doppler factor Γ ≃ 10 toward the observer.³ The jet is assumed to have an opening angle of θ ∼ Γ⁻¹ ∼ 0.1, so that the typical dissipation radius can be estimated from a simple geometry argument to be ${r}_{\mathrm{diss}}\approx {\rm{\Gamma }}{r}_{\mathrm{blob}}^{{\prime} }=3\times {10}^{17}\,\mathrm{cm}$.

The dissipation radius will not come into play until Section 6, when we consider external radiation fields characteristic of some FSRQs, namely, thermal radiation from the accretion disk and the DT, and molecular emission from the BLR. In sources where the BLR radius is larger than the jet dissipation radius, r_BLR > r_diss, the radiation zone of the jet will be situated within the BLR (see middle panel of Figure 1). In this case, we will assume that the component originating from these external fields is boosted into the jet blob. In addition, the CRs that escape the jet will have to propagate through the external fields before escaping into extragalactic space from the blazar environment. In cases where the blob dissipation occurs outside the BLR (r_diss > r_BLR; see top panel of Figure 1), CRs only experience the nonthermal radiation inside the blob, and those that escape the jet are effectively emitted by the source.

Zoom In Zoom Out Reset image size

We assume the radius of the BLR (Ghisellini & Tavecchio 2008) and the DT (Cleary et al. 2007; Hönig & Beckert 2007; Ghisellini & Tavecchio 2008; Kishimoto et al. 2011a, 2011b) to scale with the square root of the accretion disk luminosity L_disk:

$$\begin{eqnarray}&&{r}_{\mathrm{BLR}}={10}^{17}{\left(\displaystyle \frac{{L}_{\mathrm{disk}}}{{10}^{45}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/2}\,\mathrm{cm},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{r}_{\mathrm{DT}}=2.5\times {10}^{18}{\left(\displaystyle \frac{{L}_{\mathrm{disk}}}{{10}^{45}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/2}\,\mathrm{cm}.\end{eqnarray} \tag{ 2 }$$

Following the previous discussion, Equation (1) implies that CR interactions with external fields are relevant only for FSRQs with ${L}_{\mathrm{disk}}\gt 9\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. These sources will be denoted high-luminosity FSRQs (HL-FSRQs), as indicated in Figure 1. For other FSRQs, as well as for BL Lacs, a single plasma blob scenario containing only nonthermal jet emission characterizes the interaction environment. In the next sections and until Section 6 we will focus on this simplified one-zone scenario. A further description of the three-zone FSRQ model will be given in Section 6.

The bolometric luminosity of the disk, L_disk, is roughly proportional to that of the jet, L_jet, based on the phenomenological relationship from Inoue & Totani (2009), Inoue et al. (2010), Abazajian et al. (2011), and Lusso et al. (2010). In this work, the blazars are represented by L_γ, the integrated luminosity of the SED above 100 MeV. For the relationships between L_disk, L_jet, and L_γ, we follow the tabulated values in Murase et al. (2014), where for FSRQs L_disk can be empirically related to L_γ as ${L}_{\gamma }\,[\mathrm{erg}\,{{\rm{s}}}^{-1}]\approx {10}^{17}\,{({L}_{\mathrm{disk}}[\mathrm{erg}{{\rm{s}}}^{-1}])}^{0.683}$.

2.2. Particle Interactions

The transport of accelerated nuclei in the source is modeled as a coupled partial differential equation (PDE) system including the colored species in Figure 2. Details may be found in Biehl et al. (2017), where the same method is used for solving the transport equations based on the NeuCosmA code (Baerwald et al. 2012). The PDE reads for the isotope i

$$\begin{eqnarray}&&\displaystyle \frac{\partial {N}_{i}^{{\prime} }}{\partial t^{\prime} }=\displaystyle \frac{\partial }{\partial E^{\prime} }(-b^{\prime} (E^{\prime} ){N}_{i}^{{\prime} }(E^{\prime} ))-\displaystyle \frac{{N}_{i}^{{\prime} }(E^{\prime} )}{{t}_{\mathrm{esc}}^{{\prime} }}+{\tilde{Q}}_{{ji}}^{{\prime} }(E),\end{eqnarray} \tag{ 3 }$$

where $b^{\prime} (E)=E^{\prime} t{{\prime} }_{\mathrm{loss}}^{-1}$ (with the energy-loss rate $t{{\prime} }_{\mathrm{loss}}^{-1}$) and $t{{\prime} }_{\mathrm{esc}}^{-1}$ is the escape rate. The PDE system is to be solved for the differential particle densities ${N}_{i}^{{\prime} }\,[{\mathrm{GeV}}^{-1}\,{\mathrm{cm}}^{-3}]$ in the blob frame. The coupling of the PDEs arises because of the injection term

$$\begin{eqnarray}&&{\tilde{Q}}_{{ji}}^{{\prime} }(E)={Q}_{i}^{{\prime} }(E)+{Q}_{j\to i}^{{\prime} }(E),\end{eqnarray} \tag{ 4 }$$

which allows for injection from an acceleration zone ${Q}_{i}^{{\prime} }$, as well as for injection from other species j with the term ${Q}_{j\to i}^{{\prime} }$, such as from photodisintegration or β^± decays. In this work, we consider pure injection compositions, which means that there is only one injection term ${Q}_{i}^{{\prime} }$ for the injected species i (such as protons, helium-4, or iron-56):

$$\begin{eqnarray}&&{Q}_{i}^{{\prime} }\propto {({E}_{i}^{{\prime} })}^{-2}\times \exp (-{({E}_{i}^{{\prime} }/{E}_{i,\max }^{{\prime} })}^{2}).\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

The maximum energy ${E}_{i,\max }^{{\prime} }$ is determined by equating the acceleration rate,

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}^{{\prime} -1}=\displaystyle \frac{\eta \,{c}^{2}\,Z\,e\,B^{\prime} }{E^{\prime} },\end{eqnarray} \tag{ 6 }$$

with the sum of synchrotron loss, photodisintegration, and photomeson production rates and the light-crossing time of the blob. Here, we choose $\eta =1$ for the acceleration efficiency,⁴ and B' is the magnetic field strength in the blob. Unless noted otherwise, the normalization of the nuclear injection spectrum is determined by a fixed ratio compared to the gamma-ray luminosity, namely, the CR loading factor f_CR = 10.

In order to set up the PDE system, automated techniques are used: Starting with a database of nuclear isotopes⁵ (white boxes in Figure 2), fast spontaneous nucleon or α-particle emitters (with a very conservative selection threshold ${\tau }_{0}\lt {10}^{-10}\,{\rm{s}}$) are integrated out from the beginning, which means that they are not explicitly considered but instead replaced by their daughters. Then, starting with the heaviest possible injection isotope iron-56, a recursive algorithm determines all possible photodisintegration and decay (and also mixed) paths, i.e., the isotopes that will be dominantly populated. Such isotopes are marked in blue/green in the figure, and for these the coupled differential equation system is automatically set up and solved. The colored isotopes marked by crosses are potentially interesting (fast enough) beta emitters, which may contribute to neutrino production in the jet—which we take into account. Included radiation processes for the nuclei are β^± decays and spontaneous nucleon (and light nuclei) emissions,⁵ photodisintegration with TALYS (Koning et al. 2007; A ≥ 12) and CRPropa2 (Kampert et al. 2013; A < 12), photomeson production (Mucke et al. 2000; Hümmer et al. 2010; Biehl et al. 2017), photopair production (Blumenthal 1970), synchrotron cooling, and particle escape (depending on the CR escape assumption).⁶

As opposed to Boncioli et al. (2017) and Biehl et al. (2017), we use a time-dependent computation for the fluences, i.e., we explicitly integrate the neutrino and CR fluxes instead of using the steady-state approach. The solver integrates the blob system up to a few times ${t}_{\mathrm{flare}}^{{\prime} }$, while the injection of nuclei is switched off at ${t}_{\mathrm{flare}}^{{\prime} }$ (and the electrons are assumed to maintain the radiation field). While this approach is computationally more expensive than a steady-state solution and leads to similar results in the jet zone, it is more convenient for the modeling of the CR transport between the jet, BLR, and DT zones (see Section 6).

The magnetic field inside the blob plays a role in the cooling of charged CRs, mesons, and muons through synchrotron radiation, in the escape rate of charged CRs if a diffusion mechanism is assumed (see Section 5), and in the acceleration rate Equation (6). The magnetic field is obtained through the relation

$$\begin{eqnarray}&&{u}_{B}^{{\prime} }=\displaystyle \frac{{\epsilon }_{B}\,{L}_{\gamma }^{{\prime} }\,{t}_{\mathrm{flare}}^{{\prime} }}{V^{\prime} }\,\propto {r}_{\mathrm{blob}}^{{\prime} -2},\end{eqnarray} \tag{ 7 }$$

where ${u}_{B}^{{\prime} }={B}^{{{\prime} }^{2}}/8\pi$ is the magnetic energy density and _B is the magnetic loading. We fix the value of _B for all sources in this study, with the brightest jet (introduced in Section 7) having a magnetic field strength of 5 G, following Murase et al. (2014). Such a constant value of the magnetic loading for all the sources, in this case ${\epsilon }_{B}=6.3\times {10}^{-3}$, leads to the magnetic field strength being proportional to ${L}_{\gamma }^{1/2}$ and inversely proportional to the size of the blob. In Appendix B we explore a different choice of magnetic field values, assuming that the SEDs are generated via leptonic SSC models.

The blazar emission power in the source frame is represented by L_γ, the isotropic-equivalent integrated luminosity of the SED above 100 MeV. From the blazar sequence modeling, the SED and hence the comoving frame photon spectrum are constructed phenomenologically as a function of L_γ (Murase et al. 2014). The γ-ray emission power in the comoving frame of the blob, ${L}_{\gamma }^{{\prime} }$, is connected to the black hole frame quantity via ${L}_{\gamma }={{\rm{\Gamma }}}^{4}{L}_{\gamma }^{{\prime} }$, where a factor Γ² comes from boosting the luminosity itself and a factor Γ² from the solid angle boost. Due to the fact that all the nuclear interactions (photomeson production, photopair production, and photodisintegration) require a certain threshold energy, the relevant target photons are mainly the soft photons located in the first hump and the middle range of the SED.

We emphasize that the blazar SEDs are taken from observations, where we assumed that they emerge self-consistently from all participating radiation processes. However, it remains important to check whether and under what conditions this assumption can hold. The feedback of hadronic photons to the SED is nontrivial and can be difficult to compute, even in the simplest proton-only case. As shown in Figure 4 of Gao et al. (2017), the feedback may show up in the X-ray range—the saddle of the two humps. In this work, we aim to study which object classes can efficiently produce neutrinos or CRs. Therefore, we leave the fully self-consistent construction of SEDs, including nuclear interactions, to future works.

In this section, we discuss two distinctive scenarios where different physics takes place for nuclei in the jet: a low-luminosity prototype that is optically thin to photonuclear interactions, representative of what we call the "nuclear survival" regime, and a high-luminosity jet that is optically thick to photonuclear interactions at the maximum energy, representing the "nuclear cascade" regime.

3.1. Nuclear Survival Regime

The example in Figure 3 illustrates the nuclear survival case for iron-56 injection, where ${L}_{\gamma }={10}^{44.6}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and the SED is shown in the top left panel. From the interaction rates in the top right panel, one can clearly see that the source is optically thin to disintegration at the maximum energy ${E}_{\max }^{{\prime} }\sim {10}^{9}\,\mathrm{GeV}$ (which is determined by the light-crossing time)—although disintegration, which strongly increases with energy, is the leading photonuclear process in the radiation zone. The initial nucleus interacts so rarely that no significant cascading into lower-mass nuclei can be observed (bottom left panel), and thus the primary nucleus survives. Since the maximum energy is determined by the light-crossing time, the primary spectrum in the source corresponds to the injection spectrum extending to the maximum injection energy and without spectral breaks (see bottom right panel). The small offset between the injection spectrum and the in-source density of iron-56 is due to adiabatic cooling (see Section 5).

Zoom In Zoom Out Reset image size

3.2. Nuclear Cascade Regime

The nuclear cascade regime is illustrated in Figure 4 by a high-luminosity jet with ${L}_{\gamma }^{{\prime} }={10}^{48.8}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ for iron-56 injection, where the corresponding SED is shown in the top left panel. Note that only interactions with the jet radiation are considered, without additional components from BLR or DT. The maximum energy is limited by photodisintegration, and the source is optically thick (τ ∼ 100) to this process (see top right panel), resulting in efficient disintegration of the primary nucleus into lower-mass nuclei (see bottom left panel). Significant fractions of energy are transferred into densities of neighboring nuclei and light fragments, such as nucleons and α particles. As opposed to high target density cases in gamma-ray bursts (GRBs; Biehl et al. 2017), the disintegration rate strongly decreases with energy. The in-source density at the highest energies (bottom right panel) is dominated by secondary nuclei from photodisintegration reprocessing the injected isotopes into a lighter mass composition. This is characteristic for the high-energy CRs in the nuclear cascade regime, whereas the primary spectrum at lower energies is not notably modified.

Zoom In Zoom Out Reset image size

Regarding the dominant radiation process (top right panel), we find that disintegration typically dominates over photomeson production for the primaries. It has been demonstrated in Anchordoqui et al. (2008), as reproduced in this work, that photomeson production can dominate over photodisintegration for special SEDs in certain energy ranges. We observe this behavior in certain rare SEDs at the most extreme luminosities, which we however do not consider to be representative. Note that in both Anchordoqui et al. (2008) and this work a superposition model for photomeson production has been used, whereas the realistic photomeson cross sections will be somewhat smaller and need further study; see discussions in Boncioli et al. (2017) and Biehl et al. (2017).

3.3. Parameter Space Study

By scanning the parameter space of ${r}_{\mathrm{blob}}^{{\prime} }$ and L_γ, we find that the two nuclear regimes, defined by the optical thickness to photonuclear interactions at the highest energy, can occur over a wide range of parameters. We show this result in Figure 5 for three different injection isotopes: protons, helium-4, and iron-56. Note that in this figure the shape of the SED is assumed to follow the blazar sequence as a function of L_γ as we will introduce later in Section 7 (but without the external radiation fields), whereas ${r}_{\mathrm{blob}}^{{\prime} }$ only affects the photon density normalization. From Figure 5 we see that by varying L_γ along a fixed line of ${r}_{\mathrm{blob}}^{{\prime} }=3.0\,\times {10}^{16}\,\mathrm{cm}$, corresponding to our prototypes (points A and B), one can effectively control the level of disintegration and, as we will see later, the balance between a more efficient neutrino source and a more efficient UHECR source. A similar effect can be achieved by changing ${r}_{\mathrm{blob}}^{{\prime} }$ for fixed L_γ. As we will later associate L_γ with certain object classes of the blazar sequence, it is important to remember that the splitting point between nuclear survival and nuclear cascade depends somewhat on the blob size. Comparing the middle and right panels of Figure 5, we see that iron-56 can populate a nuclear cascade at lower jet luminosities compared to helium-4, since iron-56 has a higher photodisintegration rate.

Zoom In Zoom Out Reset image size

Figure 5 contains contours of the maximum energy, and it can be clearly seen that larger accelerators tend to allow for higher acceleration energies. For the case of protons (left panel), these are determined by the available time for acceleration in the case of lowest luminosities, while at the highest luminosities, photomeson interactions are responsible for limiting the maximum energy. For the case of nuclei, photodisintegration limits the maximum energy at the highest luminosities and determines the transition to the nuclear cascade regime. For low luminosities the maximum energy does not depend on the size of the blob, neither for protons nor for the nuclei. That is because the dynamical timescale is directly proportional to the blob size, while the acceleration rate is proportional to the magnetic field, which scales inversely with the blob size. Above a certain luminosity the shape of the contours is more complicated owing to the interplay between the acceleration and photodisintegration rates. The middle and right panels of Figure 5 reveal that the maximum energy is obtained in a narrow region where the transition between the two regimes occurs. In this region the acceleration is strong enough to reach high energies, while energy losses from photodisintegration are not too restrictive. The small features in the contours originate from the variation of SEDs along the blazar sequence.

The ejected neutrino fluence for both the nuclear survival and cascade regimes is shown in Figure 6, where the contributions from interactions of the primary, the secondary nuclei, the nucleons, and β decay are shown separately (where the secondary nuclei and nucleons are produced in the nuclear cascade). In the nuclear survival regime neutrinos are produced off the injection isotope, with very small contributions from secondary nuclei and nucleons. On the other hand, for sources with a nuclear cascade, most of the neutrinos at high energies are produced off secondary isotopes and nucleons.

Zoom In Zoom Out Reset image size

Contrary to the GRB case, discussed in Biehl et al. (2017), where three different regimes were identified,⁷ we only find two distinctive scenarios for blazars. The main difference is that for GRBs the interaction rate is about constant beyond the break in the SED, whereas it strongly increases for blazars in the relevant energy range (see Figures 3 and 4, top right panels). While the nuclei in GRBs continue to disintegrate efficiently along the cascade, in blazars disintegration ceases much faster because the disintegration rate quickly drops with decreasing energy. We therefore do not find a case where most energy is dumped into nucleons, which would be comparable to the optically thick regime of the GRB models (Biehl et al. 2017). Rather, neutrino production in the nuclear cascade regime of blazars is typically dominated by the secondary nuclei produced in the nuclear cascade. This again highlights the importance of better photomeson production models off nuclei pointed out in Biehl et al. (2017), as our results are based on a simple superposition model where the total interaction and inclusive pion production cross sections scale with A.

The dependence of the neutrino spectrum on the injection composition is shown in Figure 7. In the nuclear survival case, the maximum energy per nucleon scales ∝Z/A ≃ 1/2 for most isotopes (except nucleons), which means that we observe a weak dependence on the composition. In the nuclear cascade regime, the neutrino spectrum basically follows the maximum energy per nucleon as well, which, however, depends in a nontrivial way on the competition of the disintegration, photomeson, and acceleration timescales. In fact, from Figure 5 one can see that the maximum primary energy E_max for point B varies by only a factor of a few with injection composition, which means that the neutrino energy roughly scales ∝E_max/A following the energy/nucleon. Consequently, the neutrino spectrum extends to higher energies for lighter injection compositions. We therefore find that loading (high-luminosity) blazars with nuclei helps reduce the tension with the nonobservation of blazar neutrinos by IceCube in terms of both spectrum and normalization.

Zoom In Zoom Out Reset image size

We have previously discussed CR densities inside the source, which are relevant for neutrino production in the jet. In order to address CR injection into the intergalactic space, one has to specify how CRs escape from the source.

The escape of particles from the source environment on the macroscopic scale is typically modeled via diffusion and advection processes. In order to fully address this effect, the detailed profile of the plasma, such as the geometry, plasma wave, and magnetic field configurations, is required to compute the diffusion and convection coefficients on a spatially resolved grid (Chen et al. 2016). However, in one-zone calculations, a global energy-independent escape rate is frequently assumed (Diltz et al. 2015; Petropoulou et al. 2015; Gao et al. 2017), as a fraction of the free-streaming rate. In this work, while keeping the simplicity and efficiency of the one-zone calculation, we use a more refined approach: we distinguish between escape by diffusion and escape by advection, two scenarios that correspond respectively to the most conservative and most aggressive escape assumptions. We expect that any realistic escape assumption will be in between these two scenarios. In either case, the maximum particle energy is computed self-consistently from equating acceleration, escape, and energy-loss rates.

Escape by diffusion/direct escape. Escape by diffusion of charged CRs can be described by an effective escape rate ${t}_{\mathrm{esc}}^{{\prime} -1}=D^{\prime} \,{t}_{\mathrm{flare}}^{{\prime} -2}/{c}^{2}$ (from ${r}_{\mathrm{blob}}^{{\prime} }=\sqrt{D^{\prime} \,{t}_{\mathrm{esc}}^{{\prime} }}$), where D' is the (spatial) diffusion coefficient.⁸ For Bohm diffusion, which may be the most extreme assumption here, $D^{\prime} \simeq {{cR}}_{L}^{{\prime} }$, where ${R}_{L}^{{\prime} }\propto E^{\prime}$ is the Larmor radius, which leads to hard escape spectra. If the Larmor radius can reach the size of the blob (implying that the maximum energy is determined by the flaring timescale), all particles that do not interact will escape. Note that this scenario conservatively implies magnetic confinement of charged nuclei, while neutrons can escape over the free-streaming timescale. We therefore include adiabatic cooling of the charged species (which is expected if the blob expands), allowing the particle densities to decay after ${t}_{\mathrm{flare}}^{{\prime} }$.

In the top panels of Figure 8 we show the rates for the diffusion-dominated escape as solid curves. In the nuclear survival case, the maximum energy is limited by the light-crossing time, which implies that the Larmor radius at the highest energy can reach the size of the blob; however, the escape rate cannot exceed the free-streaming escape rate. Up to the maximum energy, the escape rate grows with the Larmor radius ∼E; it therefore acts as a "high-pass filter" on the ejected CRs and results in hard emission spectral indices (see solid curve in the left panel in the second row). In the nuclear cascade case (see top right panel), the maximum energy is limited by photodisintegration, which means that the Larmor radius cannot reach the size of the blob at the maximum energy, and only a fraction of the charged particles can escape. Neutrons, however, which are abundantly produced in the nuclear cascade, can freely escape. The corresponding escape spectra are shown as solid curves in the right panel in the second row.

Zoom In Zoom Out Reset image size

A similar escape mechanism has been assumed by Baerwald et al. (2013), Globus et al. (2015), and Zhang et al. (2017) for GRBs and tidal disruption events; for example, the assumption that only a fraction close to the border of the region proportional to the Larmor radius can escape ("direct escape") leads mathematically to the same results. Another motivation for the preference of harder spectra comes from a recent result by the Auger Observatory, which suggests that source spectra are preferentially hard (Aab et al. 2017b) when fitting a generic class of UHECR accelerators.

Escape by advection. CRs transported along the jet effectively escape if the magnetic fields and densities decay quickly enough and cooling is inefficient. We assume free-streaming escape as the most aggressive escape assumption in this case. The corresponding escape rates are shown as dashed curves in the top panels of Figure 8.

Note that our advective escape assumption is similar to Murase et al. (2014) in the nuclear survival case, where the effective escape has been described by $t{{\prime} }_{\mathrm{esc}}^{-1}={f}_{\mathrm{esc}}\,t{{\prime} }_{\mathrm{flare}}^{-1}$, with the escape fraction ${f}_{\mathrm{esc}}=(1-\min [1,{t}_{\mathrm{flare}}^{{\prime} }/{t}_{\mathrm{cool}}^{{\prime} }])$ and ${t}_{\mathrm{cool}}^{{\prime} }$ the dominant cooling timescale. However, compared to Murase et al. (2014), the effect of photohadronic cooling is already included in our self-consistent nuclear cascade computation (the in-source spectra are suppressed if, e.g., photodisintegration dominates), which is why we do not take into account an additional cooling-driven suppression (which is relevant in the nuclear cascade case). Note that it is a prerequisite for advective escape that the CRs must not cool over a timescale that is much longer than ${t}_{\mathrm{flare}}^{{\prime} }$ if the interactions in the BLR are to be taken into account. In particular, any adiabatic expansion of the blob would suppress this escape component, which may be evaded if the jet is collimated. Consequently, we do not include adiabatic cooling in this scenario, and after ${t}_{\mathrm{flare}}^{{\prime} }$ densities in the source simply decay by escape and interactions.

The second and third rows in Figure 8 compare the effect of the two escape hypotheses on the ejected CR spectra. The total escape spectra in the advective case resemble the E⁻² injection spectrum over a wide range of energies since the escape rate is flat in energy, i.e., they are generally softer than in the diffusive case. Note that in the nuclear survival case (left panels), near the spectral cutoff, the fluence is higher for the advective escape assumption simply because of the absence of an adiabatic cooling term. The diffusion assumption, on the other hand, hardens the spectra of charged particles by one power (second row). Since neutrons are not magnetically confined, their spectra are identical for both escape assumptions. As mentioned previously, in reality we expect the escape spectrum to lie in between these two extreme cases.

The ejected mass composition in the advection scenario is orthogonal to what is expected from UHECR observations (see bottom right panel). Since in this study we are mostly interested in UHECRs, we conclude that in the nuclear cascade case a diffusive escape mechanism is potentially more effective, at the highest energies, in producing the right trend of the ejection composition as a function of energy within the source itself, whereas in the nuclear survival case, it allows for hard ejection spectra compatible with the Auger fit (Aab et al. 2017b).

For FSRQs, radiation from outside the jet is observed as well, such as the thermal radiation from the DT and the line emissions from the so-called BLR. In our model, only HL-FSRQs may have BLR and DT regions large enough to physically contain the blob (see middle panel of Figure 1). The photons from the accretion disk can be back-scattered in the BLR and have physical importance for interactions in the blob. At the same time, messengers produced in the blob may interact in the BLR and the DT, producing neutrinos and reducing the total CR output.

For the parameterization of the additional zones, we follow the model from Murase et al. (2014); see our Figure 1 (bottom panel). Note that the classification of FSRQs into high- and low-luminosity types is performed in order to identify sources that require additional radiation zones that impact secondary production. The classes of FSRQs with lower L_γ (i.e., blazars where broad lines are present but have negligible impact on secondary production) are treated with the jet model only, equal to BL Lacs with similar L_γ. In greater detail, the additional ingredients for the HL-FSRQs are as follows:

• 1.  
  External radiation from the accretion disk is reprocessed and isotropized in the BLR and DT and boosted into the blob, where it contributes to the target photon spectrum in addition to the nonthermal radiation (as illustrated in the middle panel of Figure 1).
• 2.  
  CRs escaping from the jet zone are boosted into the black hole frame and reprocessed in the external photon fields as they travel through the BLR, where again they may interact or escape. Those that escape propagate through the DT region, where they encounter thermal radiation from the warm dust. In both zones, the free-streaming timescale, which gives the escape rate of all CR species, corresponds to the size of the region.
• 3.  
  In the BLR and DT, additional photomeson production may occur. The neutrinos produced contribute to the total neutrino fluence of the source.

In order to illustrate the effect of the external radiation fields, we use as an example an FSRQ with luminosity ${L}_{\gamma }={10}^{48.8}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, with the same jet properties as introduced in Section 2.2 (see Figure 4). While the same magnetic field strength is considered in the jet, we assume the magnetic field strength in the BLR and DT to be much smaller, i.e., ${B}^{{\prime} }\ll 1\,\mathrm{mG}$, a valid assumption if these regions consist of gas and clouds analogous to supernova remnants. Consequently, we neglect synchrotron emission and magnetic confinement in the BLR and DT.

We illustrate the target photon densities in the three zones in Figure 9, each in the respective reference frame. Unlike BL Lacs and lower-luminosity FSRQs, the jet zone receives additional target photons from the external radiation produced in the BLR and the DT. Starting with the accretion disk radiation, we assume that a fraction τ_sc is isotropized in the BLR through Thomson scattering and that this component is Lorentz-boosted into the blob. This fraction is given by the optical thickness of the BLR to Thomson scattering, assumed to be τ_sc = 0.01 (Blandford & Levinson 1995; Tavecchio et al. 2012). In the middle panel of Figure 9 we show in green the density spectrum of photons from the disk that are scattered in the BLR, based on Elvis et al. (1994). The total energy density of scattered photons u_disk corresponds to the integral of this spectrum and is normalized according to the total disk luminosity using

$$\begin{eqnarray}&&{u}_{\mathrm{disk}}=\displaystyle \frac{{\tau }_{\mathrm{sc}}}{4\pi \,{r}_{\mathrm{BLR}}^{2}\,c}\,{L}_{\mathrm{disk}},\end{eqnarray} \tag{ 8 }$$

where L_disk is the disk luminosity, given by the empirical relation discussed in Section 2. Like all external components, this radiation is assumed to be isotropic in the BLR. In the jet rest frame, external photons appear boosted owing to the relativistic motion of the blob, thus contributing to the target photon spectrum for CR interactions in the jet zone. The photon energy in the jet rest frame is then given by ${E}_{\gamma }^{{\prime} }={\rm{\Gamma }}{E}_{\gamma }$, and the energy density receives a factor Γ² compared to the BLR (see left panel of Figure 9). We further consider the H i and He ii Lyα lines emitted by the BLR gas molecules, which are excited by the disk radiation (peaks in Figure 9). In our simplified approach, the photon density of the hydrogen line in the BLR is given by

$$\begin{eqnarray}&&{u}_{{\rm{H}}{\rm{I}}}=\displaystyle \frac{{f}_{\mathrm{cov},{\rm{H}}}}{4\pi \,{r}_{\mathrm{BLR}}^{2}\,c}\,{L}_{\mathrm{disk}},\end{eqnarray} \tag{ 9 }$$

where f_cov,H is the hydrogen covering factor of the BLR, assumed to be 0.1 (Liu & Bai 2006; Ghisellini & Tavecchio 2008), which represents the fraction of the disk luminosity re-emitted as the hydrogen broad line. The helium line density is calculated similarly but using a covering factor twice as small, so that u_{H i} = 2u_{He ii} (Tavecchio & Ghisellini 2008). Finally, the thermal photon density from the DT is shown as blue curves in Figure 9. This is a blackbody distribution of temperature 500 K (in accordance with typical DT temperatures of 100–1000 K; Cleary et al. 2007; Malmrose et al. 2011). Besides being present in the DT zone itself, this component also contributes to the target photon spectrum in the BLR as well as the jet, since the DT is assumed to surround these zones (see bottom panel of Figure 1). The total IR photon density is given in the black hole frame by

$$\begin{eqnarray}&&{u}_{\mathrm{IR}}=\displaystyle \frac{{f}_{\mathrm{cov},\mathrm{DT}}}{4\pi \,{r}_{\mathrm{DT}}^{2}\,c}\,{L}_{\mathrm{disk}},\end{eqnarray} \tag{ 10 }$$

where ${f}_{\mathrm{cov},\mathrm{DT}}=0.5$ (Murase et al. 2014) is the value assumed for the DT covering factor.

Zoom In Zoom Out Reset image size

Note that in this model the photon density of the external fields (disk, broad lines, and IR radiation from the DT) does not vary for sources with different disk luminosity or BLR size, since

$$\begin{eqnarray}&&{u}_{\mathrm{disk};\mathrm{BL};\mathrm{IR}}\,\propto \displaystyle \frac{{L}_{\mathrm{disk}}}{{r}_{\mathrm{BLR}}^{2}}\simeq \mathrm{const},\end{eqnarray} \tag{ 11 }$$

because of the scaling expressed in Equation (1).

The impact of external radiation fields on the interaction rates is shown in Figure 10. While photodisintegration (of iron-56 in this case) increases in the jet zone and the source becomes optically thick to photodisintegration above PeV energies, protons are not significantly affected, because the source is still optically thin to photomeson production. The nuclei that escape from the jet zone into the BLR and DT then undergo additional disintegration at the highest energies, reducing the efficiency of HL-FSRQs as UHECR sources. In cases where the BLR radiation is present but the DT is too small, the additional transport through the BLR would not significantly affect the UHECRs (dashed curves in the middle panel).

Zoom In Zoom Out Reset image size

We couple the BLR and DT zones in the numerical calculation by re-injecting the CRs that escape from the jet into the BLR, and then subsequently into the DT (see Figure 1 and Section 2). The solver starts by integrating the blob system up to $\sim 50\,{t}_{\mathrm{flare}}^{{\prime} }$, while the injection of nuclei is switched off at ${t}_{\mathrm{flare}}^{{\prime} }$ (and the electrons are assumed to maintain the radiation field). The CRs that escaped from the blob during the integration time are injected into a new set of equations describing the radiation environment of the BLR. There we solve an initial value problem in the black hole frame, starting with the CR densities ${N}_{i}(E,t={t}_{\mathrm{flare}})\ne 0$ ejected from the previous zone, using a homogeneous set of equations without primary nuclei injection. The DT region is treated similarly, starting with the BLR ejecta. The CR densities from the jet have to be boosted into the black hole frame and rescaled to the new volume of the region. Since no further injection takes place in the outer regions, the CR densities decrease over the light-crossing timescale of the respective zone through free-streaming escape and interactions with the external photon fields.

The instantaneous luminosity of the emitted CR densities (integrated over energy) is shown in Figure 11 (for diffusive escape from the blob). Note that the time axis is logarithmic and the timescales for the outer zones are larger by several orders of magnitude. At the transition between jet and BLR zone one can clearly see the impact of the decay or outflow of in-source densities into the BLR. Although the interaction rates are smaller in the outer zones, they contribute sizably to the final result owing to the much longer timescales.

Zoom In Zoom Out Reset image size

The result does not change if CR emission is beamed at the blob–BLR transition, where the emission would be distributed over a smaller volume but experience the same target density. However, the beaming may be important at the BLR–DT transition; here the geometry of the DT has to be such that the jet direction traverses the radiation field of the torus (which is not limited to the torus, but drops gradually outside its boundaries). We will show later that the DT contribution does not qualitatively affect our results, although neutrino production in the DT might be relevant in some cases.

The neutrino spectra emitted by an HL-FSRQ are shown in Figure 12 for the two escape assumptions discussed earlier. The neutrino fluences from each zone are shown separately and have to be added to obtain the total. The corresponding CR spectra are discussed in Appendix A. Since in the advection case more lower-energy CRs can escape the jet (having softer spectra), the BLR is populated with a higher initial density at $t={t}_{\mathrm{flare}}^{{\prime} }$. Therefore, a higher number of primaries are available for neutrino production in the outer zones compared to the diffusion case (darker blue shaded regions in Figure 12). For the pure proton injection (top right panel) and advective escape, the outer zones (in particular the DT) notably broaden the peak at tens of PeV and increase the flux by a factor of ∼2 compared to diffusive escape. This result is found to be consistent with previous studies (Murase et al. 2014).

Zoom In Zoom Out Reset image size

The injection of helium-4 (middle panels) results in a similar behavior, except that the jet contribution is smaller as already discussed in Section 4. For the same reason as in the proton case, the advection case produces more neutrinos. Since the BLR is optically thick for heavier nuclei, most of the secondaries from iron-56 injection (bottom panels) are absorbed in the BLR, such that neutrino emission from the DT is lower compared to the other two cases. Since HL-FSRQs clearly lie inside the nuclear cascade regime, they emit mostly nucleons and photodisintegrated secondary nuclei.

Independent of the injection composition, the escape through advection produces more neutrinos since more CRs are injected in the BLR, which contributes at the same level as the production in the jet. Although the IR photon density is equal in the BLR and the DT (Figure 9), the latter zone is more relevant for neutrino production for proton or helium-4 injection because of its much larger size—which increases its optical depth to photopion production compared to the BLR (see Figure 10).

It has long been proposed that the peak frequencies of the blazar SEDs are anticorrelated with their bolometric luminosities, which is referred to as "blazar sequence" (Fossati et al. 1998). This trend becomes more prominent from the recent analysis (Ghisellini et al. 2017) of a sample of 747 blazars with measured redshifts from the Fermi-3LAC catalog, compared to the small sample of merely 33 EGRET blazars in Fossati et al. (1998). Despite the debate on whether this is influenced by a selection bias due to the limitation of the sensitivity of the observation instruments, this phenomenological prescription does adequately describe at least the sample of the observed blazars. Therefore, it provides us with simplified and representative SEDs in each luminosity bin. We list the used luminosity bins in Table 1.

Table 1.  Models and Parameters Considered in This Work

ID	Type	N	Lγ	${B}_{\mathrm{jet}}^{{\prime} }$	rBLR	rDT	Lν/Lγ	Lν/Lγ
		in 3LAC	(log10, erg s−1)	(mG)	(pc)	(pc)	(p)	(Fe)
11	HL-FSRQ	49	50.3	5000	2.54	64	1.7 × 10−1	8.7 × 10−2
10			49.6	2300	0.73	16	2.4 × 10−1	4.7 × 10−2
9 (B)			48.8	900	0.18	4.0	1.5 × 10−1	2.5 × 10−2
8	LBL & FSRQ	384	47.7	260	⋯	⋯	1.2 × 10−2	1.2 × 10−4
7			46.1	40	⋯	⋯	4.2 × 10−4	1.8 × 10−5
6	IBL & FSRQ	285	45.7	28	⋯	⋯	2.1 × 10−4	2.1 × 10−5
5			45.5	21	⋯	⋯	8.6 × 10−5	1.8 × 10−5
4 (A)			44.6	7.1	⋯	⋯	3.3 × 10−6	1.2 × 10−6
3	HBL	29	43.5	2.2	⋯	⋯	1.1 × 10−7	4.1 × 10−8
2			42.5	0.69	⋯	⋯	3.9 × 10−9	1.4 × 10−9
1			41.6	0.23	⋯	⋯	1.4 × 10−10	4.7 × 10−11

Note. The third column shows the number of observed blazars with measured redshift in the Fermi-3LAC catalog (Ackermann et al. 2015). Note that among model IDs 4 to 8, there are large overlaps between the FSRQ and BL Lac populations. From model IDs 1 to 8, we ignored the BLR and DT zones, since either they are negligible in all BL Lacs or the blob in the jet lies outside the regions in those FSRQs. The last columns show the results from the simulations (obtained neutrino-to-gamma-ray luminosity ratio for a CR loading f_CR = 10 and considering diffusive CR escape).

Download table as:  ASCIITypeset image

The higher-luminosity bins are mostly populated by FSRQs from observations, as in models 7–11 of Table 1. In our scenario, an emitting blob of a fixed radius and distance to the central black hole is assumed, but all the physical quantities in the table increase with luminosity. Therefore, the blob can be located inside (see Table 1, models 7–8) or outside the BLR and dusty torus (Table 1, models 9–11; see discussion in Section 2.1). The corresponding target photon SEDs in the blob frame are shown in the top five curves of Figure 13 for models 7–11, respectively. In particular, for the top three curves (models 9–11), the contributions from the BLR and DT are clearly seen, on top of the smooth two-hump structure (photons produced from the blob internally). For models 7 and 8, the contributions of the external photon fields are (by construction) not present in the interactions. The lower-luminosity bins are mainly populated by BL Lacs, corresponding to models 1–6 of Table 1 and the bottom six curves of Figure 13, where the BLR and dusty torus contributions are absent. To account for the subthreshold low-luminosity blazars, we roughly extrapolated the blazar sequence down to ${L}_{\gamma }\sim {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (models 1–4).

Zoom In Zoom Out Reset image size

We now apply our model to the sources of the blazar sequence. For low-luminosity sources only the radiative processes in the jet will be considered, given the geometry arguments presented before, while the simulation of HL-FSRQs will follow the three-zone model introduced in Section 6.

Figure 14 shows the maximum CR energy ejected by each source, for different injected isotopes and for either of the escape assumptions discussed in Section 5. For HBLs, where photonuclear/photohadronic interactions are subdominant up to the highest injection energies, the maximum energy is determined from the Hillas condition in the blob. As a result, we have ${E}_{\max }\propto {B}^{{\prime} }\propto {L}_{\gamma }^{1/2}$ along the blazar sequence, since a constant energy partition ${\epsilon }_{{\rm{B}}}$ is assumed for the magnetic fields. In the range indicated by the gray band, the maximum energy is proportional to the charge of the injected isotope, which corresponds to the requirement that cosmic accelerators follow the Peters cycle (Peters 1961; rigidity-dependent maximum energy), an often implied relationship in textbook physics such as the Hillas plot and in recent phenomenological studies (e.g., Aab et al. 2017b). This relationship breaks down at the transition between the nuclear survival and nuclear cascade regimes at ${L}_{\gamma }\gtrsim {10}^{45}\mbox{--}{10}^{46}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, where photohadronic interactions become efficient and limit the maximum energy of ejected CRs. A further increase of the luminosity results in a reduction of the maximum energy for the reasons discussed in Section 3. The effect is more pronounced for heavier injection isotopes, as we can see from the three cases displayed in Figure 14. Note that the maximum energy is not significantly affected by the escape scenario.

Zoom In Zoom Out Reset image size

We now turn to the production efficiency of UHECRs and neutrinos across the blazar sequence. We define the efficiency in the production of UHECRs as the ratio of the total energy of the time-integrated spectrum of all ejected CR species above 10⁹ GeV in the black hole frame to the total injected energy in CRs above the same energy as

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{CR}}=\displaystyle \frac{{\displaystyle \sum }_{\mathrm{all}\mathrm{CRs}}{\int }_{0}^{\infty }{\int }_{{10}^{9}\mathrm{GeV}}{Q}_{\mathrm{esc}}^{{\prime} }({E}^{{\prime} },{t}^{{\prime} })\,{{dE}}^{{\prime} }\,{{dt}}^{{\prime} }}{{\int }_{0}^{\infty }{\int }_{{10}^{9}\mathrm{GeV}}{Q}_{\mathrm{inj}}^{{\prime} }({E}^{{\prime} },{t}^{{\prime} })\,{{dE}}^{{\prime} }\,{{dt}}^{{\prime} }},\end{eqnarray} \tag{ 12 }$$

where ${Q}_{\mathrm{inj}}^{{\prime} }$ is the injection rate from the acceleration zone of the injected isotope, as defined in Equation (5). Similarly, the neutrino production efficiency ${\epsilon }_{\nu }$ is obtained by performing the sum in the numerator over all neutrino species, and both integrals over all energies.

This quantity does not depend on the CR loading factor f_CR of the sources, which cancels out in Equation (12) (as opposed to the numbers listed in the last columns of Table 1, which are computed for f_CR = 10).

The neutrino production efficiencies are shown in the left panel of Figure 15. Clearly, HL-FSRQs can be identified to be the best neutrino emitters of the blazar sequence, with ${\epsilon }_{\nu }\sim {10}^{-3}\mbox{--}{10}^{-1}$. For less luminous blazars the photon field densities are too low for efficient neutrino production, with HBL Lacs being the least efficient neutrino sources. From the same plot we can see that HBL Lacs follow the expected relation ${L}_{\nu }\sim {\epsilon }_{\nu }{L}_{\gamma }\sim {L}_{\gamma }^{2};$ above ${L}_{\gamma }={10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, this scaling becomes weaker, and beyond ${L}_{\gamma }={10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, it becomes L_ν ∼ L_γ. There is hardly any dependence on the UHECR escape mechanism, apart from the ∼2 times more efficient production for the advective case, as discussed in Section 6.

Zoom In Zoom Out Reset image size

In the right panel of Figure 15 we plot the UHECR transfer efficiency. The overall tendency is the opposite of the neutrino efficiency: HL-FSRQs cannot efficiently emit UHECRs owing to strong photohadronic interactions and also to the additional radiation zones that suppress the UHECR part of the spectrum, while lower-luminosity objects such as LL-FSRQs and BL Lacs may allow for more efficient CR survival and escape. HBL Lacs are expected to be efficient sources of UHECR regardless of the injected isotope, since 100% of the injected UHECRs are emitted. In the advective escape scenario (dashed curves in Figure 15) blazars become efficient UHECR source candidates. Note, however, that the spectra of emitted CRs are generally softer. For diffusion-dominated CR escape, the production efficiency is low. The UHECR transfer efficiency should be interpreted together with the maximum achievable energy (see Figure 14): while, for example, at low energies a large fraction of the injected energy is transferred, the magnetic field is too low to allow for very high UHECR energies of light isotopes. Therefore, each injection isotope has a "sweet spot" where UHECR can be accelerated to sufficiently high energies and escape from the source environment. For example, for helium injection, this sweet spot seems to be around 10⁴⁶ erg s⁻¹. If, in addition, neutrinos are to be efficiently produced, higher luminosities are preferred; L ∼ 10⁴⁸ erg s⁻¹ seems to be a good trade-off between neutrino and UHECR production.

It is worth noting that in Ghisellini et al. (2017), with the inclusion of 747 blazars from the Fermi-3LAC catalog, the blazar sequence becomes more prominent and the SEDs for each luminosity bin can be constructed with more refined features. In comparison with the blazar sequence considered in this work, the main differences appear in the gamma-ray slopes and in the luminosity ratio of the second to the first hump (Compton dominance), in particular for low-luminosity blazars. Since the overall shapes of the SEDs are similar, as demonstrated in Figure 9 of Ghisellini et al. (2017), we do not expect significant changes of the neutrino and CR spectra obtained in this work.

We have studied blazars as sources of UHECRs and neutrinos including the injection of isotopes heavier than hydrogen into the jet. We have identified two important regimes, depending on jet luminosity and size of the blob: In the nuclear survival regime, corresponding to low luminosities or large blob sizes, the source is found to be optically thin to photonuclear interactions. Neutrinos are mostly produced as a result of photomeson production off the primary nuclei, and the maximum CR energies behave as expected from Peters cycle, i.e., E_max ∼ Z. In the nuclear cascade regime, corresponding to high luminosities or small blob sizes, the source is optically thick to disintegration of the injected nuclei at the highest energies, and the nuclear cascade becomes more populated. Neutrinos are copiously produced in photomeson processes off the secondary nuclei generated in the nuclear cascade, and the maximum CR energies saturate owing to the photonuclear processes. Compared to GRBs, the nuclear cascade ceases faster with decreasing mass number because the disintegration rate depends more strongly on energy. In addition, the neutrino peak flux and its energy depend on the injected nuclear composition, where heavy injection compositions are more in favor of observed data.

The evolution of neutrino and CR production efficiencies has been studied over the blazar sequence, including the corresponding SEDs of target photon spectra based on observations. In order to have a more refined model for HL-FSRQs, for which the jet production region may lie within the BLR, we have included the external radiation fields in the jet and the propagation of CRs through the BLR and DT radiation fields. Our main conclusion is that the neutrino production efficiency strongly increases with luminosity from an extremely low value for HBLs, over LBLs/FSRQs, and then saturates for HL-FSRQs. For HBLs, we observe a scaling ${L}_{\nu }\propto {L}_{\gamma }^{2}$, whereas for HL-FSRQs, we find L_ν ∝ L_γ. On the contrary, the UHECR transfer efficiency is high for HBLs and then decreases for LBLs/FSRQs, and it is low for HL-FSRQs. This opposite behavior is expected since the strong radiation fields lead to photonuclear disintegration and effective neutrino production at the same time. If, in addition, the maximum achievable CR energy is taken into account, it seems that the optimal luminosity for UHECR acceleration and escape is around ${L}_{\gamma }\sim {10}^{46}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ for isotopes heavier than helium. The precise level of disintegration can be somewhat controlled with the blob size, which is typically a parameter coming from self-consistent radiation models.

Note that while our observations hold for individual sources, the contribution of HBLs versus LBLs versus FSRQs to the diffuse CR and neutrino fluxes depends on the luminosity and redshift distribution functions of the sources. However, for individual sources there may be a "sweet spot" in the LBL/FSRQ range around ${L}_{\gamma }\sim {10}^{48}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ (for helium injection), where both UHECRs and neutrinos are efficiently produced—which corresponds to high-luminosity BL Lacs and intermediate-luminosity FSRQs. In order to test the origin of CRs with neutrinos, this may be the region to search for correlations. Furthermore, note that we use a one-zone model for the jet, whereas multizone models may lead to an enhanced neutrino production for BL Lacs by enhancing the radiation fields from the relative motion of the zones (Tavecchio et al. 2014; Tavecchio & Ghisellini 2015). However, in these cases the level of photonuclear disintegration will be enhanced as well, which means that the anticorrelation between neutrino and CR production continues to hold.

On a side note, we have also discussed the impact of the UHECR escape mechanism from the jet, where we have focused on a diffusion scenario, leading to hard ejection spectra motivated by recent Auger fits, and an advection scenario, leading to unmodified escape at lower energies (compared to the acceleration spectrum) and motivated by the physics of collimated jets without adiabatic cooling. While we have not found qualitative differences in our observations for neutrinos versus CRs over the blazar sequence, the ejected CR spectrum and composition, as well as the transfer efficiency to UHECRs, can be quantitatively very different. Both the spectrum and composition of the diffusion scenario agree better with recent Auger observations, whereas the advection scenario leads to more efficient UHECR ejection. The overall neutrino production is hardly affected by the escape assumption, but we observe a significant contribution of neutrino production in the BLR of HL-FSRQs for advective escape, which is not present for diffusive escape. Note that since the assumed diffusive escape mechanism is rather conservative, and the advective escape assumption quite aggressive, we expect that any realistic scenario should in fact lie in between.

We conclude that HL-FSRQs are very efficient neutrino emitters, whereas UHECRs from blazars will efficiently come from low- and intermediate-luminosity objects. These conclusions hold in the one-zone model for the jet zone with target photons derived from observed SEDs. These observations may be especially relevant for individual source correlation studies and catalog stacking searches, whereas the contributions to the diffuse flux require the input from the luminosity distribution. For example, in spite of their low neutrino luminosity, a sufficiently large number of LBLs that are not detectable otherwise can still dominate the diffuse neutrino flux. Another variable is the CR loading, which may vary for different blazar families or evolve over the blazar sequence. Further studies will reveal whether the UHECR and neutrino fluxes can be powered at the same time and how the recent stacking results for blazars need to be interpreted.

We would like to thank Kohta Murase for fruitful discussions.

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant No. 646623).

In this section we discuss the effect of the external zones of an HL-FSRQ and the escape assumption on the ejected CR spectrum. In Figure 16 the boarders of the shaded regions represent the total CR fluence transferred between the successive zones. Contrary to the behavior in neutrinos (compare with Figure 12), the CR fluence is reduced owing to photonuclear losses in the outer zones. As already discussed in Figure 15, the diffusive escape assumption results in an approximately 10-fold lower CR output from the jet zone, which is even further suppressed at high energies during propagation through the BLR and the DT (compare left and right panels). The net impact of the outer zones is therefore smaller for diffusion, in agreement with the behavior of the ejected neutrino spectra. In the advection case this impact is more pronounced for heavier injection composition, since the optical thickness to photonuclear interactions is higher for nuclei. The outer zones do not affect the general conclusion about the spectral index, which remains softer for advection compared to diffusion.

Zoom In Zoom Out Reset image size

We have so far assumed that the energy budget of the magnetic field corresponds to a constant fraction of the total energy of the jet. It is the simplest assumption since we do not explicitly compute the components of the SED, as explained in the last paragraph of Section 2.2. In this appendix we explore parameter ranges that are closer to those considered in SSC models that are able to generate the SEDs (e.g., Tavecchio et al. 2010; Dermer et al. 2014, where the hadronic components are subdominant). In such models, the ratio of the energy densities u_B to u_syn can be immediately derived from the SSC parameter Y_SSC, which roughly corresponds to the luminosity ratio of the second hump to the first hump in the observed SED. Then, the magnetic field strength is seen to lie in the range 0.1–1 G for BL Lacs and B ∼ a few G for FSRQs. In addition, studies of diffusive shock acceleration (Inoue & Tanaka 2016) suggest a lower value for the acceleration efficiency compared to our choice of η = 1 (see Equation (6)).

The values of B derived from SSC models can be well described by a power-law relation $B^{\prime} \sim {L}_{\gamma }^{1/5}$ along the blazar sequence (in contrast to ${L}_{\gamma }^{1/2}$ as in the main text, which yields much lower magnetic fields for BL Lacs). This results in a running value of the energy partition factor ${\epsilon }_{B}={\epsilon }_{B}({L}_{\gamma })$ instead of a constant one. As a consequence, the dimmest HBL has a magnetic field strength of ${\boldsymbol{B}}^{\prime} =0.09\,{\rm{G}}$, while the brightest FSRQ has B' = 5 G, both consistent with the results from Tavecchio et al. (2010) and Dermer et al. (2014). Furthermore, we now decrease the acceleration efficiency to η = 0.1.

In Figure 17 we show the maximum injection energy of three isotopes under these conditions. For low-luminosity BL Lacs, in spite of the lower acceleration efficiency, the high magnetic fields lead to a higher maximum injection energy, regardless of the injected element (see Figure 5). The separation between the nuclear survival and nuclear cascade regimes does not change noticeably. The increase in maximum energy is also observed in the ejected neutrino and CR spectra for low-luminosity blazars, as shown in Figure 18. The increased maximum CR energies in the jet lead to a 10-fold increase in ejected neutrinos. Recall, however, that blazars in the nuclear survival case are generally not very efficient neutrino sources (see Figure 15). On the other hand, for the ejected CRs the different parameter set leads to a flux decrease by one order of magnitude in the case of diffusive escape (see middle panel of Figure 18). That is because the Larmor radius cannot reach the size of the region, even at the highest energies, and therefore CRs are always affected by confinement. For advective escape, the normalization of ejected CR spectrum is not affected, since this escape mechanism does not depend on the magnetic field strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have also studied the effect of this parameter set on the maximum ejected energy across the entire blazar sequence, as shown in Figure 19. We note that the qualitative conclusions drawn in Section 7 are unchanged: there is an intermediate-luminosity range where the ejected CR energy is maximal. The maximum ejection energy in the low-luminosity range is affected by the softer scaling of ${E}_{\max }\propto {B}^{{\prime} }\propto {L}_{\gamma }^{1/5}$. As a result, CRs ejected by HBLs can reach energies that are 10 times higher compared to the result in Section 7 (see Figure 14). For iron-56 injection, the luminosity range in which the ejected CRs reach the highest energies moves from LBL/FSRQ to HBL/LBL. This aspect, however, is not observed for helium-4 and protons.

Zoom In Zoom Out Reset image size

Finally, in Figure 20 we show the impact on the transfer efficiency (see Figure 15). While the conclusions are qualitatively unaffected, the modified parameter set results in a higher neutrino production efficiency of HBLs and an order of magnitude lower UHECR transfer efficiency of BL Lacs and low-luminosity FSRQs. On the other hand, under these conditions all BL Lacs down to the lowest luminosities are capable of ejecting CRs above 10⁹ GeV, as opposed to the parameter set discussed previously, where the low magnetic fields prevented HBLs from ejecting UHECRs in the case of a diffusive escape mechanism.

Zoom In Zoom Out Reset image size