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

MAGIC (Aleksić et al. 2016) is a system of two 17 m ground-based IACTs, operating in stereoscopic mode since fall 2009 at the Roque de los Muchachos Observatory (288 N, 178 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 2017 September 24 (MJD 58020) until 2017 October 4 (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 Aartsen et al. (2018). In this Letter, we present additional MAGIC data on TXS 0506+056 collected until 2017 November 2 (MJD 58059), adding up to a total exposure of ∼47 hr. After data-quality cuts based on the atmospheric conditions, telescope response, and stereo event rate, about 41 hr 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.

Zoom In Zoom Out Reset image size

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 (2017 October 3–4, already reported in Aartsen et al. 2018), and a second one on MJD 58057 (2017 October 31). Detection of TXS 0506+056 with MAGIC was also achieved (at the level of 5.7σ) by stacking all of the data with a similarly low VHE gamma-ray flux level (F(E > 90 GeV) <4.5 × 10⁻¹¹ cm⁻² s⁻¹) measured outside of 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)).

Zoom In Zoom Out Reset image size

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–400 GeV by a simple power law ${dN}/{dE}\approx {E}^{\gamma }$, 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 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.

Comprehensive MWL coverage of TXS 0506+056 during MJD 58019−30 is reported in Aartsen et al. (2018). For this Letter, 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), Ultraviolet and Optical Telescope (UVOT) and X-Ray Telescope (XRT) on board 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 hr 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 (Aartsen et al. 2018) 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 (Tavecchio et al. 2010).

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 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 electron-positron pair production interactions with low-energy photons ($\gamma \gamma \to {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 escapes, 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 Bethe–Heitler pair production ($p\gamma \to {{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, as 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 Letter.

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 Domínguez et al. (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₁ = 2 and n₂) 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}_{\mathrm{jet}}=\pi {R}^{2}{{\rm{\Gamma }}}^{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 (Ghisellini et al. 2010). For the two states, we assume the same values for the jet bulk Lorentz factor Γ_j = 22, viewing angle θ_view = 08 (implying δ ≃ 40), sheath bulk Lorentz factor Γ_s = 2.2, and the radius of the emission region R = 10¹⁶ cm (corresponding to a variability timescale of one day). Note that we take θ_view < 1/Γ_j, where the EC emission appears more luminous relative to the synchrotron, SSC, and hadronic cascade emission than at θ_view ∼ 1/Γ_j, which is 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}_{{\rm{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}_{{\rm{p}},\max }$ suitable for reproducing the observations. An upper limit of ${E}_{{\rm{p}},\max }\simeq {10}^{18}$ eV is expected as the theoretical maximum allowed by plausible mechanisms of particle acceleration (see Section 4). A lower limit of ${E}_{{\rm{p}},\max }\simeq {10}^{14}$ eV in the comoving frame is imposed to allow production of a neutrino with observed energy ∼290 TeV, as lower values of ${E}_{{\rm{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 Bethe–Heitler cascade emission, as well as jet energetics arguments.

Table 3.  Parameters for the Jet-sheath Model for ${E}_{{\rm{p}},\max }={10}^{16}$ (eV)

State	MJD 58029−30	Lower VHE
B (G)	2.6	2.6
Emin (eV)	3.2 × 108	2.0 × 108
${E}_{\mathrm{br}}$ (eV)	7.0 × 108	9.0 × 108
Emax (eV)	8 × 1011	8 × 1011
n1	2	2
n2	3.9	4.4
Ue (erg cm−3)	4.4 × 10−4	3.6 × 10−4
UB (erg cm−3)	0.27	0.27
Up (erg cm−3)	1.8	0.7
Pe (erg s−1)	2 × 1042	1.6 × 1042
Ep (erg s−1)	8 × 1045	3 × 1045
PB (erg s−1)	1.2 × 1045	1.2 × 1045

Download table as:  ASCIITypeset image

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}_{{\rm{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 day-timescale 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⁻⁴. This implies that for the target photons in such pγ interactions, the γγ optical depth is ∼0.1 at the corresponding threshold ${E}_{\gamma \gamma }\simeq {m}_{e}^{2}{c}^{4}/\epsilon \simeq (20{m}_{e}^{2}/{m}_{\pi }{m}_{p}){E}_{\nu }\,\simeq 12$ GeV (E_ν/290 TeV). As the density of target photons is roughly inversely proportional to (Figure 2), significant γγ absorption 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}_{{\rm{p}},\max }={10}^{16}$ eV, which also yields the highest predicted neutrino event rates. Both higher and lower values of ${E}_{{\rm{p}},\max }$, in conjunction with X-ray and gamma-ray constraints, still allow acceptable solutions, with lower neutrino rates, as shown in Figures 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 Aartsen et al. (2018), we predict 90% confidence level lower and upper limits to the muon neutrino energy of 206 TeV and 6.3 PeV for ${E}_{{\rm{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; Aartsen et al. 2018). 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 gamma-ray emission (Aartsen et al. 2018). 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. Note that 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 gamma-rays. 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, and electron density) are in the range commonly derived for BL Lac objects similar to TXS 0506+056 based on purely leptonic models (Ghisellini et al. 2010). On the other hand, the photo-meson cascade and Bethe–Heitler 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⁴⁵ − 4 × 10⁴⁶ 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¹⁴–10¹⁸ eV in the comoving frame.