The Blazar TXS 0506+056 Associated with a High-energy Neutrino: Insights into Extragalactic Jets and Cosmic-Ray Acceleration

A neutrino with energy ∼290 TeV, IceCube-170922A, was detected in coincidence with the BL Lac object TXS 0506+056 during enhanced gamma-ray activity, with chance coincidence being rejected at ∼3σ level. We monitored the object in the very-high-energy (VHE) band with the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes for ∼41 hr from 1.3 to 40.4 days after the neutrino detection. Day-timescale variability is clearly resolved. We interpret the quasi-simultaneous neutrino and broadband electromagnetic observations with a novel one-zone lepto-hadronic model, based on interactions of electrons and protons co-accelerated in the jet with external photons originating from a slow-moving plasma sheath surrounding the faster jet spine. We can reproduce the multiwavelength spectra of TXS 0506+056 with neutrino rate and energy compatible with IceCube-170922A, and with plausible values for the jet power of . The steep spectrum observed by MAGIC is concordant with internal γγ absorption above ∼100 GeV entailed by photohadronic production of a ∼290 TeV neutrino, corroborating a genuine connection between the multi-messenger signals. In contrast to previous predictions of predominantly hadronic emission from neutrino sources, the gamma-rays can be mostly ascribed to inverse Compton upscattering of external photons by accelerated electrons. The X-ray and VHE bands provide crucial constraints on the emission from both accelerated electrons and protons. We infer that the maximum energy of protons in the jet comoving frame can be in the range ∼1014 – 1018 eV.


INTRODUCTION
The birthplace of ultra-high-energy cosmic rays (UHECRs), the most energetic particles known in the Universe with energies exceeding 10 18 eV, is a longstanding mystery (Dawson et al. 2017). The observed distribution of their arrival directions in the sky favor a predominantly extragalactic origin (Aab et al. 2017). Among the numerous candidate sources that have been proposed (Hillas 1984), one of the most promising are active galactic nuclei (AGN) that eject powerful jets of magnetized plasma at relativistic velocities (Biermann & Strittmatter 1987). Blazars are AGN with their jets oriented close to the line of sight of the observer. Their predominantly non-thermal spectral energy dis-to the energy range of UHECRs (Bykov et al. 2012). These hadrons can interact with ambient matter or lowenergy photons to generate gamma-rays, as well as neutrinos, as envisaged in hadronic and lepto-hadronic scenarios of blazar emission (for reviews see Halzen 2016;Meszaros 2017). Neutrinos can be considered the "smoking gun" of hadron acceleration. In contrast to cosmic rays that may be deflected by intervening magnetic fields while propagating to the observer, photons and neutrinos are expected to point back to their origin in the sky, providing critical insights into the sources of cosmic rays.
In September 2017, the IceCube neutrino observatory revealed an event designated IceCube-170922A with 56.5% probability of being a truly astrophysical neutrino. The best-fit reconstructed direction was at 0.1 • from the sky position of the BL Lac object TXS 0506+056 1 . The most probable energy was found to be 290 TeV (311 TeV), with the 90% C.L. lower and upper limits being 183 TeV (200 TeV) and 4.3 PeV (7.5 PeV) respectively, assuming a spectral index of -2.13 (2.0) for the diffuse astrophysical muon neutrino spectrum (Aartsen et al. 2014). No additional excess of neutrinos with lower energy was found from the direction of TXS 0506+056 near the time of the alert (Paper 1).
Extensive follow-up observations revealed TXS 0506+056 to be active in all electromagnetic (EM) bands, most notably in GeV gamma-rays monitored by the Fermi-LAT (Large Area Telescope; Y. T. Tanaka 2017), and in veryhigh-energy (VHE) gamma-rays above 100 GeV, detected for the first time from this object by the MAGIC (Major Atmospheric Gamma-ray Imaging Cherenkov) telescopes (Mirzoyan 2017). The redshift has been recently measured to be z = 0.3365 ± 0.0010 (Paiano et al. 2018).
The probability for a high-energy neutrino to be detected by chance coincidence with a flaring blazar from Fermi-LAT catalogs was found to be disfavored at a 3σ confidence level, mostly due to the precise determination of the direction of IceCube-170922A (Paper 1).
These measurements offer a unique opportunity to explore the interplay between energetic photons, neutrinos and cosmic rays. 2 Hereafter we interpret within a coherent scenario the available multi-messenger data pre-sented in Paper 1, together with additional MAGIC data for TXS 0506+056, under the assumption that the association of IceCube-170922A and the blazar in an active state is genuine, and that the neutrino and EM emission arise from the same region in the object. Given the appreciable probability of IceCube-170922A being of atmospheric origin, and that chance coincidence is rejected only at ∼ 3σ level, we note that the validity of the association is still an open question.
SED data strictly simultaneous with the neutrino event is not available for all wavelengths, and the duration of the active state responsible for the neutrino emission is uncertain. On the other hand, the data from Fermi reveals enhanced gamma-ray activity of TXS 0506+056 lasting several months. We consider physical models of the source that are compatible with the detection of one neutrino by IceCube during 0.5 years in the energy range reported in Paper 1. By invoking a dense field of low-energy photons originating outside the jet as targets for photo-hadronic interactions, the measured neutrino event can be interpreted consistently with the EM observations. Furthermore, while the detection of the ∼ 290 TeV neutrino alone indicates acceleration of protons in the jet of this object to energies of at least several times 10 15 eV, the combination with the EM data allow us to probe the maximum energy that they can attain.

VHE GAMMA-RAY AND BROADBAND EMISSION OF TXS 0506+056
The rapid dissemination of the sky position of IceCube-170922A triggered an extensive multi-wavelength (MWL) campaign by many telescopes on ground and in space, from radio frequencies up to VHE gammarays. In the VHE band, TXS 0506+056 was followed up by several Imaging Atmospheric Cherenkov Telescopes (IACTs), with the earliest one starting a few hours after the neutrino alert (Paper 1).

MAGIC VHE gamma-ray Observations and Results
MAGIC (Aleksić et al. 2016) is a system of two 17m ground-based IACTs, operating in stereoscopic mode since fall 2009 at the Roque de los Muchachos Observatory (28.8 • N, 17.8 • W, 2200 m a.s.l.), on the island of La Palma, Canary Islands (Spain). The telescopes record Cherenkov light from extended air shower events starting from 30 GeV up to ∼100 TeV within a field of view of ∼10 square degrees.
MAGIC observations of the sky position of IceCube-170922A from September 24 (MJD 58020) till October 4, 2017 (MJD 58030) yielded the first detection of a VHE signal from an astrophysical object in a direction compatible with that of the high-energy neutrino, as described in Paper 1. In this work, we present additional MAGIC data on TXS 0506+056 collected until November 2, 2017 (MJD 58059), adding up to a total exposure of ∼47 hours. After data quality cuts based on the atmospheric conditions, telescope response and stereo event rate, about 41 hours of data were selected for further analysis, employing the standard MAGIC analysis framework MARS. Parameter cuts optimized for low energies and tuned with data from the Crab Nebula were used (Aleksić et al. 2016). Results for gamma-rays above 90 GeV are given in Figure 1 and Table 1.
The VHE gamma-ray flux is variable, increasing by a factor of up to ∼6 within one day. The probability of a constant flux is less than 0.3%. Two periods of enhanced VHE gamma-ray emission can be clearly distinguished: one on MJD 58029-58030 (Oct 3-4, 2017, already reported in Paper 1) and a second one on MJD 58057 (Oct 31, 2017). Detection of TXS 0506+056 with MAGIC was also achieved (at the level of 5.7σ) by stacking all data with a similarly low VHE gamma-ray flux level (F(E>90 GeV) < 4.5 × 10 −11 cm −2 s −1 ) measured outside these two periods. This allowed us to distinguish between the states with high and low VHE gamma-ray emission. Note that, the low state is unlikely to exemplify the typical quiescent state of the source according to the historical data in Figure 2 c).
We could not reveal any difference in the spectral shape in periods with different flux levels or data samples. The measured differential photon spectrum can be described over the energy range of 80 GeV to 400 GeV by a simple power law dN/dE ≈ E γ , with a spectral index ranging from (-4.0 ± 0.3) to (-3.5 ± 0.4), see Table  2. This is significantly steeper than the spectrum measured by Fermi-LAT quasi-simultaneously, with spectral index (-2.0 ± 0.2) in the energy range 100 MeV to 300 GeV. The systematic uncertainties of MAGIC spectral measurements can be divided into: < 15% in energy scale, 11-18% in flux normalization and ± 0.15 for the energy spectrum power-law slope (Aleksić et al. 2016). In a view of the lack of evidence for spectral variability within the data sets reported here, we adopt the spectral index derived from the stacked data as a common value for the lower VHE state.

Multi-wavelength SED
Comprehensive MWL coverage of TXS 0506+056 during MJD 58019-30 is reported in Paper 1. For this work, we analyzed additional open access MWL data with public tools during the period from MJD 58019 to MJD 58059 from the following instruments: KVA (Lindfors et al. 2016), UVOT and XRT onboard the Neil Gehrels Swift observatory (Swift) (Roming et al. 2005), NuSTAR (Harrison et al. 2013) and Fermi-LAT (Atwood et al. 2009). For each waveband considered, data within 24 hours around the MAGIC observations were combined and considered as quasi-simultaneous for the modeling. Figure 2 a) shows the SED restricted to the first period of VHE gamma-ray enhanced emission observed from MJD 58029 to 58030. Good coverage of simultaneous MWL data is also found for the MAGIC observations at the lowest flux in Figure 2 b). The second period of enhanced VHE gamma-ray flux (MJD 58057) has limited MWL coverage and is not used for the modeling below.
Comparison of the light curve in Figure 1 with that measured by Fermi-LAT (Paper 1) suggests the low VHE state to be more representative of the VHE average emission during the six months of enhanced GeV flux. On the other hand, the VHE emission of the flare state may be considered an upper bound. By compiling broadband SEDs quasi-simultaneous with the MAGIC observations, we extrapolate such considerations to all other wavelengths. The obtained SEDs appear typical of BL Lac objects with similar luminosities ).

INTERPRETING THE BROADBAND AND NEUTRINO EMISSION
Production of high-energy neutrinos in astrophysical environments is expected to occur mainly through the decay of charged pions produced in inelastic collisions between high-energy protons and ambient target particles (Halzen 2016;Meszaros 2017), which can be either matter (pp interactions) or low-energy photons (pγ interactions). In AGN jets that typically consist of lowdensity plasma, the latter is generally the favored channel (e.g. Mannheim 1995). Photo-hadronic production of neutrinos with energy E ν ∼ 290 TeV requires interactions between protons with energy E p 6 PeV and target photons with energy exceeding the corresponding threshold, m π m p c 4 /E p 440 eV, that is, in the UV to soft X-ray band.
FSRQs are generally considered to be promising sources of high-energy neutrinos, on account of the inferred high density of target photons for the pγ channel (Atoyan & Dermer 2001;Murase et al. 2014). In contrast, BL Lac objects such as TXS 0506+056 have generally been thought to be inefficient neutrino emitters (Murase et al. 2014), due to the relatively low density of UV to soft X-ray synchrotron photons expected inside the jet. Thus, the association of a high-energy neutrino with TXS 0506+056 is not trivial to interpret. Mod-els in which hadronic emission components are more prominent for the gamma-rays may allow more neutrino emission, at the expense of relatively large values for the power in accelerated protons that may not be so favorable from an energetics perspective (see e.g. Cerruti et al., in preparation, for the specific case of TXS 0506+056).
Higher efficiency of neutrino production without invoking high values of the proton power may still be achievable by considering target photon fields external to the jet. One such scenario in which external photons naturally occur for BL Lac objects is the so called spinelayer or jet-sheath scenario, as suggested by Tavecchio et al. (2014) (see also Righi et al. 2017). This model postulates structured jets, consisting of a faster core (or spine) surrounded by a slower sheath (or layer). Such a velocity structure is supported by observational and phenomenological evidence, as well as numerical simulations (Tavecchio et al. 2014). Due to significant relative motion between the two structures, the density of photons originating from the slower sheath appears highly boosted in the frame comoving with the faster jet (Ghisellini et al. 2005), exceeding that of the synchrotron photons arising within the jet, and thus enhancing the rate of pγ reactions and consequent neutrino emission. These external photons serve not only as targets for pγ interactions, but also as seed photons for IC upscattering, so that EC emission may also play an important role.
As long as neutrinos are produced via the pγ channel, an important consequence is the associated absorption of gamma-rays via γγ interactions (Waxman & Bahcall 1999;Dermer et al. 2007). Since the cross section for γγ interactions is on average ∼ 10 3 times larger than for pγ interactions, the presence of target photons that provide a reasonable production efficiency of neutrinos at a given energy inevitably implies strong gamma-ray absorption above a corresponding threshold. Detection of such a spectral feature constitutes a critical test of photohadronic neutrino production as well as the neutrinoblazar association itself.

Model description
Hereafter we model the broadband EM and neutrino emission from TXS 0506+056, building on the jet-sheath scenario of Ghisellini et al. (2005); Tavecchio et al. (2014). These earlier studies focused on the leptonic emission and/or the neutrino emission and lacked explicit treatment of hadronic radiative processes. Here we present the first extension of the jet-sheath model to account for all relevant hadronic processes, as well as a first comparison of the model with multi-messenger observations of a particular object.
The radiative output is calculated in the reference frame comoving with the jet flow, assuming isotropic distributions for both electrons and protons. This is then transformed to the observer frame via the Doppler factor δ, which accounts for effects due to the jet's relativistic motion (Madejski & Sikora 2016).
The calculation of the leptonic emission processes, which include synchrotron, SSC and EC emission, follows Ghisellini et al. (2005). We note that in EC models such as the jet-sheath scenario, the angular distribution of IC seed photons in the jet comoving frame will be highly anisotropic due to relativistic boosting. This leads to a more narrowly peaked beaming pattern for the EC emission compared to the synchrotron and SSC emission, (Dermer 1995;Ghisellini et al. 2005).
We account for all relevant hadronic radiative processes, including photo-meson-induced cascade emission (Mannheim 1993), Bethe-Heitler (BH) pair cascade emission (Petropoulou & Mastichiadis 2015), and synchrotron radiation from protons (Aharonian 2000) and muons. The former is the most important component in our model. For protons of sufficiently high energy, the gamma-rays resulting from photo-meson reactions can be energetic enough to undergo electronpositron pair production interactions with low-energy photons (γγ → e + e − ), thereby triggering secondary pair cascades (Mannheim 1993). This process redistributes the energy of photons down to lower energies until it falls below the threshold for γγ interactions and escape, generally resulting in a spectrum with roughly equal power over a broad energy range (from optical to VHE for the cases shown below). From the spectral properties of the external radiation fields considered here, γγ transparency is expected below a few hundreds of GeV, the energy band best accessible to MAGIC. The direct products of pγ reactions including neutrinos and photons are described using the analytical treatment of Kelner & Aharonian (2008). The secondary pair cascading processes were implemented following the formalism of Boettcher et al. (2013) and extending it to include IC in addition to synchrotron processes for the pairs. An additional hadronic component arises from BH pair production (pγ → pe + e − ), which can be observationally relevant in some cases. Gamma-ray emission via synchrotron radiation of protons or muons is mostly unimportant for the range of magnetic fields and proton energies considered here, except for a limited region of parameter space. The results presented here have been cross-checked with the codes developed in Cerruti et al. (2015) and Zech et al. (2017).
As mentioned above for EC emission, in the jet comoving frame, target photons from the sheath should appear highly anisotropic. So far, there has been little discussion in the literature concerning photohadronic neutrino and cascade emission in anisotropic photon fields. The anisotropy may induce an effect analogous to EC for the neutrino emission, but the non-trivial propagation of charged pions and muons in magnetic fields before their decay complicates the picture. On the other hand, the accompanying cascade emission is expected to be mostly isotropic, since the secondary e ± pairs are likely isotropized by magnetic fields as they degrade in energy. Thus, the neutrino emission may be more narrowly beamed compared to the cascade emission. For the viewing angles and magnetic fields considered here, we estimate that this effect may enhance the observed neutrino flux by a factor of ∼2-3 relative to the case assuming isotropy. However, this is not explicitly included in our calculations, as a full investigation of the effect requires Monte Carlo simulations that are beyond the scope of this paper.
The main parameters of the model specifying the jet environment (magnetic field, Doppler factor, size of acceleration/emission region) and the particle populations (density and energy distribution of electrons and protons) are adjusted to consistently reproduce the MWL data and the observed neutrino event rate. The absorption of gamma-rays during their propagation to the observer caused by γγ interactions with the extragalactic background light is modeled following Dominguez (2011), and is expected to be minor at the measured redshift of TXS 0506+056 at energies below 400 GeV (∼10% at 400 GeV and ∼50% at 4 TeV). Figure 2 shows results for models corresponding to the low and high states (see details in the caption). Table  3 lists their parameters (all in the jet comoving frame): jet magnetic field B, parameters specifying the electron distribution (a broken power-law with minimum energy E min , break energy E br and maximum energy E max , and spectral indices n 1 = 2 and n 2 ), and the energy densities in relativistic electrons U e , magnetic fields U B , and relativistic protons U p . The jet power is evaluated as P jet = πR 2 Γ 2 c(U e + U p + U B ), for which the individual contributions from relativistic electrons P e , magnetic fields P B , and relativistic protons P p are also given in Table 3  . For the two states, we assume the same values for the jet bulk Lorentz factor Γ j = 22, viewing angle θ view = 0.8 • (implying δ 40), sheath bulk Lorentz factor Γ s = 2.2, and the radius of the emission region R = 10 16 cm (corresponding to a variability timescale of one day). Note that we take θ view < 1/Γ j , where the EC emission appears more lu-minous relative to the synchrotron, SSC and hadronic cascade emission than at θ view ∼ 1/Γ j , the viewing angle typically assumed in the literature, by virtue of the former's more narrowly peaked beaming pattern (Dermer 1995). The energy distribution of relativistic protons is assumed to be a simple power-law with spectral index 2 and an exponential cutoff at maximum energy E p,max . We limit B to a range such that muons and pions produced in hadronic interactions do not suffer significant energy losses before decaying. We explore a range of E p,max suitable for reproducing the observations. An upper limit of E p,max 10 18 eV is expected as the theoretical maximum allowed by plausible mechanisms of particle acceleration (see Section 4). A lower limit of E p,max 10 14 eV in the comoving frame is imposed to allow production of a neutrino with observed energy ∼290 TeV, since lower values of E p,max imply significantly less pγ interactions above the corresponding energy threshold. Compensating this with higher proton power is not feasible, as the latter is limited by X-ray constraints on BH cascade emission, as well as jet energetics arguments.

Results
Despite the large number of model parameters, the combined constraints from the EM and neutrino channels are quite powerful in assessing acceptable regions of the parameter space. Particularly strong constraints on the broadband hadronic cascade emission, and therefore on E p,max in conjunction with the primary proton luminosity, come from the combination of the X-ray and gamma-ray bands, which also constrain the spectrum of primary leptons. The low state SED in Figure 2 b) is more tightly constraining, thanks to the additional data from the NuSTAR satellite, which cover the range in between the two main SED components.
VHE gamma-ray data from MAGIC provide several important constraints for the modeling. First, the size of the emission region is limited by the observed daytimescale variability, which is unresolvable in the Fermi-LAT data. Second, the VHE spectrum depends sensitively on the energy distribution of primary electrons, in particular its break and maximum energies. Finally and most crucially, strong spectral steepening due to internal γγ absorption is robustly predicted in the VHE band, as an inevitable byproduct of photohadronic production of a ∼290 TeV neutrino. The fraction of multi-PeV parent protons that undergo pγ interactions to produce ∼290 TeV neutrinos is ∼ 10 −4 . This implies that for the target photons in such pγ interactions, the γγ optical depth is ∼ 0.1 at the corresponding threshold E γγ m 2 e c 4 / (20m 2 e /m π m p )E ν 12 GeV (E ν / 290 TeV). As the density of target photons is roughly inversely proportional to (Fig. 2), significant γγ ab-sorption is expected above ∼ 100 GeV. The steep spectrum measured by MAGIC relative to Fermi-LAT indeed matches this expectation, providing a unique confirmation of the one-zone lepto-hadronic interpretation, the pγ production channel, as well as the association between the neutrino and the blazar.
A good fit to the data is found for an intermediate value of E p,max = 10 16 eV, which also yields the highest predicted neutrino event rates. Both higher and lower values of E p,max , in conjunction with X-ray and gamma-ray constraints, still allow acceptable solutions, with lower neutrino rates, as shown in Figure 2 d) and e).
The muon neutrino fluxes predicted in our model are shown in Figure 2, considering neutrino oscillations into equal fractions among the three flavors during their propagation (Aartsen et al. 2015). Convolved with the neutrino effective area reported in Paper 1, we predict 90% confidence level lower and upper limits to the muon neutrino energy of 206 TeV and 6.3 PeV for E p,max = 10 16 eV. This predicted range is in good agreement with the observed one, namely 183 TeV (200 TeV) and 4.3 PeV (7.5 PeV) for a spectral index of -2.13 (2.0) (Paper 1). The expected detection rate of muon neutrinos in the above energy interval for the models in Table  3 are about 0.17 events in 0.5 years for the higher VHE emission state, and 0.06 events in 0.5 years for the lower VHE emission state. These numbers are conservative in that they do not account for a potential contribution to the event rates due to interactions of tau-neutrinos that induce muons with a branching ratio of 17.7%. Both rates are in agreement with the detection of the single neutrino during the period of enhanced GeV gammaray emission (Paper 1). As mentioned in Section 3.1, the neutrino fluxes may be higher by a factor of ∼ 2 − 3 if the anisotropy of the target photon field for pγ interactions is properly taken into account.
The neutrino flux is also constrained by limits on the hadronic emission from the EM observations, most effectively from the X-ray and VHE bands.
In agreement with many previous studies of the broadband emission of blazars, the SEDs here comprise mostly leptonic emission. In particular, the second SED component peaking in the GeV band is largely accounted for by EC emission, for which the enhanced luminosity relative to the other emission components at θ view < 1/Γ is important for achieving acceptable fits. This is in contrast to some earlier studies of neutrino emission from blazars that predicted predominantly hadronic gammarays. Our models show that this is not necessarily the case, and neutrino emission can be commensurate with the widely accepted view of primarily leptonic gamma-ray emission from blazars. This implies that the values of the parameters that suitably reproduce the SEDs (in particular magnetic field, Doppler factor, electron density) are in the range commonly derived for BL Lac objects similar to TXS 0506+056 based on purely leptonic models ). On the other hand, the photo-meson cascade and BH components must generally be subdominant.
From the comparison of the model parameters for the low and high states, the enhanced emission can be attributed to changes in the energy distribution of the electrons in the jet (harder high-energy slope and larger average energy) and increased power of the proton population. Such variability is quite consistent with the behavior commonly displayed by BL Lac objects during active states. In both cases, the power of the jet is in the range 10 45 to 4 × 10 46 erg/s. These values are somewhat larger than that derived for BL Lac objects through purely leptonic models that do not account for accelerated protons (Ghisellini et al. 2009). The ratio of energy density in protons to electrons is ∼ 3600 and ∼ 1700 for the high and low states respectively, similar to that inferred in various astrophysical environments.
From the combined EM and neutrino data, we infer that the maximum energy of protons in the blazar emission region is in the range 10 14 −10 18 eV in the comoving frame.

IMPLICATIONS FOR UHECR ACCELERATION IN AGN JETS
As shown above, the multi-messenger observations of TXS 0506+056 can be consistently interpreted by considering the acceleration of electrons and protons in the same region of the jet, which interact with a conceivable source of underlying external photons. The detection of the single neutrino may provide the first compelling evidence for proton acceleration in AGN jets, while the MWL observations constrain key parameters such as the magnetic field and the Doppler factor of the relevant region. With minimal additional considerations from plausible theories of particle acceleration, we are in a unique position to infer the spectral distribution of protons in AGN jets, and discuss their viability as sources of UHECRs.
Various mechanisms have been proposed for accelerating UHECRs in AGN jets, involving strong shocks, magnetic reconnection, shear flows, etc. (Bykov et al. 2012;Sironi et al. 2015). Independent of the specific mechanism, as viewed in the comoving frame of the jet, the timescale for accelerating particles up to energy E in a relativistic flow can be expressed as t acc = 10ηE/(ZeBc), where e is the elementary charge, c is the velocity of light, and η is a parameter that characterizes the properties of magnetic disturbances responsible for the acceleration. For mildly relativistic shock velocities as expected in AGN jets, η 1 corresponds to the so-called Bohm limit whereby the turbulence is fully developed and allows the fastest possible acceleration. To be compared are timescales for loss processes due to adiabatic expansion t ad = 2R/c, synchrotron radiation t syn (E), photo-meson production t γπ (E), and BH pair production t γe ± (E) (e.g. Cerruti et al. 2015). The latter two depend on the spectra of the target photon fields, and can be evaluated separately for the internal photons (synchrotron from accelerated electrons) and external photons considered here following Petropoulou & Mastichiadis (2015). The maximum attainable proton energy E p,max can be estimated by balancing t acc (E) with Z = 1 and the shortest loss timescale, t loss = min[t ad , t syn , t γπ , t γe ± ]. Moreover the particles gyro-radius r g = E/ZeB should not exceed R so that it remains confined in the acceleration region. Figure 3 shows the comparison of the timescales and the gyro-radius constraint for the model of the high state discussed above. For η 1, we see that protons can be accelerated up to E p,max 10 18 eV, limited by adiabatic and photo-meson losses on external photons. All values of E p,max 10 18 eV discussed above can be entirely compatible, with lower values of E p,max corresponding to the balance of t ad and t acc with larger values of η that may reflect lower levels of turbulence (e.g. E p,max 10 16 eV for η = 100; see Fig.3).
To contribute to the observed UHECRs, particles must escape from the acceleration region into ambient intergalactic space. The energy of escaping protons E p,max as seen by an observer will appear boosted by the relativistic bulk motion of the jet so that E p,max = Γ j f el E p,max on average, where f el ≤ 1 is a factor that accounts for energy loss during the escape process. For our models with Γ j = 22 and E p,max 10 14 − 10 18 eV, E p,max 2 × 10 15 − 2 × 10 19 eV, as long as f el ∼ 1. The higher end of this range would be well within the domain of UHECRs. If nuclei heavier than protons are also accelerated in the same conditions, even higher energies may be possible.

DISCUSSION AND CONCLUSIONS
The combined MWL and neutrino observations of TXS 0506+056 can be interpreted consistently in terms of electrons and protons that are accelerated in the same region of the jet and interact with external low-energy photons, originating from the slower sheath surrounding the faster jet spine. Such jet-sheath structure is independently well motivated. The enhanced luminosity for EC emission relative to the other emission components at small viewing angles is important for adequately reproducing the MWL SEDs. Most notably, the prominent spectral steepening observed at a few tens of GeV by MAGIC confirms the internal γγ absorption that is robustly expected as a consequence of pγ production of a ∼290 TeV neutrino. Thus, our interpretation reinforces the association between the multi-messenger signals. While the bulk of the gamma-rays are inferred to be external Compton emission from electrons, a nonnegligible contribution can arise from cascade emission induced by protons, most notably in the hard X-ray and VHE gamma-ray bands.
Jet structure as envisaged in the jet-sheath scenario can be conclusively tested with future very long baseline interferometry (VLBI) measurements with high angular resolution. We note that an alternative source of external photons in BL Lac objects may be provided by radiatively inefficient accretion flows around the central supermassive black holes (Righi et al, in preparation).
An interesting question concerns the potential contribution of the entire BL Lac population to the diffuse neutrino emission detected by IceCube (Aartsen et al. 2017). In principle, BL Lac objects and FSRQs could account for a sizable fraction of the observed intensity (Tavecchio et al. 2014;Righi et al. 2017). The association of IceCube-170922A with the flaring phase of TXS 0506+056 is consistent with this scenario, although the relative rarity of these sources may imply additional constraints on their contribution (Murase & Waxman 2016).
The inferred maximum energy of protons in the blazar emission region may be consistent with an important contribution to UHECRs from protons and/or heavy nuclei accelerated in the blazar region. Although stronger constraints are not possible with a single observed neutrino event, future multi-messenger observations of blazars will offer a more critical probe of UHECR acceleration in the inner regions of AGN jets.    Figure 2. Spectral energy distribution for the enhanced VHE gamma-ray emission state (a), MJD 58029 to 58030) and the lower VHE gamma-ray emission state (LS, b)) modeled with the jet-sheath scenario with Ep,max = 10 16 eV. Symbols corresponding to data-points from different facilities and observation epochs are described in the legend. The curves represent individual emission components while the thick black curve shows the total predicted emission. The leptonic emission from the jet includes synchrotron (dashed blue), synchrotron-self-Compton (SSC, red dashed dotted), and external Compton (EC) emission (dark red dotted). Synchrotron emission from the sheath is denoted by the light green dashed line. The hadronic emission components are photo-meson-induced cascade (purple dotted), Bethe-Heitler pair cascade (yellow dash-dotted) and muon-synchrotron (yellow dotted). Predicted (anti-)neutrino spectra are marked by (light-)magenta (dashed) solid lines, the blue vertical line shows the energy ∼290 TeV of the observed neutrino. A comparison of the two solutions is also shown with the archival data from ASDC (c)). Results for different values of Ep,max are compared for the enhanced VHE gamma-ray emission state (d), MJD 58029 to 58030) and the lower VHE gamma-ray emission state (Low state, e)).  . Timescales for proton acceleration (tacc with Z = 1) when η = 1 and η = 100, compared with those for adiabatic loss t ad , synchrotron loss tsyn, photo-meson loss tγπ, and BH pair production loss t γe ± , as functions of the proton energy Ep in the comoving frame, for the high state. Also shown is the energy Ep,g = eBR where the proton gyroradius rg = R.