© The Authors 2024

1 Introduction

The origin of the cosmic rays and the astrophysical neutrinos is currently one of the central issues in high-energy astrophysics. The IceCube experiment, located at the South Pole, observes a diffuse flux of astrophysical neutrinos with energies between 100 TeV and a few petaelectronvolts (Aartsen et al. 2013a, 2014, 2015a,b, 2013b), thought to be emitted when cosmic rays interact with ambient radiation (pγ) or matter (pp). Crucially, although these hadronic processes should be accompanied by the emission of high-energy gamma rays, so far the cosmic neutrino sky does not seem to significantly correlate with gamma-rayemitting sources. This suggests that the photons emitted in cosmic ray interactions are being attenuated or otherwise losing their energy, either in the astrophysical source or en route to Earth (Murase et al. 2016; Hooper 2016; Halzen et al. 2019; Murase 2020; Halzen & Kheirandish 2020).

One of the most promising candidate source classes is that of blazars, which are active galactic nuclei (AGN) whose relativistic jet is pointing in a direction close to our line of sight (Urry & Padovani 1995). Because of this, the radiation that is emitted by the jet is relativistically boosted and so are the neutrinos potentially produced in the jet.

At the population level, gamma-ray-bright blazars have been excluded as a major contributor to the IceCube diffuse flux and limits have been set on their collectively emitted flux of teraelectronvolt-petaelectronvolt neutrinos (Aartsen et al. 2017b; Abbasi et al. 2023), although these limits depend on the average spectral shape of the emitted neutrino energy flux. Palladino et al. (2019) have shown that for the blazar population to contribute significantly to the diffuse sub-petaelectronvolt IceCube flux, a numerous population of very low-luminosity sources would have to be very highly loaded in cosmic rays, while the sources that are brightest in gamma rays must necessarily be poor cosmic-ray accelerators, given the overall lack of associations mentioned above. Such a trend has also been suggested by multiepoch blazar modeling (Petropoulou et al. 2020a). Furthermore, Petropoulou et al. (2022) have recently proposed a model of magnetically loaded jets that may provide a natural explanation for this inverse relation between gamma-ray luminosity and cosmic-ray loading.

In spite of their subdominant contribution to the diffuse IceCube astrophysical flux, some observations seem to point increasingly toward some level of neutrino emission from blazars. From one side, this comes in the form of statistical correlations with certain blazar subclasses. For instance, Giommi et al. (2020a) have demonstrated an association at the 3.2σ level between IceCube tracks and the positions of high- and intermediate-synchrotron peaked BL Lacs (HBLs and IBLs respectively, see also Padovani & Resconi 2014; Padovani et al. 2016). More recently, Buson et al. (2022), 2023) have highlighted a significant spatial association between a blazar catalog and the IceCube sky map with a p value ~10⁻⁷, while subsequently Bellenghi et al. (2023) obtained no significant correlation when utilizing an updated dataset.

Another important source of evidence are sporadic spatial associations between high-energy IceCube events and the positions of individual blazars. The blazar with the highest significance thus far is TXS 0506+056, a BL Lac object for which neutrino emission has been established at the 3.5σ level. This includes two observations: (1) a 290 TeV neutrino was detected by IceCube in 2017 in spatial coincidence with the source and in temporal coincidence with a six-month gamma-ray flare (Aartsen et al. 2018); and (2) archival IceCube data show an excess of neutrinos with energies of tens of teraelectronvolts around the source’s position during 2014 and 2015, although this coincided with a low state in gamma rays (Aartsen et al. 2018). Furthermore, after 10 yr of observations, the second hottest spot in the IceCube sky was coincident with the position of TXS 0506+056 (Aartsen et al. 2020), which can be regarded as additional evidence in support of this blazar as a neutrino source candidate. The hottest spot from the same analysis is consistent with the position of the Seyfert galaxy NGC 1068 (Abbasi et al. 2022), which also hosts an active black hole but is not a blazar.

Beside TXS 0506+056, IceCube has also detected high- energy events in spatial coincidence with other individual blazars of different classes, among which are PKS 142441 (Kadler et al. 2016; Gao et al. 2017), GB6 J1040+0617 (Garrappa et al. 2019), 3HSP J095507.9+355101 (Giommi et al. 2020b; Paliya et al. 2020), PKS 1502+106 (Franckowiak et al. 2020), and PKS 0735+178 (Sahakyan et al. 2022). These individual associations have been subsequently dissected by theoretical works by means of numerical leptohadronic modeling. For the 2017 flare of TXS 0506+056, some models suggest that the high-energy gamma rays co-emitted with the neutrinos interact through pair production with lower-energy photons in the jet, generating electromagnetic cascades whose energy is ultimately emitted by the source in the X-ray range (e.g., Cerruti et al. 2019; Gao et al. 2019; Keivani et al. 2018; Sahu et al. 2020). It is worth noting that while the very-high-energy gamma rays produced through hadronic processes cascade down to the X-ray regime, in the gigaelectronvolt gamma-ray range the source is optically thin in these models, which allows them to explain the bright state in gigaelectronvolt gamma rays observed during neutrino emission. The peak of this gigaelectronvolt emission as well as the peak of the optical emission are typically dominated by leptonic processes, while a potential hadronic signature is limited to the X-ray and, potentially, teraelectronvolt ranges.

As shown by Padovani et al. (2019), data show that TXS 0506+056 possesses a broad line region (BLR) surrounding the active black hole, which is a feature typical of flat-spectrum radio quasars (FSRQs). This leads to the reclassification of TXS 0506+056 as a masquerading BL Lac. This was also shown by Sahakyan et al. (2022) to be the case for the recent IceCube candidate PKS 0735+178, which motivated the authors to include a BLR when modeling the source. More recent modeling results by Acharyya et al. (2023) also support the existence of a BLR photon field interacting with the jet emission region in this source.

The presence of external photon fields form the BLR can lead both to enhanced neutrino production and more extensive electromagnetic cascades, thus spreading the hadronic gamma rays in energy down to the megaelectronvolt or even X-ray regime. This external field model has been suggested as a viable scenario to explain the 2017 neutrino event from TXS 0506+056 and the simultaneous multiwavelength data (Keivani et al. 2018; Petropoulou et al. 2020a; Zhang et al. 2020) and has since been successfully applied to other candidate sources (e.g., Padovani et al. 2019; Rodrigues et al. 2021a; Petropoulou et al. 2020b; Oikonomou et al. 2021; Gasparyan et al. 2021; Sahakyan et al. 2022). This model generally predicts a connection between neutrino emission and the fluxes in either X-rays or megaelectronvolt gamma rays, depending on the spectral shape of the external fields assumed and the bulk Lorentz factor of the jet.

Under certain conditions, theoretical works have shown that electromagnetic cascades can be extremely efficient in attenuating gigaelectronvolt gamma rays, to the point where neutrino emission could theoretically be correlated with a flare in the megaelectronvolt band while a low flux could simultaneously be observed by the Fermi Large Area Telescope (LAT) in gigaelectronvolt gamma rays. This scenario has been tested on the 20142015 emission from TXS 0506+056, by invoking high-density external photon fields originating in a BLR (Reimer et al. 2019; Rodrigues et al. 2019; Petropoulou et al. 2020a) or in the black hole corona, with the neutrino emission zone lying in the close vicinity of the black hole (e.g., Xue et al. 2021). Petropoulou et al. (2020b) have also tested a similar corona model for the extreme blazar 3HSP J095507.9+355101, leading to very efficient neutrino production accompanied by bright emission of megaelectronvolt photons. On the one hand, this type of model could in principle help explain the general lack of associations between IceCube events and Fermi-LAT sources, which inspired the idea of gamma-ray-suppressed neutrino flares (Kun et al. 2023). On the other hand, such high optical thicknesses generally require an environment with extremely high photon densities or Doppler boosts, which are possible only in very specific regions of parameter space. Under less extreme conditions, we should more likely expect some level of gigaelectronvolt emission during neutrino flares, as shown by the more “conventional” one-zone models cited previously. Furthermore, during hadronic flares there should be an overall increase in gamma-ray emission in the jet and it is unclear whether the additional attenuation is always sufficient to counterbalance this effect.

In this work, we perform individual modeling of 324 gammaray blazars from the Candidate Gamma-ray Blazar Survey (CGRaBS) catalog (Healey et al. 2008), 237 of which are FSRQs. This corresponds to 34% of all 619 FSRQs in the latest Fermi-LAT catalog by Ajello et al. (2020).

The sources in this sample have relatively high gamma-ray luminosities, typically above 10⁴⁵ erg s⁻¹. On the one hand, the majority of the models cited above predict these sources to be efficient neutrino emitters (e.g., Murase et al. 2014); on the other hand, as discussed above, bright blazars are the most constrained by IceCube stacking analyses (Aartsen et al. 2017b; Abbasi et al. 2023), which can place strong limits on their overall emission. The objective of this work, besides characterizing the physical parameters of each source in the sample through model fitting, is to provide physics-driven predictions of their potential neutrino emission, both at the individual level and as a population.

We model each source in the sample by numerically simulating the radiative interactions of cosmic-ray electrons and protons accelerated in the relativistic jet, in a one-zone framework. We then compare the predicted multiwavelength emission from each source to multiwavelength observations based on public data ranging from the radio band up to gamma rays, compiled by Paliya et al. (2017). We fit the data at frequencies above 300 GHz and thus constrain the source parameters.

In those cases where a proton contribution helps us explain the multiwavelength data, we present the emitted neutrino spectrum calculated self-consistently; in those cases where the multiwavelength data are well explained solely by electron emission, we set upper limits on the neutrino emission from the source. In the case of FSRQs, we additionally consider external photon fields from thermal and broad line emission present in the BLR. We then connect the individual results of the modeling of each source to draw conclusions on the sample as a whole, with an emphasis on potential proton acceleration and neutrino emission. Compared to previous leptohadronic studies of blazar samples (see e.g., Tavecchio et al. 2010; Boettcher et al. 2013; Petropoulou et al. 2015; Cerruti et al. 2015; Oikonomou et al. 2019; Liodakis & Petropoulou 2020), this work presents the largest self-consistently modeled sample of blazar AGN to date.

Some of the IceCube candidates mentioned above are included in this sample and we discuss their potential as neutrino emitters in the context of this model, while also comparing our results with the literature. We place blazar TXS 0506+056, which is not included in our sample, in the context of the results obtained. Finally, we show how several blazars in this sample may be promising source candidates for future neutrino detectors such as IceCube-Gen2 (Aartsen et al. 2021).

The paper is organized as follows: in Sect. 2 we describe the sample, multiwavelength data, physical model, and fitting method. In Sect. 3, we present the best-fit parameters, the predicted cosmic-ray acceleration by each source, and the respective flux of emitted neutrinos. In Sect. 4, we discuss the impact of our results in the context of future multi-messenger searches and suggest potential candidates in the sample for neutrino source searches that may help support or exclude the leptohadronic paradigm for high-luminosity blazars. In Sect. 5, we discuss some limitations of the model and we summarize our main conclusions in Sect. 6.

2 Methods

2.1 Source sample

The CGRaBS catalog (Healey et al. 2008) is a flux-limited (F_{8.4 GHz} > 65 mJy) sample of 1625 radio-loud AGN thought as potential gamma-ray emitting quasars to be detected with the Fermi-LAT. From this catalog, Paliya et al. (2017) selected 324 sources that were present in any of the Fermi-LAT catalogs available at the time and which also have multifrequency observations. From these observations, they put together a multiwavelength SED for each source. This includes near-infrared data from ESO’s GROND (Greiner et al. 2008), optical-ultraviolet data from the Swift Ultra-Violet Optical Telescope (UVOT, Roming et al. 2005), X-ray data from NuSTAR (Harrison et al. 2013), XMM-Newton (Matthews et al. 2001), Chandra (Weisskopf et al. 2000), and the Swift X-ray Telescope (XRT; Burrows et al. 2005) and Burst Alert Telescope (BAT; Krimm et al. 2013), as well as gamma-ray data from the Fermi- LAT (Atwood et al. 2009). They have also relied on archival data provided by the ASI Data Center (ASDC). Regarding Fermi- LAT data, 310 sources are included in the 3LAC (Acero et al. 2015), which was the most up-to-date at the time of the analysis by Paliya et al. (2017), 8 in the 2LAC (Ackermann et al. 2011), and two in the 1LAC (Abdo 2010). We refer the reader to Paliya et al. (2017) for more details on the multiwavelength data.

The upper panel of Fig. 1 shows a sky map with the positions of the objects in the sample. This amounts to a total of 324 blazars, including 237 FSRQs, 38 low-frequency-peaked BL Lacs (LBLs), 22 intermediate-frequency-peaked BL Lacs (IBLs), 17 high-frequency-peaked BL Lacs (HBLs), 3 narrowline Seyfert 1 galaxies (NLS1), and 3 sources classified as blazar candidates of uncertain type (BCUs) in the Fermi 3LAC or 4LAC (Acero et al. 2015; Ajello et al. 2020). Regarding the physical model, we apply to NLS1 galaxies the same geometric assumptions as used to describe FSRQs (see e.g., Paliya et al. 2019), while for BCUs we make the conservative assumption that they are described by the simpler BL Lac model (cf. Sect. 2.2 where these differences in the modeling are described).

In the lower panels of Fig. 1, we show the distribution of sources as a function of redshift and gamma-ray luminosity. The sources in the sample are those shown as dark markers (FSRQs in the left and BL Lacs on the right). For comparison, we show as light markers a source distribution sampled from the luminosity function of FSRQs on the left (Ajello et al. 2012) and BL Lacs on the right (Ajello et al. 2014). Above the Fermi-LAT sensitivity, represented here as a line corresponding to a flux of 4 × 10⁻¹² erg cm⁻² s⁻¹, this distribution of light markers follows the population of detected sources. Under that line, it is an extrapolation obtained by Monte Carlo sampling of the luminosity functions derived by Ajello et al. (2012), 2014). As we can see, the current sample covers a range of luminosities and redshifts comparable to that of the source population detected by the Fermi-LAT. At the same time, as we can see in the histograms alongside the margins of the lower plots, the present sample (black histograms) does not follow exactly the population of Fermi-LAT-detected sources (green and purple histograms). In particular, our FSRQ sample somewhat over-represents sources with L_γ ~ 10⁴⁷ erg s⁻¹ and under-represents the source population for L_γ ≳ 3 × 10⁴⁷ erg s⁻¹. Regarding BL Lacs, we see that sources lying at z ≳ 1 as well as L_γ ≲ × 10⁴⁵ erg s⁻¹ are underrepresented in our sample. These considerations will play a role in the extrapolation of our results to the Fermi-LAT population in Sect. 4.3.

The redshift of each source is adopted from the respective Fermi-LAT catalog, while the black hole mass and accretion disk luminosity values are those reported by Paliya et al. (2017). In each case, these were deduced following one of the following three approaches: 1) by modeling the big blue bump, if detected at infrared to ultraviolet frequencies, with the standard optically thick, geometrically thin accretion disk model (Shakura & Sunyaev 1973); 2) from the optical emission line and continuum luminosities by adopting the virial technique (cf. Shen et al. 2011; Paliya et al. 2021), or 3) when both of the above-mentioned methods were not possible, the disk luminosity was estimated assuming $L_{\text{disk }}~40L_{\gamma }^{0.93}$ (Sbarrato et al. 2012). In such cases, a black hole mass was assumed M_BH = 5 × 10⁸ M_⊙.

In this work, we perform individual leptohadronic modeling of each of the 324 blazars in this sample. The steps of this procedure are detailed in the following sections.

Fig. 1 Characterization of the sample. Upper panel: sky map with the positions of the sources studied in this work, which include 237 FSRQs, 88 BL Lacs, three narrow-line Seyferts, and three objects classified as BCUs in the Fermi catalogs. Lower panels: distribution of FSRQs (left) and BL Lacs (right) in the sample as a function of redshift and gamma-ray luminosity (dark points). For comparison, we show as lighter points to the distribution of the respective overall population, based on the luminosity functions by Ajello et al. (2012, 2014). The gray area marks the subthreshold region for the Fermi-LAT. The black distributions alongside the margins refer to the sample; the green (left) and purple (right) distributions refer to the population of sources detected by the Fermi-LAT (i.e., excluding the putative sources below the sensitivity threshold).

2.2 One-zone radiation model

We self-consistently model the emission of photons and secondary particles by a population of nonthermal electrons and protons interacting in the blob, using the time-dependent numerical framework AM³ (Gao et al. 2017).

We simulate the radiative zone in the jet as a spherical blob of radius¹ ${R^{\prime }}_{\text{b}}$ permeated by a uniform magnetic field of strength B′ and moving along the jet with a bulk Lorentz factor of Γ_b, as shown in Fig. 2. We assume the jet is observed at an angle θ_obs = 1/Γ_b relative to its axis, resulting in a Doppler factor δ_D = Γ_b. We refer to the distance of the blob to the supermassive black hole as the dissipation radius, R_diss.

Electrons (protons) are assumed to be accelerated to simple power-law distributions, $\text{d}{N}_{\text{e}(\text{p})}/\text{d}{\gamma }_{\text{e}}^{\prime }\propto {\gamma }^{\prime -{\alpha }_{\text{e}(\text{p})}}$, from a Lorentz factor ${\gamma }_{\text{e}(\text{p})}^{\prime \min }$ up to ${\gamma }_{\text{e}(\text{p})}^{\prime \max }$. The total power injected in electrons (protons) into the radiation zone is given by an injection luminosity $L_{\text{e}(\text{p})}^{\prime }$. These and other parameters of the model are listed in Table 1. We also assume that all particles, both charged and neutral, escape the radiation zone advectively, that is, the escape timescale is the light-crossing time of the blob, $t_{\text{esc}}^{\prime }=R_{\text{b}}^{\prime }/c$.

In every source simulation the time-dependent solver is evolved until a steady-state is reached, that is, an equilibrium between the injection and escape (or cooling) of electrons and protons. The emitted photon spectrum then becomes constant in time. At high enough energies, these photons interact efficiently with the EBL, leading to an attenuation of the gamma-ray fluxes at very high energies. We adopt the EBL model by Dominguez et al. (2011) and calculate this energy-dependent flux attenuation as a final step after each the source simulation. For this we use the energy-dependent attenuation length values tabulated in Fermipy (Wood et al. 2018).

Fig. 2 Schematic representation of two possible geometries of the one-zone FSRQ model. The broad line region is represented by the purple clouds. In the case of the blob represented on the left, its proximity to the BLR implies a large relative Doppler factor, making the external fields appear highly boosted in the jet frame (compare left-hand plot, where the external fields are shown in the rest frame of the central engine, with the middle plot, where they are shown as energy density spectra in the jet frame). In the case represented on the right, due to the high dissipation radius Rdiss the BLR photons impinge more from behind and the external fields from the BLR are less boosted, as shown in the lower right plot. The relativistic boosting of these external radiation fields follows the treatment by Ghisellini & Tavecchio (2009). We note that on the left plot, the broad lines are isotropic in the black hole frame, while the disk emission is anisotropic except for the ~1% that gets isotropized through Thomson scattering. This explains why in the middle plot the broad lines get more highly boosted than the disk emission.

2.3 External radiation fields

Out of the 324 sources in our sample 241 are FSRQs or narrowline Seyferts, which typically possess a bright accretion disk and a BLR surrounding the supermassive black hole. As an estimate of the disk luminosity Ldisk and black hole mass M_BH of each source we rely on the values derived by Paliya et al. (2017), as discussed in Sect. 2.1. Given the disk luminosity, we model the photon spectrum emitted by the disk as a multi-temperature Shakura-Sunyaev distribution, whose profile depends on the value of M_BH estimated for that same source (Shakura & Sunyaev 1976).

Following Ghisellini & Tavecchio (2009), we model the BLR as a thin shell located at a radius R_BLR = 10¹⁷(L_disk/10⁴⁵ erg s⁻¹)^1/2 cm around the supermassive black hole. We also consider a dusty torus of radius R_torus = 2.5 × 10¹⁸(L_disk/10⁴⁵ erg s⁻¹)^1/2 cm.

We assume that 1% of the luminosity of the thermal emission from the accretion disk is Thomson-scattered in the BLR, leading to an isotropic component in the BLR volume that has the same spectral shape as the disk emission and a total luminosity of 0.01 L_disk. A more significant 10% of the disk luminosity is isotropized in the volume of the BLR through broad line emission (Greene & Ho 2005). Following previous studies (e.g., Murase et al. 2014; Rodrigues et al. 2018, 2019), we consider the two lines that are typically strongest, namely the α lines of the Lyman series of hydrogen (with peak energy 10.2 eV and total luminosity 0.1 L_disk) and helium (40.8 eV and 0.05 L_disk). We model their spectra in the rest frame of the black hole as Gaussian distributions with a half-width of 5% of their respective peak energies. For the infrared emission from the torus we consider a fixed temperature of 500 K and a covering factor of 0.3. As a conservative assumption we ignore the contribution of potential corona emission, which in any instance would only contribute significantly to the blob system if it lied extremely close to the black hole.

The luminosity of the BLR photons, L_BLR , is converted into an energy density in the blob rest frame $u_{\text{BLR}}^{\prime }$ through the following relation: $$u_{\text{BLR}}^{\prime }\left(L_{\text{BLR}},{\Gamma }_{\text{b}},R_{\text{diss}},R_{\text{BLR}}\right)=\frac{{\delta }_{\text{BLR}}^{\text{rel}}{}^2\left({\Gamma }_{\text{b}},R_{\text{diss}},R_{\text{BLR}}\right)L_{\text{BLR}}}{4\pi R_{\text{ext}}^{2}c},$$(1)

where ${\delta }_{\text{BLR}}^{\text{rel}}$ is the relative Doppler factor with which the radiation is boosted into the blob, which depends on the blob’s relative position to the BLR. As illustrated in the diagram of Fig. 2, in the cases where R_diss ≤ R_blr we have ${\delta }_{\text{BLR}}^{\text{rel}}={\Gamma }_{\text{b}}$. For R_diss > R_blr, the external photons impinge on the blob less frontally, leading to lower effective Doppler factors. For example, for a value of Γ_b = 10 and R_diss = 4R_BLR, we have ${\delta }_{\text{BLR}}^{{\text{rel}}^2}\approx 0.1$ and the external fields from the BLR are unlikely to play a significant role in the blob emission. For the treatment of ${\delta }_{\text{BLR}}^{\text{rel}}$ we follow Ghisellini & Tavecchio (2009). Note that the energies of the BLR photons in the blob frame are also shifted by a factor $E^{\prime }={\delta }_{\text{BLR}}^{\text{rel }}$ E and therefore also depend on both Γ_b and R_diss.

For the infrared emission from the dust torus, Eq. (1) also applies, with the substitutions R_BLR → R_torus and ${\delta }_{\text{BLR}}^{\text{rel}}\to {\delta }_{\text{torus }}^{\text{rel }}$. For all the cases tested here, the blob lies within the radius of the torus, so that ${\delta }_{\text{torus }}^{\text{rel }}={\Gamma }_{\text{b}}$.

Fig. 3 Flowchart of the fitting procedure. The leptohadronic fit of the data from each source is performed as a three-step process. The initial step, where the optical and gamma ray broadband features are explained with a leptonic model, requires the most simulations because the vast parameter space is initially unconstrained. The number of CPU-hours quoted refers to the entire sample and is based on the code’s performance of 1.6 s for a single leptonic simulation and ~20 s for a leptohadronic simulation.

2.4 Parameter search

We search the parameter space of the leptohadronic model using a three-step method, represented schematically in Fig. 3. First, for each source the data in the infrared, optical, ultraviolet and gamma-ray bands are fit using a purely leptonic model. This is motivated by the fact that the optical and GeV gamma-ray broadband features of blazars can typically be well explained with leptonic models, even if protons are also assumed to be co-accelerated in the source, as shown by several leptohadronic modeling efforts of IceCube source candidates (e.g., Gao et al. 2019; Rodrigues et al. 2019; Sahakyan et al. 2022). At this stage we treat X-ray fluxes merely as upper limits for electron emission, although for some sources the X-ray data may actually be fit well by this emission. Radio fluxes (below 300 GHz) are considered as upper limits for the fit, because the relatively small radiation zone responsible for the emission of high-energy photons and neutrinos is necessarily optically thick to radio emission due to synchrotron self-absorption. These fluxes must therefore originate from larger emission regions, presumably located downstream in the relativistic jet.

By fitting the optical and GeV gamma-ray peaks with a leptonic model, we find the best-fit values for the four parameters describing the accelerated electron population and the four parameters describing the source properties. These are the first eight parameters listed in Table 1. We search this parameter space using a genetic algorithm (e.g., Goldberg 1989), a metaheuristic whose stochastic nature is adequate for searching large and poorly constrained parameter spaces (see Rodrigues et al. 2019, 2021a, where a similar method was used). Rather than requiring an initial guess, this method relies on a population of parameter sets that is originally randomly distributed within certain parameter boundaries. The values of these boundaries are given in Table 1 for each parameter. In every iteration, the subsequent generations are optimized through the genetic heuristic.

The goodness of fit is judged by means of a chi-squared test between the predicted spectrum and the multiwavelength data in the frequency bands we wish to fit. For data points falling in the radio, X-rays, and gamma rays above 20 GeV, the model only adds to the chi-squared value if the predicted flux overshoots the observations. This allows us to negatively score leptonic solutions that overshoot the observed frequencies without actually fitting data in these bands: $$x_{\text{partial}}^2=\frac{1}{N-N_{\text{par}}^{\text{lep}}+1}\left[{\displaystyle \sum _i^{\text{microw\hspace{0.17em}\hspace{0.17em}to\hspace{0.17em}\hspace{0.17em}UV},\sim GeV}\frac{{\left(F_i-{\overline{F}}_i\right)}^2}{{\sigma }_i^2}+{\displaystyle \sum _i^{rodio,X-ray,>20GeV}\frac{{\left(F_i-{\overline{F}}_i\right)}^2}{{\sigma }_i^2}\theta \left(F_i-\overline{F}\right),}}\right]$$(2)

where N is the number of data points to fit, $N_{\text{par}}^{\text{lep}}$ is the number of free parameters of the leptonic model and F_i and ${\overline{F}}_i$ are the predicted and observed photon flux values, respectively.

Having constrained the source geometry and the accelerated electrons, the only free variables left are those describing the accelerated protons. We fix the minimum proton Lorentz factor to ${\gamma }_{\text{p}}^{\min }={10}^2$ and the proton spectral index to α_P = 1.0. The choice of fixing ${\gamma }_{\text{p}}^{\min }$ is motivated by the fact that this parameter plays a minor role in the emission of secondary particles and radiation. In the case of α_P, this should play a role in the shape of the secondary photon emission; however, the peak proton energy is the leading factor in terms of neutrino emission.

We search the two remaining parameters, ${{\gamma }^{\prime }}_{\text{p}}^{\max }$ and ${L^{\prime }}_{\text{p}}$, using a grid scan of 35 ☓ 70 points, while fixing the remaining parameters to their best-fit values obtained in the previous step. The value ranges spanned by the grid scan are given in Table 1. As we can see, we scan values of ${10}^5<{{\gamma }^{\prime }}_{\text{p}}^{\max }<{10}^7$. The motivation for this search range is as follows: for proton Lorentz factors below this range, hadronic interactions should be inefficient in this model; above this range, the emitted neutrinos have energies that are too high compared to the sensitivities of current experiments (cf. Sect. 5).

In this second step, we wish to fit the X-ray and TeV gammaray data as well. The gamma rays emitted by proton interactions typically trigger electromagnetic cascades that lead to emission at lower frequencies, typically in the X-ray (Gao et al. 2019; Oikonomou et al. 2021; Rodrigues et al. 2021a; Sahakyan et al. 2022) or MeV gamma-ray band (Rodrigues et al. 2019; Reimer et al. 2019; Petropoulou et al. 2020a). The best fit should therefore minimize the chi-squared value including X-rays and TeV gamma rays: $${\chi }_{\text{LepHad}}^2=\frac{1}{N-N_{\text{par}}^{\text{LepHad}}+1}{\displaystyle \sum _i^{v>300GHz}\frac{{\left(F_i-{\overline{F}}_i\right)}^2}{{\sigma }_i^2},}$$(3)

where $N_{\text{par}}^{\text{LepHad}}=10$ is the number of free parameters of the full leptohadronic model.

Finally, we perform a local minimization in order to improve the fit (right panel of Fig. 3). Here, we consider all 10 best- fit leptohadronic parameter values and vary them continuously until the local minimum of ${\chi }_{\text{LepHad}}^2$ is reached. The minimization is performed using iminuit (Dembinski et al. 2020). While the goal of the previous steps was to constrain the vast parameter space of each source, the local minimization step is instead more confined, yielding a more precise result.

We consider the result with the minimum value of ${\chi }_{\text{LepHad}}^2$ to be the best-fit result, since it best describes the multiwavelength dataset relative to other results. At the same time, its absolute value is not a direct indicator of the overall quality of the fit. This is because, as discussed in Sect. 2.1, the datasets are not simultaneous and often contain points spanning a wide range of flux values within the same frequency band, reflecting the source’s intrinsic variability. For sources with high variability, the lowest possible value of ${\chi }_{\text{LepHad}}^2$ may therefore lie considerably above unity.

Once we find the best-fit result for each source following the criteria described above, we then vary the proton luminosity parameter ${L^{\prime }}_p$ in order to generate error estimates for the respective neutrino emission. We start with the best-fit value of ${L^{\prime }}_p$ and vary it both upward and downward in logarithmic increments of log $10\left({L^{\prime }}_p\right)=\pm 0.1$, until we reach a value ${\chi }_{\text{LepHad}}^2-{\chi }_{\text{LepHad }}^{2\text{ best }-\text{ fit }}=1$. This way we generate a 1_σ uncertainty region for ${L^{\prime }}_p$, which is the parameter to which the neutrino flux is most sensitive.

In the cases where the best fit is leptonic, that is, ${L^{\prime }}_p=10$ we cannot define a best-fit maximum proton Lorentz factor ${{\gamma }^{\prime }}_{\text{p}}^{\max }$. In those cases, we test the injection of protons with ${{\gamma }^{\prime }}_{\text{p}}^{\max }={10}^6$, which lies logarithmically in the middle of the range tested, and we determine the maximum value of ${L^{\prime }}_p$ for which ${\chi }_{\text{LepHad}}^2-{\chi }_{\text{LepHad}}^{2\text{best}-\text{fit}}\le 1$ by following the procedure described above. This represents the upper limit on the proton luminosity and the respective neutrino flux.

3 Results

3.1 Spectral fits

Using the methods described in the previous section, we self- consistently calculate the predicted multiwavelength and neutrino spectra from the sample. The best-fit parameter values of each source are listed in full in Table B.1 in descending order of the predicted neutrino flux. This includes the uncertainties on ${L^{\prime }}_p$ , obtained with the method described above, as well as the resulting uncertainties on the emitted neutrino flux $F_{v_{\mu }}$. We list first the sources for which the best-fit proton luminosity is incompatible with zero within 1_σ, then those where it is nonzero but compatible with zero within 1_σ, and finally those where only an upper limit can be established for the proton luminosity and the respective neutrino flux.

To support the discussion of the multiwavelength fits, we show in Fig. 4 the best-fit results for eight different sources. The total photon fluxes are shown as a thin black curve and are broken down into components of different physical origin, according to the legend.

In the two upper panels we can see on the left the results for the quasar 3C 273 and on the right for the HBL PKS 2155–304. In these cases, multiwavelength observations can be described purely by the emission from accelerated electrons, which is represented by the orange curve. Synchrotron emission from these primary electrons typically explains the nonthermal optical spectrum, together with thermal emission from the accretion disk and dust torus (gray curves). The Fermi-LAT fluxes are typically explained by a dominant leptonic component. For these two sources, that leptonic component alone can actually explain entirely these observations, without the need for a proton component. In the case of FSRQs, this is often dominated by external Compton scattering of the BLR photons, while for BL Lacs it is always dominated by synchrotron-self Compton (SSC) in this model. For a total of 218 sources, or 66% of the sample, leptonic emission can explain the fluxes across the multiwavelength spectrum. As indicated by the gray shaded area, observations below 300 GHz are considered as upper limits, as explained in Sect. 2.4; this is due to the fact that the compact region that we simulate is often optically thick to radio emission, which means these fluxes should originate in a larger zone that is not modeled here.

For these two sources, we can see that a subdominant component from proton interactions is allowed, but not required, to describe the data. In these cases, only an upper limit can be set for the proton luminosity, as determined by the 1_σ criterion explained in Sect. 2.4. This maximum contribution to the emission is shown as blue, green, yellow, and orange bands (see caption on the bottom of the figure), while the respective maximum neutrino emission is shown as a pink band at ~PeV energies.

In the case of the source PKS 2155-304 (upper right panel), we can see that at around 100 keV a potential contribution from secondaries coming from proton interactions (region shaded in blue) is tolerated within the 1_σ level, but also here it is not necessary to describe the data. We therefore can only set a constraint on the maximum proton loading of L_p/L_e < 10³ and a corresponding maximum neutrino flux level, shown as a pink band. It is also worth noting the extreme variability of this source, which results in a large dispersion of the multiwavelength fluxes. This situation becomes more common the brighter the source, since it can be more easily detected by the multiwavelength experiments in different activity states. In these cases, the best-fit model reflects only one of the possible activity states of the source.

We can see that for PKS 2155-304 the synchrotron emission from primary electrons undershoots the observed fluxes in the infrared, optical, and ultraviolet, which does not happen for a source like 3C 273. This is due to the use of a single chi-squared criterion to optimize the results over the entire multiwavelength spectrum, as described in Sect. 2.4. Since PKS 2155-304 displays large variability in the optical, ultraviolet, and X-ray ranges, the fit becomes less sensitive to those wavelengths compared to gamma rays, where the spectrum has a well-defined shape, small error bars, and a narrower overall spread in flux. The same applies to the source PKS 0735+17, shown in the right plot on the second row of Fig. 4, discussed in the following paragraph. In general, a more refined description of the low- energy bump would require a time-dependent analysis of the data using quasi-simultaneous datasets. Another possibility would be a more sophisticated algorithm that optimized simultaneously different individual features of the multiwavelength spectrum, such as the synchrotron peak frequency and Compton dominance. With the state-of-the-art methods, fitting a large number of sources requires a robust algorithm, and a more sophisticated fitting procedure lies therefore outside the scope of this work (see also Sect. 5 where other methodological limitations are discussed in depth).

Turning now to the second row of Fig. 4 (BL Lacs 3C 371 and PKS 0735+178), we see that in these cases Bethe-Heitler pair production by protons (pγ →pe⁺e^-, shown in blue) describes the X-ray fluxes, filling the “gap” between the two broadband features. These X-ray fluxes originate in synchrotron emission by the secondary pairs, while inverse Compton (the high-energy bump in the same blue curve) does not play a role in explaining observations in these two cases. The source PKS 0735+178 has recently been studied by Sahakyan et al. (2022), although in that case the authors used a dataset selected in the time domain. We compare these results in Appendix A.

We see that in both these sources the best-fit neutrino spectrum is incompatible with zero because the proton contribution dominated the emission in the X-ray band. This neutrino best fit is shown as a pink curve, accompanied by the respective error band from varying ${L^{\prime }}_p$ following the procedure discussed in Sect. 2.4. The accompanying gamma-ray emission (gray dotted curve visible in the case of 3C 371) is subsequently attenuated by the EBL and therefore does not appear in the final SED prediction, shown as a thin black curve.

In other cases, the proton contribution helps describe the data not only in the X-ray band but also in the gamma-ray band above 100 GeV. This is exemplified by the two BL Lacs on the third row of Fig. 4. In the left-hand plot, pairs from the Bethe-Heitler process help explain both the X-ray flux and a hardening of the Fermi-LAT spectrum. On the right-hand plot, pair production by high-energy gamma rays (green curve) also helps explain the X-ray emission as well as the high-energy part of the Fermi-LAT spectrum.

Overall, for a total of 62 sources, or 20% of the sample, the best-fit scenario involves proton emission to a level that is incompatible with zero within the 1a region, either because it is necessary to explain X-ray data or high-energy gamma rays.

Finally, on the last row of Fig. 4, we show two cases where the best-fit scenario includes a proton contribution but it is compatible with zero within the 1a region. As we can see in these two cases, pair production by protons and by the gamma rays emitted from photo-pion production helps explain the spectral shape in X-rays and gamma rays above ~10 GeV, but only marginally. The best-fit results of 44 sources fall in this category, which represents 14% of the sample. In this case we can still estimate a best-fit neutrino flux, but this estimate is not as as significant as for the group of sources discussed above. We therefore list these sources in Table B.1 separately from the above. We also do not take into account the best-fit proton luminosities when discussing the neutrino emission sample properties.

Fig. 4 Examples of best-fit results for 8 sources in the sample, given in the observer’s frame. The total modeled photon emission (including attenuation on the EBL) is shown as a thin black curve, while the color curves represent the different components of that emission, as detailed in the caption. We note that the individual SED components are given without the effect of EBL attenuation, while for the total flux the EBL attenuation is taken into account. Additionally, we show as gray curves the thermal emission from the dust torus and the accretion disk. On the right-hand side of each plot we show either the best-fit neutrino flux (pink curve) or the upper limit for the predicted emission (pink shade). The black points represent the multiwavelength flux observations that are fitted in each case, while the gray inverted triangles represent upper limits. Finally, the gray shaded area represents the frequency range for which data are not fitted but rather considered as upper limits to the model, as explained in Sect. 2.4. As we can see, the origin of the X-ray and VHE gamma-ray fluxes may in some cases be dominated by leptonic emission (upper row), but can also be dominated by electromagnetic cascades triggered by proton interactions. The corresponding plot for all modeled sources, as well as the numerical data, can be found in the online repository https://github.com/xrod/lephad-blazars.

3.2 Estimates on cosmic-ray loading and neutrino emission

We now turn to the estimated properties of the proton population and subsequent predictions for neutrino emission. In Fig. 5, we show the best-fit values of the physical luminosity, $L_{\text{e,p}}^{\text{phys}}={L^{\prime }}_{\text{e,p}}{\Gamma }^2/2$ on the left for electrons and on the right for protons, plotted against the disk luminosity of each source. For comparison, we show as a dotted line the relation $L_{}^{\text{phys}}=L_{\text{disk}}$ and as a solid line the best-fit power law relation for the entire sample. In the lower panels, the physical luminosity is plotted in units of the source’s Eddington luminosity L_Edd, estimated for each source based on M_BH .

As found by previous leptonic and leptohadronic modeling studies (e.g., Boettcher et al. 2013; Paliya et al. 2017), we can see that the physical luminosity in electrons loaded onto the jet follows closely the luminosity radiated by the disk (upper left plot). By fitting a power law we see that in this case it scales specifically with the square root of L_γ. In terms of the Eddington luminosity (lower left plot), this translates to a sub-Eddington electron injection rate of between 10^-3 and 10^-1 L_Edd for most sources.

In the right-hand panels of Fig. 5 we show the same values for protons. Here, we show as solid points only those results where the proton luminosity is incompatible with zero within 1a from the best fit. Best-fit results where the proton luminosity is either zero or compatible with zero within 1_σ are shown as inverted triangles. We see that the maximum allowed proton powers are almost exclusively above L_disk (upper right plot). In terms of each source’s Eddington luminosity, we can see in the lower right plot that there is an overall negative correlation of $L_{\text{p}}^{\text{phys }}/L_{\text{Edd}}~L_{\text{disk }}^{-0.3}$ Taking as an example a disk luminosity of L_disk = 10⁴⁵ erg s^-1, the model predicts an average proton luminosity of $L_{\text{p}}^{\text{phys }}=2.2L_{\text{Edd}}$. This average behavior was extrapolated for the subset of sources for which a nonzero proton contribution is necessary in this model (solid points). Most of these sources have L_γ > 10⁴⁵ erg s⁻¹ and have therefore, for the most part, proton powers $L_{\text{p}}^{\text{phys }}\lesssim L_{\text{Edd}}$ . In a few cases, the best fit lies above the Eddington limit, reaching up to $L_{\text{p}}^{\text{phys }}\sim {10}^3{L}_{\text{Edd}}$ in three sources. While jet loading at super-Eddington rates may be possible for short periods of flaring activity, it would be unrealistic to expect this during longer, steady quiescent states. A time-dependent study of the multiwavelength data could help one assert whether the datasets are selecting for flaring states in these cases.

In Fig. 6, we show the results for the ten sources with the highest predicted neutrino fluxes and for which the proton contribution is incompatible with zero. These correspond to the ten first sources in Table B.1. The first aspect worth noting is the variety of possible values of baryonic loading, defined here as ${L^{\prime }}_{\text{p}}/{L^{\prime }}_{\text{e}}$. In this subsample of ten sources, which all have comparable levels of neutrino emission, the baryonic loading ranges from 30 (PKS 2201+171) to 9 ☓ 104 (PKS 1954-388). This is due to the drastically different neutrino production efficiencies, which depend on the best-fit source geometry parameters, as discussed later in this section.

The second aspect is the different size of the 1a error bands. For example, in the cases of BL Lac GB6 J0654+5042 and FSRQ B2 2234+28A, the error band covers a region that seems to strongly violate X-ray data. This is due to the large amount of gamma-ray data points in these cases, which implies that the chi- squared value is dominated by deviations in the gamma-ray band, making the fit less sensitive to deviations in the optical and X-ray bands. Still, what we see in all cases is that the variation in the predicted flux of neutrinos generally maps onto a variation of the predicted flux in X-rays and gamma rays above ~ 10 GeV, due to the reasons discussed in the Sect. 3.1.

We now discuss the general trends of baryonic loading, neutrino production efficiency, and predicted neutrino luminosity. In the upper right panel of Fig. 7 we show the estimated baryonic loading of each source as a function of its luminosity L_γin the range from 100 MeV to 100 GeV. The solid line shows the best-fit power law considering only the sources for which the proton luminosity is incompatible with zero (solid points). We observe that the baryonic loading of the sample scales with L_γ^−0.6. For comparison, we show as a dashed curve the result of the phenomenological study by Palladino et al. (2019), who estimated the behavior of the baryonic loading necessary to explain the diffuse IceCube flux without violating constraints from lack of associations with high-luminosity sources. We see that the current model is compatible with that suggested by Palladino et al. (2019) for sources with L_γ ≈ 3 × 10⁴⁷ erg s⁻¹, with a baryonic loading using the predictions of 100. For lower luminosity sources that work predicts higher baryonic loading values, which would be necessary to reach the IceCube diffuse flux level. Since this work focuses on sources with L_γ ≲ 10⁴⁵ erg s⁻¹, we cannot extrapolate with confidence to that low-luminosity population. On the other hand, the work by Petropoulou et al. (2022) predicts a lower baryonic loading that scales somewhat similarly with L_γ. The difference in absolute value comes from a fundamental difference in the physical model: in this work we consider magnetic field strengths mainly below the Gauss level and the presence of external fields as the main driver of neutrino production; in comparison, Petropoulou et al. (2022) consider higher magnetic field strengths in the jet, assume magnetic reconnection as the main driver of cosmic ray acceleration, and do not model a specific source catalog.

In the upper right panel of Fig. 7, we show the modeled neutrino production efficiency of the population, which is shown to scale with <iqe/>. This trend is remarkably similar to that derived by Rodrigues et al. (2018). On the one hand this is unsurprising, since the underlying geometric model and treatment of external fields was similar to the one applied here (cf. also Murase et al. 2014). On the other hand, Rodrigues et al. (2018) did not focus on any specific source sample instead using only generic principles of blazar emission based on average population properties. That work also did not include a computation of the photons emitted by the cosmic-ray interactions, which opened the possibility that for high enough values of baryonic loading the high-energy emission might overshoot the assumed gamma-ray luminosity. Here, by simulating a specific sample with real data, we can estimate the effect of hadronic interactions on the SED, as discussed in Sect. 3.1.

These opposite trends of the baryonic loading, which scales as $\left(L_{\text{p}}^{\prime }/L_{\text{e}}^{\prime }\right)~L_{\gamma }^{-0.6}$, and the neutrino efficiency, which scales as $\left(L_v^{\prime }/L_{\text{p}}^{\prime }\right)~L_{\gamma }^{0.6}$, emerge independently as a collective behavior of the sample from each individual source fit. From this we can already infer the behavior of L_v with L_γ: $$L_{\gamma }^{\prime }~\left(\frac{L_{\text{p}}^{\prime }}{L_{\text{e}}^{\prime }}\right)\left(\frac{L_v^{\prime }}{L_{\text{p}}^{\prime }}\right)\left(\frac{L_{\text{e}}^{\prime }}{L_{\text{disk }}}\right)\left(\frac{L_{\text{disk }}}{L_{\gamma }}\right)L_{\gamma }=L_{\gamma }^{0.8},$$(4)

where the index 0.8 comes from introducing the relation $L_{\text{e}}^{\prime }~L_{\text{disk }}^{0.5}$, as shown in the upper left panel of Fig. 5, and the additional relation $L_{\gamma }~L_{\text{disk }}^{0.7}$ found for this sample (as can be derived directly from the values of L_disk and F_γ provided in Table B.1).

This relation between observed gamma-ray luminosity and predicted neutrino luminosity is demonstrated in the lower left plot of Fig. 7, where we show the predicted neutrino fluxes from the sample. For comparison, we show as a cross the case of TXS 0506+056 based on previous modeling results, since it is not contained in this sample (see Sect. 1). Other recent neutrino candidates that are in the sample are highlighted with black circles (see also Appendix A where we contrast the results with the literature on these sources). Finally, on the lower right plot of Fig. 7 we show that $L_v~L_{\text{MeV}}^{0.6}$, where L_MeV is the flux of gamma rays in the range between 0.1 and 100 MeV that is predicted by the model. Although we do not currently have data on that energy range, the model predicts that the neutrino luminosity scales with the emitted luminosity in MeV gamma rays, which strengthens the case for future MeV gamma-ray missions as neutrino-counterpart search machines (e.g., De Angelis et al. 2021; Caputo et al. 2022).

Fig. 5 Best-fit values of the physical luminosity injected in electrons (left) and protons (right), as a function of each source’s disk luminosity. The physical luminosities in the rest frame of the black hole are given by $L_{\text{e,p,}B}^{\text{phys}}={\Gamma }_{\text{b}}^2{L^{\prime }}_{\text{e,p,}B}/2$ In the upper plots the physical luminosity is given in erg s-1 and the dotted lines indicate Lphys=Ldisk ; in the lower panels it is given in units of the Eddington luminosity LEdd and the dotted lines indicate Lphys = LEdd . The best-fit power laws (excluding the upper limits) are shown as solid black lines. For protons, these relations are calculated using only the best-fit results that are incompatible with ${L^{\prime }}_p=0$. These are shown as solid points, while the inverted triangles represent results compatible with ${L^{\prime }}_p=0$within 1 σ. We also cite the result of the Pearson correlation coefficient r, applied to the logarithm of the corresponding variables, and its respective p value. As we can see, all variables present a statistically meaningful moderate power-law correlation with the disk luminosity, except the proton power in the upper right plot.

Fig. 6 Best-fit results of the 10 sources in the sample with the highest predicted neutrino fluxes. The plots are ordered by the value of the total predicted neutrino flux, as in Table B.1, here from left to right and top to bottom. The blue curves show the multiwavelength emission predicted by the full leptohadronic model. The best-fit baryonic loading value is provided in the caption and can also be found in Table B.1). The orange curves on the right-hand side of each panel show the best-fit emitted all-flavor neutrino fluxes, as discussed in Sect. 3.2. The corresponding plot for all modeled sources, as well as the numerical data, can be found in the online repository https://github.com/xrod/lephad-blazars.

Fig. 7 Characterization of the best-fit baryonic loading, neutrino efficiency and predicted neutrino fluxes. Upper panels: distribution of the best- fit baryonic loading (left) and neutrino production efficiency (right) as a function of the source’s gamma-ray luminosity between 100 MeV and 100 GeV. For comparison, in the left-hand-side plots we show the behavior predicted by the magnetically loaded jet model by Petropoulou et al. (2022) and the phenomenological study by Palladino et al. (2019). Lower panels: predicted muon neutrino luminosity for each source as a function of the Fermi LAT luminosity (left) and of the predicted MeV gamma-ray luminosity (right). The color code follows the previous pictures in that FSRQs are shown in green and BL Lacs in purple. The cases where a purely leptonic model lies within 1a of the best-fit result are shown as inverted triangles. The points therefore represent sources for which the data favors leptohadronic model. The solid black lines show the best-fit power law relations. For reference, we show four recent IceCube candidate sources. Of these, only TXS 0506+056 is not contained in this sample, for which we adopt the values from previous studies. For the three other blazars marked, the values of neutrino and MeV gamma-ray luminosity are those predicted by this work.

Fig. 8 Predicted muon neutrino flux from each source. Here we show only the cases where the neutrino flux is incompatible with zero within the 1σ region from the best fit. Previous neutrino candidates that are part of the sample are highlighted with a black circle and TXS 0506+056, which is not in the sample, is shown as a black cross for comparison. In both plots, we show in blue the IceCube point source sensitivity (Aartsen et al. 2017a); on the right plot, we show additionally the sensitivity band for the future IceCube-Gen2 (Aartsen et al. 2021) and the names of the associated sources with the highest predicted fluxes (cf. Fig. 6).

4 Neutrino flux predictions and implications for multi-messenger searches

4.1 Expected neutrino fluxes

Using the predictions for the neutrino luminosity emitted by each source, we can deduce the corresponding flux of muon neutrinos at Earth and compare it to the sensitivities of present and future experiments. In Fig. 8, we show the predicted total muon neutrino flux observed at Earth as a function of each source’s Fermi-LAT flux. In both plots, we show in blue the IceCube sensitivity based on the 7-yr point source analysis (Aartsen et al. 2017a), which is plotted as a band due to its dependence on declination.

In the right-hand plot, we zoom into the highest fluxes and show additionally the names of the associated source for the objects with highest flux. We also show the sensitivity band of the future IceCube-Gen2 experiment, approximated as 3-5 times better compared to the current IceCube (Aartsen et al. 2021). As we can see, the model suggests that seven sources in the sample may be within reach of the current IceCube, while more than twenty would be observable by the future IceCube-Gen2. Interestingly, both PKS 0735+178 and GB6 J040+0617, both of which have an associated high-energy IceCube alert (see Sahakyan et al. 2022; Garrappa et al. 2019), appear here at a subthreshold level for detection by IceCube. This apparent discrepancy might be due, for example, to occasional flaring activity of these sources, with a corresponding temporary increase in neutrino fluxes that is not grasped by this model due to the lack of time dependence of the data (cf. discussion in Appendix A, in particular the upper right panel of Fig. A.1 where this point is clarified for PKS 0735+178).

Fig. 9 Characterization of the neutrino spectra predicted by the model. Left: distribution of the predicted peak neutrino energy, given in the observer’s frame. Right: distribution of the predicted neutrino energy in the observer’s frame as a function of the source’s gamma-ray luminosity Lγ. The best-fit power law is shown as a black line. This best-fit function includes only the results (shown as dots) and not the upper limits (shown as inverted triangles).

4.2 Predicted neutrino event rates

Additionally to the total neutrino flux, the number of expected neutrino detections depends on the energy of the emitted neutrinos and on the declination of each source. The former is summarized in Fig. 9. On the left-hand panel, we show the distribution of the predicted peak energy of the neutrino energy flux spectrum observed at Earth, $E_v^{\text{obs2 }}\text{d}{F}_v \text{d}{E}_v^{\text{obs }}$. We can see that FSRQs have a broader distribution of peak energies that are centered around $E_v^{\text{obs }}=1\text{PeV}$ but extend up to 10 PeV. BL Lacs, on the other hand, have a slightly lower dispersion and are centered around 630 TeV. Compared to the distribution of maximum proton energies in the sample (cf. lower right panel of Fig. B.1), these distributions are much more peaked. This is due to the fact that the neutrino spectrum emerges from a convolution between the spectrum of cosmic-ray protons and the target ambient photons. In the case of FSRQs, those target photons can be the synchrotron photons in the jet or external photons from the BLR; this broader variety is the reason for the broader distribution of peak energies compared to BL Lacs, for which external photons are not available as targets in this model.

In the right-hand panel of Fig. 9, we show the distribution as a function of gamma-ray luminosity, showing a positive correlation of $L_v^{\text{obs}}~L_{\gamma }^{0.2}$, as expected from the arguments discussed in Sect. 3.2. We note that in both the Gaussian fit on the left and the power-law fit on the right, we account for all sources that have a best-fit result involving protons, even those compatible with zero baryonic loading at the 1σ level. This is because given the physical properties of the source, which are determined by the best-fit estimate, the peak energy of the emitted neutrinos is determined regardless of the actual baryonic loading. Additionally to these, we show as inverted triangles the sources for which the best fit is leptonic.

Using the information on the neutrino peak energy and the declination of each source, we can now estimate the number of events in IceCube. For each source, we adopt the declination- and energy-dependent effective area from the 7-yr point source analysis (Aartsen et al. 2017a) and convolve it with the predicted muon neutrino spectrum from the source, taking into account its specific declination.

The results of that analysis are shown in the left panel of Fig. 10, Like in Fig. 8, we only show those results where the baryonic loading is incompatible with zero within the 1σ region from the best fit. As we can see, for a few sources the rate of events in IceCube approaches one per decade. These low rates reflect an estimate based on the average emission from the sources, while potential periods of hadronic flaring should lead to temporary neutrino outbursts that may explain sporadic associations.

We then perform the same analysis for IceCube-Gen2. In that case, we consider the number of events in IceCube predicted for each source and scale it up by the ratio between the geometric effective areas of IceCube-Gen2 and the current IceCube (Aartsen et al. 2021). Given the different geometric setup of IceCube-Gen2, the relative improvement in effective area compared to IceCube is declination-dependent and ranges from a factor of 3 at the horizon up to a maximum factor of about 7 (cf. dashed and dotted curves in shown in Fig. 24 by Aartsen et al. 2021). This up-scaling therefore depends on the declination of each source.

These results are shown in the right panel of Fig. 10. We can see that several sources have a rate close or superior to about an event per decade, which suggests that IceCube-Gen2 will be able to resolve this subdominant population, even ignoring occasional flaring events.

Despite the positive correlation between gamma-ray and neutrino luminosity shown in Fig. 7, we see that there is a considerable number of relatively dim blazars with $L_{\gamma }\lesssim {10}^{-11}{\text{ergcm}}^{-2} {\text{s}}^{-1}$ erg cm⁻² s⁻¹ for which the model suggests a likely future detection.

Fig. 10 Number of events per year in IceCube (left) and in the future IceCube-Gen2 (right) predicted by the model as function of each source’s gamma-ray flux. Here we include only those sources for which the best-fit result is incompatible with zero baryonic loading within the 1σ region. The name of each associated source is shown for reference (cf. Table B.1). Green points denote FSRQs and purple points denote BL Lacs.

Fig. 11 All-flavor neutrino flux expected from the entire sample (solid black curve), separated into FSRQs (solid purple curve) and BL Lacs (solid green curve). Following the sampling method described in the main text, we then extrapolate the result to the entire Fermi-LAT blazar sample and obtain the diffuse spectrum shown as a dotted curve. For visual reference, the individual contributions are represented as thin curves (see also Table B.1). The black data points show IceCube data from the analysis by Aartsen et al. (2015a) and the blue curve shows the 9 yr IceCube sensitivity (Aartsen et al. 2018). In pink, we show the stacking limits derived for blazars by Aartsen et al. (2017b).

4.3 Diffuse neutrino flux estimates

Finally, we combine the neutrino flux predictions from the individual sources in our sample to estimate their contribution to the IceCube diffuse flux. This total contribution is shown as solid curves in Fig. 11, in green for FSRQs and purple for BL Lacs. The sum of these two components is shown as a black solid curve. As we can see, in this model, our sample contributes to the diffuse flux at a level of about 5%, reaching a peak flux of about 9 × 10⁻¹⁰ GeV cm⁻² s⁻¹ sr⁻¹ at energies of about 1 PeV. This diffuse flux extends, however, from 100 TeV to about 10 PeV. This follows in part from the range of proton energies tested in this model (cf. Table 1), since in these sources the neutrino energy follows closely that of the primary cosmic rays. We cannot exclude the possibility that cosmic rays are accelerated to ultra-high energies at least in some blazars, since this scenario is not tested in this model (see for example the discussion in Rodrigues et al. 2021b, and the connection to the ultra-high-energy cosmic ray issue), leading potentially to a diffuse neutrino flux at energies higher than dozens of PeV.

We then estimate the potential contribution of the entire population of Fermi-LAT blazars by extrapolating these results. Since our sample does not exactly follow the luminosity distribution of the Fermi-LAT population (cf. Fig. 1), rather than simply multiplying the flux with a source count factor, we estimate this flux by performing Monte Carlo sampling of 1) values of (L_γ, z) from the luminosity function by Ajello et al. (2012, 2014), as shown in Fig. 1, and considering only those generated sources whose flux lies above the Fermi-LAT sensitivity; 2) values of $L_{v_{\mu }}$ from the distribution of $L_{v_{\mu }}\left(L_{\gamma }\right)$ predicted by the model (lower left panel of Fig. 7); 3) $E_v^{\text{peak}}$ from the distribution shown in the histogram of Fig. 9; 4) a neutrino spectral shape from the pool of the model results, depending on $L_{v_{\mu }}$. All these sampling steps are done independently for FSRQs and BL Lacs. We then sum over the spectra obtained and multiply the total spectrum by a factor given by the fraction of sources in the sample for which the result is incompatible with a leptonic solution, which is 23% (8%) for FSRQs (BL Lacs). This incorporates into the result the fact that most sources in the generated population should also be compatible with a purely leptonic description and should therefore not contribute to the diffuse neutrino flux.

The result of this extrapolation is shown as a dotted curve in Fig. 11. Compared to the sample itself, we can see that it is broader in energy, which arises from the fact that the distribution with L_γ of Fermi-LAT sources population (green and purple curves along the vertical axes of the bottom panels of Fig. 1) are slightly broader compared to our sample (black curves).

As we can see, the model predicts a diffuse flux from the Fermi-LAT population at the level of 20% of the observed IceCube flux. This is in agreement with current limits set by the IceCube collaboration for gamma-ray blazars through stacking analyses (Aartsen et al. 2017b). We can see this by comparing the dotted curve with the IceCube stacking limit, shown as a pink band with a downward arrow.

This supports a scenario where the bulk of the diffuse IceCube flux currently detected originates in a population of more abundant sources, potentially undetected in gamma rays. The sporadic associations we observe between blazars and high-energy IceCube events may then arise from hadronic flares with temporarily enhanced neutrino emission.

At the same time, these results support the possibility that continued observations, as well as more sensitive future experiments, should help detect an increasingly significant diffuse contribution from gamma-ray blazars. In particular, it is worth noting that the diffuse neutrino flux predicted by the model has a considerable component up to hundreds of PeV dominated mainly by FSRQs, a regime that can be more easily probed by future experiments like IceCube-Gen2.

5 Limitations

Finally, we discuss a few caveats of our method and their implications on the interpretation of the results. The two first limitations arise from the fact that the multiwavelength datasets are not time-selected. Firstly, this leads to the fact that some of the datasets contain points with a vast dispersion in flux within the same wavelength band due to variability. This is especially the case for high-luminosity sources (cf. e.g., upper right plot of Fig. 4). On the one hand, regardless of the variability in the data, the χ² minimization method guarantees that the best fit describes an “average” activity state, which in accordance with the purpose of the study. On the other hand, is is difficult to quantify how this affects the uncertainties on the neutrino predictions in Table B.1; firstly, because the flux variability is different in different wavelength bands; secondly, because the potential contribution from the hadronic emission does not scale in a linear fashion across wavelengths. For example, if a source displays gamma-ray fluxes ten times higher than those captured by the best-fit SED (as we can see in Fig. 4 for blazar PKS 2155-304), this does not necessarily mean we would expect a neutrino flux ten times higher during such a flare (cf. e.g., Oikonomou et al. 2019; Rodrigues et al. 2021a; Gasparyan et al. 2021).

The second important implication of this lack of temporal selection of the data is that is is possible that in some cases the data in different bands represents different activity states of the source. For example, Fermi-LAT data are obtained from regular surveying over several years and is therefore more likely to represent the different states of activity of a source; in contrast, an instrument like the Swift-BAT is triggered more likely by flares. This issue could be improved, for example, by taking a dataset from pointing observations, which is typically the case for X-rays, and then selecting quasi-simultaneous observations form surveying experiments in other wavelengths. However, this would considerably limit the amount of multiwavelength data available to the point where a fit would be challenging for most sources.

The third limitation in the method is the lack of a dedicated analysis of the flux variability in the different wavelength bands. This is of particular importance for the infrared, optical, and ultraviolet, since in those bands the emission can be dominated either by thermal emission from the disk and the torus, or by nonthermal emission from the jet. In most cases, our fitting method is able to discriminate between these two scenarios based purely on the spectral shape. For example, in the fit to source OJ 535 (lower right plot of Fig. 4) we can see that the thermal components can fit well the data in this band, while in other cases the spectral shape is better fitted with a synchrotron spectrum. However, in cases like the quasar 3C 273 (upper left panel of Fig. 4), the spectral shape is ambiguous. In our model, the microwave fluxes are fit with a synchrotron component, which should lead to some variability in that band; however, the fluxes actually display little variability, within a factor of up to a few (e.g., Fernandes et al. 2020), which seems to favor instead thermal emission from the core. Resolving these ambiguous cases would require a more detailed analysis of the light curves, which would be prohibitive for such a large catalog with the methods employed here. It is therefore important to keep this ambiguity in mind when interpreting the results, since in some cases a different fit could be possible.

Another point worth mentioning is that we treat all BL Lac sources with the same model, assuming that the source does not contain a BLR. Given that two distinct IceCube candidates have been determined to be masquerading BL Lacs (Padovani et al. 2019; Sahakyan et al. 2022), it is likely that this source subclass is more representative in terms of neutrino emission. It is a fact that most of the BL Lacs in our sample are LBLs, for which it is less likely that the jet emission could hide their broad line emission and therefore more likely that they are true BL Lac objects; however, the sample also includes 22 IBLs and 17 HBLs. For those sources, a masquerading BL Lac model including a BLR could eventually also be tested, likely leading to higher neutrino fluxes due to the additional target photons available for hadronic interactions. In that sense, the current model can be seen as a conservative approach where neutrino production in IBLs and HBLs is not enhanced by any external photon fields. It is particularly important to keep this fact in mind when interpreting the extrapolation shown in Fig. 11.

Another important aspect to mention is the vastness of the parameter space of the one-zone model. First, this means that the range of parameters tested (cf. Table 1) is limited. For example, as discussed in Sect. 4.3 in the context of the diffuse flux predictions, the ultra-high-energy cosmic ray scenario is not tested in this work, since the maximum proton energy that was tested was 10 PeV. Another important factor is the magnetic field strength, which we limit to the 10 G level. As we know, in theoretical scenarios where the high-energy emission originates in proton synchrotron emission, the magnetic fields required are typically higher (Liodakis & Petropoulou 2020), as are scenarios where the cosmic rays are accelerated through magnetic reconnection in magnetically dominated jets (Petropoulou et al. 2022). In that sense, the predictions made by this work are strictly limited to the paradigm provided by the conventional one-zone leptohadronic framework.

The large size of the parameter space of the model also leads to degenerate solutions with comparable levels of goodness-of-fit, making it challenging to explore the parameter space in a complete way. The use of a genetic algorithm followed by a local minimization provides a robust and highly parallelizable method for finding a minimum in a highly non-linear space of many parameters; however, this method comes with two important caveats:

• A)

  It does not guarantee the uniqueness of the result. In fact, as we discuss in Appendix A, other works that have previously modeled sources in this sample have described multiwavelength data with different parameters, leading to different predictions on neutrino emission. At the same time, this is often also due to different multiwavelength datasets, obtained through different selections in the time domain, or, in some cases, to differences in numerical simulation method itself. We therefore do not claim completeness in the way we explored the parameter space of the model, nor the uniqueness of our results within the leptohadronic framework.

• B)

  As discussed in Sect. 2.4, we start by scanning for solutions where the X-ray fluxes are not necessarily described by purely leptonic emission but the gamma-ray peak is (cf. Fig. 3). This derives directly from the hypothesis that hadronic cascades may in some cases dominate the X-ray emission, which is informed by previous models in the literature, as discussed in Sect. 1. In that sense, this criterion may be seen as a prior underlying the parameter search. That means that the solutions obtained for which hadronic cascades dominate the X-ray emission (see e.g., second and third rows in Fig. 4) should be interpreted as being within the leptohadronic hypothesis and purely leptonic solutions may also be possible. Paliya et al. (2017) have in fact provided purely leptonic solutions for the entire sample; however, their model does not treat radiative energy losses in a time-dependent and self-consistent way as is done here.

Finally, in order to limit the search to a more manageable parameter space we have also fixed two parameters of the model, as described in Sect. 2.2: the jet viewing angle and the spectral index of the proton population.

Fixing the viewing angle to θ_obs = 1/Γ_b is motivated by the effect of relativistic beaming, which compresses the majority of the jet’s emission into a solid angle of $1/{\Gamma }_{\text{b}}^2$. In the BL Lac model, the Doppler and Lorentz factors are in fact degenerate, since they impact only the energy and flux transformation from the rest frame of the jet into the observer's frame. However, it is important to note two caveats: firstly, in FSRQs the external fields transform into the jet frame with the Lorentz factor Γ_b, while the transformation of the jet emission into the observer’s frame depends on both Γ_b and θ_obs. These two parameters are therefore not degenerate when external fields are present, which means that by fixing the observing angle we are effectively excluding a portion of the parameter space of the jet geometry. Secondly, the effect of the relativistic beaming on the boosting of the external fields into the jet frame and of the neutrino fluxes into the observer’s frame should dramatically increase for jets observed on-axis (Boettcher 2023). Given that our approach does not include this possibility, it can be described as conservative regarding both the effects of external fields on neutrino production and electromagnetic cascades, as well as the predicted levels of neutrino emission.

Regarding the spectral index of the accelerated proton spectrum, this was fixed to α_p = 1. This somewhat aggressive choice means that a relatively large fraction of the accelerated protons will be injected into the system at high energies, leading to the development of electromagnetic cascades and efficient neutrino emission. A softer injection spectrum, such as α_p = 2 as expected from first-order Fermi acceleration, would require a higher overall proton injection power in order to produce the same photon flux from proton interactions. This would then result in higher estimates of the best-fit baryonic loading in those sources for which proton-triggered cascades dominate the X-ray flux. At the same time, an energy-dependent escape mechanism such as Bohm-like diffusion, which was not explored in this work, can lead to a harder steady-state proton spectrum compared to the injected one (see e.g., Fig. 8 of Rodrigues et al. 2018). With such an escape mechanism, a soft accelerated proton spectrum would lead to a harder steady-state spectrum, assuming the numerical system is evolved long enough to achieve that steady state.

6 Conclusion

We have performed a source-by-source analysis of a sample of 324 blazars from the CGRaBS catalog, all detected in gigaelectronvolt gamma rays. Of these, 237 are FSRQs. From the physical perspective, we have adopted an external field model, which has been successfully deployed in several studies in recent years to describe individual associations between IceCube events and individual blazars. The modeling was performed through numerical, self-consistent, cosmic-ray simulations where we calculated the time-dependent interactions of cosmic-ray electrons and protons in the relativistic jet. We estimated the best-fit model parameters by fitting the predicted multiwavelength emission to a set of public data from each source. We then predicted the emitted neutrino spectrum, estimated self-consistently from the same simulation.

We have described the multiwavelength observations of the 324 blazars from the infrared up to gamma rays and derived the best-fit parameters. For 106 sources, or 33% of the sample, radiative emission from proton interactions either dominates the observed fluxes in X-rays or at least improves the fit compared to a purely leptonic scenario. In some cases, the description of high-energy gamma rays above tens of gigaelectronvolt also benefits from a hadronic contribution. For the remaining 218 sources, we have constrained the leptonic model parameters and established upper limits on the baryonic loading of the jet and the emitted neutrino spectrum. Overall, we can say that the leptohadronic paradigm, which was previously applied to individual IceCube candidate blazars, has shown to provide a solution for the entire sample. At the same time, we cannot guarantee the uniqueness of these solutions due to the vast parameter space of the model.

For the sources that have a nonzero proton contribution, we have estimated the best-fit baryonic loading, neutrino production efficiency, emitted neutrino flux, and peak energy of the emitted neutrinos and we have studied possible correlations. We have found that the required proton powers are mostly sub-Eddington, with a few exceptions that require super-Eddington accretion. We have also demonstrated that the best-fit baryonic loading has a decreasing trend as a function of the blazar’s gamma-ray luminosity. This conclusion is consistent with phenomenological expectations, given the stacking limits on the brightest sources derived by IceCube. Importantly, this conclusion was obtained here from a purely physics-driven approach through source-by-source modeling without these phenomenological considerations.

On average, the modeled neutrino luminosity has shown to depend positively on the source’s gamma-ray luminosity in the Fermi-LAT range, with a power-law relation $L_v~L_{\gamma }^{0.8}$, as well as on the predicted luminosity in megaelectronvolt gamma rays. Having tested maximum energy values of the primary protons in the range from 100 TeV up to 10 PeV, we deduced neutrino spectra that peak mostly around 1 PeV for FSRQs, but with a considerable number of sources peaking as high as 10 PeV. For BL Lacs, the neutrino spectra peak mostly around 630 TeV. For about 30 sources in the sample, we predict a neutrino flux above 10⁻¹³ erg cm⁻² s⁻¹, which lie above the sensitivity of the future IceCube-Gen2 experiment. Our results predict that within the first decade of operation, IceCube-Gen2 should have detected the ten brightest sources in this sample, even considering only an average quiescent state emission. If we add the effect of sporadic hadronic flares, the actual number of associations should be considerably higher.

We also conclude that this sample may contribute to the diffuse IceCube flux at the ~5% level. Furthermore, by extrapolating to the population of Fermi-LAT sources, we derived a contribution of 20% of the diffuse IceCube flux. This contribution is in agreement with upper limits derived in stacking analyses (Aartsen et al. 2017b). While IceCube currently seems to detect sporadic neutrino outbursts from hadronic blazar flares, our result suggests that next-generation neutrino experiments such as IceCube-Gen2 can potentially probe the steady-state emission from this elusive population of cosmic-ray accelerators.

Acknowledgements

We thank P. Padovani, P. Giommi, and S. Kiehlmann for insightful comments on the manuscript. We acknowledge funding from the German Science Foundation DFG, via the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source To Signal”.

Appendix A Comparison with recent studies of individual sources

Multiple sources in this sample have also been independently modeled with numerical, time-dependent frameworks in previous studies. In this section we compare our results with some of those in the literature, a comparison that is summarized in Fig. A.1.

On the upper row we see two results by Boettcher et al. (2013), who applied leptonic and leptohadronic models to a group of six FSRQs, four LBLs and two IBLs. Of these 12 sources, nine are also included in the sample studied in this work. On the left we see the fits to the IBL object W Comae. The dashed magenta curve shows the best-fit SED by Boettcher et al. (2013) and the solid blue curve the one from this work. The datasets that are fitted are shown as magenta and black data points, respectively. As we can see, the authors fit a time-selected dataset that captures a state of relatively enhanced activity compared to the solution found in this work, which corresponds to a quiescent state. The high gamma-ray fluxes above 10 GeV are better described by the model by Boettcher et al. (2013). It is possible that the low state in optical captured by our model should correspond to a lower state in high-energy gamma rays, which is the case in our best-fit result. However, further analysis of the quiescent state data would be necessary in order to confirm this difference in the gamma-ray peak between the two models.

In the right plot we see the two results of the modeling of S5 0716+71. In this case, the synchrotron peaks predicted by the two models correspond again to different activity states, but the gamma-ray emission is at a comparable level. In both these cases the current model predicts a lower hadronic contribution, with only an upper limit on the proton luminosity. In comparison, the results by Boettcher et al. (2013) suggest a large proton loading of ~ 10⁴, although the peak of the gamma-ray emission is dominated by inverse Compton scattering like in our case.

In the left panel on the second row of Fig. A.1 we compare our result for the HBL object PKS 2155-304 with that by Petropoulou (2014). This is a good example of a strongly variable source from optical, X-rays, and gamma rays. As we can see, the gamma-ray data points reach as high as 10⁻⁹ erg cm⁻² s⁻¹. In our model, as was the case for W Comae, we see that our result reproduces a state of both low synchrotron flux and low gamma-ray flux. On the other hand, Petropoulou (2014) have modeled quasi-simultaneous data form the source (magenta points) and have shown that those fluxes cannot be reproduced with a single-zone leptonic model. Instead, it requires either two distinct emission zones or a proton population co-accelerated with the electrons, which is the case shown here as a dashed magenta curve. In our approach we do not distinguish between data from different epochs (see discussion in Sect. 5), so our method cannot produce this kind of statement regarding any specific epoch. Instead, the only statement that can be derived from our method for this source is that the SED shown in blue is the best fit to the overall dataset from the source (in black) that can be achieved with a one-zone model, which is achieved without any necessary proton component.

On the same row, we show on the right the result of our leptohadronic modeling of the IBL object AP Librae as well as the result by Petropoulou et al. (2017). With their model, Petropoulou et al. (2017) have shown that TeV emission from AP Librae can be well described with a component from hadronic interactions. As we can see by comparing the dashed magenta and blue curves in the plot, our result supports this finding. The best-fit proton populations have a higher luminosity in our case, but a similar maximum energy within a factor of 5. Although the predicted neutrino spectrum is not shown explicitly by Petropoulou et al. (2017), the predicted flux should lie within the same order of magnitude of that predicted by this model, given the level of cascade emission seen in the right plot of Fig. 1 of that work.

The major difference in terms of multiwavelength predictions lies in the nonthermal optical emission. This is because, as shown also in Fig. 6 (e.g., for PKS 0426-380 and B2 2234+28A), the best-fit SED corresponds to a low state in optical emission. However, the predicted synchrotron peak frequency is compatible with the observed one and in this case also with that predicted by Petropoulou et al. (2017), at v_syn ≈ 3 × 10¹⁴ Hz. This means that in spite of the different flux level, the model still captures the IBL nature of the source.

In the third row of Fig. A.1, we see on the left the photon spectrum of PKS 1424-41 predicted by Gao et al. (2017), together with the respective neutrino flux as a dashed green curve. This was in the context of the detection of a high-energy IceCube event in 2016 (Kadler et al. 2016) from the direction of the source. In terms of the multiwavelength fluxes, the first difference between the two models lies in the flux level of the synchrotron emission, since the search algorithm used in this work has selected a low state compared to that by Gao et al. (2017). This is because 1) we treat radio fluxes as upper limits, while in this case they can in fact be described by a single-zone model as shown by Gao et al. (2017); 2) the large variability in the ultraviolet range of the dataset used in this work, which in this case disfavors solutions with high synchrotron fluxes. In contrast, in the previous study the data was time-selected, as represented by the magenta points, which means that a particular activity state is targeted. The second major difference is the large neutrino flux predicted by Gao et al. (2017), while here the best-fit result is compatible with a purely leptonic solution. This is because 1) Gao et al. (2017) were explaining an IceCube association, for which high neutrino fluxes are necessary; and 2) the same study derived a baryonic loading value of 10⁷, which lies considerably above the range tested in this work (cf. Table 1).

In the case of GB6 J1040+0617, spatially coincident with an IceCube alert in 2014 (Garrappa et al. 2019), we see that the neutrino flux predicted by Banik et al. (2020) has a drastically different spectral shape compared to this work (third row of Fig. A.1, right panel). This difference is due to the fact that that model includes proton interactions with high-density clouds surrounding the relativistic jet. This leads to a soft neutrino spectrum typical of proton-proton models and a gamma-ray emission that is harder at high energies compared to that found in this work.

In the left plot on the fourth row of Fig. A.1 we compare our result for the FSRQ PKS 1502+106 with that found by (Oikonomou et al. 2021) (see also Rodrigues et al. 2021a, where the same model used in this work is applied to time-selected datasets from this source). We can see that in spite of a baryonic loading only a factor of a few higher, Oikonomou et al. (2021) predict a much stronger neutrino emission, which is possible due to protons with energies above 10 PeV, which interact efficiently with infrared emission from a dust torus. In comparison, as listed in Table B.1, we set a limit on the neutrino flux of at most 3 × 10⁻¹⁶ erg cm⁻²s⁻¹, a factor of 10⁴ lower. Two aspects are worth noting in this case: firstly, a solution like that obtained by Oikonomou et al. (2021) should in principle be possible in the current framework and within the parameter set explored in this work, with a jet emission zone outside the BLR but within the radiation field of the dust torus. The reason why the genetic algorithm converged on this parameter set may be because it provides a larger fraction of parameter space with acceptable chi-squared values, therefore gathering more solutions over time and being therefore statistically preferred. The second aspect is that in the solution by Oikonomou et al. (2021) the emission that originates in hadronic interactions is actually subdominant in the X-ray range compared to the leptonic component and dominates the observations only for gamma rays above tens of GeV, near the cut off of the SED (cf. Fig. 8 of that reference). In that sense, this proton component may easily be neglected by a large-scale parameter search like the one used in the current work, when no requirements on neutrino emission are enforced.

Fig. A.1 Comparison between the results obtained in this work and ten selected results from recent multiwavelength studies, in chronological order of publication (Boettcher et al. 2013; Petropoulou 2014; Gao et al. 2017; Petropoulou et al. 2017; Banik et al. 2020; Oikonomou et al. 2021; Gasparyan et al. 2021; Sahakyan et al. 2022, 2023). In blue and orange we show respectively the best-fit photon and neutrino fluxes from the current model; the dashed magenta and green curves show the photon and neutrino spectra (when available) predicted by previous works. The magenta points represent the data fitted in the respective study, which in some cases varies significantly from the current work because of time-domain selection. The best-fit baryonic loading values are given in each caption.

On the right we show the result for 3C 279, an extremely variable quasar present in our sample that was also modeled recently by Gasparyan et al. (2021). Comparing the dashed green curve with the orange shaded area, we see that the neutrino fluxes predicted by the two models are vastly different. The first factor for this difference is the fact that Gasparyan et al. (2021) model a state of extremely high gamma ray fluxes corresponding to a flaring event in 2015. This allows for strong gamma-ray emission from efficient hadronic interactions, compared to the dataset modeled in this work. The second interesting fact is that this extremely high neutrino flux is achieved with a baryon loading of $L_{\text{p}}^{\prime }/L_{\text{e}}^{\prime }=300$, while our best-fit result has a higher value of $L_{\text{p}}^{\prime }/L_{\text{e}}^{\prime }=500$ and yet places only an upper limit on neutrino emission. This is due to the different nature of the model by Gasparyan et al. (2021), particularly the high values of maximum proton energy (about 200 PeV) and magnetic field strength (70 gauss), both of which are above our search range. This illustrates the point discussed in Sect. 5 on the existence of alternative scenarios that lie outside the parameter space considered in this work that can lead to different predictions in the neutrino sector.

In the case of PKS 0735+178 (lower left plot), Sahakyan et al. (2022) describe multiwavelength observations surrounding the detection of a high-energy IceCube event from the direction of the source in December 2021, during which the source was undergoing a multiwavelength flare of unprecedented magnitude. In contrast, the dataset modeled in this work is not time-selected and is an order of magnitude lower in flux. Additionally to this difference, the authors also use a masquerading BL Lac model, based on observations that suggest that in spite of its classification as BL Lac, the source actually possesses a BLR whose emission is outshone by the synchrotron peak (see also Padovani et al. 2019). In contrast, in this work we do not consider the effect of a BLR for BL Lacs, which leads to a slightly lower Compton dominance and a different spectral shape of the predicted X-ray and gamma-ray spectra. Furthermore, in the current model the only target photons for photo-pion production are those from synchrotron emission; compared to an external BLR field, these synchrotron photons have lower densities, leading to a slightly lower ratio L_v/L_γ, and lower frequencies, leading to a higher energy of the interacting protons and the emitted neutrinos.

Finally, in the recent study by Sahakyan et al. (2023), the authors describe time-selected multiwavelength datasets from the source PKS 0537 by with a purely leptonic model. One of their fits to multiwavelength dataset is shown in the lower right panel of Fig. A.1, compared to that obtained in this work. As we can see, the X-ray and gamma-ray spectral shapes are comparable, but the additional component from hadronic interactions can potentially help describe the slightly harder gamma-ray fluxes at the highest energies. At the same time, as we can see in the blue shaded band, the SED by Sahakyan et al. (2023) is consistent with the result that we find in the purely leptonic limit, which is contained within the 1σ uncertainty band, shown in blue.

Appendix B Best-fit parameters

In Fig. B.1 we show the distribution of the best-fit parameters for FSRQs (green) and BL Lacs (purple). The parameters are those listed in Table 1.

In the case of R_diss, only FSRQs are included in the histogram since the BL Lac model is not sensitive to this parameter. The last four histograms relate to hadronic variables; for these cases, we only include the sources whose best fit has a nonzero baryonic loading.

In Table B.1 we provide the best-fit parameter values of the leptohadronic model for the entire sample. Additionally, for each source we also provide the predicted energy-integrated muon neutrino flux, $F_{v_{\mu }}$, the predicted number of events per year in IceCube, $N_{v_{\mu }}$, and the peak energy of the neutrino spectrum in the observer’s frame, $E_v^{\text{peak}}$.

As explained in the main text, the errors in the neutrino flux $F_{v_{\mu }}$ were obtained by varying the proton injection luminosity until the change in the emission violates the observations by 1σ compared to the best-fit result. The corresponding uncertainties in $N_{v_{\mu }}\text{ and }{L}_{\text{p}}^{\prime }$ can be extrapolated from those in $F_{v_{\mu }}$, since these variables are directly proportional to each other.

The sources in Table B.1 are given in descending order of the predicted neutrino flux, starting with those for which the flux is incompatible with zero within the uncertainty. Then follow those for which the best-fit neutrino flux is nonzero but compatible with zero within the uncertainties, also provided in descending order of flux. Finally, we list those sources for which the best-fit solution is purely leptonic. These are sorted by the upper limit on the neutrino flux.

Fig. B.1 Binned distribution of the best-fit parameter values for the entire sample. In green we show the results for the FSRQs in the sample and in purple the BL Lacs. The best-fit parameter values are also listed in Table B.1, as well as in the online repository https://github.com/xrod/lephad-blazars.

References

All Tables

All Figures

Fig. 1 Characterization of the sample. Upper panel: sky map with the positions of the sources studied in this work, which include 237 FSRQs, 88 BL Lacs, three narrow-line Seyferts, and three objects classified as BCUs in the Fermi catalogs. Lower panels: distribution of FSRQs (left) and BL Lacs (right) in the sample as a function of redshift and gamma-ray luminosity (dark points). For comparison, we show as lighter points to the distribution of the respective overall population, based on the luminosity functions by Ajello et al. (2012, 2014). The gray area marks the subthreshold region for the Fermi-LAT. The black distributions alongside the margins refer to the sample; the green (left) and purple (right) distributions refer to the population of sources detected by the Fermi-LAT (i.e., excluding the putative sources below the sensitivity threshold).

In the text

Fig. 2 Schematic representation of two possible geometries of the one-zone FSRQ model. The broad line region is represented by the purple clouds. In the case of the blob represented on the left, its proximity to the BLR implies a large relative Doppler factor, making the external fields appear highly boosted in the jet frame (compare left-hand plot, where the external fields are shown in the rest frame of the central engine, with the middle plot, where they are shown as energy density spectra in the jet frame). In the case represented on the right, due to the high dissipation radius Rdiss the BLR photons impinge more from behind and the external fields from the BLR are less boosted, as shown in the lower right plot. The relativistic boosting of these external radiation fields follows the treatment by Ghisellini & Tavecchio (2009). We note that on the left plot, the broad lines are isotropic in the black hole frame, while the disk emission is anisotropic except for the ~1% that gets isotropized through Thomson scattering. This explains why in the middle plot the broad lines get more highly boosted than the disk emission.

In the text

Fig. 3 Flowchart of the fitting procedure. The leptohadronic fit of the data from each source is performed as a three-step process. The initial step, where the optical and gamma ray broadband features are explained with a leptonic model, requires the most simulations because the vast parameter space is initially unconstrained. The number of CPU-hours quoted refers to the entire sample and is based on the code’s performance of 1.6 s for a single leptonic simulation and ~20 s for a leptohadronic simulation.

In the text

Fig. 4 Examples of best-fit results for 8 sources in the sample, given in the observer’s frame. The total modeled photon emission (including attenuation on the EBL) is shown as a thin black curve, while the color curves represent the different components of that emission, as detailed in the caption. We note that the individual SED components are given without the effect of EBL attenuation, while for the total flux the EBL attenuation is taken into account. Additionally, we show as gray curves the thermal emission from the dust torus and the accretion disk. On the right-hand side of each plot we show either the best-fit neutrino flux (pink curve) or the upper limit for the predicted emission (pink shade). The black points represent the multiwavelength flux observations that are fitted in each case, while the gray inverted triangles represent upper limits. Finally, the gray shaded area represents the frequency range for which data are not fitted but rather considered as upper limits to the model, as explained in Sect. 2.4. As we can see, the origin of the X-ray and VHE gamma-ray fluxes may in some cases be dominated by leptonic emission (upper row), but can also be dominated by electromagnetic cascades triggered by proton interactions. The corresponding plot for all modeled sources, as well as the numerical data, can be found in the online repository https://github.com/xrod/lephad-blazars.

In the text

Fig. 5 Best-fit values of the physical luminosity injected in electrons (left) and protons (right), as a function of each source’s disk luminosity. The physical luminosities in the rest frame of the black hole are given by $L_{\text{e,p,}B}^{\text{phys}}={\Gamma }_{\text{b}}^2{L^{\prime }}_{\text{e,p,}B}/2$ In the upper plots the physical luminosity is given in erg s-1 and the dotted lines indicate Lphys=Ldisk ; in the lower panels it is given in units of the Eddington luminosity LEdd and the dotted lines indicate Lphys = LEdd . The best-fit power laws (excluding the upper limits) are shown as solid black lines. For protons, these relations are calculated using only the best-fit results that are incompatible with ${L^{\prime }}_p=0$. These are shown as solid points, while the inverted triangles represent results compatible with ${L^{\prime }}_p=0$within 1 σ. We also cite the result of the Pearson correlation coefficient r, applied to the logarithm of the corresponding variables, and its respective p value. As we can see, all variables present a statistically meaningful moderate power-law correlation with the disk luminosity, except the proton power in the upper right plot.

In the text

Fig. 6 Best-fit results of the 10 sources in the sample with the highest predicted neutrino fluxes. The plots are ordered by the value of the total predicted neutrino flux, as in Table B.1, here from left to right and top to bottom. The blue curves show the multiwavelength emission predicted by the full leptohadronic model. The best-fit baryonic loading value is provided in the caption and can also be found in Table B.1). The orange curves on the right-hand side of each panel show the best-fit emitted all-flavor neutrino fluxes, as discussed in Sect. 3.2. The corresponding plot for all modeled sources, as well as the numerical data, can be found in the online repository https://github.com/xrod/lephad-blazars.

In the text

Fig. 7 Characterization of the best-fit baryonic loading, neutrino efficiency and predicted neutrino fluxes. Upper panels: distribution of the best- fit baryonic loading (left) and neutrino production efficiency (right) as a function of the source’s gamma-ray luminosity between 100 MeV and 100 GeV. For comparison, in the left-hand-side plots we show the behavior predicted by the magnetically loaded jet model by Petropoulou et al. (2022) and the phenomenological study by Palladino et al. (2019). Lower panels: predicted muon neutrino luminosity for each source as a function of the Fermi LAT luminosity (left) and of the predicted MeV gamma-ray luminosity (right). The color code follows the previous pictures in that FSRQs are shown in green and BL Lacs in purple. The cases where a purely leptonic model lies within 1a of the best-fit result are shown as inverted triangles. The points therefore represent sources for which the data favors leptohadronic model. The solid black lines show the best-fit power law relations. For reference, we show four recent IceCube candidate sources. Of these, only TXS 0506+056 is not contained in this sample, for which we adopt the values from previous studies. For the three other blazars marked, the values of neutrino and MeV gamma-ray luminosity are those predicted by this work.

In the text

Fig. 8 Predicted muon neutrino flux from each source. Here we show only the cases where the neutrino flux is incompatible with zero within the 1σ region from the best fit. Previous neutrino candidates that are part of the sample are highlighted with a black circle and TXS 0506+056, which is not in the sample, is shown as a black cross for comparison. In both plots, we show in blue the IceCube point source sensitivity (Aartsen et al. 2017a); on the right plot, we show additionally the sensitivity band for the future IceCube-Gen2 (Aartsen et al. 2021) and the names of the associated sources with the highest predicted fluxes (cf. Fig. 6).

In the text

Fig. 9 Characterization of the neutrino spectra predicted by the model. Left: distribution of the predicted peak neutrino energy, given in the observer’s frame. Right: distribution of the predicted neutrino energy in the observer’s frame as a function of the source’s gamma-ray luminosity Lγ. The best-fit power law is shown as a black line. This best-fit function includes only the results (shown as dots) and not the upper limits (shown as inverted triangles).

In the text

Fig. 10 Number of events per year in IceCube (left) and in the future IceCube-Gen2 (right) predicted by the model as function of each source’s gamma-ray flux. Here we include only those sources for which the best-fit result is incompatible with zero baryonic loading within the 1σ region. The name of each associated source is shown for reference (cf. Table B.1). Green points denote FSRQs and purple points denote BL Lacs.

In the text

Fig. 11 All-flavor neutrino flux expected from the entire sample (solid black curve), separated into FSRQs (solid purple curve) and BL Lacs (solid green curve). Following the sampling method described in the main text, we then extrapolate the result to the entire Fermi-LAT blazar sample and obtain the diffuse spectrum shown as a dotted curve. For visual reference, the individual contributions are represented as thin curves (see also Table B.1). The black data points show IceCube data from the analysis by Aartsen et al. (2015a) and the blue curve shows the 9 yr IceCube sensitivity (Aartsen et al. 2018). In pink, we show the stacking limits derived for blazars by Aartsen et al. (2017b).

In the text

Fig. A.1 Comparison between the results obtained in this work and ten selected results from recent multiwavelength studies, in chronological order of publication (Boettcher et al. 2013; Petropoulou 2014; Gao et al. 2017; Petropoulou et al. 2017; Banik et al. 2020; Oikonomou et al. 2021; Gasparyan et al. 2021; Sahakyan et al. 2022, 2023). In blue and orange we show respectively the best-fit photon and neutrino fluxes from the current model; the dashed magenta and green curves show the photon and neutrino spectra (when available) predicted by previous works. The magenta points represent the data fitted in the respective study, which in some cases varies significantly from the current work because of time-domain selection. The best-fit baryonic loading values are given in each caption.

In the text

	Fig. B.1 Binned distribution of the best-fit parameter values for the entire sample. In green we show the results for the FSRQs in the sample and in purple the BL Lacs. The best-fit parameter values are also listed in Table B.1, as well as in the online repository https://github.com/xrod/lephad-blazars.
In the text