© S. Das et al. 2022

1. Introduction

Blazars are a subclass of radio-loud active galactic nuclei (AGNs) with a highly relativistic jet collimated toward the observer’s line of sight. They are considered to be prominent candidates for the origin of IceCube-detected diffuse astrophysical neutrino flux beyond ∼10 TeV (IceCube Collaboration 2013; Eichler 1979; Sikora et al. 1987; Petropoulou et al. 2015; Murase et al. 2018; Yuan et al. 2020) and may also contribute in the PeV-EeV energy range (Kalashev et al. 2013; Kochocki et al. 2021; Das et al. 2021). For the first time in September 2017, a high-energy muon-neutrino event IC-170922A (E_ν ∼ 0.3 PeV) was associated with the γ-ray flaring blazar TXS 0506+056 at 3σ significance (IceCube Collaboration 2018a). Subsequently, other events having positional coincidence with Fermi-LAT detected blazars are also observed with lower statistical significance (Garrappa et al. 2019; Franckowiak et al. 2020). Although there are several studies on this object, it is crucial to revisit the spectral properties in order to predict the multi-messenger signals from similar sources.

The explanation of the neutrino event requires synchrotron (SYN) and synchrotron self-Compton (SSC) photons as the target for pγ interactions (Gao et al. 2019; Cerruti et al. 2019; Sahu et al. 2020). Whereas, other models require an external photon field, resulting in external inverse-Compton (IC) emission (Reimer et al. 2019; Keivani et al. 2018; Rodrigues et al. 2019; Petropoulou et al. 2020). In the hadronuclear interpretation via pp interaction, the shock-accelerated protons may interact with gas clouds in the vicinity of the supermassive black hole (Liu et al. 2019) or cold protons in the jet (Banik & Bhadra 2019). Some studies invoke neutrino production from the interaction of relativistic neutron beams in the jet, originating in pγ interactions with external photons (see, e.g., Zhang et al. 2020). In Fraija et al. (2020), pγ interactions occur with seed photons produced by the annihilation of electron-positron pairs from the accretion disk.

The cosmic-ray-induced cascade from Bethe-Heitler (BH) pair production contributes near the X-ray energies in the pγ scenario, thus also limiting the hadronic component in GeV-TeV γ-rays. As a result, this constrains the astrophysical neutrinos, in many cases, to a flux level lower than that predicted by IceCube observations. Thus, an additional photon field of energy ϵ ≃ m_πm_pc⁴/20E_ν ≃ 440 eV (i.e., in the UV to soft X-ray energy band) is required. However, an “orphan” neutrino flare from this source from September 2014 to March 2015 is revealed from the analysis of archival data, at 3.5σ statistical significance (IceCube Collaboration 2018b) with 13 ± 5 signal events above the atmospheric background. This flare was not accompanied by increased activity in γ-rays, indicating that different astrophysical processes may dominate for neutrino and γ-ray flares. Often, a two-zone model is employed, considering a high opacity for GeV γ-rays in the neutrino production region (Sahakyan 2018; Xue et al. 2019, 2021).

The MAGIC collaboration recently modeled the spectrum of TXS 0506+056, observed during a multiwavelength campaign lasting 16 months from November 2017 to February 2019, covering the radio band, optical/UV, high-energy, and very high-energy (VHE, E >  100 GeV) γ-rays (Acciari et al. 2022). A γ-ray flaring activity was observed by MAGIC during December 2018. Fermi-LAT observed several short flares on timescales of days to weeks, unlike the long-term flare of 2017. At lower energies, no significant variability was observed. The observed flare was not associated with any neutrino event. Their model infers the neutrino luminosity to be lower than the detection threshold of currently operating instruments.

We analyzed the multiwavelength spectral energy distribution (SED) using a one-zone leptohadronic model. The low-energy peak results from the SYN radiation of relativistic electrons. The high-energy peak is produced by SSC and IC scattering of external photons, also called external Compton (EC) emission. These external photons may originate from the broad-line region (BLR). Although broad-line emission in the optical spectrum is not detected, TXS 0506+056 can be a masquerading BL Lac, which is intrinsically a flat spectrum radio quasar (FSRQ) with hidden broad lines and a standard accretion disk (Padovani et al. 2019). The interaction of the cosmic-ray protons with the leptonic radiation and the external photon field produces a characteristic photon spectrum resulting from the electromagnetic cascade of secondary electrons, which is constrained by the X-ray data.

Blazars are also suitable candidates for ultra-high-energy cosmic-ray (UHECR; E ≳ 10¹⁸ eV) acceleration. They can escape the jet and interact with cosmic background photons. In an earlier work, the neutrino flux originating in extragalactic propagation of UHECRs from TXS 0506+056, and a few other blazars was analyzed (Das et al. 2022). We assumed a correlation between the cosmic-ray luminosity and IceCube-detected neutrino luminosity inside the jet to scale the cosmogenic fluxes. Here we constrain the UHECR luminosity by SED modeling, consistently with the allowed flux of cosmogenic γ-rays. The extragalactic magnetic field (EGMF) is crucial for determining this line-of-sight resolved γ-ray component. Using the latest spectrum data, we coherently explain the multiwavelength SED, and predict the corresponding neutrino flux from the jet emission region and the plausibility of cosmogenic γ-ray contribution to the SED. Finally, the luminosity constraint can be used to predict the resulting cosmogenic neutrino flux at EeV energies.

The observed multiwavelength SED is difficult to resolve into components coming from interactions inside the jet or line-of-sight resolved UHECR interactions or multiple zones, for example. Hence all photons are treated on equal footing by the detectors. Thus, if UHECRs escape from the source, they produce cosmogenic gamma rays and the line-of-sight resolved component of that flux can contribute at the VHE range of the multiwavelength SED. We invoke three components: (i) purely leptonic emission; (ii) hadronic emission inside the jet; and (iii) a line-of-sight component of the hadronic emission during extragalactic propagation. This is the reason the cosmogenic gamma-ray spectrum is also used to fit the observed SED.

Our study is essentially a one-zone model, with all the jet parameters constrained by the multiwavelength SED. These jet parameters are used to calculate the luminosity in cosmic rays required inside the jet for SED and neutrino modeling. However, the parameter space is degenerate and hence adjusted to maximize the neutrino production inside the jet. Using the same normalization of the proton spectrum required for this luminosity, we also calculated the luminosity in escaping UHECRs. Hence, it is essentially the same proton spectrum, with the same normalization. However, we assume at energies ≳0.1 EeV that they escape the source because the observed SED does not allow for UHECR interactions inside the source and the simultaneous explanation of quiescent state neutrino flux.

We present the methods of our one-zone modeling in Sect. 2. In Sect. 3.1 we present the results for leptohadronic emissions inside the jet. In Sect. 3.2 we show the contribution of line-of-sight resolved cosmogenic γ-ray contribution to the SED and the subsequent constraints on the EGMF strength. We discuss our results and draw our conclusions in Sect. 4.

2. Radiative modeling

We consider the emission region in the jet to be a spherical blob of radius R′, consisting of a relativistic plasma of electrons and protons moving through a uniform magnetic field B. The bulk Lorentz factor of the jet is Γ and the doppler factor is given by δ_D = [Γ(1 − βcosθ)]⁻¹, where βc is the velocity of the emitting plasma and θ is the viewing angle. For θ ≲ 1/Γ, Γ ≈ δ_D. We inject electrons in the blob with a spectrum

$$\begin{array}{c}\hfill Q_{\mathrm{e}}^{\prime }({\gamma }_{\mathrm{e}}^{\prime })=A_{\mathrm{e}}{({\gamma }_{\mathrm{e}}^{\prime }/{\gamma }_0)}^{-\alpha -\beta {log}_{10}({\gamma }_{\mathrm{e}}^{\prime }/{\gamma }_0)}\text{for}{\gamma }_{\mathrm{e},\min }^{\prime }<{\gamma }_{\mathrm{e}}^{\prime }<{\gamma }_{\mathrm{e},\max }^{\prime }\end{array}$$(1)

to fit the observed broadband SED. The normalization of the spectrum A_e depends on the luminosity of injected electrons, and γ₀m_ec² is a reference energy fixed at 500 MeV. A quasi-steady state is reached when the injection is balanced by radiative cooling and/or escape. Empirically, the steady-state electron density distribution is given as $N_{\mathrm{e}}^{\mathrm{\prime }}({\gamma }_{\mathrm{e}}^{\mathrm{\prime }})=Q_{\mathrm{e}}^{\mathrm{\prime }}({\gamma }_{\mathrm{e}}^{\mathrm{\prime }})t_{\mathrm{e}}^{\mathrm{\prime }}$, where $t_{\mathrm{e}}^{\mathrm{\prime }}=\mathrm{m}\mathrm{i}\mathrm{n}\{t_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }},t_{\mathrm{e}\mathrm{s}\mathrm{c}}^{\mathrm{\prime }}\}$. We consider the escape timescale $t_{\mathrm{e}\mathrm{s}\mathrm{c}}^{\mathrm{\prime }}\simeq t_{\mathrm{d}\mathrm{y}\mathrm{n}}^{\mathrm{\prime }}=2R^{\mathrm{\prime }}/c$. The radiative cooling timescale is given as

$$\begin{array}{c}\hfill t_{\mathrm{cool}}^{\prime }={\displaystyle \frac{3m_{\mathrm{e}}c}{4(u_{\mathrm{B}}^{\prime }+{\kappa }_{\mathrm{KN}}{u}_{\mathrm{ph}}^{\prime }){\sigma }_T{\gamma }_{\mathrm{e}}^{\prime }}},\end{array}$$(2)

where $u_{\mathrm{B}}^{\mathrm{\prime }}=B^2/8\pi$ is the energy density in magnetic field, $u_{\mathrm{p}\mathrm{h}}^{\mathrm{\prime }}$ is the energy density of soft photons, σ_T is the Thomson scattering cross section, and κ_KN accounts for the suppression of IC emission due to the Klein-Nishina effect. We use the open-source code GAMERA to solve the transport equation for obtaining the injection spectrum at time t′ (Hahn 2016)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial N_{\mathrm{e}}^{\prime }}{\partial t}}=Q_{\mathrm{e}}^{\prime }({\gamma }_{\mathrm{e}}^{\prime },t^{\prime })-{\displaystyle \frac{\partial }{\partial {\gamma }_{\mathrm{e}}^{\prime }}}(b{N}_{\mathrm{e}}^{\prime })-{\displaystyle \frac{N_{\mathrm{e}}^{\prime }}{t_{\mathrm{esc}}}},\end{array}$$(3)

where $b=b({\gamma }_{\mathrm{e}}^{\mathrm{\prime }},t^{\mathrm{\prime }})$ is the energy loss rate of electrons.

The steady-state electron spectrum yields the SYN and SSC emission. In addition, we consider an external photon field, which is Compton upscattered by the same electrons. It is considered to be a blackbody with temperature T′ and energy density $u_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{\prime }}=(4/3){\mathrm{\Gamma }}^2{u}_{\mathrm{e}\mathrm{x}\mathrm{t}}$ in the jet frame, where the energy density in the AGN frame is $u_{\mathrm{e}\mathrm{x}\mathrm{t}}={\eta }_{\mathrm{e}\mathrm{x}\mathrm{t}}{L}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}/4\pi R_{\mathrm{e}\mathrm{x}\mathrm{t}}^{2}c$ and η_ext is the fraction of the disk luminosity. Here R_ext is the radius of the region containing the external photons. The emission blob is assumed to be at this distance along the axis of the jet. These photons can enter the relativistic jet and become doppler boosted in the comoving frame.

The steady-state proton injection spectrum is given by a power law $N_{\mathrm{p}}^{\mathrm{\prime }}({\gamma }_{\mathrm{p}}^{\mathrm{\prime }})=A_{\mathrm{p}}{\gamma }_{\mathrm{p}}^{\mathrm{\prime }-\alpha }$. The main energy loss processes of the protons are pion production (pγ → p + π⁰ or n + π⁺) and BH process (pγ → p + e⁺e⁻). The seed photons are the leptonic emission and external photons. The charged pions decay to produce neutrinos. The timescale of these interactions can be expressed as

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{t_{\mathrm{p}\gamma }^{\prime }}}={\displaystyle \frac{c}{2{\gamma }_{\mathrm{p}}^{\prime 2}}}{\displaystyle {\int }_{{ϵ}_{\mathrm{th}}/2{\gamma }_{\mathrm{p}}}^{\infty }}\mathrm{d}{ϵ}_{\gamma }^{\prime }{\displaystyle \frac{n({ϵ}_{\gamma }^{\prime })}{{ϵ}_{\gamma }^{\prime 2}}}{\displaystyle {\int }_{{ϵ}_{\mathrm{th}}}^{2ϵ{\gamma }_{\mathrm{p}}}}\mathrm{d}{ϵ}_r\sigma ({ϵ}_r)K({ϵ}_r){ϵ}_r,\end{array}$$(4)

where σ(ϵ_r) and K(ϵ_r) are respectively the cross section and inelasticity of photopion production or BH pair production as a function of photon energy ϵ_r in the proton rest frame, and $n({ϵ}_{\gamma }^{\mathrm{\prime }})$ is the target photon number density (Stecker 1968; Chodorowski et al. 1992; Mücke et al. 2000). The interaction timescale of protons inside the jet is many orders of magnitude higher than the dynamical timescale, below tens of PeV energies. Hence, in order to increase the efficiency of pγ interactions required for appreciable neutrino production, it is compelling to ignore their escape if the Eddington luminosity budget is maintained. Otherwise, the production of the same neutrino flux will require higher kinetic power in protons if the escape rate is comparable to the pγ interaction rate. Hence, an escape timescale greater than the pγ interaction timescale is assumed. The normalization A_p of the proton spectrum is calculated from the luminosity requirement arising from the in-jet hadronic contribution to the SED.

The spectrum of γ-rays from the decay of neutral pions and the spectrum of e⁺e⁻ due to the BH process are calculated using the parameterization by Kelner & Aharonian (2008). When calculating the pion decay gamma rays and electron spectra, the input proton spectrum is weighted by the rate of the corresponding process; for example, in the case of electron spectrum from the BH process, we inject $N_{\mathrm{p}}^{\mathrm{\prime }}({\gamma }_{\mathrm{p}}^{\mathrm{\prime }})\ast R_{\mathrm{B}\mathrm{H}}/R_{\mathrm{t}\mathrm{o}\mathrm{t}}$, where R_BH is the BH interaction rate (calculated using Eq. (4)). The quantity R_tot is the total interaction rate considering photopion production and BH pair production. The high-energy γ-rays are absorbed by γγ → e^± pair production with the leptonic radiation and also with the external blackbody radiation, leading to the attenuation of TeV γ-rays. The escaping γ-ray flux is given as

$$\begin{array}{c}\hfill Q_{\gamma ,\mathrm{e}sc}^{\prime }({ϵ}_{\gamma }^{\prime })=Q_{\gamma ,\pi }^{\prime }({ϵ}_{\gamma }^{\prime })({\displaystyle \frac{1-exp(-{\tau }_{\gamma \gamma })}{{\tau }_{\gamma \gamma }}}).\end{array}$$(5)

We calculate τ_γγ using the formalism given by Gould & Schréder (1967) to calculate the absorption probability per unit path length for an isotropic photon field

$$\begin{array}{c}\hfill l_{\gamma \gamma }^{-1}({ϵ}_{\gamma }^{\prime })={\displaystyle \frac{1}{2}}{\displaystyle \int \int n}({ϵ}_k^{\prime }){\sigma }_{\gamma \gamma }({ϵ}_{\gamma }^{\prime },{ϵ}_k^{\prime },\theta )(1-cos\theta )sin\theta \mathrm{d}\theta \mathrm{d}{ϵ}_k^{\prime },\end{array}$$(6)

where σ_γγ is the full pair-production cross section and the $n({ϵ}_k^{\mathrm{\prime }})$ is the combined density of soft photons and external radiation.

The high-energy electrons and positrons produced in γγ pair production ($Q_{\mathrm{e},\gamma \gamma }^{\mathrm{\prime }}$), charged pion decay ($Q_{\mathrm{e},\pi }^{\mathrm{\prime }}$), and BH process (Q′_e,BH) can initiate cascade radiation from the jet. We solve the steady-state spectrum of secondary electrons $N_{\mathrm{e},\mathrm{s}}^{\mathrm{\prime }}({\gamma }_{\mathrm{e}})$ in the jet frame using the analytical approach of Boettcher et al. (2013), including ${Q\prime }_{\mathrm{e},\mathrm{BH}}$ in the source term and the escape term to be the same as primary electrons. In a synchrotron-dominated cascade, emission from secondary electrons is given by

$$\begin{array}{c}\hfill Q_{\mathrm{s}}^{\prime }({ϵ}_{\mathrm{s}}^{\prime })=A_0{ϵ}_{\mathrm{s}}^{\prime -3/2}{\displaystyle {\int }_1^{\infty }}\mathrm{d}{\gamma }_{\mathrm{e}}^{\prime }{N}_{\mathrm{e},\mathrm{s}}^{\prime }({\gamma }_{\mathrm{e}}^{\prime }){\gamma }_{\mathrm{e}}^{\prime -2/3}{e}^{-{ϵ}_{\mathrm{s}}^{\prime }/b{\gamma }_{\mathrm{e}}^{\prime 2}}\end{array}$$(7)

with A₀ = cσ_TB′²/[6πm_ec²Γ(4/3)b^4/3] being a normalization constant, where b = B′/B_crit and B_crit = 4.4 × 10¹³ G.

Our results in Sect. 3.1 suggest that the proton spectrum inside the source is required to be cut off beyond a specific energy to explain the multiwavelength SED. The resulting neutrino flux is thus limited by this value of E′_p,max. We model the protons to escape beyond this energy if accelerated inside their source. An energy-independent escape timescale of the order ∼R/c is sufficient for the escape to dominate over photohadronic interactions inside the jet. However, considering a diffusion faster than ∝E¹ leads to negligible interaction efficiency beyond tens of PeV energies in the comoving jet frame. This is reasonable when a quasi-ballistic propagation is assumed instead of diffusive propagation inside the jet emission region.

The resulting muon neutrino flux from pγ interactions is calculated as

$$\begin{array}{c}\hfill E_{\nu }^2{J}_{\nu }={\displaystyle \frac{1}{3}}{\displaystyle \frac{V^{\prime }{\delta }_D^2{\mathrm{\Gamma }}^2}{4\pi d_L^2}}{E}_{\nu }^{\prime 2}{Q}_{\nu ,p\gamma }^{\prime },\end{array}$$(8)

where the factor 1/3 corresponds to neutrino oscillation and $Q_{\nu ,p\gamma }^{\mathrm{\prime }}$ is the total electron and muon neutrino flux from charged pion decay in the comoving frame.

3. Results

3.1. Leptohadronic emission inside the jet

During the multiwavelength campaign, the source was monitored using the Swift-XRT, Swift-UVOT, and NuSTAR, maximizing the simultaneity of observation with the MAGIC telescope (Acciari et al. 2022). From November 2017 to February 2019, a total of ∼79 h of good-quality data were collected by MAGIC. During most of this period (∼74 h), the source was found to be in a low state with average photon flux F(> 90 GeV)=(2.7 ± 2.1) × 10⁻¹¹ erg cm⁻² s⁻¹. Fermi-LAT data was also used in the spectral analysis. For the flaring state of December 2018, only simultaneous data is obtained by Fermi-LAT and MAGIC in the GeV and VHE γ-rays, and by ASAS-SN in the optical. The integral photon flux observed by MAGIC rose by an order of magnitude compared to the low state. The most significant variability was observed in the GeV band, while the X-ray variability was found to be at a lower level. The radio, optical, and UV bands showed moderate variability (Acciari et al. 2022).

The γ-ray variability timescale of TXS 0506+056 observed in October 2017 was shown to be t_var ≤ 10⁵ s (Keivani et al. 2018). The γ-ray flare of December 2018 was found to be very similar. The size of the emission region inferred from the variability is R′≲δ_Dct_var/(1 + z)≃6.75 × 10¹⁶(δ_D/30)(t_var/10⁵s) cm. We assume the radius of the emission region to be 10¹⁶ cm and Γ ≈ δ_D during the high and the low states. The value of the magnetic field was fine-tuned by fitting the optical and gamma-ray data. It is also assumed to be the same in the two states, B′=0.28 G. The muons and pions produced in hadronic interactions do not suffer significant energy losses before decaying, for this value of B′.

The luminosity distance of TXS 0506+056 is d_L  ≈  1837 Mpc, with a redshift of z = 0.3365. The total kinetic power of the jet in the AGN frame is calculated as $L_{\mathrm{k}\mathrm{i}\mathrm{n}}=L_{\mathrm{e}}+L_{\mathrm{p}}+L_{\mathrm{B}}=\pi R^{\mathrm{\prime }2}{\mathrm{\Gamma }}^{2}c(u_{\mathrm{e}}^{\mathrm{\prime }}+u_{\mathrm{p}}^{\mathrm{\prime }}+u_{\mathrm{B}}^{\mathrm{\prime }})$, where $u_{\mathrm{e}}^{\mathrm{\prime }}$, $u_{\mathrm{p}}^{\mathrm{\prime }}$, and $u_{\mathrm{B}}^{\mathrm{\prime }}$ are the energy densities of electrons, protons, and magnetic field, respectively. The maximum electron energy changes from ${\gamma }_{\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}=2\times {10}^4$ in the low state to ${\gamma }_{\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}=5\times {10}^4$ in the high state to account for the spectral variability. We vary the maximum proton energy ($E_{\mathrm{p},\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{\prime }}$) in the comoving jet frame over a wide range to find the best-fit value of 6.3 PeV, fixed for both the low and high states. The cascade emission from the steady-state secondary electron spectrum $N_{\mathrm{e},\mathrm{s}}^{\mathrm{\prime }}$ is shown by the red lines in Fig. 1. The low-energy peak of the cascade emission originates from the secondary emission of e^± pair produced in BH process, which is severely constrained by the X-ray data. As a result, the pion decay cascade at higher energies is also limited, and the contribution to the high-energy peak is not significant.

Fig. 1. Multiwavelength SED of TXS 0506+056 during the 2017–2019 campaign. Left: low-state (average) flux during the observation period. Black data points show the Fermi-LAT average spectrum and MAGIC upper limits, and Swift and NuSTAR spectra on 16 October 2018. Gray data points show the whole range of optical to X-ray spectra. Right: high-state observation in VHE γ-rays during December 2018. Black data points show the Fermi-LAT, MAGIC, and ASAS-SN contemporaneous spectra. Gray data points show the nearest observations in radio, optical, UV, and X-rays on 8 December 2018. The green, blue, and purple curves respectively correspond to SYN, SSC, and EC-BLR emissions. The red curve is the cascade emission from secondary electrons. The magenta curve is the predicted neutrino flux in both cases. The black solid and dashed PeV data points represent the IceCube flux upper limit for one detection in 0.5 yr and 7.5 yr respectively for the 2017 detection event. The gray line is the cosmogenic γ-ray flux from the extragalactic propagation of UHE protons. The orange and blue dashed lines are CTA and LHAASO point source sensitivity. See text for details.

We obtain T′=2 × 10⁵ K and $u_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{\prime }}=0.01$ erg cm⁻³ for the external photon field from fitting the SED, which is the most important target of pγ interaction for neutrino production and for IC scattering, crucial for explaining the VHE spectrum. It is also vital for γγ absorption in the jet, beyond a few hundred GeV. For a typical disk luminosity L_disk ≈ 10⁴⁶ erg s⁻¹ and the scattered disk emission to be a fraction η_ext ∼ 0.01 of the disk photon energy density, R_ext is found to be a few times 10¹⁸ cm.

The VHE flare of December 2018 does not have simultaneous observations at lower energies, and thus the constraints on the theoretical model are moderate. Nevertheless, to reduce the uncertainties, only the electron primary distribution and its corresponding luminosity is changed with respect to the low state. The parameter values used in the modeling are given in Table 1.

We fit the low-state spectrum first, and optimize the parameters δ_D, B′, and spectral indices, using the leptonic emission alone. The SYN spectrum peaks at the optical band and a log-parabola injection spectrum of electrons explains the data well. The hadronic component is then added and the power and maximum proton energy are varied to fit the X-ray data with a BH cascade. The VHE photon flux upper limits constrain the contribution from the pion-decay cascade. We also consider the absorption of VHE γ-rays in the extragalactic background light (EBL) using the Gilmore et al. model (Gilmore et al. 2012). It can be seen from the left panel of Fig. 1 that the neutrino flux is comparable to the 7.5-year averaged flux prediction from this source by IceCube. We find the neutrino flux to be roughly unchanged in modeling the low- and high-state SEDs obtained by the MAGIC multiwavelength campaign. This corresponds to a flux that produces on average one neutrino detection (as for IC-170922A) over a period of 7.5 years and explains the non-observation of neutrino events during December 2018 flare in the MAGIC waveband. A further increase in the neutrino flux leads to the violation of the X-ray data.

3.2. Cosmogenic γ-rays from UHECRs

The maximum proton energy for photohadronic interactions in the AGN frame is $E_{\mathrm{p},\max }=\mathrm{\Gamma }{E\prime }_{\mathrm{p},\max }\approx 0.17$ EeV from modeling multiwavelength SEDs of TXS 0506+056. The proton spectrum has to be cut off at 𝒪 ∼ 0.1 EeV inside the jet to satisfy the constraints from the X-ray data and to simultaneously produce neutrinos with a flux in the PeV range inferred from the detection of one event in 7.5 yr of IceCube operation. According to the Hillas condition, the maximum acceleration energy $E_{\mathrm{p},\max }^{\mathrm{acc}}\sim 2\beta cZeBR\mathrm{\Gamma }$, where the bulk Lorentz factor Γ takes into account the frame transformation from comoving jet frame to AGN frame. The gyration radius of 10²⁰ eV protons from this simplistic expression can be calculated as r_L ≈ 2.13 × 10¹⁶ cm, which is comparable to the blob radius in our modeling. Thus, it is possible to produce protons of energy higher than 0.17 EeV in the same blob for the magnetic field and the length scale considered. The pγ opacity in the jet is higher than 10⁻³ for protons with energy > 6 × 10¹⁷ eV (see, e.g., Das et al. 2022). The jet is opaque to photons beyond ∼1 TeV; therefore, if UHE protons interact inside the jet, additional gamma rays from the π⁰, π^± cascade contribute at lower energies, violating the X-ray data. Hence, we assume that UHE protons beyond this energy escape the jet. A faster escape at higher rigidities has been parameterized in earlier studies (see, e.g., Eq. (20) in Harari et al. 2014) and applied to the case inside the source in Muzio et al. (2022). Thus, it is inherently the same population of protons that interacts below ∼0.1 EeV and escapes at higher energy. This assumption allows us to use the same normalization of the proton spectrum inside and outside the jet, with the same injection spectral index.

Assuming that protons escape efficiently beyond ∼0.1 EeV and up to 10²⁰ eV as UHECRs, we calculate the cosmogenic γ-ray spectrum resulting from their propagation in the EGMF and interactions with the cosmic microwave background and EBL photons. Secondaries from these interactions initiate an electromagnetic cascade in the extragalactic medium, undergoing synchrotron radiation in the EGMF, pair production with EBL, and inverse IC scattering of background photons. The resulting spectrum is shown as the gray line in Fig. 1. We note that the cascade is sufficiently developed for E_p ≳ 40 EeV, the Greisen-Zatsepin-Kuzmin (GZK) energy. Hence, the spectrum depends more on the propagation effects than on the source parameters. The upcoming γ-ray detector Cherenkov Telescope Array (CTA) will detect γ-rays in the range from 20 GeV up to several hundred TeV, with unprecedented sensitivity. The orange dashed curve in Fig. 1 shows the differential point source sensitivity of CTA, assuming 50h of observation time and pointing to 20 degrees zenith (Gueta 2021). CTA observation of a hard and non-variable multi-TeV γ-ray spectrum will indicate the presence of UHECR acceleration inside this source. The blue dashed line is the LHAASO 1 yr sensitivity to Crab-like γ-ray point sources (Vernetto 2016).

We define the line-of-sight resolved component of the cosmogenic γ-ray spectrum as the fraction of UHECRs (ξ_B) that survives within 1° of the initial propagation direction, from the source to the observer, after the deflection in the EGMF. For neutrinos of energy ∼30 TeV, the angular resolution of IceCube for track-like events is ∼0.5°, whereas Fermi-LAT has a resolution of 3.5° to photons of energy < 100 MeV, and ∼0.15° beyond 10 GeV. We use CRPropa-3 to propagate UHECRs in the extragalactic space (Alves Batista et al. 2016). The value of ξ_B is calculated from the arrival direction of protons in a 3D simulation on the surface of a sphere of radius 100 kpc, centered at the observer. We propagate a E⁻² proton spectrum in the energy range between 0.1 and 100 EeV in a random turbulent magnetic field with a Kolmogorov power spectrum. The coherence length of the field is adjusted to 100 kpc. Thus, the secondary flux for a given luminosity is multiplied by ξ_B to obtain the line-of-sight component at the position of the observer (Das et al. 2020). The survival fraction ξ_B as a function of RMS field strength is shown in Fig. 2.

Fig. 2. Survival fraction of UHECRs along the line of sight, i.e., within 0.1° of the initial direction of propagation, as a function of the RMS strength of EGMF.

The value of B_rms can be constrained from the required luminosity in cosmogenic γ-rays by the following expression, under isotropic approximation

$$\begin{array}{c}\hfill {\displaystyle \frac{L_{\mathrm{UHE}p}}{4\pi d_L^2}}={\displaystyle \frac{1}{{\xi }_{\mathrm{B}}{f}_{\gamma ,\mathrm{p}}}}{\displaystyle {\int }_{{ϵ}_{\gamma ,\mathrm{m}in}}^{{ϵ}_{\gamma ,\mathrm{m}ax}}}{ϵ}_{\gamma }{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}{ϵ}_{\gamma }\mathrm{d}A\mathrm{d}t}}\mathrm{d}{ϵ}_{\gamma },\end{array}$$(9)

where dn/dϵ_γdAdt is the differential flux of cosmogenic γ-rays constrained by the SED. We calculate the electromagnetic cascade using the external code DINT integrated with CRPropa-3 (Heiter et al. 2018). The factor f_γ, p takes into account the fraction of injected UHECR power that goes into cosmogenic γ-rays and is fairly constant with the variation of B_rms. For TXS 0506+056 we find f_γ, p ≈ 0.156. The integrated flux of cosmogenic photons, along the line of sight of the observer, allowed by the observed SED is ∼1 × 10⁻¹¹ erg cm⁻² s⁻¹, as seen in Fig. 1. The luminosity of protons interacting inside the jet is found to be L_p = 1.6 × 10⁴⁸ erg s⁻¹ (Table 1). Using the same normalization for the escaping proton spectrum beyond 0.1 EeV, the luminosity in UHECR protons (i.e., in the range 0.1 − 100 EeV) is L_UHEp ≈ 8 × 10⁴⁷ erg s⁻¹. According to Eq. (9), for the allowed flux of cosmogenic γ-rays, this implies ξ_B ≲ 0.05. This indicates an EGMF with an RMS value higher than a few times 10⁻⁵ nG, as shown in Fig. 2. There may be no cosmogenic component along the line of sight for magnetic field strength much greater than this, otherwise, the luminosity budget is violated. It should be noted that this result is an order-of-magnitude estimate. The precise value is sensitive to the angular resolution to detect high-energy γ-rays, the coherence length of EGMF, the angular spread of jet emission, and the numerical precision of cascade calculation.

We calculate the flux of cosmogenic neutrinos produced simultaneously during the propagation of UHECRs and use the same normalization as obtained for the allowed flux of cosmogenic γ-ray spectrum. The cosmogenic neutrino spectrum peaks at 2.8 × 10¹⁸ eV with peak flux ∼4.5 × 10⁻¹³ erg cm⁻² s⁻¹, shown in Fig. 3. The IceCube-Gen2 detector will be capable of detecting neutrinos from TeV to EeV energies, with sensitivity five times higher than the currently operating IceCube experiment. We also show the 5σ sensitivity for detection of muon neutrino flux from TXS 0506+056, using the IceCube-Gen2 detector (IceCube Collaboration 2021). The black and red lines correspond to 100 days and 10 years of observation and indicate the sensitivity for neutrino flares and the time-averaged neutrino emission, respectively. The neutrino flux is lower than the detection threshold, however. Thus the γ-rays provide more stringent limits on UHECR acceleration in TXS 0506+056.

Fig. 3. All-flavor cosmogenic neutrino flux from TXS 0506+056 due to UHECR propagation along the line of sight, i.e., using the same normalization as the cosmogenic γ-ray spectrum in Fig. 1. The black and red curves correspond to 100 days and 10 years of observations of muon neutrino fluxes by IceCube-Gen2 and indicate the sensitivity for neutrino flares and the time-averaged neutrino emission, respectively.

4. Summary and conclusions

The flux variability observed by MAGIC in December 2018 was very similar to that seen in 2017. No neutrino event was detected in 2018, however, in contrast to the 2017 flare. In our one-zone modeling of the multiwavelength data from November 2017 to February 2019, we see that increased γ-ray activity does not yield increased neutrino flux, and the latter is comparable to the 7.5 yr flux prediction by IceCube, the same as that obtained for the low state. Hence, the correlation between γ-ray and neutrino activity is subjective to the model undertaken. In our one-zone modeling, we include three components: the leptonic emission inside the jet, secondary emission due to hadronic cascade induced by protons (E_p ≲ 0.1 EeV), and the line-of-sight resolved cosmogenic γ-rays due to propagation of UHECR protons (E_p ≳ 0.1 EeV) from the blazar to the Earth. The emission from the cascade is found to be subdominant in the GeV range, but important in the X-ray and VHE bands. If the same physical process is responsible for the neutrino production, then the X-ray data constrain the neutrino flux to be consistent with IceCube prediction for one event in ΔT = 7.5 yr. Hence, an additional hidden sector must be invoked to explain a higher neutrino flux, for example, the neutrino flare in 2017 and ΔT = 0.5 yr flux. This is consistent with no observation of γ-ray activity during the orphan neutrino flare during 2014–15.

We highlight our modeling as an important step in the course of study of TXS 0506+056: (i) a neutrino flare (archival data) was observed in 2014–2015, but no gamma-ray activity; (ii) increased gamma-ray activity and simultaneous detection of a neutrino event were observed in 2017; and (iii) there was no neutrino detection, but increased gamma-ray activity was observed in 2018. From the study of phase (iii), we show that indeed there is little correlation between gamma-ray and neutrino flares. Petropoulou et al. (2020) put an upper limit of (0.4 − 2) ν_μ events in ten years of IceCube operation, using multi-epoch modeling. They argue that the IC-170922A can be explained as an upward fluctuation from the average neutrino rate expected from the source, but in strong tension with the 2014–2015 neutrino flare. Rodrigues et al. (2019) shows only 2 − 5 neutrino events during the 2014-15 flare, which can be explained consistently with the X-ray constraints or high-energy γ-ray flux measured by Fermi-LAT. For all the cases presented in the two-zone model of Xue et al. (2019), the neutrino flux prediction is comparable to the 7.5 yr IceCube upper limit. In our study the neutrino event rate is ${\mathcal{N}}_{{\nu }_{\mu }+{\overline{\nu }}_{\mu }}=1\times (\mathrm{\Delta }T/7.5\mathrm{yr})$. Hence, along the lines of the one-zone leptohadronic model adopted here, our results lead to similar conclusions using yet another epoch of the blazar and thus add to the literature.

The lack of simultaneous X-ray data in the high state is a drawback to the SED modeling, although the data points shown are for the nearest observation on 2018 December 8. The parameters are varied minimally from the low state to account for this. Interestingly, our modeling does not predict any significant flare of the optical flux, where the SYN spectrum peaks, but the UV and soft X-ray fluxes are expected to change moderately. We note that a higher X-ray flux can allow for an increased neutrino flux; however, explaining the ΔT = 0.5 yr neutrino flux remains difficult, due to excess X-ray production. Many plausible alternatives exist in the literature, such as neutrino production near the supermassive black hole of the AGN, in the accretion disk, or the corona (Stecker 2013; IceCube Collaboration 2022), or multiple emission zones with increased γγ opacity in the neutrino production zone (Xue et al. 2021). Nevertheless, production of one neutrino event in 0.5 yr is achieved in Cerruti et al. (2019) and Fraija et al. (2020), among others, but the neutrino flux peak is shifted compared to the mean energy of the observed event.

The origin of external photons is a question of fundamental importance in the modeling of TXS-like blazars. In their modeling, MAGIC collaboration used an external field originating from the spine layer or the jet sheath (Ansoldi et al. 2018; Acciari et al. 2022). In our analysis, we consider it to originate from the BLR. It provides a substantial target for neutrino production by pγ processes and also inverse-Compton scattering by electrons. For this to be true, the radius of the emission region must be smaller than the radius of the BLR. In our analysis, R_ext is a few times ∼10¹⁸ cm, which is large compared to usual estimates for the BLR, $R_{\mathrm{BLR}}={10}^{17}{L}_{\mathrm{disk},45}^{0.5}$ cm (see Eq. (4) in Ghisellini & Tavecchio 2008). One possibility is that the blob lies at the edge or even outside the BLR leading to a decrease in the effective BLR photon density (Tavecchio & Ghisellini 2008). The typical energy of the photons in the AGN frame is ϵ_ext ∼ 3k_BT′/Γ ≈ 17 eV. This is comparable to that obtained in other studies (Keivani et al. 2018) and can also be considered as scattered emission from the disk. The contribution from the disk photon itself is negligible. We consider a log-parabola spectrum for the injection of electrons, to improve the fit to the observed SED. Other assumptions are often made in the literature, such as a broken power-law spectrum Xue et al. (2019).

In our modeling, for photohadronic interaction rate to dominate over escape, we need an escape timescale greater than 10⁶(R/c) at ∼1 PeV inside the emission region, and even higher at lower energies, assuming an energy-independent escape. In an energy-dependent parameterization, the escape timescale can be expressed as t_esc ≫ 10⁶(R/c)(E/10³ TeV)⁻¹ for the proton energy range interacting inside the jet. In the one-zone model, the efficient escape of UHECRs requires a rigidity-dependent diffusion rate, for example D(E)∝E² at higher energies (Globus et al. 2008; Harari et al. 2014; Muzio et al. 2022). As an alternative to the step function for the escape, as we had assumed, one can also assume a separate emission zone for the acceleration of UHE protons with lower photohadronic opacity. We do not present the analytical estimates of such an astrophysical scenario in this paper but assume a single proton population. In our analysis, the cosmogenic γ-ray spectrum remains fixed for both the low and high states. A change in the primary proton distribution will not affect the cosmogenic flux significantly because the spectrum is driven greatly by parameters guiding the extragalactic propagation. Since UHECRs are delayed in the EGMF, any observed variability in the VHE regime occurs, most likely due to increased activity inside the jet. The required luminosity in UHE protons can also be translated into a resulting flux of neutrinos at EeV energies. The cosmogenic neutrino flux predicted here, from constraints on γ-ray flux, is found to be lower than that in our earlier study (Das et al. 2022). Thus, detection of cosmogenic neutrinos from TXS 0506+056 seems unlikely with the next generation upgrade of IceCube, leaving ground-based γ-ray detectors such as CTA to test UHECR signature in the SED of blazars.

Acknowledgments

S.D. thanks Konstancja Satalecka (MAGIC Collaboration) for the correspondence regarding the multiwavelength SED data. The work of S.D. was supported by JSPS KAKENHI Grant Number 20H05852. Numerical computation in this work was carried out at the Yukawa Institute Computer Facility. S.R. was supported by a grant from NITheCS and the University of Johannesburg URC. We thank the anonymous referee for useful comments and suggestions.

References

All Tables

All Figures

Fig. 1. Multiwavelength SED of TXS 0506+056 during the 2017–2019 campaign. Left: low-state (average) flux during the observation period. Black data points show the Fermi-LAT average spectrum and MAGIC upper limits, and Swift and NuSTAR spectra on 16 October 2018. Gray data points show the whole range of optical to X-ray spectra. Right: high-state observation in VHE γ-rays during December 2018. Black data points show the Fermi-LAT, MAGIC, and ASAS-SN contemporaneous spectra. Gray data points show the nearest observations in radio, optical, UV, and X-rays on 8 December 2018. The green, blue, and purple curves respectively correspond to SYN, SSC, and EC-BLR emissions. The red curve is the cascade emission from secondary electrons. The magenta curve is the predicted neutrino flux in both cases. The black solid and dashed PeV data points represent the IceCube flux upper limit for one detection in 0.5 yr and 7.5 yr respectively for the 2017 detection event. The gray line is the cosmogenic γ-ray flux from the extragalactic propagation of UHE protons. The orange and blue dashed lines are CTA and LHAASO point source sensitivity. See text for details.

In the text

	Fig. 2. Survival fraction of UHECRs along the line of sight, i.e., within 0.1° of the initial direction of propagation, as a function of the RMS strength of EGMF.
In the text

Fig. 3. All-flavor cosmogenic neutrino flux from TXS 0506+056 due to UHECR propagation along the line of sight, i.e., using the same normalization as the cosmogenic γ-ray spectrum in Fig. 1. The black and red curves correspond to 100 days and 10 years of observations of muon neutrino fluxes by IceCube-Gen2 and indicate the sensitivity for neutrino flares and the time-averaged neutrino emission, respectively.

In the text