Leptohadronic blazar models applied to the 2014-15 flare of TXS 0506+056

TL;DR

This study assesses leptohadronic models to explain neutrino emissions from blazar TXS 0506+056 during 2014-15, finding that models can account for a few neutrino events but struggle to match higher event counts without conflicting with electromagnetic observations.

Contribution

The paper compares different leptohadronic models to explain neutrino emissions from TXS 0506+056, highlighting limitations in simultaneously fitting neutrino and electromagnetic data.

Findings

Models can explain 1-5 neutrino events without conflicting with X-ray data.

Higher neutrino event counts conflict with gamma-ray observations.

No single model can fully reconcile neutrino and electromagnetic observations.

Abstract

We investigate whether the emission of neutrinos observed in 2014-15 from the direction of the blazar TXS 0506+056 can be accommodated with leptohadronic multi-wavelength models of the source commonly adopted for the 2017 flare. While multi-wavelength data during the neutrino flare are sparse, the large number of neutrino events ($13\pm5$) challenges the missing activity in gamma rays. We illustrate that two to five neutrino events during the flare can be explained with leptohadronic models of different categories: a one-zone model, a compact core model, and an external radiation field model. If, however, significantly more events were to be accommodated, the predicted multi-wavelength emission levels would be in conflict with observational X-ray constraints, or with the high-energy gamma ray fluxes observed by the Fermi LAT, depending on the model. For example, while the external radiation field model can predict up to five neutrino events without violating X-ray constraints, the absorption of high-energy gamma rays is in minor tension with data. We therefore do not find any model that can simultaneously explain the high event number quoted by IceCube and the (sparse) electromagnetic data during the neutrino flare.

Figures

Tables

Table 1. Table 1: Selected parameter sets for each model, predicted number of neutrino events N ν subscript 𝑁 𝜈 N_{\nu} and SED quality compared to data. The corresponding neutrino rate is given by N ν / T subscript 𝑁 𝜈 𝑇 N_{\nu}/T , where T = 158 ​ days 𝑇 158 days T=158~{}\mathrm{days} is the duration of the neutrino flare. The values for 1-zone (a) and (b) are given for two representative curves from Fig.  1 (red curve from the left panel and green curve from the right panel). The physical luminosities L e , phys obs superscript subscript 𝐿 𝑒 phys obs L_{e,\mathrm{phys}}^{\mathrm{obs}} and L p , phys obs superscript subscript 𝐿 𝑝 phys obs L_{p,\mathrm{phys}}^{\mathrm{obs}} carried by electrons and protons are given by L phys obs = L iso / Γ 2 superscript subscript 𝐿 phys obs subscript 𝐿 iso superscript Γ 2 L_{\mathrm{phys}}^{\mathrm{obs}}=L_{\mathrm{iso}}/\Gamma^{2} , where L iso subscript 𝐿 iso L_{\mathrm{iso}} is the isotropic-equivalent luminosity. The physical luminosities can be compared to the Eddington luminosity, which is 4 × 10 46 ​ erg / s 4 superscript 10 46 erg s 4\times 10^{46}~{}\mathrm{erg/s} for a black hole mass of 3 × 10 8 ​ M ⊙ 3 superscript 10 8 subscript 𝑀 direct-product 3\times 10^{8}~{}M_{\odot} estimated by Padovani et al. ( 2019 ) ; note, however, that this value can be temporarily exceeded during flares.

	Quality criteria		Parameters
	$N_{\nu }$	SED	$B^{\prime }$	$R_{\mathrm{blob}}^{\prime }$	${\Gamma }_b$	${\gamma }_{p,\mathrm{obs}}^{\max }$	$L_{e,\mathrm{phys}}^{\mathrm{obs}}$	$L_{p,\mathrm{phys}}^{\mathrm{obs}}$	$f_{\mathrm{esc}}$	$T_{\mathrm{disk}}$	$L_{\mathrm{disk}}$	$R_{\mathrm{BLR}}$
Model			[G]	[cm]			[$\mathrm{erg}/s$]	[$\mathrm{erg}/s$]		[K]	[$\mathrm{erg}/s$]	[cm]
1-zone (a)	$1.8$	Compatible	$1.0$	${10}^{17}$	$9.0$	${10}^{5.8}$	${10}^{44.2}$	${10}^{49.6}$	${10}^{-2.5}$	–	–	–
1-zone (b)	$13.2$	Overshoot	$0.001$	${10}^{15}$	$7.0$	${10}^{5.7}$	${10}^{44.8}$	${10}^{50.7}$	${10}^{-2.5}$	–	–	–
C. core (blob)	$0.0$	Compatible	$0.01$	${10}^{18.7}$	$10.0$	–	${10}^{44.8}$	–	${10}^{-2.5}$	–	–	–
C. core (core)	$1.9$			${10}^{15}$		${10}^{6.1}$	${10}^{43.7}$	${10}^{49.5}$	${10}^{-2.5}$	–	–	–
Ext. field (a)	4.9	Compatible	0.6	${10}^{15.8}$	49.1	${10}^{5.9}$	${10}^{43.6}$	${10}^{48.7}$	${10}^{-4.8}$	${10}^{5.7}$	${10}^{44.7}$	${10}^{17.8}$
Ext. field (b)	4.0	Cutoff, 10 GeV	0.9	${10}^{16.3}$	48.0	${10}^{6.3}$	${10}^{42.9}$	${10}^{48.4}$	${10}^{-3.1}$	${10}^{5.3}$	${10}^{44.7}$	${10}^{17.3}$

Full text

Multi-messenger, blazar, astrophysical neutrinos, TXS 0506+056

1 Introduction

The object TXS 0506+056 is an Active Galactic Nucleus (AGN) of the blazar type, classified as a BL Lac object, with a measured redshift of $z=0.3365$ (Paiano et al., 2018). In September 2017, a muon neutrino with a reconstructed energy of about 290 TeV was observed by IceCube from a position compatible with this source in coincidence with a period of flaring in multiple wavelengths (Aartsen et al., 2018a) at a significance level of $3\sigma$. This event has enticed the multi-messenger community to explore the potential of TXS 0506+056 as a source of astrophysical neutrinos. The connection between neutrino production and the electromagnetic flare has been described by several leptohadronic ($p\gamma$) production models (Gao et al., 2018; Cerruti et al., 2018; Zhang et al., 2018; Keivani et al., 2018; Ansoldi et al., 2018; Sahakyan, 2018; Gokus et al., 2018; Righi et al., 2018), as well as hadronic ($pp$) production models (Liu et al., 2018; Sahakyan, 2018). For example, it has been concluded by Gao et al. (2018) that conventional one-zone models describing the spectral energy distribution and the neutrino event suffer from too low neutrino rates in combination with excessively high neutrino energies or sustained super-Eddington injection luminosities. The current theoretical consensus is that the geometry of the radiation zone must be more complex, involving a compact radiation core with high photohadronic interaction rates (Gao et al., 2018), or external radiation fields boosted into the jet frame, either thermal (Keivani et al., 2018) or non-thermal (Ansoldi et al., 2018).

Triggered by the multi-messenger discovery of the 2017 flare, IceCube searched their archival data for an excess from the direction of TXS 0506+056 (Aartsen et al., 2018b) for the entire duration of IceCube’s data taking. In the period between October 2014 and March 2015, a temporal clustering has been detected of 64 events in total within $3^{\circ}$ of the direction of the same source. By using a likelihood function, in which the atmospheric background is taken into account and the signal is assumed to be distributed as a power law, a 3.5$\sigma$ excess over the atmospheric background was found, with an estimated number of signal events of $13\pm 5$ (henceforth the “historical neutrino flare”). The most energetic event has a deposited energy of 20 TeV in IceCube, while most events have energies around $\sim$ 10 TeV. Interestingly, this signal was not accompanied by any significant increase in electromagnetic emission. In contrast to the 2017 flare, the multi-wavelength data from this period are very sparse and the only constraints on the spectral energy distribution (SED) can be derived from gamma-ray flux measurements by the Fermi LAT (Garrappa et al., 2019), as well as radio and optical monitoring data compiled by Padovani et al. (2018). Additionally, the Swift Burst Alert Telescope (BAT) that monitors X-ray transients and performs regular sky surveys, was not triggered during the period of the neutrino flare and did not detect TXS 0506+056 in the 15–50 keV band, implying that its flux during the neutrino flare was significantly less than 3 mCrab, or $7.2\times 10^{-11}\text{erg}~{}\text{cm}^{-2}~{}\text{s}^{-1}$ (Krimm et al., 2013). Based on Fermi data, Padovani et al. (2018) have speculated that there may be a hardening in the SED of the source above 2 GeV during the neutrino flare, although this feature may in fact not be significant (Garrappa et al., 2019). Theoretical models for the historical flare are sparse, facing the challenge that the high neutrino flux has to be accommodated with the inconspicuous SED activity (Murase et al., 2018). A possible way out could be jet-cloud/star ($pp$) interactions (Bednarek & Protheroe, 1997; Barkov et al., 2010; Wang et al., 2018), whereas the photohadronic model in Halzen et al. (2018) does not contain a self-consistent SED computation.

In this letter, we present a theoretical analysis of the neutrino and electromagnetic emission during the historical flare of TXS 0506+056. We focus on the available observational evidence during the neutrino emission period only – which was presented above. Due to the lack of enhanced gamma-ray activity, we treat the historical flare independently from the 2017 event. Motivated by the limited constraints from long-term multi-wavelength data, we do not attempt to derive a time-dependent model that explains the transitions between the neutrino bright and quiet states. The multi-messenger SED is computed using the self-consistent numerical code AM3 (Gao et al., 2017), that has been successfully applied in the interpretation of the 2017 flare (Gao et al., 2018), and that has been extended by the inclusion of external radiation fields. Apart from a conventional one-zone model, we test two other classes, namely an inverse-Compton dominated compact-core model and a model involving an external radiation field from accretion disk radiation isotropized in a broad-line region (BLR)– considering a scenario in which this source possesses features typical of Flat-Spectrum Radio Quasars (FSRQs), that are out-shined by the radiation from the jet.

2 Methods

We construct the models for the simultaneous emission of neutrinos and photons using the leptohadronic code AM3, solving self-consistently the time-dependent kinetic equations for non-thermal electrons, positrons, protons, neutrons, photons and neutrinos produced in the relativistic jet. The production region is simulated as a spherical blob of radius $R^{\prime}_{\mathrm{blob}}$ in its rest frame111Primed quantities refer to the blob rest frame, unprimed quantities to the observer’s frame, moving along the blazar jet with a bulk Lorentz factor $\Gamma_{\mathrm{b}}$. The assumption is that protons and electrons are accelerated to a power-law spectrum $dN/dE^{\prime}\propto E^{\prime-2}$, up to certain maximum Lorentz factors $\gamma_{e,\mathrm{max}}^{\prime}$ and $\gamma_{p,\mathrm{max}}^{\prime}$, and are then injected isotropically into a radiation zone of the jet. They interact with the target photons according to Hümmer et al. (2010), producing charged and neutral pions that ultimately decay into neutrinos and secondary gamma rays, electrons and positrons. These particles will feed into the electromagnetic cascade and potentially lead to signatures in the SED. The other interactions included in the model are electron synchrotron emission and synchrotron self-absorption, inverse Compton (IC) scattering by both electrons and protons, photon pair production and annihilation, and Bethe-Heitler pair production, $p\,\gamma\,\rightarrow\,p\,e^{+}\,e^{-}$. The electron synchrotron emission depends on the strength of the turbulent magnetic field in the radiation zone, $B^{\prime}$, considered randomly oriented.

The critical quantity of interest is the number of neutrinos predicted by the model, which we compute by folding the emitted fluence with the effective area given by Aartsen et al. (2018b) for the IC86b data period. For the single neutrino observed during the 2017 flare, the expected number of detectable neutrinos is likely smaller than one for different reasons (Eddington bias in Strotjohann et al. (2018) or too many associations expected in Palladino et al. (2018)). These arguments do not apply to the historical flare, where the predicted event number needs to be significantly larger than one to be compatible with observations. Since for blazars the number density of target photons (X-rays for this particular source) is low compared to other compact objects such as gamma-ray bursts, the optical depth to $p\gamma$ interactions is typically much lower than unity – which needs to be compensated by a large proton loading in order to become a significant neutrino source. The photohadronic interaction rate can be enhanced by assuming a smaller production region or a higher density in X-rays. The latter can be achieved by an external radiation field boosted into the blob frame, such as the X-ray photons initially produced by the accretion disk and scattered by the dust or cloud surrounding the jet.

We perform extensive scans within physically plausible ranges for the parameter space (for $10^{15.0}<R^{\prime}_{\mathrm{blob}}/\mathrm{cm}<10^{17.0}$, $10^{-3}<B^{\prime}/G<10$, $5<\Gamma_{\mathrm{b}}<50$, $0.06<E_{\mathrm{p,inj}}^{\prime}/\mathrm{PeV}<15$ and $10^{2.8}<L_{\mathrm{p,inj}}/L_{\mathrm{e,inj}}<10^{6.3}$). Yet, we cannot claim completeness of our scans because of the complexity of the problem. The parameter space was searched using two methods: a grid-based parameter scan and a genetic algorithm (Goldberg, 1989). The goodness of fit for each parameter set is defined according to a simple $\chi^{2}$-criterion in $\nu F_{\nu}$ between the simulated SED, and the optical, gamma-ray and X-ray constraints. Since a rigorous minimization is not feasible due to the sparsity of the data, we choose the “neutrino-loudest” areas of the parameter space.222A lepto-hadronic blazar model involves $N\sim 10$ parameters. While some parameters are correlated, a $\chi^{2}$ goodness of fit estimator produces a highly degenerate likelihood space given the underconstrained nature of the problem (only seven data points). The emitted neutrino spectra, which significantly differ from power-laws, are convolved with the effective area of IceCube at the declination of the source (IC86b data period, Aartsen et al., 2018b) to obtain the predicted number of muon track events, which is then compared to the observed signal.

3 Results

We first test a conventional one-zone model, where the radiation zone consists of a single spherical blob. The neutrinos in the jet escape the blob over the free-streaming timescale $t^{\prime}_{\mathrm{FS}}=R^{\prime}_{\mathrm{blob}}/c$ (and likewise for the photons and neutrons that survive the interactions). Charged particles escape the blob at a slower rate due to the magnetic confinement. For simplicity, we implement an energy-independent escape rate for charged particles of $t^{\prime}_{\mathrm{esc}}=f_{\mathrm{esc}}\,t^{\prime}_{\mathrm{FS}}$, with $f_{\mathrm{esc}}>1$.333We also tested alternative scenarios; (a) $f_{\mathrm{esc}}=1$ and (b) a harder proton injection spectrum $dN/dE\propto E^{-1.5}$. Both yield similar SEDs to that illustrated by the red curve in the left panel of Fig. 1 if $L_{\mathrm{p,inj}}$ is increased by a factor of 45 in case (a), and decreased by a factor of 3.8 in case (b). Fig. 1 clearly demonstrates that the models compatible with the SED (left panel) produce too few neutrinos, where at most 1.8 events are expected during the duration of the neutrino flare (red curve, parameters listed in Tab. 1). This number is limited by the X-ray constraint on the SED, which we derive from the non-detection by Swift BAT. The two bumps around the X-ray limit come from synchrotron and IC emission off e± that originate from $\gamma\gamma$ annihilation at higher energies, and from Bethe-Heitler pair production. This example demonstrates the importance of electromagnetic data across the entire spectrum to constrain theoretical models, since the electromagnetic cascade accompanying the neutrino production can be hidden in unconstrained energy ranges (such as MeV).

On the other hand, a compatible neutrino flux level implies an SED that is in tension with observations (right panel). These neutrino-compatible SEDs belong to a class of models with a strong hadronic cascade. Note that the self-consistently computed SED is very different from the ad hoc assumption in Halzen et al. (2018), and peaks at lower energies. We also find a cluster of strongly IC-dominated solutions, due to a compact emission region and low magnetic field strength, supporting a high $p\gamma$ efficiency and hence higher neutrino fluxes. These solutions cannot, however, explain the emission outside the Fermi LAT range. On the other hand, the subset of models with sufficient synchrotron emission fail to simultaneously comply with the X-ray and gamma-ray constraints. We conclude that one-zone models are in tension with observations of the historical flare. The corresponding model parameters are listed in Tab. 1.

As suggested earlier, a smaller emission region can enhance neutrino production. Following the model in Gao et al. (2018) for the 2017 flare, we speculate that during the neutrino flare a compact core is formed inside the larger emission region, sharing its Doppler factor. This scenario can be regarded as a (spatially) structured jet model where the resulting emission of photons and neutrinos originates from the superposition of both radiation zones. In the present case, the core is a highly p$\gamma$-efficient region that simultaneously explains the gamma-ray and neutrino emission with a suppressed synchrotron cascade; most electromagnetic radiation at lower energies originates from the larger blob region. The so-called spine-sheath model assumes in addition a velocity structure that allows for finer control over the multi-messenger emission at the cost of a higher number of free parameters (Ansoldi et al., 2018).

The results for the compact core model are shown in Fig. 2 for one optimized set of parameters; see also Tab. 1. Emission from the blob describes the data from radio to soft gamma rays, while its contribution above GeV energies is low. The large volume of the blob translates into small target photon densities for photo-hadronic interactions, leading to inefficient neutrino production and a dim hadronic cascade (effectively leading to a leptonic model). The higher radiation densities in the core create an IC-dominated hadronic cascade, which accounts for a gamma-ray emission that hardens above 10 GeV. The suppression of synchrotron emission, in combination with small hadronic and Bethe-Heitler cascades, can suppress X-ray emission to a minimum. This example represents a neutrino-efficient model that is not strongly constrained by X-ray emission. Parameter sets can also be found that yield more neutrino events, but are in our view unphysical, *e.g.*stronger magnetic fields in the blob than in the core.

The compact core model has also been previously applied to the 2017 flare (Gao et al., 2018). However, it requires additional fine-tuning of the core and the blob parameters to explain the temporal correlation among the optical, X-ray and gamma-ray flares. As a side effect, the compact core model can also reproduce a fast gamma-ray variability (for instance through a modulation of the core size that would have no effect on the radio to X-ray bands). However, the transition from a neutrino-quiescent to a neutrino-loud state without signature in gamma rays would require fine-tuning in the temporal evolution of the different parameters. The neutrino flux emitted by the core translates to 1.9 observed muon tracks, which is slightly higher than in SED-compatible one-zone models; however, it is still in tension with the IceCube result.

Finally, we consider the impact of an external thermal field, similarly to what has been assumed by Keivani et al. (2018) for the 2017 flare. While this source has been identified as a BL Lac due to the lack of significant broad line emission, Padovani et al. (2019) argue that TXS 0506+056 is a masquerading BL Lac that includes broad lines and thermal emission from an accretion disk as for FSRQs, which are outshined by the beamed non-thermal emission from the jet.

In that case, a fraction of the accretion disk radiation is isotropized through Thomson scattering in a BLR surrounding the disk (we fix this fraction to 1%) and a more significant fraction (around 10%, Greene & Ho, 2005) is re-emitted as atomic broad lines. If the emitting region lies within the BLR, these components will appear boosted in the jet frame and interact with the non-thermal particles (see e.g.  Rodrigues et al., 2017, regarding the relativistic transformations). The thermal continuum is modeled as a blackbody emission of temperature $T_{\mathrm{disk}}\sim 10^{5}~{}\mathrm{K}$ (Bonning et al., 2007); broad line emission is represented by the H$\alpha$ line, which is typically the brightest. The maximum proton energy is adjusted such that interactions with the thermal continuum result in neutrinos with energy $E_{\nu}\sim 100~{}\mathrm{TeV}(T/10^{5}~{}\rm{K})^{-1}$. Due to photon annihilation, the external fields also attenuate gamma rays from the jet as they cross the BLR, with maximum attenuation at $E_{\gamma}\sim 10~{}\mathrm{GeV}(E_{\nu}/100~{}\mathrm{TeV})^{-1}$. At this redshift, photons with energy $E>300$ GeV suffer additional loss due to the interaction with the extra-galactic background light (EBL, modeled in Domínguez et al. (2011)), hence the steep cutoff of the SED shown in the figures.

The results of the external field model are presented in Fig. 3. The thermal field (orange) is out-shined by the highly beamed jet radiation (red) and is invisible to the observer. The solid red curve represents parameter set (a) in Tab. 1 and leads to 4.9 neutrino events during the flare. In this model, the hump observed in the optical band originates from synchrotron emission by $e^{\pm}$ pairs from Bethe-Heitler production, while the hump in the MeV-GeV range is emitted by $e^{\pm}$ pairs from the annihilation of hadronic photons. The emission from primary electrons is therefore sub-dominant across the spectrum, which implies that it is difficult to identify such a model in simplified (such as analytical) approaches. However, as mentioned above, the high neutrino production efficiency implies a softening and a suppression of the gamma-ray spectrum above 10 GeV. A spectral softening is in tension with Fermi observations (Garrappa et al., 2019; Padovani et al., 2018). The X-ray bound is almost saturated, as well. One of the specifics of this model is the anti-correlation between VHE gamma-ray and neutrino emission, which has been previously discussed in Murase et al. (2016) (see also Xue et al. (2018)). The parameter set (b) in Fig. 3 yields 4 neutrino events; it has a slightly lower X-ray component at the cost of a higher tension with the last Fermi data point. The attenuation of high-energy gamma rays between models depends not only on the disk luminosity but also on the assumed radius of the BLR (cf. Tab. 1). Note that given the disk luminosity of the source, phenomenological relationships would suggest a BLR radius of around $3\times 10^{17}~{}\mathrm{cm}$ (Kaspi et al., 2000). In the examples shown, the values of $R_{\mathrm{BLR}}$ differ from that reference value by a factor of two or less, which is within the statistical spread of the AGN sample reported by Kaspi et al. (2000).

4 Summary and conclusions

We have tested the compatibility of the historical 2014-15 neutrino flare from the blazar TXS 0506+056 with leptohadronic (photohadronic) multi-messenger source models. Within the constraints of the sparse observations in optical, X-rays and gamma rays during the neutrino flare, we scanned the parameter space using several distinct assumptions about the geometry and environment of the blazar. In addition to conventional one-zone models, we have considered scenarios involving a compact core emission region (corresponding to a spatially structured jet) and external radiation fields, such as that from an accretion disk.

We have demonstrated that at most two to five neutrino events during the period of the flare can be expected from any of the three models in compatibility with multi-wavelength constraints. While the one-zone model saturates the available X-ray bound, the SED at gamma-ray energies can be reasonably reproduced. The electromagnetic cascade from charged and neutral pion decays can be hidden as a prominent hump at MeV energies, where no data is available.

A compact core model yields a similar expectation for the neutrino rate as the most optimistic one-zone model, accompanied by a spectral hardening in gamma rays above 10 GeV. The radiation at highest energies and the neutrinos both originate from the inverse-Compton-dominated core, whereas the X-ray data are generated by a larger emission region via a leptonic synchrotron self-Compton (SSC) process. A natural feature is a faster gamma-ray variability from the compact core compared to the slower radio-to-X-ray variability. However, a transition between neutrino-quiescent and flaring states would imply a fine-tuned correlation in the evolution of the parameters.

The external radiation field model yields SEDs with more than two neutrino events during the flare; however, the high neutrino production efficiency implies a higher optical thickness to $\gamma\gamma$ annihilation at VHEs, softening the expected gamma-ray spectrum – in minor tension with Fermi observations. Such a softening or cutoff can serve as an electromagnetic signature for an orphan neutrino flare.

While we do not claim completeness, our study demonstrates the obstacles involved in the simultaneous description of the electromagnetic SED and neutrino observations during the historical flare. Since all present models are challenged by the high neutrino event rate, we have not even attempted to describe the temporal evolution, i.e. the transition between the neutrino-quiescent and flaring states, or to achieve a unified description between the 2014-15 and 2017 flares. However, from our modeling and extensive parameter scans, we conclude that a) obtaining more than two to five neutrino events during the flare implies violating multi-wavelength constraints, particularly in gamma rays, and that b) a transition between the neutrino-quiescent and flaring states without distinctive electromagnetic activity is unlikely within the photohadronic framework. Let us finally remark that the $13\pm 5$ signal events quoted by IceCube are obtained under the assumption of a power-law spectrum with a spectral index $-2.2$. For a harder spectrum or even a different spectral shape, the signal emerging over the atmospheric background may be significantly smaller.

The authors would like to thank an anonymous reviewer for their useful comments. This work has been supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 646623).